Superradiant and lasing states in driven-dissipative Dicke models

When large numbers of emitters couple to light, there can be constructive interference effects which enhance both emission and coupling to light. Such enhancement can arise either from the bosonic stimulation due to occupation of the photon mode, or from constructive interference between emission pathways in the atoms, or both. These effects have consequences both for the dynamics of an initially excited state, and for the possibility to maintain steady states with coherent matter and light. The term superradiance is used in the literature to describe a variety of these phenomena, which are related, but distinct. The original use of the term superradiance by Dicke [1] described the enhanced transient emission from a cloud of initially excited atoms when the size of the cloud is smaller than the wavelength of the emitted light. Later, it was noted that a related model, describing many atoms placed in a single mode cavity, can undergo a ground state phase transition from a normal state at weak light–matter coupling to a superradiant (SR) state with a macroscopically occupied cavity mode [2–4]. However, the existence of this phase transition in the ground state of real systems is under debate due to the effects of diamagnetic terms in the minimal coupling Hamiltonian [5–11].

In spite of the debate regarding the ground state phase transition, a related phase transition certainly is achievable using a Raman driven scheme [12, 13], which has been realised in recent experiments with ultracold atoms in optical cavities [14–17]. In such a scheme, two low-lying states of the atom can be considered to form an effective two-level system. These levels are then connected by Raman transitions, via virtually populated excited levels of the atoms. The Raman scheme involves both light in the cavity and an external pump. This means that the effective coupling between atoms and the cavity mode can be controlled by the strength of external pump, circumventing the issues with the ground state transition. This driven dissipative system is able to show a transition very similar to the ground state transition described above. These experiments prompted much theoretical investigation [18–27] into the nature of the phase transition.

There are many similarities between the steady state superradiance transition in the Raman driven context described above and the transition seen in simple models of a two-level laser [28]. Both involve a transition from a state with an empty cavity to a macroscopically occupied photon mode, and both operate due to external driving balancing cavity loss. Microscopically, the main difference between them is in the mechanism by which the driving occurs. In a two-level description of a laser the upper level is populated via an incoherent pumping process, while in the Raman driven Dicke model the driving appears via 'counter-rotating' terms in the Hamiltonian, which provide a coherent pump. It is worth noting in this context that the SR laser discussed in [29–31] corresponds to an incoherent pumping process [31], and is distinguished from a standard laser by operating deep in the bad cavity limit of cavity QED.

In this paper we look at the ways in which these coherent and incoherent driving processes can interact and lead to a rich phase diagram with regions showing differing behaviours. By studying a model which is able to show both steady-state superradiance and standard two-level lasing, we may study whether these limiting behaviours can be continuously connected, and whether other forms of lasing or superradiance exist when both driving processes are combined. We find that the phase diagram as a function of both coherent and incoherent pumping shows two distinct phases where the photon mode is macroscopically occupied. These are continuously connected to the SR phase and the lasing phase. The lasing phase can only occur in the region where the spins are inverted and the SR phase below inversion. We go on to explore the behaviour in these two regions using using a second order cumulant expansion which is valid at large (but finite) N, backed up by exact numerics based on permutation symmetry at intermediate N. We find 3 distinct parameter regions: when the coherent pumping term is small we find that the system transitions to a lasing state at large pumping. When the coherent pumping is large there is a SR state which exists only when the incoherent pump is weak. Between these two limits a crossover is seen where the system goes SR-normal-lasing as the incoherent pump power is increased.

The structure of this paper is as follows. In section 2 we present the model we use and give some motivation for how it could be experimentally realised. In section 3 we give the mean field equations and show how these can be used to calculate the stability of the normal state. Then in section 4 we go beyond these mean field equations to derive the second order cumulant equations which allow us to examine second order moments of the distribution. We then go on, in section 5 to show how these results compare to exact solutions using permutation symmetric methods, these also allow us to calculate higher moments such as the Fano factor and probability distribution for occupation of the photon mode through both the SR and lasing transitions. In section 6 we calculate the emission spectrum of the model, finding a distinct qualitative difference between the behaviours in the lasing and SR regions of the phase diagram. In section 7 we show that a deeper understanding of the relation between lasing and SR states can be reached by considering the phase diagram as a function of the frequency of the photon mode. This reveals that the model actually has four distinct phases: the normal SR and lasing phase which typically occur for positive photon frequency but also an inverted version of each state which occur when the photon is inverted. This also leads us to uncover some interesting chaotic dynamics which occur in the mean field and cumulant equations. Finally, in section 8 we present our conclusions.

We examine a model which can show both Dicke superradiance and a standard lasing transition. The model has the Hamiltonian

