Few emitters in a cavity: From cooperative emission to individualization

We study the temporal correlations of the field emitted by an electromagnetic resonator coupled to a mesoscopic number of two-level emitters that are incoherently pumped by a weak external drive. We solve the master equation of the system for increasing number of emitters and as a function of the cavity quality factor, and we identify three main regimes characterized by well distinguished statistical properties of the emitted radiation. For small cavity decay rate, the emission events are uncorrelated and the number of photons in the emitted field becomes larger than one, resembling the build-up of a laser field inside the cavity. At intermediate decay rates (as compared to the emitter-cavity coupling) and for few emitters, the statistics of the emitted radiation is bunched and strikingly dependent on the parity of the number of emitters. The latter property is related to the cooperativity of the emitters mediated by their coupling to the cavity mode, and its connection with steady state subradiance is discussed. Finally, in the bad cavity regime the typical situation of emission from a collection of individual emitters is recovered. We also analyze how the cooperative behavior evolves as a function of pure dephasing, which allows to recover the case of a classical source made of an ensemble of independent emitters, similar to what is obtained for a very leaky cavity. State-of-art techniques of Q-switch of resonant cavities, allied with the recent capability to tune single emitters in and out of resonance, suggest this system as a versatile source of different quantum states of light.


I. INTRODUCTION
The interaction of an ensemble of identical two-level emitters with the electromagnetic field has been the subject of intense theoretical and experimental studies since the pioneering work of Dicke [1]. Among the others, superradiance (or superfluorescence) is the most celebrated effect; it is manifested as a transient effect by the delayed emission of an intense burst of light, after the atomic population has been totally inverted [2]. In essence, superradiance is characterized by the emergence of collective coherence, which leads to a strongly enhanced emission rate as compared to the case of uncorrelated emitters. Such striking consequences of the quantum interference effect in radiation-matter coupled systems has been the subject of an intense experimental interest both for atomic [2,3] and solid state [4] ensembles. Early studies on the statistical properties of the emitted light have evidenced the role of quantum fluctuations in the initial stages of the superradiant emission [5][6][7]. Whereas superradiance usually appears in the transient regime, recent theoretical studies have evidenced other signatures of collective behavior in the steady state regime, e.g. when the two-level emitters are continuously excited by incoherent light [8][9][10][11]. In these references, cooperative emission is obtained by coupling the emitters to a low quality (Q) cavity mode, which behaves as a collective channel of losses. As a matter of fact, strong confinement enhances the interaction of each emitter with the cavity field, eventually making it predominant with respect to individual relaxation channels. Collective dynamics can be induced this way, provided the atoms identically interact with the cavity mode, a condition much less stringent than in free space. Tuning the pumping power allows one to explore different regimes, leading to different intensities and statistics for the emitted radiation [10,11]. In particular, it has been shown that when the pumping strength is too low to saturate a single emitter of the ensemble, the collection of two-level systems emits a highly bunched field [9,11]. The strong temporal correlation has been attributed to the relaxation of the system in the subradiant states of the ensemble of emitters. These states behave as sources of photon pairs triggered by the driving pulses, hence the bunching effect in the correlation properties.
In this paper, we analyze the subradiant regime and go beyond previous works by keeping the cavity mode as a relevant degree of freedom. This situation is new with respect to previous theoretical studies, where the radiation mode was commonly eliminated by an adiabatic approximation. We show that the bunching degree decays with the number of emitters and oscillates with its parity, being larger when the number of emitters is even, upon converging to a single value as the number of emitters grows. We give an interpretation of this effect, based on the structure of Hilbert space of the system. Furthermore, keeping the cavity mode in the description allows to explore two other limits, depending on the cavity finesse. If this parameter is significantly increased, the strong coupling conditions lead to the emission of a field showing no temporal correlations, an effect that we iden-FIG. 1: Schematic of the system under consideration in this work: N two-level emitters incoherently pumped by an external drive, coupled to a single resonator mode. Correlations in the emitted field can be detected through a HBT set-up either from the cavity or from the emitters' loss channel, respectively. A typical incoherent pumping scheme is shown, in which high-energy excitations relax down to a well defined ground state transition (e.g. through phonon scattering) via photon emission.
tify with the emission of Poissonian radiation similar to a lasing transition. Recently, such regime has been investigated from the spectral point of view [12,13], while here we address its statistical properties. On the contrary, for decreasing cavity finesse one progressively decreases the coupling with the mode, allowing to continuously explore the boundary between collective and individual behavior of such radiation sources. When all the emitters become uncorrelated, increasing their number one by one up to large values coincides essentially with building a classical source with chaotic statistical properties. Dipole correlation can also be destroyed with pure dephasing, a situation that we eventually explore as well.
The present work becomes particularly relevant when considering that the subradiant regime is reachable if and only if each single emitter of the ensemble is predominantly coupled to the cavity mode. In particular, this allows to pump each emitter very slowly, but still faster than any individual relaxation mechanism. This is nowadays possible thanks to impressive progress in nano-fabrication technology on single emitter-cavity coupling [14][15][16][17][18][19][20][21][22][23][24][25][26][27]. Small assemblies of artificial quantum emitters, like quantum dots [14] and superconducting qubits [15] coupled to a high-Q resonator mode have recently become available, strongly motivating to investigate the new physics provided by these systems. Moreover, beyond introducing an external parameter to control the steady state collective regime, this model is of great interest also for the study of lasing transition, as well as to investigate the connections between superradi-ance and lasing [28].
The paper is organized as follows. We first introduce the model and the methods employed to solve it, both numerically and analytically (Sec. II). Then, we present results as a function of cavity finesse, and we give a physical interpretation for the different regimes found (Secs. III, IV, and V). Before summarizing the main conclusions, we give a brief description of possible realistic implementations to explore the physics investigated in this work (Sec. VI).

