A quantum model of lasing without inversion

The textbook quantum description of a laser involves a gain medium of an ensemble of atomic three-level systems (2LS) coupled to a single mode photonic cavity [1, 2]. The population inversion of the atomic ensemble required for lasing is achieved via an incoherent driving process. The microscopic origin of this inversion is usually via coherent pumping of transitions to other levels which, when adiabatically eliminated, give rise to the standard two-level model. However, population inversion is not a necessary condition for lasing and lasing without inversion can be achieved through a degree of coherent control over three-level systems [3–5].

The various configurations available to three-level atoms allow a variety of other non-trivial optical phenomena such as coherent population trapping [6, 7], stimulated Raman adiabatic passage [8, 9], and electromagnetically induced transparency [10, 11]. However, when compared to two-level systems, the many-body physics of three-level atomic systems is relatively unexplored. Such theory is not only interesting but fundamental to the development of future quantum technologies: For instance, atomic clocks [12–14], the next generation of which incorporate lasers with many-body systems in order to engineer clock schemes with quantum noise-limited fractional stability [15]. Beyond this, three-level systems have proven to be a powerful tool in generating quantum memory for photons [16, 17], atomic cooling [18–20], and quantum computing [21–23].

Of particular relevance to this work is nanoscale laser [24] and maser [25] technology. As a consequence of their size, nanolasers exhibit a number of unique physical properties the most obvious of which is the low power consumption that comes as a result of the small gain medium [26, 27] as well as the thresholdless operation that comes with being far from the thermodynamic limit [28–30]. Such devices can achieve single-photon emission and thereby allow for the creation of photonic qubits; the encoding of quantum information in the photon state [31]. There are also possibilities to observe strong light–matter coupling effects [32, 33] and squeezing [34]. The physics of these lasing systems, from conventional laser to nanolaser, from macroscopic to microscopic is a vibrant and diverse area of research with a clear range of applications. In this paper we explore the lasing behaviour across these regimes by modelling the inversionless lasing of a cavity QED model of three-level atoms from the few- to the many-body limit.

The lasing phase transition of three-level systems became the subject of intense inquiry following the discovery of lasing without population inversion (LWI) [35–37]. Initial research into this area focused on gain increases created by utilising the recoil frequency shift experienced by an atom following the emission of a photon [38]. In two-level systems, this mechanism allows the amplification of radiation to occur without a larger population in the upper level than the lower level i.e. without population inversion [39].

To identify population inversion in multi-level systems it is usual to group sets of levels into the ground and excited manifolds of an effective two-level system. For instance, Λ-type three-level systems are structured with a single upper energy level and two lower energy levels that are closely spaced (non-degenerate) or identical (degenerate). In this case population inversion is defined as the collective population of the excited states being higher than that of the lower. Moving beyond the recoil mechanism, research into LWI in the context of three- and four-level systems was significantly developed by Scully, Zhu, and collaborators in a series of seminal works [3–5, 40]. Two distinct approaches to producing LWI in Λ-type three-level systems were proposed: The first utilises coherent population trapping (CPT) in degenerate and non-degenerate systems. The system is prepared with initial quantum coherence between the two lower levels, which causes quantum interference effects that cancel absorption. Stimulated emission increasingly becomes the only pathway available for the system and therefore results in an increased gain. When the transitions from lower to upper state are completely excluded then any population in the higher energy level leads to lasing. Thus, lasing is achieved without population inversion [41].

In the second approach the two lower levels of a Λ-system are coupled by a coherent field [40]. This set-up is known as a quantum beat-laser [42]. This early model considers a beam of these atoms passing through a cavity field. The coherent field splits and mixes the two atomic sublevels, inducing an interaction of the cavity with two resulting dressed states. The coherences generated between the dressed states and the upper level produce similar interference effects to those observed in the case of CPT, suppressing absorption.

Following its theoretical conception, LWI of Λ-type three-level systems has been experimentally investigated and subsequently observed in the case of CPT [43, 44] and the quantum beat laser [40]. Furthermore, recent work has shown how multiple V-type three-level systems can generate many-body dark and nearly-dark states through dissipative stabilization of excited states [45]. Beyond this, recent studies have identified additional second-order phase transitions of three-level systems separate to the laser threshold, revealing their relationships with the fundamental U(1) symmetry [46]. There is also interest in LWI of non-Hermitian systems, which can occur without any initial coherence [47, 48].

As the number of components of a system increases, the dimension of the objects, such as the Liouvillian, that describes the many-body physics grows exponentially. This makes calculating the dynamics and steady states challenging as large dimensions result in long computation times and impractical memory usage. By using the permutation symmetry of the density matrix the size of the Liouvillian can be reduced [49–52]. This approach effectively increases the number of atoms for which a solution can be calculated. As their number increases the dynamics of the atomic systems (and cavity) approach the thermodynamic limit wherein the lasing phase transition becomes increasingly sharp. Using this method allows us to examine the properties of the system at intermediate system sizes accessing this crossover region.

The structure of this paper is as follows: in section 2, we construct a microscopic model for a Λ-type three-level system, including a semi-classical microwave drive and all of the necessary incoherent processes. In section 3, we examine the behaviour in the thermodynamic limit by analysing the mean-field equations. Linear stability analysis of these equations gives us access to the phase diagram. This allows to calculate the location of the phase transition as a function of the phase of the microwave drive revealing the necessary parameters for LWI. Section 4 moves beyond the mean-field equations by deriving the second order cumulant equations. These give 1/N corrections to the mean-field equations which we are able to compare these results with exact solutions calculated using permutation symmetric methods. This demonstrates how the behaviour of these systems approaches the thermodynamic limit. Finally, in section 5, we present our conclusions.

We begin our analysis by constructing a microscopic model, as illustrated in figure 1(a), that is able to capture the physics of a quantum beat laser [40]. To this end, we consider a collection of degenerate Λ-type atomic three-level systems where the transition between ground states is coherently driven by an applied microwave field which we treat semi-classically. The upper level and the two ground levels of each 'atom' are labelled as |e⟩, |1⟩, and |2⟩ respectively.

Zoom In Zoom Out Reset image size

The Hamiltonian for N such systems interacting with a single cavity mode is given by:

$$\begin{equation}H={H}_{\mathrm{F}}+\sum\limits _{i}^{N}{H}_{\mathrm{A}\mathrm{F}}^{(i)}+{H}_{\mathrm{M}}^{(i)},\end{equation} \tag{ 1 }$$

where H_F is the field Hamiltonian and ${H}_{\mathrm{A}\mathrm{F}}^{(i)}$ describes the ith atom and its interaction with the cavity and ${H}_{\mathrm{M}}^{(i)}$ is the interaction with the semi-classical microwave drive. These are given by:

$$\begin{equation}{H}_{\mathrm{F}}=\nu {a}^{{\dagger}}a,\end{equation} \tag{ 2 }$$

$$\begin{equation}{H}_{\mathrm{A}\mathrm{F}}^{(i)}={\omega }_{\mathrm{e}}{\rho }_{\mathrm{e}\mathrm{e}}^{(i)}+\frac{g}{\sqrt{N}}\left[a({\rho }_{\mathrm{e}1}^{(i)}+{\rho }_{\mathrm{e}2}^{(i)})+\text{h.c.}\right],\end{equation} \tag{ 3 }$$

$$\begin{equation}{H}_{\mathrm{M}}^{(i)}={\Omega}\left({\text{e}}^{\text{i}\phi }{\rho }_{21}^{(i)}+\text{h.c.}\right).\end{equation} \tag{ 4 }$$