$$\begin{eqnarray}&&H=\omega {a}^{\dagger }a+\displaystyle \sum _{i}^{N}\displaystyle \frac{{\omega }_{0}}{2}{\sigma }_{i}^{z}+g({\sigma }_{i}^{+}a+{\sigma }_{i}^{-}{a}^{\dagger })+g^{\prime} ({\sigma }_{i}^{-}a+{\sigma }_{i}^{+}{a}^{\dagger }),\end{eqnarray} \tag{ 1 }$$

which describes a single photon mode (annihilation operator a) interacting with an ensemble of two level atoms (described by the Pauli operators ${\sigma }_{i}$). A schematic diagram showing a possible realisation of this model in a cold atom setting is shown in figure 1 along with the four level scheme which can produce this model after adiabatic elimination of the upper (detuned) levels [12].

Zoom In Zoom Out Reset image size

This scheme allows us to separate the co- and counter-rotating terms of the light–matter interaction [12] to have prefactors g and $g^{\prime}$ respectively. The master equation for the density matrix, ρ, of the system is given by:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\rho }{{\rm{d}}{t}}=-{\rm{i}}[H,\rho ]+\kappa { \mathcal D }[a]+\displaystyle \sum _{{i}}^{N}{{\rm{\Gamma }}}_{\downarrow }{ \mathcal D }[{\sigma }_{{i}}^{-}]+{{\rm{\Gamma }}}_{\uparrow }{ \mathcal D }[{\sigma }_{{i}}^{+}],\end{eqnarray} \tag{ 2 }$$

which includes photon loss at rate κ, incoherent loss of atomic excitations at rate ${{\rm{\Gamma }}}_{\downarrow }$ and an incoherent pumping process at rate ${{\rm{\Gamma }}}_{\uparrow }$ as standard Lindblad terms with ${ \mathcal D }[x]=x\rho {x}^{\dagger }-1/2\{{x}^{\dagger }x,\rho \}$. We note that both the atomic pumping and dissipation terms could in principle be engineered using the same Raman scheme as illustrated in figure 1, but not involving the cavity—i.e. spontaneous Raman processes involving real emission from excited states of the atoms. In most experiments [14–17], such terms are deliberately suppressed, by using the fact that the strength of this incoherent Raman process scales as the inverse square of the laser detuning from the atomic resonance, while the coherent terms scale only as the inverse of the detuning. In general, allowing real decay also induces additional dephasing terms ${ \mathcal D }[{\sigma }_{i}^{z}]$ which have been discussed elsewhere [25, 26]. Here, for simplicity, we neglect additional dephasing, and focus on the competition between incoherent and coherent driving processes.

It is clear that when $g=g^{\prime}$, ${{\rm{\Gamma }}}_{\uparrow }=0$ this is the regular open Dicke model which we have previously studied [26]. This model has a phase transition to a SR state which spontaneously breaks a ${{\mathbb{Z}}}_{2}$ symmetry which exists since the model is invariant under the exchange a→ − a, ${\sigma }_{i}^{x}\to -{\sigma }_{i}^{x}$. The presence of the two loss terms κ and ${{\rm{\Gamma }}}_{\downarrow }$ modify the phase transition by slightly shifting the critical point to [25–27]

$$\begin{eqnarray}&&{g}_{c}^{2}N=\displaystyle \frac{1}{4{\omega }_{0}\omega }\left({\omega }_{0}^{2}+\displaystyle \frac{{{\rm{\Gamma }}}_{\downarrow }^{2}}{4}\right)\left({\omega }^{2}+\displaystyle \frac{{\kappa }^{2}}{4}\right).\end{eqnarray} \tag{ 3 }$$

In contrast, if $g^{\prime} =0$ then this shows a regular lasing transition at a particular value of ${{\rm{\Gamma }}}_{\uparrow }$ which for small κ is when ${{\rm{\Gamma }}}_{\uparrow }\simeq {{\rm{\Gamma }}}_{\downarrow }$. This transition breaks a U(1) symmetry in which the equations are invariant under the transformation $a\to a{{\rm{e}}}^{{\rm{i}}\phi }$, ${\sigma }_{i}^{-}\to {\sigma }_{i}^{-}{{\rm{e}}}^{{\rm{i}}\phi }$.

By varying $g^{\prime}$ and ${{\rm{\Gamma }}}_{\uparrow }$, the behaviour can be continuously varied between the two cases described above. This allows us to understand how the physics crosses over between these two limits, and to identify probes to experimentally distinguish between the lasing and superradiance that both occur in this model. Our aim in the remainder of this paper will be to study this model with increasing levels of sophistication, in order to identify the phases that occur, and to characterise the behaviour in each phase.

