%~Mouliné par MaN_auto v.0.27.3 2023-03-10 10:46:10
\documentclass[CRPHYS, Unicode, XML, screen, HideArticleType, published]{cedram}
\makeatletter \cdr@french@issue@titletrue \makeatother



\usepackage{amssymb}
\usepackage{mathtools}
\usepackage{graphicx}
\usepackage{array}
\makeatletter\cdr@showtranslationtrue\makeatother 


\newcommand{\rmin}{\mathrm{in}}
\newcommand{\rout}{\mathrm{out}}
\newcommand{\rmest}{\mathrm{est}}
\newcommand{\rmfeed}{\mathrm{feed}}
\newcommand{\rmbare}{\mathrm{bare}}
\newcommand{\rmeff}{\mathrm{eff}}
\newcommand{\Sr}{\mathrm{Sr}}
\newcommand{\rmY}{\mathrm{Y}}

%\newcommand{\mathcolorbox}[2]{\colorbox{#1}{$\displaystyle #2$}}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
\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}}}
\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}
\let\oldexists\exists
\renewcommand*{\exists}{\mathrel{\oldexists}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\RubricEF{CNRS Gold Medal Jean Dalibard}{Médaille d'or du CNRS Jean Dalibard}


\title{Ultra narrow linewidth frequency reference via measurement and feedback}
\alttitle{Étalon de fréquence à largeur de raie étroite par mesure et rétroaction}

\keywords{Quantum optics, Superradiance, Cavity quantum electrodynamics, Feedback, Lasers}
\altkeywords{Optique quantique, Superradiance, Électrodynamique quantique en cavité, Rétroaction, Lasers}

\author{\firstname{Diego} \lastname{Barberena}\IsCorresp}
\address{JILA, NIST, Department of Physics, University of Colorado, Boulder, CO 80309, USA}
\address{Center for Theory of Quantum Matter, University of Colorado, Boulder, CO 80309, USA}
\email[D. Barberena]{diego.barberena@colorado.edu}
\thanks{This work is supported by the AFOSR Grant No. FA9550-18-1-0319, by the DARPA (funded via ARO) Grant No. W911NF-16-1-0576, the ARO single investigator Grant No. W911NF-19-1-0210, the NSF PHY1820885, NSF JILA-PFC PHY-1734006 and NSF QLCI-2016244 grants, by the DOE Quantum Systems Accelerator (QSA) grant and by NIST}

\author{\firstname{Robert} \middlename{J.} \lastname{Lewis-Swan}}
\address{Homer L. Dodge Department of Physics and Astronomy, The University of Oklahoma, Norman, OK 73019, USA}
\address{Center for Quantum Research and Technology, The University of Oklahoma, Norman, OK 73019, USA}
\email[R. J. Lewis-Swan]{lewisswan@ou.edu}

\author{\firstname{Ana Maria} \lastname{Rey}}
\address[1]{JILA, NIST, Department of Physics, University of Colorado, Boulder, CO 80309, USA}
\address[2]{Center for Theory of Quantum Matter, University of Colorado, Boulder, CO 80309, USA}
\email[A. M. Rey]{arey@jilau1.colorado.edu }

\author{\firstname{James} \middlename{K.} \lastname{Thompson}}
\address[1]{JILA, NIST, Department of Physics, University of Colorado, Boulder, CO 80309, USA}
\email[J. K. Thompson]{jkt@jila.colorado.edu}


\begin{abstract}
The generation of very narrow linewidth light sources is of great importance in modern science. One such source is the superradiant laser, which relies on collectively interacting ultra long lived dipoles driven by incoherent light. Here we discuss a different way of generating spectrally pure light by coherently driving such dipoles inside an optical QED cavity. The light exiting cavity carries information about the detuning between the driving light and the atomic transition, but is also affected by the noise originating from all the decoherence processes that act on the combined atom-cavity system. We calculate these effects to obtain fundamental limits for frequency estimation and stabilization across a range of values of input light intensities and atom-light interaction strengths, estimate these limits in state-of-the-art cavity experiments with alkaline-earth atoms and identify favorable operating conditions. We find that the achievable linewidths are comparable to those of the superradiant laser.
\end{abstract}

\begin{altabstract}
L'obtention de sources lumineuses à largeur de raie étroite revêt une grande importance dans la science moderne. Une de ces sources est le laser superradiant, qui met en jeu des dipôles à très longue durée de vie interagissant collectivement et forcés par un champ lumineux incohérent. Nous discutons ici d'une autre manière d'obtenir une source spectralement pure, par forçage cohérent de tels dipôles dans une cavité optique QED (en régime de couplage fort). Le champ qui sort de la cavité est porteur d'informations sur le désaccord en fréquence entre le champ de forçage et la transition atomique, mais il est également affecté par le bruit dû aux processus de décohérence à l'oeuvre dans le système combiné atomes-cavité. Nous tenons compte de ces effets pour déterminer les limites fondamentales sur la mesure de la fréquence et sa stabilisation, en fonction des intensités lumineuses d'entrée et des amplitudes de couplage atomes-champ, puis nous estimons ces limites dans l'état de l'art des expériences en cavité sur des atomes alcalino-terreux et nous déterminons les régimes de fonctionnement favorables. Nous trouvons que les largeurs de raie accessibles sont comparables à celles du laser superradiant.
\end{altabstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\dateposted{2023-06-13}
\begin{document}
\maketitle

\section{Introduction}

Sources of coherent light with high frequency stability are crucial components of modern day technologies. In line with these requirements, the proposal for a continuous superradiant laser\cite{Meiser2009,CHE09,Bohnet2012,Norcia2015b,norcia2016superradiance,Schaffer2020,Zhang2021,Kazakov,Cline2022} could provide an improved active frequency reference that is largely insensitive to cavity frequency noise, which limits modern day most stable lasers~\cite{Kessler2012}. Such a device comprises an electromagnetic mode of an optical cavity that is made to interact collectively with an ensemble of atoms that host a long lived two-level transition (ground and excited states). The atoms are incoherently driven by exciting the ground state to an auxiliary third level that decays quickly onto the excited state, and this is done in a regime where the cavity linewidth is much larger than the bare atomic linewidth. Above a critical value of the incoherent pump rate, the interplay between the incoherent drive and the atom-light interactions synchronize the atomic array forcing it to collectively emit coherent light with a linewidth set by the system cooperativity times the atomic linewidth. Although the superradiant laser can open great opportunities for quantum metrology applications, the need of a strong incoherent pump can significantly heat the atoms making the experimental implementation of continuous superradiance in optical transitions challenging. Up to date, superradiance has been realized in a pulsed way, in metrologically relevant transitions~\cite{Norcia2018Super,Laske2019}, or quasi-continuously, in transitions with broader linewidths~\cite{Norcia2016c}. In this article, we analyze this same system in a complementary situation, where the long lived transition is excited directly with a laser. The light coming outside of the cavity can then be used to estimate the detuning between the driving light and the atomic transition, possibly allowing for frequency stabilization through a feedback scheme.

This system has been extensively analyzed in the past under the name of cooperative fluorescence, with studies focusing on its steady state properties\cite{Carmichael_1980,Drummondo1980,Walls1980}, dynamics\cite{Mielke:97,DRUMMOND1978}, bistability and hysteresis\cite{BONIFACIO1976,Gripp1996}, correlation functions of the output light\cite{Foster2000}, etc; and has seen a resurgence of interest from recent experiments\cite{Ferioli2022} that probe the superradiant phase transition\cite{Carmichael_1980,Walls1980} in free space. There has also been previous work in frequency stabilization at large input fields\cite{Martin2011}, where the atomic ensemble is strongly saturated, or in situations where the atomic transition is affected by motion~\cite{Westergaard2015,Tieri2015}. Here we give a comprehensive description of this system's ability to provide information on the frequency of the driving light across many values of input light intensity and effective atom-light interaction. Correctly assessing this requires two pieces: (1) the spectral response of the system to frequency fluctuations in the input light, and (2) the noise properties of the output light, which is affected by the presence of the atoms and is in general much noisier than that of a purely coherent light source.

Our results allow us to evaluate the advantages and disadvantages (for frequency stabilization) of operating the system in different parameter regimes. We determine the spectral purity of the emitted light after feedback stabilization and observe that it can be compatible with the one generated by a superradiant laser. Our analysis also leads us to the conclusion that artificially increasing the effective spontaneous emission rate by means of an extra depumping process may sometimes be advantageous in terms of frequency stability.

\section{Ideal Model}

In this section we will analyze the interaction of an ensemble of atoms with a cavity QED mode in the absence of single particle decoherence and in the presence of a constant laser drive [see Fig.~\ref{fig:Schematic}$\MK$(a)]. We discuss the ideal case first because formulas are less cumbersome and the logic behind the steps can be explained more smoothly. The inclusion of the extra decoherence in section~3 will lead to technical complications, but the spirit behind the calculations will be the same. The ideal case will also constitute a standard with which to compare later results.

The analysis we will pursue here consists of four parts. First, we compute the non-equilibrium steady state of the system within the mean field approximation. Second, we calculate the linear response to two kinds of perturbations: (i) fluctuations in the input light frequency and (ii) quantum noise of the input light. Third, using the response of the output light to (i) we construct an estimator that allows us to infer the frequency of the driving light, while (ii) gives us the quantum noise in this estimator and provides a measure of its efficacy. We finalize by using (i) within a closed feedback loop to calculate the effective linewidth of the frequency stabilized input light source.

The system we are studying consists of an ensemble of $N$ two-level atoms with atomic transition $\omega_a$ collectively coupled to a cavity QED mode with resonance frequency $\omega_c$ and power decay linewidth $\kappa$[see Fig.~\ref{fig:Schematic}$\MK$(a)]. The Hamiltonian describing this system is (in units where $\hbar=1$):
\begin{equation}
\hat{H}=\sum_{i=1}^N\frac{\omega_a}{2}\left(1+\sigma_i^z\right) +\omega_c\hat{a}^\dagger\hat{a}+g \sum_{i=1}^N\left(\hat{a}\hat{\sigma}_i^++\hat{a}^\dagger\hat{\sigma}_i^-\right),
\end{equation}
where $\hat{\sigma}_i^{x,y,z}$ are Pauli matrices describing the two-level system of atom $i$, $\hat{a}$ ($\hat{a}^\dagger$) is an annihilation (creation) operator describing the cavity mode and $2g$ is the single photon Rabi frequency. We further add a laser drive with instantaneous phase $\phi(t)$ [i.e. instantaneous frequency $\omega_d(t)=\dot{\phi}(t)$] and flux of $\alpha_{\rmin}^2$ photons per second. To account for the properties of the light outside the cavity, we describe the evolution of the system in terms of Heisenberg--Langevin equations\cite{Gardiner1985}, written in the lab frame as:
\begin{equation}
\begin{split}
\partial_t\hat{a}&=-\left(i\omega_c+\frac{\kappa}{2}\right)\hat{a}-ig\hat{S}^-+\sqrt{\kappa}\left[\alpha_{\rmin}e^{i\phi(t)}+\delta\hat{A}_{\rmin}(t)\right]\\
\partial_t\hat{S}^-&=-i\omega_a\hat{S}^-+2ig\hat{S}_z\hat{a}\\
\partial_t\hat{S}_z&=ig\left(\hat{a}^{\dagger}\hat{S}^--\hat{a}\hat{S}^{+}\right),
\end{split}
\end{equation}
where $\hat{S}_{x,y,z}=\sum_{i=1}^N \hat{\sigma}_i^{x,y,z}/2$ are collective spin operators and $\delta\hat{A}_{\rmin}(t)$ accounts for the quantum fluctuations of the input field. Formally, $\delta\hat{A}_{\rmin}$ acts on the Hilbert space associated to the continuum of modes outside the cavity and is constructed of Heisenberg operators in the infinite past\cite{Gardiner1985}. Given that the input state is coherent, it is entirely characterized by: $\langle\delta\hat{A}_{\rmin}(t)\rangle=0$, and $\langle\delta\hat{A}_{\rmin}(t)\delta\hat{A}_{\rmin}^\dagger(t')\rangle=\delta(t-t')$, with all other second order expectations equal to 0. The output field, which describes the properties of the light outside of the cavity in the infinite future, can then be calculated using the input-output relation $\hat{A}_{\rout}(t)=\alpha_{\rmin}e^{i\phi(t)}+\delta\hat{A}_{\rmin}(t)-\sqrt{\kappa}\hat{a}$.

In what follows we assume that the cavity is locked to the drive [$\omega_c=\omega_d(t)=\dot\phi(t)$] and close to resonance with the atomic transition [see Fig.~\ref{fig:Schematic}$\MK$(c)]. In the rotating frame of the drive/cavity, the equations of motion take the simpler form:
\begin{equation}\label{eqn:Heis-LangeIdeal}
\begin{split}
\partial_t\hat{a}&=-\frac{\kappa}{2}\hat{a}-ig\hat{S}^-+\sqrt{\kappa}\left[\alpha_{\rmin}+\delta\hat{A}_{\rmin}(t)\right]\\
\partial_t\hat{S}^-&=i\Delta_t\hat{S}^-+2ig\hat{S}_z\hat{a}\\
\partial_t\hat{S}_z&=ig\left(\hat{a}^{\dagger}\hat{S}^--\hat{a}\hat{S}^{+}\right),
\end{split}
\end{equation}
where $\Delta_t\equiv\dot{\phi}(t)-\omega_a=\omega_d(t)-\omega_a$ is the atom-drive detuning, which we assume to be small. Locking the drive to the atoms corresponds to having $\Delta_t=0$, but in a real laser $\Delta_t$ is a fluctuating quantity. We thus need to estimate it and correct it towards 0 [see Fig.~\ref{fig:Schematic}$\MK$(d)].

\begin{figure}
\centering
\includegraphics[width=0.9\textwidth]{Figure1_v6.pdf}
\caption{(a) Schematic of the system: the cavity is driven by coherent light ($\alpha_{\rmin}$) and interacts with the atoms with strength $g$. The atoms are subject to dephasing ($\gamma_d$), spontaneous emission into free space ($\gamma$) and an additional depumping process [$\gamma_p$, see panel (b)]. The light coming out of the cavity ($\alpha_{\rout}$) is measured and this information is used to apply feedback on the input light source. (b) The depumping process can be engineered by means of a rapidly decaying auxiliary third level. (c) The cavity is frequency locked to the input light [$\omega_c=\omega_d(t)$], whose frequency fluctuates in time. They are detuned by $\Delta_t$ from the atomic transition. (d) Schematic of the feedback scheme: when $\Delta_t=0$ the system reaches a steady state. Small $\Delta_t$ perturbs this condition and this is manifested in atomic observables and the output light quadrature $Y_{\rout}$. (e) In the ideal system (without single particle decoherence), the steady state is characterized by a perfect cancellation between the torques induced by the input field and the self-radiated field of the atoms.}\label{fig:Schematic}
\end{figure}


\subsection{Mean field analysis}

The mean field equations of motion are obtained by replacing the operators ($\hat{a},\hat{S}^-,\hat{S}_z,\hat{A}_{\rout}$) by c-numbers ($\alpha,J,Z,\alpha_{\rout}$) and omitting $\delta\hat{A}_{\rmin}(t)$ in Eq.~\eqref{eqn:Heis-LangeIdeal}. The only stable steady state solution at $\Delta_t=0$ occurs below a drive threshold, $\alpha_{\rmin}\leq \alpha_{\rmin}^c=gN/(2\sqrt{\kappa})$, and is characterized by:
\begin{equation}
\alpha=0\quad J=-\frac{i N}{2}\sin\theta\quad Z=-\frac{N}{2}\cos\theta
\end{equation}
where $\sin\theta=\alpha_{\rmin}/\alpha_{\rmin}^c$. Physically, at nonzero inversion ($Z$) the atoms radiate a field into the cavity that precisely cancels the field generated by the input light. This leads to zero intracavity photons ($\alpha=0$) and a collective Bloch vector pointing at an angle ($\theta$) on the southern hemisphere, as shown in Fig.~\ref{fig:Schematic}$\MK$(e). Using the input-output relation, we find that $\alpha_{\rout}=\alpha_{\rmin}$, i.e. all the light is reflected.

Above the threshold $\alpha_{\rmin}^c$ there are no stable steady state solutions and mean field theory predicts persistent dynamical oscillations. The relaxation towards the steady state is caused by quantum-fluctuation-induced diffusion between classical trajectories, and leads, at very long times, to a very highly mixed steady state with exactly zero $Z$~\cite{Carmichael_1980}. This quantity is continuous but non-analytic at the threshold point, reminiscent of second order phase transitions, though notions of symmetry-breaking are related to other observables instead~\cite{Hannukainen2018,Link2019}.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{Fig2_v6.pdf}
\caption{(a) Response of the system ($\hat{Y}$ quadrature of output light) to a static atom-drive detuning ($\Delta_0$) for fixed $\alpha_{\rmin}=\alpha_{\rmin}^c/\sqrt{2}$. (b) Representative spectral density of frequency fluctuations after feedback stabilization. (c) Representative response of the system to frequency fluctuations. (d) Steady state atomic inversion $Z/N$ as a function of input light intensity in the presence of single particle decoherence for three different values of effective cooperativity $NC_{\rmeff}=4,10,100$ (red, blue and purple, respectively). Solid lines are stable solutions while dashed lines are unstable. Regions of bistability are shaded.}\label{fig:IdealSystem}
\end{figure}

\subsection{Fluctuations}
We now linearize Eq.~\eqref{eqn:Heis-LangeIdeal} with respect to the mean field steady state solution. The equations are more easily written in terms of fluctuations of the intracavity field quadratures $\hat{X}=(\hat{a}+\hat{a}^\dagger)/2$, $\hat{Y}=(\hat{a}-\hat{a}^\dagger)/(2i)$, of input field quadratures $\delta\hat{X}_{\rmin}=(\delta\hat{A}_{\rmin}+\delta\hat{A}_{\rmin}^\dagger)/2$, $\delta\hat{Y}_{\rmin}=(\delta\hat{A}_{\rmin}-\delta\hat{A}_{\rmin}^\dagger)/(2i)$ and of spin fluctuations perpendicular to the mean field Bloch vector: $\hat{S}_x$ and $\hat{S}_\perp=\hat{S}_y\cos\theta+\hat{S}_z\sin\theta$. This results in the following sets of decoupled linear equations:
\begin{equation}
\partial_t
\begin{pmatrix}
\hat{S}_x\\ \hat{Y}
\end{pmatrix}=
\begin{pmatrix}
0& gN\cos\theta\\
-g & -\frac{\kappa}{2}
\end{pmatrix}
\begin{pmatrix}
\hat{S}_x\\
\hat{Y}
\end{pmatrix}+
\begin{pmatrix}
\frac{N}{2}\Delta_t \sin\theta\\ \sqrt{\kappa}\delta\hat{Y}_{\rmin}
\end{pmatrix},\quad \partial_t
\begin{pmatrix}
\hat{S}_\perp\\
\hat{X}
\end{pmatrix}=
\begin{pmatrix}
0& gN\\
-g\cos\theta & -\frac{\kappa}{2}
\end{pmatrix}
\begin{pmatrix}
\hat{S}_\perp\\
\hat{X}
\end{pmatrix}+
\begin{pmatrix}
0\\
\sqrt{\kappa}\delta\hat{X}_{\rmin}
\end{pmatrix}.
\end{equation}
The direction parallel to the mean field Bloch vector, $\hat{S}_{||}=\hat{S}_z\cos\theta-\hat{S}_y\sin\theta$, is only sensitive to fluctuations at second order due to the curvature of the Bloch sphere [$\delta\hat{S}_{||}=1/2-(\hat{S}_x^2+\hat{S}_{\perp}^2)/N$]. From the linearized equations and the input-output relation we can obtain expressions for the Fourier modes of the output field, defined as $\hat{A}_{\rout}(\omega)=\int e^{i\omega t}\hat{A}_{\rout}(t)\,dt$, with similar definitions for $\Delta(\omega),\,\delta\hat{A}_{\rmin}(\omega),\,\hat{a}(\omega)$ and $\hat{S}_{x,y,z}(\omega)$. Omitting transient effects\footnote{The expression for the output signal $\hat{Y}_{\rout}(\omega)$ is seemingly missing an independent contribution coming from atomic projection noise, but this is an incorrect assertion. Atomic projection noise in the steady state is entirely determined by dynamical noise ($\delta\hat{A}_{\rmin}$), a fact that is further established by computing the steady state value of $\hat{S}_x(\omega)$ and noting that the correct variance of $\hat{S}_x$ (in agreement with~\cite{Carmichael_1980,Drummondo1980,Walls1980}) is recovered using correlation functions of $\delta\hat{A}_{\rmin}$.}:
\begin{equation}\label{eqn:IdealFluctuations}
\begin{split}
\hat{Y}_{\rout}(\omega)&=-\frac{\big(gN\sin\theta\sqrt{\kappa}/2\big)}{\big(i\omega-\lambda_+\big)\big(i\omega-\lambda_-\big)}\Delta(\omega)+\frac{\big(i\omega+\lambda_+\big)\big(i\omega+\lambda_-\big)}{\big(i\omega-\lambda_+\big)\big(i\omega-\lambda_-\big)}\delta\hat{Y}_{\rmin}(\omega),
\end{split}
\end{equation}
where $\lambda_{\pm}=(\kappa/4)\pm\sqrt{(\kappa/4)^2-g^2N\cos\theta}$. The output light quadrature $\hat{Y}_{\rout}$ has two contributions: one proportional to $\Delta(\omega)=\int e^{i\omega t}\Delta_t\,dt$ (signal) and another one proportional to $\delta\hat{Y}_{\rmin}(\omega)$ (noise).

The real parts of $\lambda_\pm$ are the decay constants of the system. The smallest one determines how long it takes the system to relax to the steady state. When $\alpha_{\rmin}\to\alpha_{\rmin}^c$ ($\theta\to\pi/2$), this timescale increases without bounds, rendering the steady state condition harder to achieve. In the limit that $g\sqrt{N}\ll \kappa$, relevant for narrow linewidth transitions, these constants become $\lambda_+\approx \kappa/2$ and $\lambda_-\approx 2g^2N\cos\theta/\kappa$ (thus $\lambda_+\gg \lambda_-$). The combination $g^2/\kappa$ appears often and is typically rewritten as $4g^2/\kappa\equiv C\gamma$, where $\gamma$ is the natural linewidth of the transition and $C=4g^2/(\kappa\gamma)$ is the cooperativity, which depends only on cavity geometry.

Equations~\eqref{eqn:IdealFluctuations} relate the fluctuations of quantum observables in the steady state to the fluctuations of the input field ($\delta\hat{X}_{\rmin},\,\delta\hat{Y}_{\rmin}$). For example, fluctuations of $\hat{Y}_{\rout}$ define $S_{Y_{\rout}}(\omega)$, the noise spectral density of the output light, according to:
\begin{equation}
\left\langle \delta\hat{Y}_{\rout}(\omega)\delta\hat{Y}_{\rout}(\omega')\right\rangle\equiv \frac{2\pi S_{Y_{\rout}}(\omega)\,\delta\left(\omega+\omega'\right)}{4}.
\end{equation}
Using $\langle\delta\hat{Y}_{\rmin}(\omega)\delta\hat{Y}_{\rmin}(\omega')\rangle=2\pi\delta(\omega+\omega')/4$ and $\langle\delta\hat{Y}_{\rmin}(\omega)\rangle=0$, it follows that $S_{Y_{\rout}}(\omega)=1$, which is the value associated to shot noise. A similar result holds for $S_{X_{\rout}}(\omega)$, which indicates that the output light is in a coherent state\cite{Somech2022}.

\subsection{Frequency estimation}
A glance at Eq.~\eqref{eqn:IdealFluctuations} indicates that the steady state output light quadrature $\hat{Y}_{\rout}$ carries information about the drive detuning $\Delta_t$. This suggests that this system can be used as a passive frequency reference, similar in spirit to Ramsey spectroscopy, where an observable carrying information about $\Delta_t$ is measured after an operation time $T$. In the presence of a static detuning $\Delta_t=\Delta_0$ [$\Delta(\omega)=2\pi\Delta_0\delta(\omega)$], the most straightforward observable that serves as an estimator of $\Delta_0$ is constructed from the average of $\hat{Y}_{\rout}$ over time. Within the linear regime:
\begin{equation}
\hat{Y}_{\rmest}\equiv\frac{1}{T}\int_0^T\hat{Y}_{\rout}(t)\,dt\longrightarrow \left\langle\hat{Y}_{\rmest}\right\rangle=\frac{\Delta_0}{2g\cot\theta/\sqrt{\kappa}}.
\end{equation}
This linear relation is valid for $\Delta_0\lesssim NC\gamma$ [condition obtained by solving the mean field equations for finite static $\Delta_t=\Delta_0$, see Fig.~\ref{fig:IdealSystem}$\MK$(a)], beyond which a typical dispersive response is observed. A real measurement of $\hat{Y}_{\rmest}$ over a finite time is subject to quantum noise coming from the fluctuations of the output light, which are related to $\delta\hat{Y}_{\rmin}(\omega)$ in Eq.~\eqref{eqn:IdealFluctuations}. The estimator variance after an integration time $T$ is calculated using Eq.~\eqref{eqn:IdealFluctuations}:
\begin{equation}
\begin{split}
\left\langle\delta\hat{Y}_{\rmest}^2\right\rangle\equiv \left\langle\left(\hat{Y}_{\rmest}-\left\langle\hat{Y}_{\rmest}\right\rangle\right)^2\right\rangle=\frac{S_{\rmY_{\rout}}(0)}{4T}=\frac{1}{4T},
\end{split}
\end{equation}
and from this we can compute the sensitivity:
\begin{equation}
\delta\Delta_0^2=\frac{\left\langle\delta\hat{Y}_{\rmest}^2\right\rangle}{\left(\partial_{\Delta_0}\left\langle\hat{Y}_{\rmest}\right\rangle\right)^2}=\frac{g^2(\cot\theta)^2}{\kappa T}=\frac{C\gamma (\cot\theta)^2}{4 T}.
\end{equation}
As an example, consider the $^1S_0\to ^3\hspace{-0.1cm}P_0$ transition in $^{87}$Sr ($\lambda_a=698$ nm) after one second integration time and operating at $\theta=45^\circ$. Using the cavity parameters in~\cite{norcia2018cavity,Norcia2016c}, this corresponds to a relative frequency resolution of $\delta\Delta_0/\omega_a \approx 10^{-17}$. Experimentally, such a measurement of $\hat{Y}_{\rout}$ is achieved by means of a homodyne detection setup where the input light field acts as the local oscillator.


\subsection{Stabilization via feedback}\label{sec:feed}
The light coming outside of the cavity can be used within a closed feedback loop to provide automatic frequency stabilization. As shown in the previous section, we can estimate the detuning between the laser drive and the atomic transition, and we can then use this information to correct the frequency of the laser drive. However, the output light is a fluctuating quantity and hence the feedback process introduces noise into the corrected frequency. The objective in this section is to estimate the size of this noise, following the discussion in Ref.~\cite{Riehle2004}.

We begin with the expression for the output quadrature $\hat{Y}_{\rout}(\omega)$ given in Eq.~\eqref{eqn:IdealFluctuations}:
\begin{equation}\label{eqn:YtoR}
\hat{Y}_{\rout}(\omega)=\underbrace{-\frac{\big(gN\sin\theta\sqrt{\kappa}/2\big)}{\big(i\omega-\lambda_+\big)\big(i\omega-\lambda_-\big)}}_{R(\omega)}\hat{\Delta}(\omega)+\underbrace{\frac{\big(i\omega+\lambda_+\big)\big(i\omega+\lambda_-\big)}{\big(i\omega-\lambda_+\big)\big(i\omega-\lambda_-\big)}\delta\hat{Y}_{\rmin}(\omega)}_{\hat{N}(\omega)},
\end{equation}
where $R(\omega)$ is the response of the signal $\hat{Y}_{\rout}(\omega)$ to the perturbation $\hat{\Delta}(\omega)$, and $\hat{N}(\omega)$ describes the quantum fluctuations of $\hat{Y}_{\rout}(\omega)$. Note also that we are now writing $\hat{\Delta}(\omega)$ as an operator because we are considering the steady state of the feedback process, where the input $\hat{\Delta}$ already includes the effects of quantum noise.

We assume the feedback signal is obtained from the output light using a linear filter $\beta(\omega)$, which gives us a second equation $\hat{\Delta}_{\rmfeed}(\omega)=\beta(\omega)\hat{Y}_{\rout}(\omega)$. Finally, we have a third relation connecting the bare detuning in the absence of feedback [$\Delta_{\rmbare}(\omega)$] to the corrected detuning: $\hat{\Delta}(\omega)=\Delta_{\rmbare}(\omega)-\hat{\Delta}_{\rmfeed}(\omega)$. From these three equations we can obtain a direct relation between $\hat{\Delta}(\omega)$ and $\Delta_{\rmbare}(\omega)$:
\begin{equation}
\hat{\Delta}(\omega)=\frac{1}{1+R(\omega)\beta(\omega)}\Delta_{\rmbare}(\omega)-\frac{\beta(\omega)}{1+R(\omega)\beta(\omega)}\hat{N}(\omega)\approx \frac{1}{R(\omega)\beta(\omega)}\Delta_{\rmbare}(\omega)-\frac{\hat{N}(\omega)}{R(\omega)},
\end{equation}
where in the last equality we have assumed $|\beta(\omega)R(\omega)|\gg 1$ for the frequencies of interest (i.e. below the unity gain frequency of the feedback loop). If the original laser drive is already noisy, then $\Delta_{\rmbare}(\omega)$ becomes stochastic. However, the feedback process suppresses this noise by the large factor $R(\omega)\beta(\omega)$, while at the same time introducing the extra noise $\delta\hat{\Delta}(\omega)=-\hat{N}(\omega)/R(\omega)$:
\begin{equation}
\delta\hat{\Delta}(\omega)\approx\frac{\big(i\omega+\lambda_+\big)\big(i\omega+\lambda_-\big)}{gN\sqrt{\kappa}\sin\theta/2}\hat{Y}_{\rmin}(\omega)
\end{equation}
with second order average:
\begin{equation}
\left\langle\delta\hat{\Delta}(\omega)\delta\hat{\Delta}(\omega')\right\rangle=2\pi\delta(\omega+\omega')\left[\frac{ \big(\omega^2+\lambda_+^2\big)\big(\omega^2+\lambda_-^2\big)}{g^2N^2\kappa\sin\theta^2}\right].
\end{equation}
The coefficient in front of $2\pi\delta(\omega+\omega')$ defines the spectral density of frequency fluctuations $S_{\Delta}(\omega)$~\cite{Riehle2004,Martin2011}. In the limit $\kappa\gg g\sqrt{N}$, this will be of the form:
\begin{equation}\label{eqn:idealLW}
S_{\Delta}(\omega)=2\pi\Delta\hspace{-0.05cm}f \left[1+\left(\frac{\omega}{\kappa/2}\right)^2\right]\left[1+\left(\frac{\omega}{\omega_S}\right)^2\right],
\end{equation}
with $\omega_S\ll \kappa/2$ [see Fig.~\ref{fig:IdealSystem}$\MK$(b)]. The value of $S_{\Delta}(\omega)$ at $\omega=0$ determines the effective linewidth $\Delta\hspace{-0.05cm}f$. In this ideal scenario, $\Delta\hspace{-0.05cm}f=C\gamma \cot\theta^2/4$, which is the same scaling as the one that appears in studies of the superradiant laser\cite{Meiser2009}. As $\theta\to\pi/2$ ($\alpha_{\rmin}\to\alpha_{\rmin}^c$), $\Delta \hspace{-0.05cm}f$ approaches 0, and thus it is ultimately limited by finite size effects. The corner frequency $\omega_S=NC\gamma \cos\theta/2$ marks the point at which $S_{\Delta}(\omega)$ switches from constant to $\sim\omega^2$ behaviour, which in a real laser are characteristic of the Schawlow-Townes noise floor and of shot noise (due to a finite laser output power) respectively.

The construction of the feedback loop [i.e. the specific choice of $\beta(\omega)$] also requires information on the response function $|R(\omega)|$ to guarantee large low-frequency gain [$R(0)\beta(0)$] and keep the loop stable. In this case, $|R(\omega)|$ behaves as [see Eq.~\eqref{eqn:YtoR}]
\begin{equation}\label{eqn:Response}
|R(\omega)|\propto\frac{1}{\big[1+(2\omega/\kappa)^2\big]^{1/2}\big[1+(\omega/\omega_R)^2\big]^{1/2}}
\end{equation}
and is schematically depicted in Fig.~\ref{fig:IdealSystem}$\MK$(c). In this ideal case, $\omega_R$ is also $NC\gamma\cos\theta/2$, though this is not generic. Furthermore, a realistic calculation of $ \Delta \hspace{-0.05cm}f$, $\omega_S$ and $\omega_R$ requires that single particle sources of decoherence be included in the analysis.


\section{Non ideal model}
In any realistic system there is dephasing and incoherent decay, and these processes will modify both the steady state of the system and its linear response to perturbations. This is encoded in the Heisenberg--Langevin equations that replace Eq.~\eqref{eqn:Heis-LangeIdeal}:
\begin{equation}\label{eqn:Heis-LangeNonIdeal}
\begin{split}
\partial_t\hat{a}&=-\frac{\kappa}{2}\hat{a}-ig\hat{S}^-+\sqrt{\kappa}\left[\alpha_{\rmin}+\delta\hat{A}_{\rmin}(t)\right]\\
\partial_t\hat{S}^-&=i\Delta_t\hat{S}^-+2ig\hat{S}_z\hat{a}-\left(\frac{\gamma+\gamma_d+\gamma_p}{2}\right)\hat{S}^-+\sqrt{N\gamma}\hat{F}_{\gamma}^-(t)+\sqrt{N\gamma_{p}}\hat{F}_{\gamma_p}^-(t)+\sqrt{N\gamma _d}\hat{F}_{\gamma_d}^-(t)\\
\partial_t\hat{S}_z&=ig(\hat{a}^{\dagger}\hat{S}^--\hat{a}\hat{S}^{\dagger})-(\gamma+\gamma_p)\left(\hat{S}_z+\frac{N}{2}\right)+\sqrt{N\gamma }\hat{F}_{\gamma}^z(t)+\sqrt{N\gamma_p}\hat{F}_{\gamma_p}^z(t),
\end{split}
\end{equation}
where $\gamma_d$ is the dephasing rate, $\gamma$ is the natural spontaneous emission rate of the transition, $\gamma_p$ is the rate of the depumping process depicted in Fig.~\ref{fig:Schematic}$\MK$(b), $\hat{F}_{\gamma}^-$ and $\hat{F}_{\gamma}^z$ are noise operators associated to spontaneous emission (notice $\hat{F}_{\gamma}^z$ is hermitian), and $\hat{F}_{\gamma_d}^-$ is a noise operator associated to dephasing. We have also included an additional depumping process, $\gamma_p$, with associated noise operators $\hat{F}_{\gamma_p}^-$ and $\hat{F}_{\gamma_p}^z$, and which can be engineered via coherent driving onto a rapidly decaying auxiliary third level [see Fig.~\ref{fig:Schematic}$\MK$(b)]. In principle the noise operators are defined in terms of the degrees of freedom which cause the associated decay processes, but they are in general difficult to access, unlike $\hat{A}_{\rout}(t)$. In practice, we are interested mostly in their effect on the system, which is encoded in their lowest order correlators. The only nonzero ones (up to hermitian conjugation) are\cite{Tieri2017,Meystre2007}:
\begin{equation}
\begin{matrix*}[l]
\left\langle \hat{F}_{\gamma}^-(t)\hat{F}_\gamma^+(t')\right\rangle&=\left\langle \hat{F}^-_{\gamma_p}(t)\hat{F}_{\gamma_p}^+(t')\right\rangle&=\delta(t-t')\\[2pt]
\left\langle \hat{F}_\gamma^-(t)\hat{F}_{\gamma}^z(t')\right\rangle&=\left\langle \hat{F}_{\gamma_p}^-(t)\hat{F}_{\gamma_p}^z(t')\right\rangle&=\frac{\left\langle\hat{S}^-\right\rangle}{N}\delta(t-t')\\[2pt]
\left\langle \hat{F}_{\gamma}^z(t)\hat{F}_{\gamma}^z(t')\right\rangle&=\langle \hat{F}_{\gamma_p}^z(t)\hat{F}_{\gamma_p}^z(t')\rangle&=\left(\frac{\langle\hat{S}_z\rangle}{N}+\frac{1}{2}\right)\delta(t-t')
\\[4pt]
\left\langle \hat{F}^-_{\gamma_d}(t)\hat{F}_{\gamma_d}^+(t')\right\rangle&=\left(\frac{1}{2}-\frac{\langle\hat{S}_z\rangle}{N}\right)\delta(t-t')\\[4pt]
\left\langle \hat{F}^+_{\gamma_d}(t)\hat{F}_{\gamma_d}^-(t')\right\rangle&=\left(\frac{1}{2}+\frac{\left\langle\hat{S}_z\right\rangle}{N}\right)\delta(t-t').
\end{matrix*}
\end{equation}
Similarly to the ideal case, this system hosts resonant cooperative fluorescence\cite{BONIFACIO1976}. As in the previous section, we first investigate the system within the mean field approximation, then compute the fluctuations of the output light with respect to the steady state and analyze how these incoherent processes modify our results for frequency estimation and feedback.


\subsection{Mean field analysis}
Once again, we replace ($\hat{a},\hat{S}^-,\hat{S}_z,\hat{A}_{\rout}$) by c-numbers ($\alpha,J,Z,\alpha_{\rout}$) and omit all noise operators. The resulting equations are nonlinear but their steady state solution can be parameterized in terms of $z=Z/N$ as:
\begin{equation}\label{eqn:NonIdealSS}
\alpha=\frac{2\alpha_{\rmin}}{\sqrt{\kappa}\left(1-2NC_{\rmeff}z\right)}\quad\frac{J}{N}=i\sqrt{\frac{\gamma+\gamma_p}{\Gamma}}\left(\frac{\alpha_\rmin}{\sqrt{16 I_0}}\right)\left(\frac{2NC_{\rmeff}z}{1-2NC_{\rmeff}z}\right)\quad \alpha_{\rout}=\alpha_{\rmin}\left(\frac{2zNC_{\rmeff}+1}{2zNC_{\rmeff}-1}\right),
\end{equation}
where $C_{\rmeff}\equiv4g^2/(\Gamma\kappa)=C \gamma/\Gamma$ is the effective cavity cooperativity, $\Gamma=\gamma+\gamma_d+\gamma_p$ is the total dephasing rate, $I_0=N^2C_{\rmeff}(\gamma+\gamma_p)/16$ and $z$ satisfies:
\begin{equation}\label{eqn:NonIdealZ}
\frac{\alpha_{\rmin}^2}{8I_0}=-\frac{(z+1/2)}{z}\left(z-\frac{1}{2NC_{\rmeff}}\right)^2.
\end{equation}
The natural scale for $\alpha_{\rmin}^2$ is set by $I_0$ and is connected to the ideal threshold field ($\alpha_{\rmin}^c$) by $I_0=(\alpha_{\rmin}^c)^2(\gamma+\gamma_p)/\Gamma$. The fate of the transition present in the ideal model depends on the value of $NC_{\rmeff}$. When $NC_{\rmeff}<8$, there is only one real solution to Eq.~\eqref{eqn:NonIdealZ} and the mean field observables are smooth functions of the input field strength $\alpha_{\rmin}$ [see Fig.~\ref{fig:IdealSystem}$\MK$(d) for $NC_{\rmeff}=4$]. Notably, there is one ``dark point'' where $\alpha_{\rout}=0$, i.e. no light comes out of the cavity [see Eq.~\eqref{eqn:NonIdealSS} with $z=-1/(2NC_{\rmeff})$]. In the context of cooperative resonance fluorescence, as discussed before in Ref.~\cite{Carmichael_1980,BONIFACIO1976,Walls1980}, the smooth behavior observed in this regime is a manifestation of the disappearance of the phase transition observed in the absence of decoherence.

When $NC_{\rmeff}>8$, there are input field amplitudes $\alpha_{\rmin}$ for which there are three real solutions for $z$. Two of them correspond to stable steady states and one is in an unstable branch [see Fig.~\ref{fig:IdealSystem}$\MK$(d) for $NC_{\rmeff}=10,100$]\cite{BONIFACIO1976}. This leads to bistability when:
\begin{equation}
\frac{1}{16}\left(1-\sqrt{1-\frac{8}{NC_{\rmeff}}}\right)\left(3+\sqrt{1-\frac{8}{NC_{\rmeff}}}\right)^3<\frac{2\alpha_{\rmin}^2}{I_0}<\frac{1}{16}\left(1+\sqrt{1-\frac{8}{NC_{\rmeff}}}\right)\left(3-\sqrt{1-\frac{8}{NC_{\rmeff}}}\right)^3
\end{equation}
In the limit $NC_{\rmeff}\gg 1$, the lower and upper bounds become $16/(NC_{\rmeff})$ and $1$~\cite{Tucker2020}, respectively. This bistable behaviour is accompanied by large discontinuous jumps between stable solutions and thus the presence of hysteresis~\cite{BONIFACIO1976} instead of the smooth second order transition observed in the ideal case~\cite{Carmichael_1980}. It is preferable to work in non-bistable regimes to avoid the large fluctuations associated with these switches.

\subsection{Fluctuations}
We linearize Eq.\eqref{eqn:Heis-LangeNonIdeal} about any of the mean field steady state solutions and look at $\hat{Y}$, $\hat{Y}_{\rout}$ and $\hat{S}_x$, which carry information about $\Delta_t$. They satisfy:
\begin{equation}
\partial_t
\begin{pmatrix}
\hat{S}_x\\
\hat{Y}
\end{pmatrix}=
\begin{pmatrix}
-\frac{\Gamma}{2}& -2gNz\\
-g & -\frac{\kappa}{2}
\end{pmatrix}
\begin{pmatrix}
\hat{S}_x\\
\hat{Y}
\end{pmatrix}+
\begin{pmatrix}
(iJ)\Delta_t+\sqrt{N\gamma} \hat{F}_{\gamma}^{x}(t)+\sqrt{N\gamma_p}\hat{F}_{\gamma_p}^{x}(t)+\sqrt{N\gamma_d}\hat{F}_{\gamma_d}^{x}(t)\\ \sqrt{\kappa}\delta\hat{Y}_{\rmin}
\end{pmatrix},
\end{equation}
where $\hat{F}_{\gamma}^x=(\hat{F}_{\gamma}^-+\hat{F}_{\gamma}^+)/2$, $\hat{F}_{\gamma_p}^x=(\hat{F}_{\gamma_p}^-+\hat{F}_{\gamma_p}^+)/2$ and $\hat{F}_{\gamma_d}^x=(\hat{F}_{\gamma_d}^-+\hat{F}_{\gamma_d}^+)/2$ are the corresponding quadrature noise operators. We can then calculate the behaviour of $\hat{Y}_{\rout}$ through $\hat{Y}_{\rout}=\delta\hat{Y}_{\rmin}-\sqrt{\kappa}\hat{Y}$. In frequency space, omitting transient behaviour:
\begin{multline}\label{eqn:NonIdealYout}
\hat{Y}_{\rout}(\omega)=-\frac{2g\sqrt{\kappa}(iJ)}{\big(i\omega-l_+\big)\big(i\omega-l_-\big)}\Delta(\omega)+\frac{\big(i\omega+m_+\big)\big(i\omega+m_-)}{\big(i\omega-l_+\big)\big(i\omega-l_-\big)}\delta\hat{Y}_{\rmin}(\omega)\\[3pt] +\frac{\big(\kappa/2\big)\sqrt{NC_{\rmeff}\Gamma}}{\big(i\omega-l_+\big)\big(i\omega-l_-\big)}\left[\sqrt{\gamma}\hat{F}_{\gamma}^x(\omega)+\sqrt{\gamma_p}\hat{F}_{\gamma_p}^x(\omega)+\sqrt{\gamma_d}\hat{F}_{\gamma_d}^x(\omega)\right],
\end{multline}
where $l_{\pm}=\frac{\kappa+\Gamma}{4}\pm\sqrt{\big(\frac{\kappa-\Gamma}{4}\big)^2+2g^2Nz}$ and $m_{\pm}=\frac{\kappa-\Gamma}{4}\pm\sqrt{(\frac{\kappa+\Gamma}{4})^2+2g^2Nz}$. The smallest decay constant, $l_-$, is modified due to the presence of the extra dissipation sources and is always finite. In the limit $g\sqrt{N},\Gamma\ll \kappa$, these constants reduce to $l_+\approx m_+\approx \kappa/2$, $l_-\approx \Gamma(1-2NC_{\rmeff}z)/2$ and $m_-\approx -\Gamma(1+2NC_{\rmeff}z)/2$.

The noise spectral density of $\hat{Y}_{\rout}$ at $\omega=0$ can be calculated from Eq.~\eqref{eqn:NonIdealYout}:
\begin{equation}\label{eqn:NonIdealEstimatorVariance}
\begin{split}
S_{Y_{\rout}}(0)=\frac{\left(1+2NC_{\rmeff}z\right)^2}{(1-2N C_{\rmeff}z)^2}+\frac{4NC_{\rmeff}}{\left(1-2NC_{\rmeff}z\right)^2}
\end{split}
\end{equation}
and is shown in Fig.~\ref{fig:NonIdealNoise}$\MK$(a) as a function of $\alpha_{\rmin}^2/I_0$. The smallest values of noise improve with increasing $NC_{\rmeff}$, but fall within bistable regions. Outside those regions, reducing $NC_{\rmeff}$ brings the noise closer to shot noise levels.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{figures/Fig4_Rescaled.pdf}
\caption{(a) Power spectral density of output quadrature $\hat{Y}_{\rout}$ as a function of input light intensity $\alpha_{\rmin}^2$ for $NC_{\rmeff}=0.5,4,10,100$ (green, red, blue and purple, respectively). Solid lines are stable solutions and dashed lines are unstable. Bistable regions are shaded. (b) Effective linewidth, measured in units of $\Delta\hspace{-0.05cm}f_0=(4\pi)^{-1}(C\gamma/4)[\Gamma/(\gamma_p+\gamma)]$. The ``dark point'' for each $NC_{\rmeff}$ is marked by the black circle (nonexistent for the green curve).}\label{fig:NonIdealNoise}
\end{figure}


\subsection{Stabilization via feedback}
Applying the formulas of Section~\ref{sec:feed} leads to the following effective linewidth:
\begin{equation}\label{eqn:NonIdealLinewidth}
2\pi\Delta \hspace{-0.05cm}f=\left(\frac{C\gamma}{4}\right)\left(\frac{\gamma+\gamma_p+\gamma_d}{\gamma+\gamma_p}\right)\left[\frac{(1+2NC_{\rmeff}z)^2+4NC_{\rmeff}}{2\left(N C_{\rmeff}+2NC_{\rmeff}z\right)\left(-2NC_{\rmeff}z\right)}\right],
\end{equation}
which is parameterized by $z$ and hence is an implicit function of the input field intensity $\alpha_{\rmin}^2$ (since $z$ depends on $\alpha_{\rmin}$). We have expressed Eq.~\eqref{eqn:NonIdealLinewidth} as a product of three different factors to make the discussion clearer. The first factor is the scaling present both in the ideal section and in the superradiant laser. The second factor $\Gamma/(\gamma+\gamma_p)$ strongly penalizes the effective linewidth when $\gamma_d\gg \gamma+\gamma_p$, which is a consequence of a reduction in the size of the radiating dipole moment $J\propto \sqrt{(\gamma+\gamma_p)/\Gamma}$ [see Eq.~\eqref{eqn:NonIdealSS}]. Under these conditions, an adjustable depumping rate provides an advantage when the natural $\gamma$ is small, as is the case in narrow linewidth transitions. If $\gamma_p\gg \gamma_d$, then $(\gamma+\gamma_p)/\Gamma \approx 1$ and there is no penalization. The third factor includes all the effects of operating at different $\alpha_{\rmin}$ and is difficult to estimate a priori. It can be of order $1$ or much bigger than $1$ depending on $\alpha_{\rmin}$, as shown in the next subsection.

In Fig.~\ref{fig:NonIdealNoise}$\MK$(b) we show $\Delta\hspace{-0.05cm}f$ as a function of input light intensity for different values of $NC_{\rmeff}$, normalized to $\Delta\hspace{-0.05cm}f_0=(4\pi)^{-1}(C\gamma/4)[\Gamma/(\gamma+\gamma_p)]$. When $NC_{\rmeff}>8$, the $\alpha_{\rmin}$ that minimizes $\Delta\hspace{-0.05cm}f$ falls within regimes of bistability. Outside of these regions, we find that it is preferable to work at more moderate values of $NC_{\rmeff}$ for fixed $\alpha_{\rmin}^2/I_0$. Reducing $NC_{\rmeff}$ too much is not convenient either because the size of the signal acquired through the output light is also reduced, negatively impacting $\Delta\hspace{-0.05cm}f$ [see Fig.~\ref{fig:NonIdealNoise}$\MK$(b) for $NC_{\rmeff}=0.5$].


\subsection{Operating regimes}
We describe here two possible operating regimes corresponding to different values of $\alpha_{\rmin}$, discuss the achievable effective linewidths $\Delta \hspace{-0.05cm}f$, and calculate the corner frequencies $\omega_S$ and $\omega_R$.


\subsubsection{Strong input field}
Here we take $\alpha_{\rmin}^2\gg I_0=N^2C_{\rmeff}(\gamma_p+\gamma)/16$, a situation previously analyzed in Ref.~\cite{Martin2011}, and which leads to $z\sim 1/\alpha_{\rmin}^2\to 0$, $\alpha_{\rout}\approx -\alpha_{\rmin}$, and $J$ (dipole moment) $\sim 1/\alpha_{\rmin}$. Since $J$ shrinks with larger $\alpha_{\rmin}$, the size of the signal necessary to estimate frequency fluctuations is reduced, while the noise in the output light does not change significantly [see Eq.~\eqref{eqn:NonIdealEstimatorVariance} at $z=0$]. The effective linewidth is:
\begin{equation}
2\pi\Delta \hspace{-0.05cm}f=\frac{C\gamma}{4}\left(\frac{\Gamma}{\gamma_p+\gamma}\right)\left(\frac{\alpha_{\rmin}^2}{4I_0}\right)\left(1+4NC_{\rmeff}\right)
\end{equation}
In this case, apart from $\Gamma/(\gamma+\gamma_p)$, there is an extra factor proportional to the input intensity that is much larger than one by assumption. This is the linewidth found in Ref.~\cite{Martin2011}, except for the factor $(1+4NC_{\rmeff})$ which includes the additional noise due to the single particle decoherence processes. Even when $NC_{\rmeff}>8$, the strong input field region is free from bistability, which is advantageous if manipulating $NC_{\rmeff}$ is not possible or easy. On the other hand, the achievable effective linewidth is then also penalized from the extra $(1+4NC_{\rmeff})$ factor.

In contrast to the ideal case, $\omega_S$ and $\omega_R$ are now different. In the limit $g\sqrt{N},\Gamma\ll \kappa$ they are $\omega_S=\Gamma\sqrt{NC_{\rmeff}+1/4}$ and $\omega_R=\Gamma/2$. It is desirable to have a large $\omega_S$, either by increasing $NC_{\rmeff}$, at the expense of a larger linewidth $\Delta\hspace{-0.05cm}f$, or by modifying the depumping rate $\gamma_p$. An important final piece of information is an estimate of the size of detunings $\Delta_t$ under which the linear approximation is valid. Solving the mean field equations in the presence of a static detuning $\Delta_0$ leads to a dispersive curve similar to Fig.~\ref{fig:IdealSystem}$\MK$(a), but where the maximal value is obtained at $\Delta_0\approx NC_{\rmeff}\Gamma\alpha_{\rmin}/\sqrt{8I_0}=NC\gamma\alpha_{\rmin}/\sqrt{8I_0}$ instead of $\Delta_0\approx NC\gamma$ (ideal case). Linearity is thus valid for a larger range of $\Delta_t$, as compared to the ideal case, at the expense of a smaller slope and hence lower sensitivity.

\subsubsection{Dark point}
This corresponds to no light coming outside the cavity, namely $\alpha_{\rout}=0$. At this operating point, the effective linewidth is:
\begin{equation}\label{eqn:darkpointlw}
2\pi\Delta \hspace{-0.05cm}f=\frac{C\gamma}{2}\left(\frac{\Gamma}{\gamma+\gamma_p}\right) \left(\frac{NC_{\rmeff}}{NC_{\rmeff}-1}\right).
\end{equation}
and only exists when $NC_{\rmeff}>1$. In this case, apart from $\Gamma/\gamma$, there is a third factor which is of order one when $NC_{\rmeff}\gtrsim 2$. Hence, the achievable linewidth can be of the same size as that of the superradiant laser ($\sim C\gamma$)~\cite{Meiser2009}, provided $\gamma+\gamma_p$ and $\Gamma$ are comparable. Consider again the $^1S_0\to\,^3P_0$ transition in $^{87}\Sr$ in the cavity of Ref.~\cite{norcia2018cavity}. Using $g=2\pi\times 4$ Hz, $\kappa=2\pi\times 160$ kHz, $\gamma_d=2\pi\times 3$ Hz, an artificially augmented spontaneous emission rate $\gamma_p+\gamma=\gamma_d$ and $N=10^5$ atoms leads to $\Delta\hspace{-0.05cm}f=0.5$ mHz at $NC_{\rmeff}\approx 6.6$.
%to $\Delta\nu=0.5$

Within the approximations $g\sqrt{N},\Gamma\ll\kappa$, the corner frequencies are now $\omega_S=\sqrt{NC_{\rmeff}\Gamma^2}=\sqrt{NC\gamma \Gamma}$ and $\omega_R=\Gamma$. Increasing $\omega_S$ while avoiding bistability ($NC_{\rmeff}<8$) requires increasing the depumping rate $\gamma_p$. Furthermore $\gamma_p$ also enhances $\omega_R$, which could be as small as $\gamma$ when $\gamma_d=\gamma_p=0$, and this relaxes technical constraints on the feedback loop gain. Note also that if we apply a depumping strong enough that $\Gamma\sim NC\gamma/2$, then $\omega_S \sim NC\gamma$, similar to the result for the superradiant laser at optimal pumping rate. Finally, the response ($\hat{Y}_{\rout}$) to large static detunings $\Delta_0$ is still dispersive, as in Fig.~\ref{fig:IdealSystem}$\MK$(a), with the maximal value attained at $\Delta_0\approx NC_{\rmeff}\Gamma=NC\gamma$ (similar to ideal case).

This ``dark point'' is representative of the properties of the system close to the minimum $\Delta\hspace{-0.05cm}f$ [see Fig.~\ref{fig:NonIdealNoise}$\MK$(b)], and could be easily identified by looking at the reflected intensity. Numerically, we find that Eq.~\eqref{eqn:darkpointlw} is at most within a factor of $2$ of the optimal $\Delta\hspace{-0.05cm}f$ when $2<NC_{\rmeff}<8$.

\section{Conclusions and outlook}
We have calculated the effective linewidth of a laser stabilized to an ensemble of atoms inside a QED cavity across various values of input light intensity and atom-light interaction strength, measured in terms of the cooperativity parameter $NC_{\rmeff}$. We find that working with moderate values of $NC_{\rmeff}\in [2,8]$ is preferable and can lead to effective linewidths comparable to those of the superradiant laser, provided the total dephasing and effective spontaneous emission rates are of the same size. These predictions open a path for a complementary implementation of ultranarrow linewidth light which does not require a strong incoherent pump but instead coherent drive and feedback and thus a more versatile and amenable avenue for the experimental generation of continuous coherent light with long lived dipoles. Furthermore, concerns about the finite lifetime of atoms in the cavity can be addressed by continuous atom loading schemes~\cite{Cline2022}.

The achievable linewidths using this configuration in state-of-the-art cavities with alkaline-earth atoms are in the mHz range. Further improvements would require reducing the ratio $g^2/\kappa$ while keeping $NC_{\rmeff}$ fixed, which could be realized using a larger number of atoms.


\begin{table}[!hp]
\caption{List of parameters I}
%\vspace*{-10pt}
\centering
\hspace*{-25pt}
\begin{tabular}{ | >{\centering\arraybackslash}m{1em} m{14em} |>{\centering\arraybackslash}m{9.5em} m{15.5em}|}
\hline $2g$ & Single photon Rabi frequency & $\omega_a$, $\omega_c$, $\omega_d$ & Atomic, cavity, drive frequencies\\

$\kappa$ & Cavity power decay rate & $\Delta_t=\omega_d-\omega_a$ & Drive-atom detuning\\

$\gamma$&Spontaneous emission rate & $\Gamma=\gamma+\gamma_p+\gamma_d$& Total dephasing rate\\

$\gamma_p$& Depumping rate &$C=4g^2/(\kappa\gamma)$& Cavity cooperativity\\

$\gamma_d$& Dephasing rate &$C_{\rmeff}=C\gamma/\Gamma$ & Effective cooperativity\\

$N$ & Atom number& $\alpha_{\rmin}^c=gN/(2\sqrt{\kappa})$ & Ideal critical input field\\

$\Delta\hspace{-0.05cm}f$ & Effective linewidth &$I_0=(\alpha_{\rmin}^c)^2(\gamma+\gamma_p)/\Gamma$& Input field scale\\
\hline
\end{tabular}
%\vspace*{25pt}
\end{table}
%\vspace*{15pt}

\begin{table}[!hp]
\caption{List of parameters and variables II ($O$ stands for any variable)}
\centering
\vspace*{-10pt}
\hspace*{-23pt}
\begin{tabular}{ | >{\centering\arraybackslash}m{10em} >{\centering\arraybackslash}m{10em} |>{\centering\arraybackslash}m{10em} >{\centering\arraybackslash}m{13em}|}
\hline $\hat{a}=\hat{X}+i\hat{Y}$ & Intracavity field & $\alpha$ & Mean field intracavity field\\

$\hat{S}_z$, $\hat{S}^-=\hat{S}_x-i\hat{S}_y$ & Atomic variables &$Z=Nz, J$ & Mean field atomic variables\\

$\hat{A}_{\rmin}=\hat{X}_{\rmin}+i\hat{Y}_{\rmin}$& Input field& $\alpha_{\rmin}$ & Mean field input field\\

$\hat{A}_{\rout}=\hat{X}_{\rout}+i\hat{Y}_{\rout}$& Output field& $\alpha_{\rout}$ & Mean field output field\\
\hline
\vspace{0.25em} $\hat{F}_{\gamma}^z$, $\hat{F}_{\gamma}^-=\hat{F}_{\gamma}^x-i\hat{F}_{\gamma}^y$& Noise due to $\gamma$& $\hat{F}_{\gamma_{d}}^-=\hat{F}_{\gamma_{d}}^x-i\hat{F}_{\gamma_{d}}^y$& Noise due to $\gamma_d$ \\
$\hat{F}_{\gamma_p}^z, \hat{F}_{\gamma_p}^-=\hat{F}_{\gamma_p}^x-i\hat{F}_{\gamma_p}^y$& Noise due to $\gamma_p$ &$\hat{S}_O(\omega)$& Noise spectral density of $O$\\

\hline
\vspace{0.25em} $R(\omega)$& \hspace{-3em}Response function Eq.~\eqref{eqn:Response}&$O(\omega)=\int e^{i\omega t} O(t)\,dt$ & Fourier transform of $O$ \\
$\omega_S$ & \hspace{-3em}Corner frequency of $S_{\Delta}(\omega)$&$\omega_R$ & Corner frequency of $R(\omega)$\\
\hline
\end{tabular}
\vspace*{50pt}
\end{table}

\section*{Conflicts of interest}
The authors have no conflict of interest to declare.

\section*{Acknowledgements}
We thank E. Y. Song and A. Ludlow for a careful reading and comments on the manuscript, and J. R. K. Cline and D. J. Young for helpful discussions.



\bibliographystyle{crunsrt}
\bibliography{crphys20221085}
\end{document}