The cavity is described by the annihilation operator a and has frequency ν. To describe the atomic degrees of freedom we use the notation ${\rho }_{ab}^{(i)}=\vert a\rangle \langle b{\vert }_{i}$ where each state |k⟩ has an associated energy ω_k (ℏ = 1). For simplicity we assume that the two lower atomic energy levels are degenerate and use this as our reference energy such that ω₁ = ω₂ = 0. Furthermore, we assume that the coupling strength, g, between the light and matter is weak compared to the bare cavity and atomic transition frequencies and so we perform the rotating-wave approximation on this term, ignoring the counter-rotating terms [2]. The factor of $\sqrt{N}$ in the light matter coupling term ensures a well-defined thermodynamic limit when taking N → ∞. With this scaling when the system size is changed the collective light mater coupling remains constant as the number of emitters is increased [52–54]. Following the approach from reference [3], we have already moved to a rotating frame and performed a rotating wave approximation ³ with respect to the microwave field which is described by the Rabi frequency Ω which has a phase of ϕ relative to the light–matter coupling. A gauge transformation which removes the phase from this term in the Hamiltonian would introduce an equivalent phase on the light–matter coupling g. We note that the model outlined here does not include any microwave induced relaxation between the two lower levels.

To include the effects of incoherent processes along with the coherent ones included in the above Hamiltonian we write a Lindblad form quantum master equation for the density matrix ρ of the combined cavity and atomic degrees of freedom. This takes the form [4, 55]:

$$\begin{align}\hfill \dot{\rho }& =-\mathrm{i}\left[H,\rho \right]+\kappa \mathcal{D}[a]\rho +{{\Gamma}}^{{\uparrow}}\sum\limits _{i}^{N}\left(\mathcal{D}[{\rho }_{\mathrm{e}1}^{(i)}]+\mathcal{D}[{\rho }_{\mathrm{e}2}^{(i)}]\right)\rho +{{\Gamma}}^{{\downarrow}}\sum\limits _{i}^{N}\left(\mathcal{D}[{\rho }_{1\mathrm{e}}^{(i)}]+\mathcal{D}[{\rho }_{2\mathrm{e}}^{(i)}]\right)\rho ,\hfill \end{align} \tag{ 5 }$$

where the usual Lindblad superoperator is $\mathcal{D}[x]\rho =x\rho {x}^{{\dagger}}-\frac{1}{2}\left\{{x}^{{\dagger}}x,\rho \right\}$. We have introduced cavity relaxation at rate κ, and a common incoherent pumping rate Γ^↑ from both lower levels to |e⟩, mirrored by an incoherent loss of atomic excitations at rate Γ^↓. These incoherent processes can be realised by using a coherent drive to an additional level, as is standard for regular two-level lasing [2], or as the result of a weak coherent drive applied directly to the relevant transition. The second case would, however, restrict the range of parameters to Γ^↑ ⩽ Γ^↓.

This model has a U(1) symmetry which corresponds to shifting the phase of the cavity and atomic coherences a → e^iθ a, ρ_e1 → e^−iθ ρ_e1, and ρ_e2 → e^−iθ ρ_e2. The lasing transition of this model occurs when this symmetry is spontaneously broken. For the purpose of comparison a parallel description for the case of two-level instead of three-level systems is provided in appendix A.

This model could be realised in a variety of physical systems. An obvious approach utilises cold atoms trapped in optical cavities. In these systems fine control over the required level structure, interaction strengths and decay rates can be achieved [56]. Other potential platforms where this kind of physics could be realised includes atomic ensembles in vapour cells [57] or rare earth ion doped solid-state ensembles [58] where Λ-like Raman level structures have been shown to enable quantum memories. Impurities in diamond [25] also hold promise if a suitable level structure can be engineered.

In the thermodynamic limit of a large ensemble of emitters, N → ∞, it is well known that this class of model is described exactly by mean-field theory [1]. This allows us to characterise the point at which the lasing phase transition occurs as a function of the parameters of the master equation. In this limit the transition is identified by the point at which the cavity population per atom ⟨a^† a⟩/N becomes finite; below the lasing transition, in the normal state, $\langle {a}^{{\dagger}}a\rangle \sim \sqrt{N}$ while above the transition ⟨a^† a⟩ ∼ N. To derive the mean-field equations we make use of the fact that all of the emitters are identical and then use the master equation, (5) to find the evolution of each of the required operators using the relation $\langle \dot{x}\rangle =\mathrm{Tr}[x\dot{\rho }]$. In the limit where mean-field theory is valid we may neglect second order and higher cumulants and write a closed set of equations by breaking second order correlations such as ⟨aρ_e1⟩ = ⟨a⟩⟨ρ_e1⟩ [1, 52]. We will return to the effect that approximation has at finite N in section 4. We hence obtain equations that describe the evolution of the expectation values of the complex photon amplitude and atomic degrees of freedom:

$$\begin{equation} \langle {\dot{\rho }}_{\mathrm{e}1}\rangle =\mathrm{i}\left[{\omega }_{\mathrm{e}}\langle {\rho }_{\mathrm{e}1}\rangle +\frac{g}{\sqrt{N}}\langle {a}^{{\dagger}}\rangle \left(\langle {\rho }_{11}\rangle -\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle +\langle {\rho }_{21}\rangle \right)-{\Omega}\,{\text{e}}^{-\text{i}\phi }\langle {\rho }_{\mathrm{e}2}\rangle \right]-{{\Gamma}}_{\phi }\langle {\rho }_{\mathrm{e}1}\rangle ,\end{equation} \tag{ 6 }$$

$$\begin{equation}\qquad \langle {\dot{\rho }}_{\mathrm{e}2}\rangle =\mathrm{i}\left[{\omega }_{\mathrm{e}}\langle {\rho }_{\mathrm{e}2}\rangle +\frac{g}{\sqrt{N}}\langle {a}^{{\dagger}}\rangle \left(\langle {\rho }_{22}\rangle -\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle +\langle {\rho }_{12}\rangle \right)-{\Omega}\,{\text{e}}^{\text{i}\phi }\langle {\rho }_{\mathrm{e}1}\rangle \right]-{{\Gamma}}_{\phi }\langle {\rho }_{\mathrm{e}2}\rangle ,\end{equation} \tag{ 7 }$$

$$\begin{equation}\qquad \langle {\dot{\rho }}_{12}\rangle =\mathrm{i}\left[-\frac{g}{\sqrt{N}}\left(\langle {a}^{{\dagger}}\rangle \langle {\rho }_{1\mathrm{e}}\rangle -\langle a\rangle \langle {\rho }_{\mathrm{e}2}\rangle \right)+{\Omega}\,{\text{e}}^{\text{i}\phi }\left(\langle {\rho }_{22}\rangle -\langle {\rho }_{11}\rangle \right)\right]-{{\Gamma}}^{{\uparrow}}\langle {\rho }_{12}\rangle ,\end{equation} \tag{ 8 }$$

$$\begin{equation}\qquad \langle {\dot{\rho }}_{\mathrm{e}\mathrm{e}}\rangle =-\frac{g}{\sqrt{N}}\left[\mathrm{i}\langle a\rangle (\langle {\rho }_{\mathrm{e}1}\rangle +\langle {\rho }_{\mathrm{e}2}\rangle )+\text{c.c.}\right]+{{\Gamma}}^{{\uparrow}}\left(\langle {\rho }_{11}\rangle +\langle {\rho }_{22}\rangle \right)-2{{\Gamma}}^{{\downarrow}}\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle ,\end{equation} \tag{ 9 }$$

$$\begin{equation}\,\langle {\dot{\rho }}_{11}\rangle =\left[\mathrm{i}\left(\frac{g}{\sqrt{N}}\langle a\rangle \langle {\rho }_{\mathrm{e}1}\rangle +{\Omega}\,{\text{e}}^{\text{i}\phi }\langle {\rho }_{21}\rangle \right)+\text{c.c.}\right]-{{\Gamma}}^{{\uparrow}}\langle {\rho }_{11}\rangle +{{\Gamma}}^{{\downarrow}}\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle ,\end{equation} \tag{ 10 }$$