We start our analysis by looking at the mean-field equations for the cavity and the spins. As demonstrated in previous work [26], even with dephasing and dissipation, the mean-field analysis gives a reasonable picture of the states that can arise. This decomposition requires us to break correlations between the different parts of the system at first order, resulting in the following set of equations,

$$\begin{eqnarray}&&{\partial }_{t}\langle a\rangle =-\left({\rm{i}}\omega +\displaystyle \frac{\kappa }{2}\right)\langle a\rangle -{\rm{i}}{N}(g\langle {\sigma }^{-}\rangle +g^{\prime} {\langle {\sigma }^{-}\rangle }^{* }),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\partial }_{t}\langle {\sigma }^{-}\rangle =-\left({\rm{i}}{\omega }_{0}+\displaystyle \frac{{{\rm{\Gamma }}}_{T}}{2}\right)\langle {\sigma }^{-}\rangle +{\rm{i}}(g\langle a\rangle +g^{\prime} {\langle a\rangle }^{* })\langle {\sigma }^{z}\rangle ,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\partial }_{t}\langle {\sigma }^{z}\rangle =4g\mathrm{Im}[\langle a\rangle {\langle {\sigma }^{-}\rangle }^{* }]-4g^{\prime} \mathrm{Im}[\langle a\rangle \langle {\sigma }^{-}\rangle ]-{{\rm{\Gamma }}}_{T}\langle {\sigma }^{z}\rangle +{{\rm{\Gamma }}}_{\uparrow }-{{\rm{\Gamma }}}_{\downarrow },\end{eqnarray} \tag{ 6 }$$

where ${{\rm{\Gamma }}}_{T}=({{\rm{\Gamma }}}_{\downarrow }+{{\rm{\Gamma }}}_{\uparrow })$. When $g^{\prime} =0$ these are just the Maxwell–Bloch equations for a laser [32].

The normal state has ${\langle a\rangle }_{{\rm{ns}}}={\langle {\sigma }^{-}\rangle }_{{\rm{ns}}}=0$ and ${\langle {\sigma }^{z}\rangle }_{{\rm{ns}}}=({{\rm{\Gamma }}}_{\uparrow }-{{\rm{\Gamma }}}_{\downarrow })/{{\rm{\Gamma }}}_{T}$. Linearising around this solution with $\delta a=\langle a\rangle -{\langle a\rangle }_{{\rm{ns}}}$, $\delta s=\langle {\sigma }^{-}\rangle -{\langle {\sigma }^{-}\rangle }_{{\rm{ns}}}$ etc gives the following matrix

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\left(\begin{array}{c}\delta a\\ \delta {a}^{* }\\ \delta s\\ \delta {s}^{* }\end{array}\right)=\left(\begin{array}{cccc}-{\rm{i}}\omega -\kappa /2 & 0 & -{\rm{i}}{g} & -{\rm{i}}{g}^{\prime} \\ 0 & {\rm{i}}\omega -\kappa /2 & {\rm{i}}{g}^{\prime} & {\rm{i}}{g}\\ {\rm{i}}{g}{\langle {\sigma }^{z}\rangle }_{{\rm{ns}}} & {\rm{i}}{g}^{\prime} {\langle {\sigma }^{z}\rangle }_{{\rm{ns}}} & -{\rm{i}}{\omega }_{0}-{{\rm{\Gamma }}}_{T}/2 & 0\\ -{\rm{i}}{g}^{\prime} {\langle {\sigma }^{z}\rangle }_{{\rm{ns}}} & -{\rm{i}}{g}{\langle {\sigma }^{z}\rangle }_{{\rm{ns}}} & 0 & {\rm{i}}{\omega }_{0}-{{\rm{\Gamma }}}_{T}/2\end{array}\right)\left(\begin{array}{c}\delta a\\ \delta {a}^{* }\\ \delta s\\ \delta {s}^{* }\end{array}\right).\end{eqnarray} \tag{ 7 }$$

This 4 × 4 matrix equation arises since $\langle a\rangle$ couples both $\langle {\sigma }^{+}\rangle$ and its conjugate $\langle {\sigma }^{-}\rangle$. We may ignore the equation of motion for $\langle {\sigma }^{z}\rangle$ since this only couples to itself in the linearised system, and so does not contribute to the stability analysis.