II. MODEL AND METHODS
The system under study consists of N two-level emitters coupled to the mode of an electromagnetic cavity of annihilation (creation) operator a (a † ) with equal coupling strength, g, and it is schematically pictured in Fig. 1. To focus on collective effects, we will assume that the resonator and all the emitters have the same transition frequency, ω. The range of validity of this ideal model in practical implementations is discussed in Sec. VI. The model can be rigorously described by the Tavis-Cummings Hamiltonian (in the rotating wave approximation) where J + = i σ i + is the collective dipole operator, and σ i − = |g i e i | is the lowering operator of the single emitter i, whose excited (ground) state is denoted by |e i (|g i ). The Hilbert space corresponding to the assembly of resonant atoms is equivalent to the tensor product of the 2-dimensional space of each single emitter, i.e. H (C 2 ) ⊗N , spanned by the basis of factorized states in the form |i 1 |i 2 ...|i n , where i n = {g n , e n } (in the following we shall refer to this basis as the uncoupled basis). When including possible dissipation channels and external driving, the Markovian dynamics of the full system is given bẏ The Lindblad dissipative terms, modeling incoherent loss and pumping of each emitter, explicitly read where the quantity γ is the individual emission rate in loss channels distinct from the cavity mode for each emitter, and P x quantifies the excitation rate. As pictured in Fig. 1, incoherent pumping can be obtained by resonantly exciting higher energy states, that further relax down to the levels of interest. The decay of the cavity mode with rate κ through imperfect mirrors is described by the following term whereas the pure dephasing of the emitters, which is unavoidable in solid-state artificial atoms, can be taken into account by a Liouvillian term of the form [29] Following standard input-output theory [30], the output field at the photodetectors is proportional to the a operator (and its conjugate). We are interested in the temporal correlations of the field radiated by the cavity mode. Experimentally, an ordinary Hanbury-Brown Twiss (HBT) setup can be used, giving access to the steady state second-order correlation at zero-time delay [31] g (2) In order to calculate the latter quantity, we numerically solve Eq. (2) either by direct integration or by Quantum Monte Carlo (QMC) simulations [32,33]. Direct integration allows one to obtain accurate results at lower computational time costs, for which we use an extension of the code already employed in a previous work [34]. However, the tensor structure of the phase space makes the computational resources needed for a direct integration rapidly diverge. On the contrary (wave function) QMC gives the solution as a stochastic average over quantum trajectories of wavefunctions unravelings, drastically reducing the amount of space needed and allowing us to push forward the previous computational limits [33]. If κ is large enough, the cavity mode can be adiabatically eliminated from the equations, giving rise to an effective dynamics described by the reduced master equationρ where The physical meaning of this equation is the following: the cavity acts as a common relaxation channel for the atoms, whose dissipation rate is quantified by the coupling strength with the collective atomic dipole J − , i.e. Γ = 4g 2 /κ. Namely, adiabatic elimination is valid when κ N Γ [10]. In this regime, the cavity field, and thus the output field is proportional to the collective dipole degree of freedom, J − . Consequently, the temporal autocorrelation function at zero time delay can be expressed as Note that the same detection properties could be obtained by registering the far field where the atomic dipoles interfere, as it is pictured in Fig. 1. Also for Eq. (9) we can employ numerical approaches either based on direct integration or on QMC, respectively. As usual, once the master equation is solved for ρ the expectation values can be calculated with the general rule Ô = Tr{Ôρ}, whereÔ identifies the operator corresponding to a given observable.
Sweeping the cavity finesse allows to continuously explore the transition from a collective behavior driven by L Γ (ρ) to an individual one, corresponding to L i (ρ). In this work we analytically study some limits of these regimes. If the dynamics is individual, it is naturally described in the uncoupled basis mentioned above. If the dynamics is collective, an order appears in the system due to the build up of correlations between dipoles belonging to different atoms [2], manifested by the occurrence of coherences in the uncoupled basis. In the latter case, a natural basis is provided by the eigenstates |J, M of the total angular momentum {J 2 , J z } (in the following, we shall talk about the coupled basis of the atomic states). Note that among them, the states satisfying J − |J, M = 0 are not coupled to the electromagnetic field and are thus defined as "subradiant".