$$\begin{equation}\langle {\dot{\rho }}_{22}\rangle =\left[\mathrm{i}\left(\frac{g}{\sqrt{N}}\langle a\rangle \langle {\rho }_{\mathrm{e}2}\rangle -{\Omega}\,{\text{e}}^{\text{i}\phi }\langle {\rho }_{21}\rangle \right)+\text{c.c.}\right]-{{\Gamma}}^{{\uparrow}}\langle {\rho }_{22}\rangle +{{\Gamma}}^{{\downarrow}}\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle ,\end{equation} \tag{ 11 }$$

$$\begin{equation} \langle \dot{a}\rangle =-\mathrm{i}\left[\nu \langle a\rangle +g\sqrt{N}\left(\langle {\rho }_{1\mathrm{e}}\rangle +\langle {\rho }_{2\mathrm{e}}\rangle \right)\right]-\frac{\kappa }{2}\langle a\rangle ,\end{equation} \tag{ 12 }$$

where Γ_ϕ = (Γ^↑ + 2Γ^↓)/2 is the total atomic dephasing rate. This is induced by the presence of the two incoherent rates, Γ^↑ and Γ^↓, which describing energy relaxation processes and hence give rise to an associated dephasing process. From these equations we can extract the conditions under which the lasing phase transition occurs. To derive an expression for the location of the phase transition we consider the linear stability of the normal state fixed point. It is straightforward to show that the coherences ⟨a⟩_ns = 0, ${\langle {\rho }_{1\mathrm{e}}\rangle }_{\mathrm{n}\mathrm{s}}=0$, ${\langle {\rho }_{2\mathrm{e}}\rangle }_{\mathrm{n}\mathrm{s}}=0$ and ${\langle {\rho }_{12}\rangle }_{\mathrm{n}\mathrm{s}}=0$ and populations ${\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle }_{\mathrm{n}\mathrm{s}}={{\Gamma}}^{{\uparrow}}/2{{\Gamma}}_{\phi }$, ${\langle {\rho }_{11}\rangle }_{\mathrm{n}\mathrm{s}}={\langle {\rho }_{22}\rangle }_{\mathrm{n}\mathrm{s}}={{\Gamma}}^{{\downarrow}}/2{{\Gamma}}_{\phi }$ are always a solution of equations (6)–(12). We can then linearise around this fixed point to find the point at which it becomes unstable and hence the system enters the lasing phase. We find that the only terms which give non-trivial contributions to this linear stability are δa = ⟨a⟩ − ⟨a⟩_ns, $\delta {\rho }_{1\mathrm{e}}=\langle {\rho }_{1\mathrm{e}}\rangle -{\langle {\rho }_{1\mathrm{e}}\rangle }_{\mathrm{n}\mathrm{s}}$ and $\delta {\rho }_{2\mathrm{e}}=\langle {\rho }_{2\mathrm{e}}\rangle -{\langle {\rho }_{2\mathrm{e}}\rangle }_{\mathrm{n}\mathrm{s}}$. ⁴ Hence we can write the linearised equation for these fluctuations

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix}\hfill \delta a\hfill \\ \hfill \delta {\rho }_{1\mathrm{e}}\hfill \\ \hfill \delta {\rho }_{2\mathrm{e}}\hfill \end{matrix}\right)=J\left(\begin{matrix}\hfill \delta a\hfill \\ \hfill \delta {\rho }_{1\mathrm{e}}\hfill \\ \hfill \delta {\rho }_{2\mathrm{e}}\hfill \end{matrix}\right),\end{equation} \tag{ 13 }$$

where the Jacobian matrix is

$$\begin{equation}J=\left(\begin{matrix}\hfill -\mathrm{i}\nu -\frac{\kappa }{2}\hfill & \hfill -\mathrm{i}g\sqrt{N}\hfill & \hfill -\mathrm{i}g\sqrt{N}\hfill \\ \hfill \mathrm{i}g\left({\langle {{\Delta}}_{\mathrm{e}1}\rangle }_{\mathrm{n}\mathrm{s}}-{\langle {\rho }_{12}\rangle }_{\mathrm{n}\mathrm{s}}\right)/\sqrt{N}\hfill & \hfill \mathrm{i}{\omega }_{\mathrm{e}}-{{\Gamma}}_{\phi }\hfill & \hfill \mathrm{i}{\Omega}\,{\text{e}}^{\text{i}\phi }\hfill \\ \hfill \mathrm{i}g\left({\langle {{\Delta}}_{\mathrm{e}2}\rangle }_{\mathrm{n}\mathrm{s}}-{\langle {\rho }_{21}\rangle }_{\mathrm{n}\mathrm{s}}\right)/\sqrt{N}\hfill & \hfill \mathrm{i}{\Omega}\,{\text{e}}^{-\text{i}\phi }\hfill & \hfill \mathrm{i}{\omega }_{\mathrm{e}}-{{\Gamma}}_{\phi }\hfill \end{matrix}\right),\end{equation} \tag{ 14 }$$

and we have introduced the operators

$$\begin{equation}{{\Delta}}_{\mathrm{e}1}={\rho }_{\mathrm{e}\mathrm{e}}-{\rho }_{11},\quad {{\Delta}}_{\mathrm{e}2}={\rho }_{\mathrm{e}\mathrm{e}}-{\rho }_{22}.\end{equation} \tag{ 15 }$$

When the real part of one of the eigenvalues of equation (14) becomes positive the normal state is unstable and the cavity is occupied, thereby allowing us to identify the phase transition from the eigenvalues of the characteristic polynomial of the Jacobian:

$$\begin{align}\hfill & {\lambda }^{3}+\left(2{{\Gamma}}_{\phi }+\mathrm{i}(\nu +2{\omega }_{\mathrm{e}})+\frac{\kappa }{2}\right){\lambda }^{2}+\left[{g}^{2}\frac{({{\Gamma}}^{{\downarrow}}-{{\Gamma}}^{{\uparrow}})}{{{\Gamma}}_{\phi }}+({{\Gamma}}_{\phi }+\mathrm{i}{\omega }_{\mathrm{e}})\left({{\Gamma}}_{\phi }+\mathrm{i}(2\nu +{\omega }_{\mathrm{e}})+\kappa \right)+{{\Omega}}^{2}\right]\lambda \hfill \\ \hfill & \hspace{25.0pt}+{g}^{2}\frac{({{\Gamma}}^{{\downarrow}}-{{\Gamma}}^{{\uparrow}})}{{{\Gamma}}_{\phi }}\left[{{\Gamma}}_{\phi }+\mathrm{i}\left({\omega }_{\mathrm{e}}+{\Omega}\,\mathrm{cos}\left(\phi \right)\right)\right]+\frac{1}{2}\left(\kappa +2\mathrm{i}\nu \right)\left[{\left({{\Gamma}}_{\phi }+\mathrm{i}{\omega }_{\mathrm{e}}\right)}^{2}+{{\Omega}}^{2}\right]=0,\hfill \end{align} \tag{ 16 }$$

where λ is an eigenvalue of J. Therefore, solving for λ corresponds to finding the roots of equation (16). The location of the phase transition as a function of the pumping rate Γ^↑ has a non-trivial dependence on the phase of the applied microwave drive, ϕ. Hence, the structure of the eigenvalues is more complicated than that of a two-level system which we present for comparison in appendix A.

3.1. Mean-field results

We can now use the mean-field equations and linear stability analysis to examine the optimal conditions for observing lasing in this model. To develop intuition for the underlying physics it is useful to first look in more detail at the structure of the Hamiltonian. The splitting of the energy levels by the microwave field effectively shifts the resonant cavity frequency. This has a significant impact on the lasing phase transition. In order to capture this behaviour and identify the shifted resonant frequency, we rewrite the Hamiltonian, equation (1), in terms of the following dressed states which diagonalises the terms which do not involve the cavity