When the real part of (at least) one of the eigenvalues of this matrix becomes positive the normal state is unstable to perturbations and the photon mode acquires a macroscopic occupation. This linear stability approach allows us to very quickly explore large regions of parameter space. In figure 2 we explore this phase diagram as a function of both the coherent pumping term $g^{\prime}$ and the incoherent pumping term ${{\rm{\Gamma }}}_{\uparrow }$. In order to show the full range of ${{\rm{\Gamma }}}_{\uparrow }/{{\rm{\Gamma }}}_{\downarrow }$ (i.e. from 0 to $\infty$) we vary the incoherent pump while keeping the total decay of the spins, ${{\rm{\Gamma }}}_{T}$, constant, so that ${{\rm{\Gamma }}}_{\uparrow }/{{\rm{\Gamma }}}_{T}=1$ corresponds to ${{\rm{\Gamma }}}_{\uparrow }/{{\rm{\Gamma }}}_{\downarrow }=\infty$.

Zoom In Zoom Out Reset image size

We see that there are two distinct regions in which the normal state is unstable, one, at small values of ${{\rm{\Gamma }}}_{\uparrow }$, is continuously connected to the SR region while the other, at small $g^{\prime}$, is connected to the lasing transition. We label these regions 'SR' and 'laser' respectively on the phase diagram. We see that, even when the photon loss rate is very small, these two regions are always separated by a (possibly small) normal region which is always present exactly at the point where the spins become inverted ${{\rm{\Gamma }}}_{\uparrow }={{\rm{\Gamma }}}_{\downarrow }$. The main change as κ is reduced is to the shape of the lasing region which we see move towards the inversion line as the losses are reduced. We note that the lasing region can never cross this line from above while the SR region never crosses it from below.

In figure 2, for the regions where the normal state is unstable, we also indicate whether there exists only a single unstable mode, or a pair of unstable modes. In the latter case, this corresponds to a complex conjugate pair of eigenvalues, and describes a mode which oscillates as well as becoming unstable. As we will discuss further below, this is exactly as expected for the laser, as the frequency of the lasing mode generally matches the cavity frequency, and so the lasing solution is expected to be time dependent (oscillatory) in the frame we are working in. This means that the instability to lasing should involve an eigenvalue with an imaginary part, corresponding to the lasing frequency. Since eigenvalues come in complex conjugate pairs, the instability in this case is a Hopf bifurcation, where a complex conjugate pair of eigenvalues simultaneously become unstable. In contrast, the 'standard' SR state is expected to be stationary, and the instability to it thus corresponds to a mode which grows without oscillating, thus a single mode is unstable at the phase boundary. The phase diagram observed is compatible with this observation, however more complicated behaviour will be seen below.

While mean-field theory is useful to determine the structure of the phase diagram, it can only capture the behaviour in the limit $N\to \infty$. As such, it is not always possible to directly compare mean-field theory results to exact results at finite N. We can go beyond mean-field theory by finding the equations for all the second moments of the photon and atomic distributions. This approach allows us to describe quantities such as the photon number accurately and without the need to explicitly break symmetries in the initial conditions. Indeed, in writing these equations we will choose to keep only those terms that respect the symmetries in the problem. This approach is directly analogous to the semiclassical approach used in laser theory [32], which considers equations of motion for the photon number, while incorporating both spontaneous and stimulated processes.

The second moments of the photons obey the equations of motion

$$\begin{eqnarray}&&{\partial }_{t}\langle {a}^{\dagger }a\rangle =-\kappa \langle {a}^{\dagger }a\rangle -2N(g\mathrm{Im}[{C}^{a+}]+g^{\prime} \mathrm{Im}[{C}^{a-}]),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\partial }_{t}\langle {aa}\rangle =-(2{\rm{i}}\omega +\kappa )\langle {aa}\rangle -2{\rm{i}}{N}({{gC}}^{a-}+g^{\prime} {C}^{a+}),\end{eqnarray} \tag{ 9 }$$

while the photon matter correlations are

$$\begin{eqnarray}{\partial }_{t}{C}^{a+} & = & -\left({\rm{i}}\omega +\displaystyle \frac{\kappa }{2}+\tilde{{\rm{\Gamma }}}\right){C}^{a+}+{\rm{i}}{\omega }_{0}{C}^{a+}\\ & & -{\rm{i}}{g}[(N-1){C}^{+-}+\langle {\sigma }^{z}\rangle \langle {a}^{\dagger }a\rangle +1/2(1+\langle {\sigma }^{z}\rangle )]\\ & & -{\rm{i}}{g}^{\prime} [(N-1){C}^{++}+\langle {\sigma }^{z}\rangle \langle {aa}\rangle ],\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}{\partial }_{t}{C}^{a-} & = & -\left({\rm{i}}\omega +\displaystyle \frac{\kappa }{2}+\tilde{{\rm{\Gamma }}}\right){C}^{a-}-{\rm{i}}{\omega }_{0}{C}^{a-}\\ & & -{\rm{i}}{g}[(N-1){({C}^{++})}^{* }-\langle {\sigma }^{z}\rangle \langle {aa}\rangle ]\\ & & -{\rm{i}}{g}^{\prime} [(N-1){C}^{+-}-\langle {\sigma }^{z}\rangle \langle {a}^{\dagger }a\rangle +1/2(1-\langle {\sigma }^{z}\rangle )],\end{eqnarray} \tag{ 11 }$$