III. SECOND-ORDER CORRELATIONS OF THE EMITTED FIELD
As mentioned above, the regime studied in the literature so far is the adiabatic elimination regime, where the cavity mode operator can be assumed as a ∝ J − . We assume a fixed and weak incoherent pumping rate, P x /g. In order to understand the collective behavior and to link it to the previous literature, we plot in Fig. 2a the average cavity mode occupation, n a = a † a , the average collective occupation of the emitters, n J = J + J − , and the total emitters population, N n x (where n x = σ i + σ i − in the steady state. The second-order correlations are plotted in Fig. 2b, and they have been calculated both from the cavity mode operator, Eq. (6), and the collective dipole operator, Eq. (9), respectively. The results are shown for N = 5 emitters in the cavity, which represents a good compromise to simultaneously have a well converged numerical solution for Eq. (6) and to evidence the effects of collective behavior. As it can be seen in Fig.  2b, the two calculated correlations don't match for most of the parameter space, showing that the cavity mode has to be taken into account to obtain the correct result when the cavity finesse is large. The two curves converge to the same value when κ/g 50, where the adiabatic approximation starts to be satisfied. We first focus on the case where adiabatic elimination is valid: two different situations emerge. The case where Γ = 4g 2 /κ > γ corresponds to a collective dynamical behavior for the atomic ensemble. In the low pumping regime under study, the emitted field is strongly bunched, g (2) (0) > 1, confirming previous theoretical works [9,11] on a similar model. In the two-atom case, this behavior has been attributed to the emission of so called "superradiant" photon pairs [9], typically separated by a delay scaling like the inverse of the pumping rate. This regime has been further explored for a larger number of emitters in [10,11], where the authors show that the atomic ensemble has mostly relaxed in the subradiant states. This is in agreement with Fig. 2a where the average intensity emitted by uncorrelated atoms, N n x is compared to the average fluorescence field emitted by the collective atomic dipole, J + J − . It can be clearly observed that in a κ/g range spanning several decades, the condition J + J − < N n x is usually satisfied, a signature of the anticorrelation of the atomic dipoles and of the efficient relaxation in the sub-radiant states. Note that the subradiant regime is reached when P x < Γ, a condition that is not fulfilled anymore for κ/g > 10 3 . In that case the steady state emission becomes super-radiant, i.e. J + J − > N n x , and correspondingly g (2) (0) approaches 1, a limit that is exactly reached only for N 1 as evidenced in previous work [11]. Finally, when κ/g is increased up to a point where γ > Γ (very bad cavity regime), emission appears as if it was coming from an independent collection of emitters. Under such conditions, the limiting values are such that J + J − → N n x and g (2) (0) → 2(1 − 1/N ), a limit that we analyze in detail below.
We focus now on the statistical properties of the system by looking at the dependence on the number of emitters, the results are shown in Fig. 3 with the same parameters used in Fig. 2. Keeping the cavity mode as a relevant degree of freedom allows one to explore the temporal correlations on the entire κ/g range, which can only be done by the full numerical solution of the master equation. At this point, we go beyond previous works. First, we notice that a transition to amplified spontaneous emission clearly occurs for very large cavity finesse, i.e. as soon as κ/g < 0.001 (P x > κ), with the temporal fluctuations becoming Poissonian independently of the number of emitters. A temporal picture for this effect is that as soon as P x > κ, the cavity stores the photons on a timescale longer that the typical delay between two emission events by the atomic assembly, blurring out any temporal correlation feature. This situation has been already evidenced for the single emitter case in [34][35][36], and it was identified as a lasing threshold signature. Here, the system becomes a source of radiation with Poissonian statistical noise, reminiscent of a few emitters laser. The case for N = 2 in the strong emitter-cavity coupling regime was also studied for a similar model in Ref. [37,38] as a function of pump power. As a matter of fact, three main effects on measurable quantities can be identified: (i) the average cavity population becomes larger than one, setting the onset of stimulated emission (see also Fig. 2a), (ii) the population of each emitter is clamped to a definite value, here, n x 0.5 (i.e. 5n x → 2.5 in Fig. 2a), and (iii) g (2) (0) → 1 (more clearly evidenced in Fig. 3 for any N = 2, 3, 4, 5).
On increasing κ/g, it is evident from Fig. 3 that in the subradiant range there is an unexpected and quite pronounced parity-dependent oscillation of the photon correlation value at zero-time delay, a behavior that has not been discussed in the literature, to our knowledge. Such parity oscillations are not limited to the adiabatic elimination region (see the inset on the left panel for κ/g = 0.45 and up to a number of emitters N = 8). The latter calculation required consistent numerical optimization and was performed on a cluster, where the full master Eq. (2) could be solved without truncation effects on the convergence. Although quite difficult to be analyzed in the dressed state picture, the oscillatory behavior of correlations as a function of the number of emitters can be roughly interpreted as a consequence of the collective coupling in Eq. (8), which dominates the dynamics over the individual relaxation channels in Eqs. (3) and (4). The former couples the states with the same quantum number J and is diagonal in the coupled basis. On the contrary, the other channels act on each emitter and are responsible for coherences between states with different values of J. We will investigate this regime in detail in Sec. IV. The bunching behavior reaches a maximum for κ g, slightly depending on the number of emitters. An interesting even/odd behavior has also been found in the spectral response [12] and it would be interesting to relate the two phenomena in future developments of the present work.
Finally, in the opposite limit of κ/g becoming large enough that γ > Γ the radiation is characterized by a field emitted from independent atoms. Oscillations tend to disappear, as it can be seen in Fig. 3 (see also inset on the right hand side). In this regime, the single emitters undergo individual dynamics, and the dipoles become uncorrelated (see also the large κ limit in Fig. 2a, where J + J − → N n x . Restricting our considerations to the correlation function, an analytical expression of g 2 (0) is derived in Sec. V, where we show that g (2) (0) = 2(1 − 1/N ). Note that in the limit of a large number of atoms, g (2) (0) → 2, which corresponds to the emission of an incoherent chaotic field.

IV. PARITY-RELATED OSCILLATIONS
In this Section we calculate analytical expressions for the bunching degree of the emitted field, in the two and three emitters cases, which provide a physical explanation for the observed parity dependence. The analytic treatment is valid in the regime of adiabatic elimination of the cavity mode dynamics. To sustain our analysis, we have also used a Quantum Monte Carlo (QMC) simulation [33]. We have calculated g (2) (0) as a function of the number of emitters N , with N being as large as 10. The results are plotted in Fig. 4a by simulating the reduced master equation, Eq. (7). As it can be seen, we alternatively recovered the result shown in Fig. 3 at a specific value of κ/g = 100, by generalizing to a larger number of emitters that would be unaccessible through direct numerical integration. Again, bunching and oscillatory behaviors are clear signatures of this collective low pumping regime, which we are going to analytically describe below.

A. Two emitters case
In the following of this Section we use the coupled basis for the emitters. Such basis of states is schematically indicated in Fig. 4b. In particular, the two-emitter case is pictured in the first panel of Fig. 4b, and we will focus on this diagram here. In this specific case, the basis |J, M is non-degenerate and the Hilbert space consists of the three symmetric (or Dicke) states |1, 1 = |e 1 , e 2 , |1, 0 = (|e 1 , g 2 + |g 1 , e 2 )/ √ 2, |1, −1 = |g 1 , g 2 , and the antisymmetric state |0, 0 = (|e 1 , g 2 − |g 1 , e 2 )/ √ 2 (i.e., the usual triplet and singlet states of a spin 1 system). The steady state (and g (2) (0) statistics) is the result of a competition between collective (specifically L Γ ) and individual decoherence channels (i.e., L Px , L γ , L γ * ). The results shown in Fig. 4a are obtained without pure dephasing and with P x = 10γ. In the steady state, low pumping equally populates the subradiant states |1, −1 and |0, 0 . As it was first discussed by Temnov and Woggon [9], one can decompose the problem in two mutually exclusive cycles. In the regime under study, the Γ-related de-excitation events are much more probable than the excitation events driven by the pumping rate P x . Hence, as soon as an excitation is stored in a state where the two mechanisms are possible, it is more likely that the excitation will undergo a collective emission than an upward pump jump. The microscopic scenario will then involve one cycle (pictured in green) consisting of the excitation and relaxation of the symmetric state |1, 0 , leading to the emission of single photons. In the other cycle (pictured in red), the excitation is first stored in |0, 0 . Then, during the time evolution, the only possibility is to excite the state |1, 1 , leading to the fast emission of two photons. The system prepared in the subradiant state |0, 0 behaves thus like a deterministic source of photon pairs. Each pair is emitted on a timescale 1/Γ, the delay between two pairs being typically 1/P x . Consequently, the emitted field is expected to show bunching, g (2) (0) ∼ Γ/P x , an order of magnitude that we exactly recover below by analytically solving the approximate master equation.
In fact, to go beyond this qualitative interpretation we have computed the expression for g (2) (0) in terms of the relevant elements in the steady state density matrix at first order in P x /Γ. The diagonal matrix elements of the density matrix are written with the row (column) index only, e.g. J, M |ρ|J, M = ρ J,M . In this two-emitter case the evolution equations for the diagonal elements do not involve off-diagonal terms and the solution is obtained by solving a system of rate equations governing the evolution of the atomic populations. With lengthy but straightforward calculations, we express Eq. (7) on the basis above, and solve it for the steady stateρ = 0.
We finally obtain an analytical estimation g (2) (0) = 16 9 which for P x /g = 10 −3 and κ/g = 100 gives g (2) (0) 17.8, a value that is very close to the numerical results (obtained from both direct numerical integration and QMC).

B. Three emitters case
In the N = 3 case, the simple picture above cannot be applied directly. The most important difference emerges immediately from the very nature of the Hilbert space. As pictured in Fig. 4b, right panel above, the J = J min = 1/2 subspaces have a 2-level system structure. As a consequence, the subradiant states corresponding to |J min , −M Jmin do not behave as deterministic sources of photon pairs as in the N = 2 case, but eventually give rise to single-photon cycles when they are pumped (light blue arrows in Fig. 4b). The presence of these extra loops spoils the bunching of the emitted field, and this gives a simple explanation for the abrupt loss of g (2) (0) bunching when going from two to three emitters.
To sustain this argument we have computed the second-order correlation g (2) (0) in terms of the density matrix elements, which for the N = 3 case reads (12) We followed the same notation as in Eq. (10) for the matrix elements. The subradiant states are written |1/2, −1/2; 1 and |1/2, −1/2; 2 [39], where the second label is a shorthand notation related to the g = 2 extra degeneracy (see App. A and the related figure). After writing the master equation in this basis, we then proceeded with a perturbative solution of the linear system for the steady-state matrix elements in powers of (P x /Γ). Differently from the N = 2 case, the equations of motion for the diagonal terms now involve off-diagonal matrix elements, preventing one from describing the dynamics with a simple set of rate equations. Nevertheless, within the same perturbative picture, we can consider these terms to be much smaller than the diagonal ones. Indeed, the coherences are brought by some residual individualization of the P x -driven channel that, for the parameters chosen, acts as a perturbation over the collective dynamics induced by the Γ events. The latter argument has been checked with the steady-state density matrix elements obtained numerically. As detailed in App. B, the final result reads We obtain the analytical value g (2) (0) 9.46, which compares very well with the numerically calculated value 9.55 reported in Fig. 4a from QMC simulations. Notice that the population of the state | 3 2 , 3 2 cannot be neglected and was also taken into account in the analytical calculation.

C. Generalization to N > 3 emitters
Considering the energy level diagrams in the generic even and odd cases, as depicted in Fig. 4b, lower panels, one can always find the two patterns above (i.e., N = 2 and N = 3) as distinctive features of each scenario. For increasing number of emitters more subradiant states become available, among which the population arising from the P x reservoir will be shared (bold levels in the figure). This spreading of population unavoidably leads to the decay of the bunching degree as N is increased, as observed in Figs. 3 and 4a, and confirmed by the general trend anticipated in [9] where a semi-classical approach was used. In the limit of large large numbers, such as N = 10000, the bunching degree of g (2) (0) converged to 2, as for classical chaotic sources. Note that in our case taking the limit of large numbers requires some caution, as one could exit the validity of the adiabatic elimination approximation, κ N Γ, and eventually enter strong collective light-matter coupling regime (see, e.g., recent works [41,42]). For the parameters chosen in the simulation in Fig. 4a, the condition for adiabatic elimination limits the maximum number of emitters to about 2000.

V. INDIVIDUALIZATION
We consider now the case where the cavity finesse is decreased down to a point where the coupling of each emitter to the mode, Γ, becomes of the same order of magnitude of the coupling to its individual decay channel, γ. As we have seen in Figs. 2b and 3, this leads to totally uncorrelated emission events from the atomic dipoles. In this regime, which we define as the "individualization," the parity-induced oscillations vanish and the bunching degree converges to g (2) (0) = 2(1 − 1/N ) (see Fig. 3, inset on the right), a behavior that we demonstrate below. Note that the same value for the parameter g (2) (0) has been computed in the strong pumping case [11]. As a matter of fact, a large pump power leads to total population inversion, a state where the dipoles necessarily become uncorrelated. Moreover, it is interesting to show in the context of the present work that an alternative way of obtaining a source of uncorrelated dipoles is to increase the pure dephasing rate, γ * . The latter effect is shown in Fig. 5, where a full numerical solution based on direct integration has been employed for κ/g = 100. It is evident that individualization is reached as soon as γ * > Γ, washing out the oscillations and leading to the same limit for the second order correlation parameter, as clearly displayed in the inset.
For the sake of completeness, we have computed the steady state second-order correlation g (2) (τ ), for any delay τ .
Adiabatic elimination is valid and this quantity can be calculated as a function of the atomic operator, namely g (2) (τ ) = J + (0)J + (τ )J − (τ )J − (0) / J + (0)J − (0) 2 . We make the assumption that the steady state of the ensemble can be factorized in the uncoupled basis, which is reasonable as each emitter behaves individually. Thus, developing (see App. C) the expression of J ± with respect to the individual atomic operator σ i ± leads to where we have introduced the intensity per emitter, the first and the second order correlation functions for the quantum emitter i, as As initially mentioned in seminal lectures [43] and papers [44], and as it clearly appears in Eq. (14), the second order correlation function involves the emission of pairs of atomic excitations. When τ = 0, this expression simplifies to confirming our numerical simulations (see right inset of Fig. 2). The bunching at small delays is due to an interference resulting from the indistinguishability of the atoms in the emitting pairs [43]. Note that surprisingly, this behavior was never observed experimentally. As a matter of fact, the interference is washed out as soon as individual first order coherence, g i (τ ), vanishes [31,[44][45][46][47]. This phenomenon takes place on the timescale of dephasing processes, leading to the following expression: eventually corresponding to a particle-like behavior. In fact, let's consider N independent sources each emitting classical identical particles one-by-one. Such a modelsource is characterized by a total emission intensity I = In the stationary regime we have I i (t) = I i (t + τ ) = I 0 , since we suppose that the N sources are characterized by the same main current. The total second order correlation signal is therefore g (2) tot (τ ) = I(t)I(t + τ ) I(t) I(t + τ ) (17) which becomes However, our N idealized sources are single particle emitters i.e. g (2) i (0) = 0. This clearly mimics the quantumlike behavior of two-level systems. Therefore, we get Eq. (16). Interestingly, from basic textbook quantum optics we obtain for a pure |n photon Fock state the expression g (2) (0) = 1 − 1/n. The result is therefore also intuitive from this point of view.
The limit given in Eq. (16) was already observed in the pioneering antibunching experiment by Kimble, Dagenais and Mandel with flying sodium atoms coherently driven [31], and it was theoretically interpreted in the following years [31,[44][45][46][47]. In these references, the source of dephasing was mainly attributed to the finite size of the atomic source with a random number (i.e. a Poisson distribution) of emitters inducing a spatial averaging out beyond the typical coherence volume of the fluorescence light [31,46]. In nowadays experiments with solid state emitters confined in nanometre-size systems, this problem could be ultimately overcome by considering a good control of the optical detection. However, notice that dephasing processes are always present in solid state quantum emitters (due for example to coupling with phonons in solid matrices), and that the measurement of g (2) (0) involves actually an average over finite values of τ near τ 0. If the typical dephasing time τ * = 1/γ * is larger than the time resolution of the HBT setup (e.g. the jitter time of single photon detectors) this may explain why experiments reported so far with different quantum emitters [48][49][50][51][52][53][54] tend to evidence a behavior that is closer to Eq. (16) rather than the present prediction described by Eq. (15). Observing the 2(1 − 1/N ) behavior, and its transition to 1 − 1/N , could be a challenging issue, e.g. requiring a control over the local temperature of the emitters.

VI. PHYSICAL REALIZATIONS
The results presented in this work suggest the intriguing possibility of studying the effects of cooperative emission with a mesoscopic number of quantum emitters in a high-finesse resonator. Moreover, depending on the specific physical implementation of the proposed model, different regimes of cavity QED might be explored and characterized in terms of the statistical properties of the emitted radiation. To this purpose, we stress the rele-vance of Fig. 3 in providing a roadmap to experimental investigations. In fact, depending on the specific system of few quantum emitters coupled to a single cavity mode, different regimes of κ/g might be realized.
However, some of the requirements that we assumed in our analysis might be difficult to be simultaneously realized in a physical implementation of the proposed model. Namely, we assumed identical quantum emitters, equal dissipation rates, low pumping of each individual emitter as compared to its coupling to the cavity mode. In order to realize the conditions γ * , γ < P x < Γ, each quantum emitter should have a high beta factor (quantifying the coupling to the cavity loss channel). Usually, this is quite challenging with solid state artificial atoms, like quantum dots, where inhomogeneous broadening (which can be related to γ * ) would inevitably mix the different J states. Atomic systems in optical resonators would certainly offer an ideal implementation for what concerns the atomic-like properties, but it would be difficult to have a mesoscopic number of such atoms at fixed spatial positions. Alternatively, atomic-like impurities in solid state resonators, e.g. semiconductors, might represent a valid scenario for this type of physics to be investigated.
As a first example, we consider semiconductor quantum dots, especially made in III-V materials, as a promising approach to coupling in a deterministic way [22,55] a few quantum emitters to the same cavity mode, e.g. using photonic crystal cavities [14]. Such systems can have g 0.1 meV [22,23] and κ 0.01 meV [56]. We can thus expect to be in the range κ/g ∼ 0.1 for such a physical implementation, where subradiant vs. lasing behaviors could be studied under weak driving conditions. As these nanostructures are usually grown randomly, and strong inhomogeneous fluctuations between different quantum emitters occur, techniques have been recently developed to post-process the sample and obtain more homogeneous ensembles [57]. Electrical control of different quantum emitters has also been shown to provide an effective tuning mechanism [58]. Moreover, growth of pyramidal dots in pre-patterned positions [59] is a very promising technique to place more than one such nanostructures in the antinodes of the cavity field in order to maximize coupling.
For active centers in solid-state resonators, the coupling to the cavity mode is usually several orders of magnitude smaller, dictated by the smaller oscillator strength of an atomic-like impurity as compared to a quantum dot. In such case, the conditions for high-κ/g values might be investigated, i.e. the crossover from subradiant to superradiant emission, and maybe to individualization. Recently, single to few Nitrogen-Vacancy centers in diamond have been implanted in pre-determined spatial positions [60]. Resonators in diamond have already been realized with quality factors ranging from few hundreds in early attempts [61,62] up to 4000 more recently [63]. Furthermore, NV-centers are not the only type of active defects in diamond, see e.g. recent advances with Cr-centers [64] and Silicon-vacancies (SiV centers) [65]. Among solid state emitters, we also include single molecules in standing wave resonators made from a Bragg mirror and a resolution limited focused beam through an optical fibre [66]. The latter approach has the advantage that the cavity can be formed at any arbitrary position, once the single molecules have been selected on the Bragg mirror surface. Such quantum emitters display radiation-limited zero-phonon line, thus showing very limited pure dephasing rate with almost radiation-limited linewidths [67].
Finally, a promising platform to realize a controlled experiment might be with circuit QED devices, where the coupling of a few qubits to a single resonator has already been realized [15]. High tunability of the qubits physical properties, such as transition energy, dissipation rates, pure dephasing and effective oscillator strength, make these systems particularly attractive to closely realize the model in Eqs. (1)(2), and to experimentally address the physical results reported in the present work.

VII. CONCLUSIONS AND PERSPECTIVES
We have presented an extensive theoretical analysis of the quantum statistical properties of a source of electromagnetic radiation made of few two-level emitters coupled a single resonator mode. We have assumed a model in which the two-level systems are incoherently pumped by a weak external drive, and dissipate their energy through spontaneous emission in free space, pure dephasing, or other loss channels that can be described within a Liouvillian formalism. In this sense, the present work should be considered as an extension of our previous report on single emitter cavity QED [34], where we have analyzed the crossover from a single-photon source (emitting anti-bunched radiation) to a single-emitter lasing regime.
However, the case of more than one emitter studied here presents a much richer phenomenology depending on the specific light-matter coupling regime of each single emitter with the common radiation field, which can be identified from the statistical properties of the emitted radiation. Throughout this work, we have assumed a tuning parameter given by the ratio of the cavity decay rate to the single emitter-cavity coupling strength, which we could numerically investigate by spanning over 10 decades. As a function of this parameter, we have clarified the crossover from lasing (low cavity decay rate as compared to emitters cavity coupling) to chaotic light emission (very large cavity decay rate), going through sub-and super-radiance in steady state. Although in the superradiant regime the emitted radiation would present a close-to-Poissonian statistics, we have explicitly shown that it would not reach the level of coherence displayed by a laser source. Between these two regimes, sub-radiant emission takes place and the emitted radiation undergoes super-Poissonian intensity fluctuations.
We have shown that the strong bunching in the degree of the second-order coherence oscillates with the parity of the number of emitters, eventually converging to a single value in the large number limit. We have clarified the origin of this novel behavior by using a suitable basis for the system of few emitters, showing that it is due to the presence of extra-loops in the radiative cascade from subradiant states when the number of emitters is odd, thus reducing the degree of bunching with respect the even case. Finally, we have also studied the conditions for reaching an incoherent source of chaotic radiation, which can be done either by coupling the ensemble of emitters to a very leaky cavity, or by enhancing the pure dephasing rate (a feature that is quite natural for solidstate quantum emitters).
As a last remark, notice that the cavity finesse has been the only parameter to be significantly changed in order to produce very different output fields. Recently, fast dynamical control of this parameter has been demonstrated in different scenarios [68,69], suggesting that a system composed of few emitters inside a resonant cavity could be a versatile source of very different quantum states of light. An immediate and important extension of this work related to its experimental and technological implications would be to investigate the effects caused by the inhomogeneous distribution of the quantum emitters, which is a typical setting for smiconductor quantum dots. A further seemingly important extension would be to explore the consequences of the ultra-strong coupling regime (available nowadays for superconducting qubits in microwave resonators) on the statistics of the emitted radiation. In the latter scenario, anti-rotating wave terms in the hamiltonian (1) must be taken into account [70] in order to properly get the emitted radiation and statistics. multaneous diagonalization of the angular momentum J (with J max = N/2) algebra with the permutation group P N of N atoms [39]. We write the states in the coupled basis as where J and M J represent the collective dipole and its z−projection labels whereas i enumerates the states belonging to the λ permutation irreducible representation.
The derivation reduces to the use of standard branching rules for Young tableaux with the theorem in Sec. V of Ref. [39] or with the tensor method [40]. Since we are dealing with a tensor product of two-dimensional Hilbert spaces, C 2 , then for 3 emitters the branching rules give three types of tableaux, out of which we must cancel those with more than two rows, as schematically outlined in Fig. 6. By means of Young symmetrizers, the coupled basis is easily obtained in terms of the uncoupled states as In this appendix we show the analytical details for calculating the autocorrelation at zero-time delay for the N = 3 case.
Once we neglect the off-diagonal terms (see text in Sec. IV B), the two degenerate states |1/2, ±1/2; 1 , |1/2, ±1/2; 2 give the same values, and here they are considered simply as |1/2, ±1/2 . The lowest matrix element is calculated from the constrain on the trace of the density matrix. It reads (B2) In powers of P x /Γ we finally obtain g (2)  .

(B3)
The same argument could have been used for the twoemitters case, for which the solution reads g 2 (0) 4 9 Γ P x + 13 27 − 2 27 In the case we are considering, namely Γ/P x = 40, the error in neglecting all the terms but the first one in Γ/P x is of about 3% for Eq. (B4), while for eq. (B3) it would have been ∼ 30 %. This formally explains why in the two-emitters case the linear approximation works so well. Indeed, setting Γ/P x = 40 in Eq. (B3) we have g (2) (0) 9.46 as given in the text.
Appendix C: Individualization and chaotic limit Developing the expression for J ± with respect to individual atomic operators leads to g (2) (τ ) = i,j,k,l σ i + (0)σ j + (τ )σ k − (τ )σ l − (0) k,l σ k We are considering the case where the atomic dipoles are independent and incoherently excited. It means that whatever i, j we have being ξ i either σ i − or σ i + . It is immediate to see that whenever i = j, k, l (and permutations) the four operators expectation value is zero. We are left with i = j = k = l and the case where the indices are equal to each other in pairs. As a consequence where I k = σ k + (0)σ k − (0) For zero delay, the denominator D involves J + J − = i,j σ i + σ j − , which simply writes i σ i + σ i − = N I, as the dipoles are independent. As for the numerator, the case i = j = k = l gives no contribution as the field emitted by a single atom is antibunched, the pairs (i = k, j = l) and (i = l, j = k) contribute with 2N (N − 1) terms equal to I 2 (by definition, we have g (1) l (0) = 1). Eventually, one finds the g (2) (0) = 2(1 − 1/N ) result mentioned in the text. It is worth noticing that in the large number limit we have g (2) (0) = 2 corresponding to a chaotic chaotic source [43,44]. It is also interesting to remark that with a distribution of N classical electromagnetic emitters we still have g (2) (τ ) = k g (2) k (τ ) + k =l g (1) k (τ )(g (1) l (τ )) * N 2 However, the classical Schwartz inequality imposes g (2) i (0) ≥ 1 and therefore we now obtain for the N emitters The value 2(1 − 1/N ) is therefore somewhere intermediate between the pure particle like behavior, 1 − 1/N (see main text), and the pure classical undulatory prediction, 2 − 1/N (in agreement with the usual Einstein interpretation of the fluctuation formula for the blackbody radiation).