$$\begin{equation}\vert +\rangle =\frac{1}{\sqrt{2}}\left({\text{e}}^{\text{i}\phi /2}\vert 2\rangle +{\text{e}}^{-\text{i}\phi /2}\vert 1\rangle \right)\end{equation} \tag{ 17 }$$

and

$$\begin{equation}\vert -\rangle =\frac{1}{\sqrt{2}}\left({\text{e}}^{\text{i}\phi /2}\vert 2\rangle -{\text{e}}^{-\text{i}\phi /2}\vert 1\rangle \right).\end{equation} \tag{ 18 }$$

Using these states we move to a frame rotating at the cavity frequency to give the Hamiltonian

$$\begin{align}\hfill {H}_{\mathrm{R}}& =\sum\limits _{i}^{N}{\Delta}\omega {\rho }_{\mathrm{e}\mathrm{e}}^{(i)}+{\omega }_{+}{\rho }_{++}^{(i)}+{\omega }_{-}{\rho }_{--}^{(i)}+\frac{g}{\sqrt{N}}\left(a\left[\mathrm{cos}\left(\frac{\phi }{2}\right){\rho }_{e+}^{(i)}-\mathrm{i}\,\mathrm{sin}\left(\frac{\phi }{2}\right){\rho }_{\mathrm{e}-}^{(i)}\right]+\text{h.c.}\right)\hfill \end{align} \tag{ 19 }$$

where Δω = ω_e − ν and ω_± = ±Ω. This is illustrated in figure 1(c).

Outside of certain parameter regimes the characteristic polynomial, equation (16), does not have straightforward solutions and so we numerically compute the eigenvalues of the Jacobian, equation (14), to find the phase diagram as a function of the atomic pumping rate Γ^↑ and atom-cavity detuning Δω. It is important to note that both Γ^↑ and Γ^↓ describe two individual incoherent processes, as such we define the total rate Γ_T = 2(Γ^↑ + Γ^↓) which allows us to normalise the rates. Examining how this phase diagram changes for different values of the phase of the microwave drive ϕ provides a complete picture of the lasing behaviour in the thermodynamic limit. Figure 2 shows these phase diagrams for different values of ϕ. The y-axis is in terms of Γ^↑/Γ_T. We vary this quantity while keeping Γ_T constant such that when Γ^↑/Γ_T = 0 there is no pumping and when 2Γ^↑/Γ_T = 0.5 there is no atomic loss. In these diagrams we highlight where Δω = ±Ω since this corresponds to the detunings where the cavity is on resonance with the microwave-dressed atomic states described above.

Zoom In Zoom Out Reset image size

In figure 2(a), for ϕ = 0 we see a single minimum in the lasing region at Δω = Ω. This can be understood as follows: As previously stated, the coherent microwave field mixes the lower two atomic energy levels with driving strength Ω. This results in the hierarchy of energy levels shown in figures 1(b) and (c) where ω_± = ±Ω. As such this forces a particular phase dependence between the state of the atoms and the emitted photons. The term describing the interaction between the cavity and the atomic transitions in Hamiltonian (3) accounts for this dependence and thereby defines a phase reference for the microwave field. As a result, applying the microwave drive with phase ϕ = 0 favours the formation of the |+⟩ state and hence we observe that lasing occurs when the cavity is tuned to the frequency of this state. In the opposite limit of ϕ = π/2, shown in figure 2(c), we see two symmetric local minima at Δω = ±Ω (with some corrections proportional to the cavity relaxation rate κ). In that case the interplay between light–matter interaction and the microwave drive favours the formation of an equal superposition of the |+⟩ and |−⟩ states and so lasing occurs at both frequencies. When the microwave phase is between these limits as in figure 2(b), where ϕ = π/4, we see a smooth crossover between these limits. In this case the minima at Δω = Ω remains largely unchanged. In all cases we see that close to the minimum required pumping the lasing state is stable while the inversion is still negative. As the pump power required is increased the observed inversion level also increases. The results shown here outline the range of possible behaviours, for completeness however, in appendix B we also show how the phase diagram depends on the total rate Γ_T and the microwave phase ϕ.

To further understand how the location of these minima behave we track their locations in terms of both their frequency and pump strength as a function of the microwave phase ϕ. We capture this behaviour in figure 3 where the location of the minimum at positive detuning is shown in red while the one at negative detuning is in blue. We see that there are large plateaus close to ϕ = 0 and ϕ = ±π where the microwave drive allows lasing from only one of the dressed states. This results in the lowest value of the incoherent pump strength at each of these points. At the centre of these plateaus only a single positively or negatively detuned minimum exists, its associated Δω corresponds to the cavity mode being exactly resonant with one of the atomic dressed states, as expected. A shift away from that to larger detuning is observed in the overlap regions, and is particularly pronounced for the respective 'higher lying' (i.e. stronger pumping to reach the phase transition) minimum. From these results we are able to find an exact description of the critical pumping required to achieve lasing for particular cases. Exactly at points ϕ = 0, ±π we are able to find a simple analytic expression for the critical pumping rate from the characteristic polynomial in equation (16):

$$\begin{equation}{{\Gamma}}_{\text{crit}}^{{\uparrow}}={{\Gamma}}_{\mathrm{T}}+8\frac{{g}^{2}}{\kappa }\left(1-\sqrt{1+\frac{3{{\Gamma}}_{\mathrm{T}}}{16}\frac{\kappa }{{g}^{2}}}\right),\end{equation} \tag{ 20 }$$

which is valid exactly on resonance. Here, we are incoherently pumping from both lower levels such that in the limit g ≫ κ we find ${{\Gamma}}_{\text{crit}}^{{\uparrow}}/{{\Gamma}}_{\mathrm{T}}=1/4$. For the same limit in the case of the two-level system, where the incoherent pumping is only from a single lower level, we find instead the equivalent result ${{\Gamma}}_{\text{crit}}^{{\uparrow}}/{{\Gamma}}_{\mathrm{T}}=1/2$.

Zoom In Zoom Out Reset image size

To explore this further it useful to use the full mean-field equations (6)–(12) to find the steady state population inversion and cavity population generated by this model as a function of incoherent pumping Γ^↑/Γ_T. We define the population inversion as the difference between the excited state population and the total population in the ground state manifold ρ_ee − (ρ₁₁ + ρ₂₂). In figure 4 we show how this changes as the phase of the microwave field is varied from ϕ = 0 to ϕ = 3π/4 for two different values of the total atomic decoherence rate Γ_T. The main result here is that we observe a strong dependence on the nature of the phase transition with changing ϕ. Lasing without inversion when the phase of the microwave drive is tuned to ϕ ≲ π/2 such that the |+⟩ state is favoured. This only occurs when Γ_T is small enough since when this is increased a larger population is required in the excited state to overcome the decay. When ϕ ≳ π/2 significant inversion is required before the phase transition occurs. This is because, for these values of ϕ, the microwave drive populates the |−⟩ state which is far detuned.

Zoom In Zoom Out Reset image size

A significant difference from the usual two-level laser model is the lack of mode-locking. As we see in appendix A above the lasing transition the population inversion locks to a fixed value and all of the additional energy from the pumps goes into the cavity field. We do observe behaviour close to mode-locking for ϕ = 3π/4 and Γ_T = Ω, whereas for Γ_T = 3Ω this occurs at ϕ = π/2. This effect originates from a different mechanism to the usual case. In the case where the laser is far detuned (for the largest value of ϕ shown) we see that the population inversion decreases above the transition threshold. When there is a large coherent field inside the cavity the atomic dynamics becomes too fast for the microwave drive to have any effect and so the dressing of the atomic states is reduced and the inversion along with it. In the limit of a very large occupation of the cavity the correct picture for the atomic system becomes the bare states.