here e.g. ${C}^{a+}=\langle a{\sigma }_{i}^{+}\rangle$. There are also matter-matter correlations given by

$$\begin{eqnarray}&&{\partial }_{t}{C}^{++}=2{\rm{i}}{\omega }_{0}{C}^{++}-2\tilde{{\rm{\Gamma }}}{C}^{++}-2{\rm{i}}\langle {\sigma }^{z}\rangle (g{({C}^{a-})}^{* }+g^{\prime} {C}^{a+}),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}{\partial }_{t}{C}^{{zz}} & = & 8\langle {\sigma }^{z}\rangle (g\mathrm{Im}[{C}^{a+}]-g^{\prime} \mathrm{Im}[{C}^{a-}])\\ & & -2{{\rm{\Gamma }}}_{\downarrow }({C}^{{zz}}+\langle {\sigma }^{z}\rangle )-2{{\rm{\Gamma }}}_{\uparrow }({C}^{{zz}}-\langle {\sigma }^{z}\rangle ),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\partial }_{t}{C}^{+-}=-2\tilde{{\rm{\Gamma }}}{C}^{+-}-2\langle {\sigma }^{z}\rangle (g\mathrm{Im}[{C}^{a+}]-g^{\prime} \mathrm{Im}[{C}^{a-}]).\end{eqnarray} \tag{ 14 }$$

where e.g. ${C}^{+-}=\langle {\sigma }_{i}^{+}{\sigma }_{j\ne i}^{-}\rangle$ gives the correlation between ${\sigma }^{+}$ at site i and ${\sigma }^{-}$ at a different site j. We also need the equation for $\langle {\sigma }^{z}\rangle$, which does not break the symmetry, and takes a form similar to that in equation (6) but without breaking the second order correlations,

$$\begin{eqnarray}&&{\partial }_{t}\langle {\sigma }^{z}\rangle =4g\mathrm{Im}[{C}^{a+}]-4g^{\prime} \mathrm{Im}[{C}^{a-}]-{{\rm{\Gamma }}}_{T}\langle {\sigma }^{z}\rangle +{{\rm{\Gamma }}}_{\uparrow }-{{\rm{\Gamma }}}_{\downarrow }.\end{eqnarray} \tag{ 15 }$$

These reproduce the equations in [26] when $g^{\prime} =g$ and ${{\rm{\Gamma }}}_{\uparrow }=0$.

In figure 3 we show how these equations can be used to examine the phase diagram which was plotted in figure 2. We may check that the region with a macroscopic photon number is in very good agreement with that predicted by mean field theory, with the SR dome at large $g^{\prime}$ and small ${{\rm{\Gamma }}}_{\uparrow }$ and the lasing region at large ${{\rm{\Gamma }}}_{\uparrow }$ and small $g^{\prime}$. By looking at the spin state we see directly that the SR region is where the spins are not inverted and the lasing region corresponds to spin inversion. The vertical white region in figure 3(b) shows exactly where this spin inversion occurs. This is always at the point where ${{\rm{\Gamma }}}_{\uparrow }={{\rm{\Gamma }}}_{\downarrow }$ (at least up to corrections of order $\kappa /N$) since here the photons have no occupation and so the dynamics are purely determined by the driven atoms.

Zoom In Zoom Out Reset image size

To find exact numerical solutions of this system we can make use of the permutation symmetry of the individual density matrix elements. This allows us to find exact solutions at intermediate N. A full description of the method we use can be found in [26] while the code can be found at [33]. Another library which implements the same algorithm has also recently been released [34]. Similar techniques have been employed to study spin ensembles [35], simple lasing models [36], coherent surface plasmons [37], the competition between collective and individual decay channels [38], the behaviour of an ensemble of Rydberg polaritons [39], equilibrium properties of a model with a larger local Hilbert space [40], subradiant states in the Dicke model [41] and to explore the effect of individual losses on transient SR emission [42]. In this section, we use these exact solutions both to check the validity of the cumulant method presented above and calculate higher order moments of e.g. the photon distribution which are neglected in the cumulant calculation. Indeed, from the exact results we can calculate not only moments, but the full probability distribution of the photon number.

Results comparing exact numerics to the cumulant expansion are shown in figure 4. The top row shows a sweep through the lasing region (below the SR dome), the middle column goes through both the SR and lasing regions while the top row shows a sweep above the lasing region. The first two columns show the reduced photon number $\langle {a}^{\dagger }a\rangle /N$ the spin inversion level $\langle {\sigma }^{z}\rangle$, which allow comparison to the mean-field and cumulant approximations. We see that the match between the cumulant equations and the exact solution is good across all parameters, even at the relatively small value of N = 25 shown here which is far from the thermodynamic limit, with only some small deviations at large pump strengths.

Zoom In Zoom Out Reset image size

The third column of figure 4 shows the Fano factor of the photon distribution,

$$\begin{eqnarray}&&F=\displaystyle \frac{\langle {a}^{\dagger }{{aa}}^{\dagger }a\rangle -{\langle {a}^{\dagger }a\rangle }^{2}}{\langle {a}^{\dagger }a\rangle }.\end{eqnarray} \tag{ 16 }$$

Calculating this quantity requires evaluating moments beyond second order, so is not accessible from our mean-field or cumulant approaches. The Fano factor is one possible measure of the dispersion (spread) of the photon number distribution. For any coherent (Possonian) state the Fano factor is 1. This means that in the thermodynamic limit of an ideal laser the Fano factor of the lasing state is unity both above and below threshold. Directly at threshold the fluctuations give an infinite value for F [43]. However, it is known that finite size effects and finite non-radiative loss (i.e. finite β factor) smear out the peak in the Fano factor so that it can become a very broad feature [43]. This is the effect we see here, where in the lasing phase the Fano factor is elevated far above unity by a very broad feature. The peak associated with the transition to the SR state is however much narrower. This can be clearly seen in the Fano factor at $g^{\prime} =g$ where the structure is made up from these two peaks. This also explains why the cumulant expansion is less accurate in the lasing phase than it is in the SR phase: the cumulant expansion only describes the behaviour of approximately Gaussian states, which the lasing phase we see does not satisfy.

In figure 5 we show the full probability distribution of the photon mode, P(n), for the same sets of parameters as in figure 4, i.e. with N = 25. This allows us to see more clearly exactly what happens to the photon distribution as we go through the various thresholds. In the lasing regime as the pump power is increased the distribution moves to larger and larger n. The Fano factor does not drop down to 1 in this regime since for the parameters chosen here the peak of P(n) only reaches n = 4 at the largest pump strength. The opposite effect is seen in the SR regime where increasing pumping kills the coherent state. In the crossover regime we see both effects: at small pumping the SR phase is present which is killed by the incoherent drive, but at larger pumping a coherent state again appears as the system undergoes a transition to a lasing state.

Zoom In Zoom Out Reset image size

The emission spectra of light–matter systems can reveal subtle features of the dynamics not available from the steady state behaviour. A well known example of this is the Mollow triplet, seen in the fluorescence from a resonantly driven two-level system [28]. By examining the spectrum of the light which escapes the cavity we may be able to see the difference between the SR and lasing phases.

The simplest quantity to look at, accessible from both our exact numerics and the cumulant expansion is the emission spectrum which is the Fourier transform of the ${g}^{(1)}(t)$ correlation function,

$$\begin{eqnarray}&&S(\nu )={\int }_{-\infty }^{\infty }\langle {a}^{\dagger }(t)a(0)\rangle {{\rm{e}}}^{{\rm{i}}\nu t}\,{\rm{d}}{t}.\end{eqnarray} \tag{ 17 }$$

We can obtain this quantity from our exact numerics using the quantum regression theorem. We find the steady state ${\rho }_{{\rm{ss}}}$ as before and then initialise using the state $a{\rho }_{{\rm{ss}}}$. The correlation function is then given by the time evolution of ${a}^{\dagger }$ from this state [28]. The spectrum is then obtained by taking a Fourier transform.

In the large N limit we may also use the second order cumulant equations and quantum regression theorem to calculate the same quantity. To do this we must solve the matrix equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\bf{C}}(t)={ \mathcal M }{\bf{C}}(t),\end{eqnarray} \tag{ 18 }$$

where the coupled set of correlation functions required are

$$\begin{eqnarray}&&{\bf{C}}(t)=\left(\begin{array}{c}\langle {a}^{\dagger }(t)a(0)\rangle \\ \langle a(t)a(0)\rangle \\ \langle {\sigma }^{+}(t)a(0)\rangle \\ \langle {\sigma }^{-}(t)a(0)\rangle \end{array}\right),\end{eqnarray} \tag{ 19 }$$

and the evolution matrix is given by