Whilst mean-field theory provides useful information regarding the phase transition for traditional lasing where N → ∞, it is not appropriate for small or intermediate numbers of three-level systems. To analyse this limit in the following we derive a set of second-order cumulant equations which provide a systematic way of calculating corrections of order 1/N to the mean-field results allowing access to the intermediate N regime. Along with this we also provide exact numerical solutions to our microscopic model for small N. These exact simulations make use of the permutation symmetry of the Liouvillian [51] allowing us to simulate around N ∼ 10 three-level systems.

To derive the second-order cumulant equations we retain correlations between pairs of operators and instead break the third order moments into products of first and second order correlations [52, 59, 60]. This procedure then produces a closed set of equations for the first and second moments which in the limit N → ∞ reproduce the mean-field equations and include 1/N corrections which vanish in this limit. The full set of equations is cumbersome and so we give full details of the derivation in appendix B. Furthermore, a parallel derivation of the cumulant equations for the equivalent model with two-level systems is provided in appendix A for comparison.

For small values of N we make use of exact simulations of the master equation (5) to find the steady state density operator. Directly simulating this scales exponentially with the number of three-level systems included and is difficult even for N ∼ 5. Instead we make use of the permutation symmetry of the Liouvillian which arises because each atom is identical and couples identically to the cavity mode. This approach is fully outlined in references [50, 51] whilst the code which our results are based on can be found at [61]. This technique has been previously applied to a variety of models with two-level systems. For example plasmons [49], spin ensembles [62], subradiant states [63, 64], superradiant phase transitions in the Dicke model [50, 65] and fully connected Ising models [66]. The generalisation of these approaches to three-level systems is straightforward. We show the compression achieved by using this method in figure 5. The dimension of the Liouville space, L_D, that describes N copies of our three-level system scales exponentially as 3^2N however using the permutation symmetry this scaling is reduced to polynominal in the system size. For example at N = 2, the permutation approach produces reduces L_D from 81 to 45 while at N = 8, the full space has L_D = 43 046 721 which is reduced down to a much more manageable L_D = 12 870 by the permutation symmetry. We note here that the full simulations need to include the Hilbert space of the cavity resulting in an increase to L_D by a factor of P², where P is the number of Fock states used to represent the cavity field.

Zoom In Zoom Out Reset image size

We begin by looking at how the results at finite N scale towards the thermodynamic limit. In figure 6 we show how the second-order cumulant results converge upon those of the mean-field as N increases for a representative set of parameters both above and below the lasing phase transition and for different values of ϕ. For each case we calculate the reduced photon number ⟨a^† a⟩/N and population inversion. Alongside these results we also show the exact numerical steady state solutions of the master equation (5) for N systems, truncating the cavity Hilbert space to P = 10 Fock states. This allows for a direct comparison with the mean-field theory discussed in the previous section and captures the importance of cross-correlations to the microscopic description of this system.

Zoom In Zoom Out Reset image size

For both ϕ = 0 and ϕ = 3π/4 we see very good agreement between the exact permutation solutions and the corresponding second-order cumulant solutions, even at relatively small values of N. The agreement is generally better for the atomic inversion but even the results for the cavity photon number show convergence as N increases. At the same time we see that the cumulant results approach the mean-field limit as N increases. We thus show that by using a combination of these techniques it is possible to access the behaviour across the whole range of system sizes from the microscopic N ∼ 1 to the macroscopic N → ∞. We note however that the parameters chosen here are fairly specific. We needed to make sure that the steady state photon number in the lasing state is small enough that the resulting coherent state could be represented with only P = 10 Fock states. It is likely that, for regimes where the photon number is larger, the required value of N to see convergence between these different approaches also increases. This is supported by the results outlined in appendix A where the two-level system shows similar agreement for an equivalent set of parameter choices.

Sweeping the pumping rate through the lasing transition, figure 7 provides a comparison between the steady state solutions for population inversion and cavity population calculated from the second order cumulant terms and permutation approach for N = 2 and N = 8. As before, mean-field theory provides comparison for these in the N → ∞ limit. The two distinct cases ϕ = 0 and ϕ = 3π/4 where the transition is below and above inversion are explored in order to demonstrate that this effect can be seen using all of our approaches. We again see that, even for small values of N the agreement between the exact solutions and the cumulant expansion is very good, only breaking down slightly when the cavity occupation is large. This is region where the atom-cavity interaction has the most significant effect meaning that higher order cumulants can be generated. In every case shown by figure 7 we see the phase transition sharpening as N increases with the discontinuity only apparent in the mean-field results. We note that, figures 7(a) and (b) include data points in which the Fock space truncation P = 10 is a limiting factor at the largest pumping strengths shown. This slightly reduces the observed photon occupation compared to the exact value, as such we have indicated when the occupation of the highest energy Fock state exceeds 1/1000. This is most apparent at high pumping with ϕ = 0 which allows for a large cavity occupation.

Zoom In Zoom Out Reset image size

To examine the nature of the light emitted by the cavity it is useful to examine higher order correlation functions of the photon distribution. These calculations are only straightforward for the exact results and so we are limited to the small numbers of emitters which are possible with the permutation approach. In figure 8 we show both the second order coherence

$$\begin{equation}{g}^{(2)}(0)=\frac{\langle {a}^{{\dagger}}{a}^{{\dagger}}aa\rangle }{\langle {a}^{{\dagger}}{a\rangle }^{2}},\end{equation} \tag{ 21 }$$

and Fano factor

$$\begin{equation}\mathcal{F}=\frac{\langle {a}^{{\dagger}}a{a}^{{\dagger}}a\rangle -{\langle {a}^{{\dagger}}a\rangle }^{2}}{\langle {a}^{{\dagger}}a\rangle }.\end{equation} \tag{ 22 }$$

As the system is driven through the lasing transition we expect both of the quantities to show signatures of the change in the state of the radiation. The second order coherence reduces from g⁽²⁾(0) = 2 which is characteristic of thermal light to g⁽²⁾(0) = 1 which is characteristic of coherent light [2]. The change between these two states becomes sharper as the thermodynamic limit is reached. The Fano factor in contrast is a measure of the scale of the fluctuations in the photon number, as such it shows a distinct peak at the transition point which corresponds to the large photon number uncertainty at this point. We observe the expected trend for both of these effects in the results for the three-level systems, with the signatures becoming clearer for N = 8. Our results for the two-level model in appendix A, where we are able to go to much larger N, are even more pronounced and clearly approach the expectations for the thermodynamic limit.

Zoom In Zoom Out Reset image size

We have constructed a Hamiltonian for a Λ-type three-level system where the introduction of a microwave field provides a fully quantum model of the lasers originally described in reference [40]. We explicitly demonstrate how the interplay between the microwave dressed states and the light–matter interaction can lead to LWI for the correct choice of phase.

We characterised the system dynamics through mean-field equations, the solutions of which give the exact behaviour in the thermodynamic limit N → ∞, and demonstrate how the location of the minimum in the lasing phase diagram is detuned by an amount Δω ≈ Ω as a result of the introduction of the microwave field. Linear stability analysis allowed us to derive a complete analytical expression for the location of the phase transition at certain values of the phase of the microwave field, ϕ. We show how the properties of the lasing transition can be optimised by adjusting this detuning along with the microwave phase.

We go beyond mean-field theory to examine the behaviour of the system at finite system sizes. To do this we use two approaches: a cumulant expansion which gives access to corrections of order 1/N and exact numerical solutions, calculated using permutation symmetry which allow us to exactly find the steady state density operator of the system at small and intermediate numbers of three-level systems. The combination of these approaches lets us characterise the properties of the lasing transition across the whole range of possible system sizes with good agreement in the expected regimes. For the parameters chosen here we were able to see agreement between the cumulant expansion and exact numerics even at N = 8.

The microwave field described here introduces a degree of quantum control over the laser phase transition of the system. Furthermore, the framework outlined would provide a genuine advantage for understanding intermediate N systems where cross-correlations accurately captured by the cumulant expansion have a significant impact upon steady state behaviour. Such systems, with a bespoke architecture of energy levels, are at the heart of quantum technologies beyond lasing such as atomic clocks. Therefore, with the quantum control demonstrated here alongside a rich analysis of finite sized systems, there is great potential to build upon this work by exploring further quantum systems which can be experimentally realised through nanophotonics.

NW and EMG acknowledge funding from the EPSRC Grant No. EP/T007214/1.

The data that support the findings of this study are openly available at the following URL/DOI: https://github.com/NickWerren/LWI_res [68].

For comparison, and to highlight the features that uniquely arise from the microwave-driven three-level Λ structure, we here provide the analogous results for a collection of incoherently pumped two-level systems coupled to a lossy cavity mode. Our approach follows previous work, specifically the textbook Dicke and Scully-Lamb laser model [1, 2]. Here we use—to the extent possible—the same parameters and notation as in the main text to aid with this comparison. We again describe the cavity using a single mode approximation

$$\begin{equation}{H}_{\mathrm{F}}=\nu {a}^{{\dagger}}a.\end{equation} \tag{ A.1 }$$

The collection of two-level atomic systems and their interaction with the cavity are described by the following Hamiltonian

$$\begin{equation}{H}_{\mathrm{A}\mathrm{F}}^{(i)}={\omega }_{\mathrm{e}}{\rho }_{\mathrm{e}\mathrm{e}}^{(i)}+\frac{g}{\sqrt{N}}\left({a}^{{\dagger}}{\rho }_{\mathrm{g}\mathrm{e}}^{(i)}+a{\rho }_{\mathrm{e}\mathrm{g}}^{(i)}\right),\end{equation} \tag{ A.2 }$$

where we again use the energy of the ground state as our reference ω_g = 0 and ${\rho }_{jk}^{(i)}=\vert j\rangle \langle k{\vert }^{(i)}$ is an operator on the ith two-level system. In the atom-field Hamiltonian, we make the rotating wave approximation, discarding the counter-rotating terms, and assume g is real. As such, the resulting quantum master equation for N systems is [cf equation (5)]:

$$\begin{equation}\dot{\rho }(t)=-\mathrm{i}\left[H,\rho (t)\right]+\kappa \mathcal{D}[a]+\sum\limits _{i}{{\Gamma}}^{{\uparrow}}\mathcal{D}[{\rho }_{\mathrm{e}\mathrm{g}}^{(i)}]+{{\Gamma}}^{{\downarrow}}\mathcal{D}[{\rho }_{\mathrm{g}\mathrm{e}}^{(i)}],\end{equation} \tag{ A.3 }$$

where

$$\begin{equation}H={H}_{\mathrm{F}}+\sum\limits _{i}{H}_{\mathrm{A}\mathrm{F}}^{(i)}.\end{equation} \tag{ A.4 }$$

This model again has a U(1) symmetry corresponding to the phase of the cavity field operator and atomic coherence which is spontaneously broken at the lasing phase transition.

As with the three-level system, we construct a set of mean-field equations for the operators of the two-level system which are as follows:

$$\begin{equation}\langle {\dot{\rho }}_{\mathrm{g}\mathrm{e}}\rangle =-\mathrm{i}\left({\omega }_{\mathrm{e}}\langle {\rho }_{\mathrm{g}\mathrm{e}}\rangle -\frac{g}{\sqrt{N}}\langle a\rangle \langle {\sigma }_{\mathrm{z}}\rangle \right)-\frac{{{\Gamma}}_{\mathrm{T}}}{2}\langle {\rho }_{\mathrm{g}\mathrm{e}}\rangle ,\end{equation} \tag{ A.5 }$$

$$\begin{equation}\langle {\dot{\sigma }}_{\mathrm{z}}\rangle =-\frac{2g}{\sqrt{N}}\left(\mathrm{i}\langle a\rangle \langle {\rho }_{\mathrm{e}\mathrm{g}}\rangle +\text{c.c.}\right)+{{\Gamma}}_{\mathrm{T}}\langle {\sigma }_{\mathrm{z}}\rangle +({{\Gamma}}^{{\uparrow}}-{{\Gamma}}^{{\downarrow}}),\end{equation} \tag{ A.6 }$$

$$\begin{equation}\langle \dot{a}\rangle =-\mathrm{i}\left(\nu \langle a\rangle +g\sqrt{N}{\rho }_{\mathrm{g}\mathrm{e}}\right)-\frac{\kappa }{2}\langle a\rangle ,\end{equation} \tag{ A.7 }$$

where Γ_T = Γ^↑ + Γ^↓, σ_z = ρ_ee − ρ_gg, and ⟨ρ_ee⟩ + ⟨ρ_gg⟩ = 1. From the mean-field equations, linear stability analysis of the normal state can be used to further understand the phase transition. The normal state is always defined such that ${\langle a\rangle }_{\mathrm{n}\mathrm{s}}={\langle {\rho }_{\mathrm{e}\mathrm{g}}\rangle }_{\mathrm{n}\mathrm{s}}=0$ and ${\langle {\sigma }_{\mathrm{z}}\rangle }_{\mathrm{n}\mathrm{s}}=({{\Gamma}}^{{\uparrow}}-{{\Gamma}}^{{\downarrow}})/{{\Gamma}}_{\mathrm{T}}$. Linearising around this solution, such that δa = ⟨a⟩ − ⟨a⟩_ns, $\delta {\rho }_{\mathrm{e}\mathrm{g}}=\langle {\rho }_{\mathrm{e}\mathrm{g}}\rangle -{\langle {\rho }_{\mathrm{e}\mathrm{g}}\rangle }_{\mathrm{n}\mathrm{s}}$ and $\delta {\sigma }_{\mathrm{z}}=\langle {\sigma }_{\mathrm{z}}\rangle -{\langle {\sigma }_{\mathrm{z}}\rangle }_{\mathrm{n}\mathrm{s}}$ the resulting matrix describes the linear stability of the normal state is,

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix}\hfill \delta a\hfill \\ \hfill \delta {\rho }_{\mathrm{g}\mathrm{e}}\hfill \\ \hfill \delta {\sigma }_{\mathrm{z}}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill -\mathrm{i}\nu -\frac{\kappa }{2}\hfill & \hfill -\mathrm{i}g\sqrt{N}\hfill & \hfill 0\hfill \\ \hfill \mathrm{i}\frac{g}{\sqrt{N}}{\langle {\sigma }_{\mathrm{z}}\rangle }_{\mathrm{n}\mathrm{s}}\hfill & \hfill -\mathrm{i}{\omega }_{\mathrm{e}}-\frac{{{\Gamma}}_{\mathrm{T}}}{2}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -{{\Gamma}}_{\mathrm{T}}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill \delta a\hfill \\ \hfill \delta {\rho }_{\mathrm{g}\mathrm{e}}\hfill \\ \hfill \delta {\sigma }_{\mathrm{z}}\hfill \end{matrix}\right).\end{equation} \tag{ A.8 }$$

We note that to find the phase transition from this matrix we do not need to consider the fluctuations of δσ_z since this decouples from the other quantities and is always stable. By calculating the characteristic polynomial of this matrix,

$$\begin{equation}{\lambda }^{2}+\left[\frac{{{\Gamma}}_{\mathrm{T}}+\kappa }{2}-\mathrm{i}({\omega }_{\mathrm{e}}+\nu )\right]\lambda +{g}^{2}\frac{\left({{\Gamma}}^{{\uparrow}}-{{\Gamma}}^{{\downarrow}}\right)}{{{\Gamma}}_{\mathrm{T}}}+\left(\frac{{{\Gamma}}_{\mathrm{T}}}{2}-\mathrm{i}{\omega }_{\mathrm{e}}\right)\left(\frac{\kappa }{2}-\mathrm{i}\nu \right)=0,\end{equation} \tag{ A.9 }$$