$$\begin{eqnarray}{ \mathcal M }=\left(\begin{array}{cccc}{\rm{i}}\omega -\kappa /2 & 0 & {\rm{i}}{g} & {\rm{i}}{g}^{\prime} \\ 0 & -{\rm{i}}\omega -\kappa /2 & -{\rm{i}}{g}^{\prime} & -{\rm{i}}{g}\\ -{\rm{i}}{g}{\langle {\sigma }^{z}\rangle }_{{\rm{ss}}} & -{\rm{i}}{g}^{\prime} {\langle {\sigma }^{z}\rangle }_{{\rm{ss}}} & {\rm{i}}{\omega }_{0}-{{\rm{\Gamma }}}_{T}/2 & 0\\ {\rm{i}}{g}^{\prime} {\langle {\sigma }^{z}\rangle }_{{\rm{ss}}} & {\rm{i}}{g}{\langle {\sigma }^{z}\rangle }_{{\rm{ss}}} & 0 & -{\rm{i}}{\omega }_{0}-{{\rm{\Gamma }}}_{T}/2\end{array}\right).\end{eqnarray} \tag{ 20 }$$

A similar approach was taken in [39] to calculate correlation functions of an ensemble of Rydberg polaritons. Here ${\langle {\sigma }^{z}\rangle }_{{\rm{ss}}}$ is the steady state value of the spin inversion and the initial condition for this problem is given by the relevant quantities from the second moment equations.

These two methods then allow us to calculate the spectrum as shown in figure 6. The simplest interpretation can be found by looking at the results of the cumulant expansion at N = 500 (the right-hand column in figure 6). When the counter-rotating term is weak (i.e. $g^{\prime} =0.1g$), and the incoherent pump is strong, the spectrum has a single peak at the cavity frequency which gets more intense as the pumping is increased. This monochromatic emission is typical of a simple laser. We may also note that the emergence of this non-zero single frequency in the lasing state is directly related to the Hopf bifurcation that leads to the lasing state. In the crossover region, $g^{\prime} =g$ at low pump powers (in the SR dome) we see a triplet like structure, reminiscent of the Mollow triplet, with the sideband at a frequency which is not simply set by the cavity. As the incoherent pumping is increased the spectrum changes to look more that seen for the case of the laser. At large $g^{\prime} =3g$ (above the lasing region in figure 2) we see a very similar structure but the residual peak at the cavity frequency is much weaker. This peak is still present since the cumulant calculations are done at finite N and the parameters chosen are relatively close to the lasing region of the phase diagram. The N = 25 results form the exact numerics in the left-hand column of figure 6 show a similar qualitative behaviour but the relatively small value of N means that the features are less easy to distinguish.

Zoom In Zoom Out Reset image size

In contrast to the steady state measurements discussed in previous sections, we see here that the emission spectrum of the cavity provides an easy way to identify the differences between the SR and lasing regimes. This difference can be particularly clearly seen in the spectrum as a function of cavity frequency ω. In figure 7 we show this calculated both in the lasing phase and the SR phase. In the laser the location of the peak is approximately linear in the cavity frequency consistent with the interpretation of this being a weak light–matter coupling effect, while in the SR phase the sideband has a more complex dependence on the detuning.

Zoom In Zoom Out Reset image size

In the model we consider, the effective cavity frequency, ω, is actually the detuning between the cavity and external pump frequency (see the level diagram in figure 1). As such, it is reasonable to ask how the behaviour changes when this detuning is negative, i.e. when the pump is blue detuned from the cavity mode. From previous work, we know that if only photon loss is present then the normal state in this parameter range is inverted [21]. In this section we therefore discuss what happens to this inverted regime when there is incoherent pumping and decay of the spins.

In figure 8 we show the phase diagram as a function of both the incoherent pumping rate and photon detuning for various values of $g^{\prime}$. If $g^{\prime} \ll g$ then the only phase that is present is the normal lasing state which appears at large ${{\rm{\Gamma }}}_{\uparrow }$ and for resonant frequencies $\omega \approx {\omega }_{0}$. In the middle column of figure 8 we see that when $g=g^{\prime}$ more phases appear, both above and below inversion. When $\omega \gt 0$ we see the normal SR and lasing phases at small and large ${{\rm{\Gamma }}}_{\uparrow }$ respectively. At negative ω we see two very similar phases but at opposite sides of the phase diagram. To understand these we note that the model has a duality under the replacements $\omega \leftrightarrow -\omega$, $g\leftrightarrow g^{\prime}$, ${{\rm{\Gamma }}}_{\uparrow }\leftrightarrow {{\rm{\Gamma }}}_{\downarrow }$: such a transformation leaves the steady states unchanged up to an inversion of ${\sigma }^{z}$. This then means that the two phases seen at negative ω are the same as those when ω is positive but with photon creation and annihilation swapped, i.e. the inverted lasing phase at small ${{\rm{\Gamma }}}_{\uparrow }$ is where the dominant process is photon and spin excitation creation via the term $g^{\prime} {a}^{\dagger }{\sigma }^{+}$ in the Hamiltonian. This state is similar to one discussed in [44, 45] where a parametric drive is used to engineer similar behaviour.