the critical pumping rate at the lasing phase transition is found to be

$$\begin{equation}{{\Gamma}}_{\text{crit}}^{{\uparrow}}={{\Gamma}}_{\mathrm{T}}\left(\frac{1}{2}+\frac{{{\Gamma}}_{\mathrm{T}}\kappa }{{g}^{2}}\left[\frac{1}{8}+\frac{{({\Delta}\omega )}^{2}}{{({{\Gamma}}_{\mathrm{T}}+\kappa )}^{2}}\right]\right),\end{equation} \tag{ A.10 }$$

where Δω = ω_e − ν is the atom-cavity detuning. It is clear that the resonance condition here is that the cavity should be at the same frequency as the atomic transition ν = ω_e. At this point and in the coherent limit of g ≫ κ this reduces to ${{\Gamma}}_{\text{crit}}^{{\uparrow}}={{\Gamma}}_{\mathrm{T}}/2={{\Gamma}}^{{\downarrow}}$. As expected the optimum lasing transition occurs exactly when gain and loss are balanced. This can be clearly seen in figures A1 and A2 where we show the phase diagram as a function of detuning and the steady state properties as we go through the lasing threshold. Notice that the gain locking shown in figures A2(a) and (b) is perfect [28, 67], in contrast to that observed in the three-level system model. We see that as the cavity decay rate is lowered the results approach the limit described above.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A.1. Cumulant equations

The second order cumulants for two-level system laser model provide a useful benchmark for the results of this paper. The relevant equations for the matter terms obey the following

$$\begin{equation}\langle \dot{{\sigma }_{\mathrm{z}}}\rangle =-\frac{2g}{\sqrt{N}}\left(\mathrm{i}\langle a{\rho }_{\mathrm{e}\mathrm{g}}\rangle +\text{c.c.}\right)+{{\Gamma}}_{\mathrm{T}}\langle {\sigma }_{\mathrm{z}}\rangle +({{\Gamma}}^{{\uparrow}}-{{\Gamma}}^{{\downarrow}}),\end{equation} \tag{ A.11 }$$

$$\begin{equation}\langle \dot{{\rho }_{\mathrm{e}\mathrm{g}}{\rho }_{\mathrm{g}\mathrm{e}}}\rangle =-{{\Gamma}}_{\mathrm{T}}\langle {\rho }_{\mathrm{e}\mathrm{g}}{\rho }_{\mathrm{g}\mathrm{e}}\rangle +\mathrm{i}\frac{g}{\sqrt{N}}\left(\langle a{\rho }_{\mathrm{e}\mathrm{g}}\rangle \langle {\sigma }_{\mathrm{z}}\rangle +\text{c.c.}\right).\end{equation} \tag{ A.12 }$$

We have broken the third order moments that arise into products of first and second order moments by setting the third cumulants to zero. We choose to write the symmetry preserving versions of these equations such that we set all terms which break the U(1) symmetry of the model to 0. This approach has been shown to give good agreement with exact results since the full quantum model always respects this symmetry for finite N [60]. Here the cross correlation terms denote operators acting on different atoms i.e. $\langle {\rho }_{\mathrm{e}\mathrm{g}}{\rho }_{\mathrm{e}\mathrm{g}}\rangle \equiv \langle {\rho }_{\mathrm{e}\mathrm{g}}^{(i)}{\rho }_{\mathrm{e}\mathrm{g}}^{(j)}\rangle$ with i ≠ j. Since, within the model we consider, all atoms are identical and hence all correlations between different pairs of atoms are identical we are able to drop the explicit labelling of these. The second moments of the field-matter terms obey the following equations of motion:

$$\begin{equation}\langle \dot{a{\rho }_{\mathrm{e}\mathrm{g}}}\rangle =-\left[\mathrm{i}\left(\nu -{\omega }_{\mathrm{e}}\right)+\frac{1}{2}\left(\kappa +{{\Gamma}}_{\mathrm{T}}\right)\right]\langle a{\rho }_{\mathrm{e}\mathrm{g}}\rangle -\mathrm{i}\frac{g}{\sqrt{N}}\left[(N-1)\langle {\rho }_{\mathrm{e}\mathrm{g}}{\rho }_{\mathrm{g}\mathrm{e}}\rangle +\langle {a}^{{\dagger}}a\rangle \langle {\sigma }_{\mathrm{z}}\rangle +\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle \right].\end{equation} \tag{ A.13 }$$

Finally, the second moments of the field-field terms obey the following equations of motion:

$$\begin{equation}\langle \dot{{a}^{{\dagger}}a}\rangle =-\kappa \langle {a}^{{\dagger}}a\rangle +g\sqrt{N}(\mathrm{i}\langle a{\rho }_{\mathrm{e}\mathrm{g}}\rangle +\text{c.c.}).\end{equation} \tag{ A.14 }$$

Together these provide a closed set of equations for the second moments which can be used to find the 1/N corrections to the mean-field equations. A comparison between the behaviour described by the second order cumulants, mean-field equations, and the exact solution to the full quantum model obtained from the permutation symmetric method are shown in figures A3 and A4. These again show how the exact solution approaches the cumulant results at relatively small N. The trend of moving towards the mean-field result at large N can clearly be seen. These results should be compared to figures 6 and 7 of the main text. We see very similar behaviour in both cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Further in figure A5 we show the second order coherence and Fano factor of the steady state photon distribution as defined in equations (21) and (22) of the main text for this two-level model. We see, much more clearly than in the three-level model, that when N is large there saturation of g⁽²⁾ ≃ 2 in the normal state which crosses over to g⁽²⁾ ≃ 1 in the lasing state. The accompanying Fano factor has a significant peak around the location of the transition.

Zoom In Zoom Out Reset image size

In figure B1 we show a more complete set of results to go alongside those in figure 2 of the main text. We see from the phase diagram that changing the phase by π switches the |+⟩ and |−⟩ states and so inverts the sign of the detuning. The effects of changing the total atomic decay rate are also shown; we see that increasing Γ_T uniformly increases the threshold pumping rate as the system needs to be more strongly driven to overcome the extra losses. The qualitative features of the phase and detuning dependence are unaffected.

Zoom In Zoom Out Reset image size

We demonstrated in figure 4 of the main text how the atomic population inversion and cavity occupation change as the system is driven through the lasing threshold for various microwave drive phases ϕ. These results were all obtained by optimising the detuning and tuning the system exactly to the correct resonance condition. Here, in figure B2, instead we show similar results, but instead keep the detuning fixed to Δω = Ω. This again quantitatively changes the results but the behaviour is qualitatively the same as seen in the main text.

Zoom In Zoom Out Reset image size

We have thus far assumed that the two lower levels are degenerate. We now lift this assumption in order to show that the results presented in the main text are robust to the inclusion of this extra parameter. To do this we show the mean-field phase diagram in figure B3 and the steady state inversion and cavity population in figure B4 both at finite detuning. We see very similar qualitative behaviour to that for the case where the two-lower levels are degenerate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For completeness, in this section we give the full set of second order cumulant equations which describe the three-level system model. The procedure behind generating these is the same as that presented in appendix A.

The first set of equations are for the first moments of the atomic operators. We again only include those terms which do not break the symmetry of the model. These are just the same as the mean-field equations presented in the main text except here we do not break the second order correlations which appear

$$\begin{equation}\langle {\dot{\rho }}_{\mathrm{e}\mathrm{e}}\rangle =-\frac{g}{\sqrt{N}}\left[\mathrm{i}(\langle a{\rho }_{\mathrm{e}1}\rangle +\langle a{\rho }_{\mathrm{e}2}\rangle )+\text{c.c.}\right]+{{\Gamma}}^{{\uparrow}}(\langle {\rho }_{11}\rangle +\langle {\rho }_{22}\rangle )-2{{\Gamma}}^{{\downarrow}}\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle ,\end{equation} \tag{ C.1 }$$