Zoom In Zoom Out Reset image size

In this region of parameter space we see some artefacts in the results found by time evolving the cumulant expansion (most visible in the spin dynamics). These artefacts reflect the fact the photon number and spin inversion do not reach a steady state, but instead one has limit cycles and chaotic dynamics. We discuss this behaviour more below. As $g^{\prime}$ is increased further we see that the phases below inversion—the normal SR and inverted lasing phases—coalesce, while the two phases at large ${{\rm{\Gamma }}}_{\uparrow }$ are suppressed. Indeed, by the point $g^{\prime} =3g$, shown in the right-hand column of figure 8, the standard lasing phase has completely vanished. Although the SR and lasing phases join continuously to each other at $g^{\prime} =3g$, one can nonetheless distinguish two distinct behaviours on crossing the phase boundary: it remains the case that one can distinguish whether a single eigenvalue or a pair of eigenvalues becomes unstable. At this largest value of $g^{\prime}$, the absence of a steady state photon number extends over a wider portion of the phase diagram at negative ω.

Careful comparison of the mean field and cumulant expansion phase diagrams reveal that close to the upper boundary of the inverted SR phase there is a region where mean-field stability analysis shows that the normal state is stable and yet the cumulant equations show a macroscopic photon population. The equations in this regime are bistable: while the normal state is stable, there is also a stable SR solution to the equations. Hence we find that in this region of the phase diagram the results of the cumulant equations depend on the initial conditions.

7.1. Evolution of chaotic attractors

To characterise the non-steady-state behaviour in more detail we look at the dynamics of the mean field equations as the frequency of the photon mode is changed. These results are shown in figure 9. We plot the dynamics of the photon mode using an initial condition with finite $\langle a\rangle$ to break the symmetry (the normal state with $\langle a\rangle =0$ is always a solution to the MF equations, but in the SR phase this solution is unstable), and show the behaviour after waiting sufficiently long that any transient behaviour has vanished. To plot the bifurcation diagram [46] in the centre of figure 9 we plot the value of $\mathrm{Im}\langle a\rangle$ each time the photon crosses the $\mathrm{Re}\langle a\rangle =0$ line from right to left.

Zoom In Zoom Out Reset image size

In the region where ω is closest to zero we see that the dynamics are regular, the photon amplitude undergoes a simple periodic limit cycle which results in a single point in the bifurcation diagram. As the frequency is decreased (away from 0) the system goes through a sequence of transitions to more complex (but still regular) orbits which eventually lead to the chaotic region at around $\omega =-2.0$. At larger negative detunings again we find that the dynamics become regular. At positive detunings, e.g. $\omega =0.33$, the limit cycle transitions towards a fixed point: as we approach this, the limit cycle distorts, and the time evolution pauses near the points that ultimately become the fixed point values. These fixed points are indicated by crosses in this panel. The occurrence of chaotic behaviour within the lasing region is not too surprising since it is known that in certain limits the Maxwell-Bloch equations can be mapped to the Lorentz equations [47].

We have examined in detail the steady state phase diagram of a model which is able to continuously cross over from standard two-level lasing to Dicke superradiance by changing the balance between coherent and incoherent pumping processes. We have seen that when the photon energy is positive these two phases are distinct, the lasing phase only exists when the spins are inverted and the SR phase when the spins are not inverted.

By using a combination of approaches: mean-field theory valid in the thermodynamic limit, a cumulant expansion which gives access to the large N asymptotics and exact numerics which reach intermediate N we have shown that while the steady state photon number is similar in the two phases the spectral properties of the emission from the cavity are very different.

The phase diagram as a function of the effective photon energy ω reveals two more phases when $\omega \lt 0$. An inverted lasing phase at weak incoherent pump strengths and an inverted SR phase at large incoherent pump. This also revealed the evolution of the dynamics through a period doubling route to a chaotic attractor.

Future studies could examine the role of the nonlinearity which arises in the Hamiltonian due to the ac Stark shift of the cavity mode [19, 21, 27]. This is known to generate limit cycles in the mean field dynamics, it would be interesting to see what effect the incoherent processes considered here have on these phases.

PK acknowledges support from EPSRC (EP/M010910/1). JK acknowledges support from EPSRC programs TOPNES (EP/I031014/1) and 'Hybrid Polaritonics' (EP/M025330/1). The research data supporting this publication can be accessed at http://dx.doi.org/10.17630/24ca8410-1db7-4f73-9bdd-105dc55eafe7.