$$\begin{equation}\langle \dot{{\rho }_{11}}\rangle =\left[\mathrm{i}\left(\frac{g}{\sqrt{N}}\langle a{\rho }_{\mathrm{e}1}\rangle -{\Omega}\,{\text{e}}^{-\text{i}\phi }\langle {\rho }_{12}\rangle \right)+\text{c.c.}\right]-{{\Gamma}}^{{\uparrow}}\langle {\rho }_{11}\rangle +{{\Gamma}}^{{\downarrow}}\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle ,\end{equation} \tag{ C.2 }$$

$$\begin{equation}\langle \dot{{\rho }_{22}}\rangle =\left[\mathrm{i}\left(\frac{g}{\sqrt{N}}\langle a{\rho }_{\mathrm{e}2}\rangle +{\Omega}\,{\text{e}}^{-\text{i}\phi }\langle {\rho }_{12}\rangle \right)+\text{c.c.}\right]-{{\Gamma}}^{{\uparrow}}\langle {\rho }_{22}\rangle +{{\Gamma}}^{{\downarrow}}\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle ,\end{equation} \tag{ C.3 }$$

$$\begin{equation}\langle \dot{{\rho }_{12}}\rangle =\mathrm{i}\left[\frac{g}{\sqrt{N}}\left(\langle a{\rho }_{\mathrm{e}2}\rangle -{\langle a{\rho }_{\mathrm{e}1}\rangle }^{\ast }\right)-{\Omega}\,{\text{e}}^{\text{i}\phi }(\langle {\rho }_{11}\rangle -\langle {\rho }_{22}\rangle )\right]-{{\Gamma}}^{{\uparrow}}\langle {\rho }_{12}\rangle .\end{equation} \tag{ C.4 }$$

The relevant second order terms which describe correlations between different atoms are given by

$$\begin{equation}\langle \dot{{\rho }_{1\mathrm{e}}{\rho }_{\mathrm{e}1}}\rangle =\left(\mathrm{i}\left[\frac{g}{\sqrt{N}}\langle a{\rho }_{e1}\rangle (\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle -\langle {\rho }_{11}\rangle -\langle {\rho }_{12}\rangle )-{\Omega}\,{\text{e}}^{-\text{i}\phi }\langle {\rho }_{\mathrm{e}2}{\rho }_{1\mathrm{e}}\rangle \right]+\text{c.c.}\right)-{{\Gamma}}_{\phi }\langle {\rho }_{1e}{\rho }_{e1}\rangle ,\end{equation} \tag{ C.5 }$$

$$\begin{equation}\langle \dot{{\rho }_{2\mathrm{e}}{\rho }_{\mathrm{e}2}}\rangle =\left(\mathrm{i}\left[\frac{g}{\sqrt{N}}\langle a{\rho }_{\mathrm{e}2}\rangle (\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle -\langle {\rho }_{22}\rangle -{\langle {\rho }_{12}\rangle }^{\ast })+{\Omega}\,{\text{e}}^{-\text{i}\phi }\langle {\rho }_{\mathrm{e}2}{\rho }_{1\mathrm{e}}\rangle \right]+\text{c.c.}\right)-{{\Gamma}}_{\phi }\langle {\rho }_{2\mathrm{e}}{\rho }_{\mathrm{e}2}\rangle ,\end{equation} \tag{ C.6 }$$

$$\begin{align}\hfill \langle \dot{{\rho }_{\mathrm{e}2}{\rho }_{1\mathrm{e}}}\rangle & =\mathrm{i}\frac{g}{\sqrt{N}}\left[\langle a{\rho }_{\mathrm{e}2}\rangle \left(\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle -\langle {\rho }_{11}\rangle -{\rho }_{12}\right)-{\langle a{\rho }_{\mathrm{e}1}\rangle }^{\ast }\left(\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle -\langle {\rho }_{22}\rangle -{\rho }_{12}\right)\right]\hfill \\ \hfill & \quad -\mathrm{i}{\Omega}\,{\text{e}}^{\text{i}\phi }\left(\langle {\rho }_{1\mathrm{e}}{\rho }_{\mathrm{e}1}\rangle -\langle {\rho }_{2\mathrm{e}}{\rho }_{\mathrm{e}2}\rangle \right)-{{\Gamma}}_{\phi }\langle {\rho }_{\mathrm{e}2}{\rho }_{1\mathrm{e}}\rangle ,\hfill \end{align} \tag{ C.7 }$$

where Γ_ϕ = (2Γ_↓ + Γ_↑)/2. Again here the two operators which appear act on different atoms and we have used the permutation symmetry of the model to realise that for each pair these must be identical. We have also broken the third order correlations which instead appear as products of first and second moments following the same procedure described in appendix A.

There is then a set of terms which describe the correlations between the atoms and the field. These are given by

$$\begin{align}\hfill \langle \dot{a{\rho }_{\mathrm{e}1}}\rangle & =-\mathrm{i}\frac{g}{\sqrt{N}}\left[(N-1)\left(\langle {\rho }_{1\mathrm{e}}{\rho }_{\mathrm{e}1}\rangle +{\langle {\rho }_{\mathrm{e}2}{\rho }_{1\mathrm{e}}\rangle }^{\ast }\right)+\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle \right]-\mathrm{i}{\Omega}\,{\text{e}}^{-\text{i}\phi }\langle a{\rho }_{\mathrm{e}2}\rangle \hfill \\ \hfill & \quad -\mathrm{i}\frac{g}{\sqrt{N}}\langle {a}^{{\dagger}}a\rangle (\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle -\langle {\rho }_{11}\rangle -{\langle {\rho }_{12}\rangle }^{\ast })-\left[\mathrm{i}\left(\nu -{\omega }_{\mathrm{e}}\right)+\frac{1}{2}\left(\kappa +{{\Gamma}}_{\phi }\right)\right]\langle a{\rho }_{\mathrm{e}1}\rangle ,\hfill \end{align} \tag{ C.8 }$$

$$\begin{align}\hfill \langle \dot{a{\rho }_{\mathrm{e}2}}\rangle & =-\mathrm{i}\frac{g}{\sqrt{N}}\left[(N-1)\left(\langle {\rho }_{2\mathrm{e}}{\rho }_{\mathrm{e}2}\rangle +\langle {\rho }_{\mathrm{e}2}{\rho }_{1\mathrm{e}}\rangle \right)+\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle \right]-\mathrm{i}{\Omega}\,{\text{e}}^{\text{i}\phi }\langle a{\rho }_{\mathrm{e}1}\rangle \hfill \\ \hfill & \quad -\mathrm{i}\frac{g}{\sqrt{N}}\langle {a}^{{\dagger}}a\rangle \left(\langle {\rho }_{\mathrm{e}\mathrm{e}}\rangle -\langle {\rho }_{22}\rangle -\langle {\rho }_{12}\rangle \right)-\left[\mathrm{i}\left(\nu -{\omega }_{\mathrm{e}}\right)+\frac{1}{2}\left(\kappa +{{\Gamma}}_{\phi }\right)\right]\langle a{\rho }_{\mathrm{e}2}\rangle .\hfill \end{align} \tag{ C.9 }$$

Finally, the only second moment of the field required is the photon number operator which has the equations of motion

$$\begin{equation}\langle \dot{{a}^{{\dagger}}a}\rangle =-\kappa \langle {a}^{{\dagger}}a\rangle +g\sqrt{N}\left[\mathrm{i}\left(\langle a{\rho }_{\mathrm{e}1}\rangle +\langle a{\rho }_{\mathrm{e}2}\rangle \right)+\text{c.c.}\right].\end{equation} \tag{ C.10 }$$

Together the solutions to this set of equations can be used to calculate the behaviour of any quantity of interest.