Quantum Langevin approach for superradiant nanolasers

A new approach for analytically solving quantum nonlinear Langevin equations is proposed and applied to calculations of spectra of superradiant lasers where collective effects play an important role. We calculate lasing spectra for arbitrary pump rates and recover well-known results such as the pump dependence of the laser linewidth across the threshold region. We predict new sideband peaks in the spectrum of superradiant lasers with large relaxation oscillations as well as new nonlinear structures in the lasing spectra for weak pump rates. Our approach sheds new light on the importance of population fluctuations in the narrowing of the laser linewidth, in the structure of the lasing spectrum, and in the transition to coherent operation.


I. INTRODUCTION
Progress in various technologies has enabled considerable size reductions of lasers. Nowadays quantum dot photonic crystal [1][2][3][4][5][6][7], micropillar [8][9][10], plasmonic [11] and other kinds of nanolasers [12] are intensively investigated. This research is motivated by fundamental questions, such as the minimum size of lasers and the role of quantum effects. The miniaturization of nanolasers is also driven by applications, in nano-electronics for example, where energy-efficient nanolasers are directly incorporated into nano-chips [13][14][15]. The high density of photon states in nanocavities leads to Purcell enhancement [16] of spontaneous emission into the nanolaser mode, large gain and to the rapid increase of laser power even at small pump rates.
Nowadays there is great interest in superradiant lasers, which are lasers that combine a large gain with a small cavity operating in the so-called bad-cavity regime [17][18][19]. In this regime, the polarization cannot necessarily be adiabatically eliminated and collective spontaneous emission into the lasing mode is significant. Superradiant lasers have been experimentally realized, for example with cold alkaline earth atoms [20][21][22][23], rubidium atoms [24], and with quantum dots [25] as the active medium. Superradiant lasers are less sensitive to cavity-length fluctuations, which is important for atomic clocks [20,21,24]. Superradiance leads to interesting collective effects, such as excitation trapping [22,24] and superthermal photon statistics [10,25,26], with possible applications in high-visibility optical imaging [27].
An analytical description of superradiant nanolasers and their spectra is complicated by the facts that their quantum noise is not a perturbation, that the equations are nonlinear, and that the polarization of the active medium cannot be adiabatically eliminated. We will address these issues in this paper, where we present an analytical approach to understanding superradiant lasers.
The quantum theory of lasers began with applications of methods of classical statistical radio-physics first for lasers comprising a cavity with high quality (Q) factor [28][29][30] and later also for low-Q cavity lasers [31,32]. In many papers the fluctuations of amplitude and phase of laser radiation are considered separately, in the frame of rate equations, where the active-medium polarization is adiabatically eliminated [33][34][35][36][37]. This approach is satisfactory for usual semiconductor lasers and leads to various analytical results, but the approximations leading to the usual rate equations are not always justified for superradiant nanolasers.
Presently nanolasers are theoretically modeled either by rate equations as in [37,38], by numerical solution of the density matrix equations as in [22,[39][40][41], or by systems of equations for correlations as in the cluster expansion [25,42] or cumulant expansion [43,44] methods. Numerical analysis of superradiant emission and lasing has recently led to new and interesting results, such as mechanical effects in photon-atom interactions [45], lasing with a millihertz linewidth and rapid emitter number fluctuations [46], Wigner functions for semiconductor heterostructures [47], transition from superradiance to regular lasing by varying the coherent and incoherent driving [44], sub-and superradiance in multimode optical waveguides [48], and photon-antibunching in the fluorescence from an optical nanofiber-tip [49]. However, complementary analytical methods to model nanolasers without adiabatic elimination of polarization, that would apply to superradiant nanolasers, are not well developed.
Here we use Heisenberg-Langevin equations, which are very convenient for the description of lasers [50], to describe systems ranging from LEDs to superradiant nanolasers. The method also describes the limit of nonsuperradiant lasers, where it will reproduce some wellknown results. We follow the approach by Lax, see for example Ref. [32]: operators are treated as stochastic variables, while quantum properties such as non-vanishing commutation relations are taken into account by diffusion coefficients.
Our first new application of the method will be the description of both the lasing field and the active-medium polarization by symmetric (S) and anti-symmetric (A) combinations of quadrature operators. Quadratures of the polarization have been used in Lamb's semiclassical laser theory [30], while quadratures of the electromagnetic field were used, for example, in the analysis of the driven Van der Pol oscillator applied to lasers in [36] and in quantum optics to squeezed states of the radiation field [50]. Symmetric and anti-symmetric combinations of the quadratures of a quantum field are also used in entanglement criteria [51,52]. Our work is an extension since, to our knowledge, quadratures as well as their S/A combinations for both the lasing field and polarization combined have not been used in laser theory before. The approach has several advantages: it does not require a quantum phase operator [53], and furthermore we find the equations for the S/A quadrature combinations less cumbersome than for the density matrix [50].
Our second new application is the linearization of the quantum Langevin equations, where fluctuations are not necessarily small compared to mean values. This goes beyond a small-signal analysis as, for example, in Refs. [36,54,55] and requires other approximations. We describe our method in detail, and it may be useful also for other physical systems with resonant nonlinear interactions of light with matter, see examples in Ref. [56].
Section II introduces the quantum Maxwell-Bloch equations with dissipation for a two-level laser and we rederive some key results of semiclassical laser theory [50,57] to be used for comparison later in the paper.
In Sec. III we demonstrate a method, used before in [58,59], for approximately solving the nonlinear quantum Maxwell-Bloch equations by use of a Fourier transform. With this method, we extend the semiclassical approach by taking into account spontaneous emission into the lasing mode below the semiclassical threshold, where population fluctuations can be neglected. In subsequent sections we do take population fluctuations into account. We obtain the expression for the laser linewidth as a function of the population inversion, as in [59], reproduce the well-known laser linewidth at small excitation, as in [57] and derive the beta-factor for bad-cavity lasers, as in [38], three results to illustrate the efficiency of our method.
Secs. IV-VI are the main parts of the paper. In Sec. IV we represent the lasing field and polarisation by the S-and A-combinations of quadrature-operators, and derive the central linear equations of our approach. We will show that in linear approximation only the Scombinations interact with population fluctuations, while the A-combinations do not. We justify the approximations made to linearize the initial nonlinear Maxwell-Bloch equations.
In Sec. V we solve the equations for the Acombinations, show that they describe the laser output power and the linewidth in the high-pump limit, and reproduce the formula for the laser linewidth in that limit.
In Sec. VI we derive expressions for lasing spectra and show examples of the analysis of spectra for superradiant and non-superradiant lasers. The final two sections contain discussions and our conclusions.
The novelty of our method is that we analyticaly describe the laser below, near and above the threshold by the same set of stochastic equations, taking into account the field, polarisation and population quantum fluctuations, spontaneous emission into the lasing mode and full laser dynamics without adiabatic elimination of polarization. Such accurate treatment is in particular important for superradiant lasers, where collective effects among the emitters need to be taken into account. With our method we reproduce well-known results and identify features, in particular in the laser spectra, that largely went unnoticed. In particular, we calculate the full spectrum of the lasing field below as well as above threshold and identify and explain the appearance of a broad spectral background above threshold and a multi-peak structure above and below threshold. In particular, we focus on the role of population fluctuations and nonlinear polarisation dynamics in superradiant lasers.

II. QUANTUM MAXWELL-BLOCH EQUATIONS. SEMICLASSICAL LASER MODEL
In order to keep the analysis simple, we consider a stationary single-mode laser, shown schematically in Fig. 1, with a large number N 0 1 of homogeneously broadened identical two-level emitters, with their transition frequency ω 0 equal to the cavity mode frequency. We Figure 1. Scheme of the two-level laser, with parameters and operators as defined in the main text. It is common to assume that γ ⊥ 2κ, γ so that the polarization can be adiabatically eliminated and the laser is well described by rate equations. This does not work for the superradiant nanolasers considered here, where γ ⊥ < 2κ.
write the quantum Maxwell-Bloch equations (MBE) for such a laser [50] in the rotating-wave approximation with carrier frequency ω 0 , Hereâ is the annihilation operator of the laser cavity mode, the operatorv = i N0 i=1 f iσi describes the polarization of the emitters,σ i describes transitions from the excited to the ground state of the i th emitter; f i characterizes the coupling of the i th emitter to the lasing mode. We also define the average coupling f = N −1 0 i f 2 i . Furthermore, Ω 0 is the vacuum Rabi frequency and the total excited-and ground-state population operatorsN e,g are defined as the sums N0 i=1n (e,g) i , wheren e i (n g i ) are operators of populations of excited (ground) states of the i th emitter. The operatorN =N e −N g is then the population inversion. In general we use 'hats' to denote operators, while mean values of operators are indicated by the absence of a hat, for example N = N . We will consider the stationary case, so mean values do not depend on time. The laser field leaves the cavity through the mirror at the cavity decay rate 2κ; γ is the population relaxation rate of the upper lasing level, γ P is the pump rate from the lower to the upper level; γ ⊥ /2 is the polarization relaxation rate (so that γ ⊥ is the width of the lasing transition). The total number of emitters is assumed preserved,N e +N g = N 0 , so we can rewritê The quantum operators in Eq. (1) are interpreted as stochastic variables, where Langevin forces and their associated correlation strengths ensure correct quantum properties. In more detail, we introduced the Langevin forcesF α , with α taken from the set {a, v, a + , v + , N e }. These describe white noise, have zero mean (i.e. F α = 0), and are delta-correlated in time: F α (t)F β (t ) = D αβ δ(t − t ), where D αβ are the diffusion constants. In the Fourier-domain the cross-correlation of these Langevin forces is given by From Eq. (1a) and its Hermitian conjugate we then obtain in terms of the stationary mean photon number n = â +â . In combination with Eq. (1c) we can eliminate the atom-field correlations and obtain the energy conservation law 2κn + γ N e = γ P N g .
Stationary solutions of Eq. (1) are readily obtained, if we neglect the Langevin forces and replace operators by c-numbers:â → a,v → v,N → N . This gives the conventional stationary equations for the stationary solutions of the semiclassical laser model [50] Eqs. (5a), (5b) have non-vanishing solutions if From the energy conservation law (4), we find the stationary number of photons in the semiclassical model to be So the semiclassical model predicts lasing, n > 0, when N 0 > N th and when the dimensionless pump rate P exceeds the semiclassical lasing threshold P th For a dimensionless pump rate P smaller than P th , lasing is absent and n = 0 in the semiclassical laser model.

III. ANALYSIS NEGLECTING POPULATION FLUCTUATIONS
Nanolasers have large beta factors, so that spontaneous emission into the lasing mode is non-negligible [2][3][4][5]. Below and close to the semiclassical threshold, nanolasers are not well described by the standard semiclassical model, which neglects spontaneous emission and predicts zero photons. One can improve upon this when γ ⊥ 2κ, γ , in which case polarization can be adiabatically eliminated. This is typically the case for semiconductor lasers. In this case, the laser can be described, for all pump rates, by quantum rate equations (QRE) [38,60,61], that do take into account spontaneous emission into the lasing mode. The intensity noise spectra can then also be found by small-signal analysis of the QRE [38]. However, the QRE cannot be applied to superradiant lasers, where γ ⊥ < 2κ and the active-medium polarization cannot be adiabatically eliminated.
In this Section we go beyond the rate-equation approach in the low-pump limit by including the dynamics of the polarization, rather than eliminating it adiabatically. Instead, following [58,59], our main approximation will be that we neglect population fluctuations. This is a good approximation for small pump rates, when the material gain is much smaller than the cavity loss, such that fluctuations of the populations, and thereby of the material gain, do not significantly change the net cavity gain [34]. We also take into account spontaneous emission into the lasing mode and introduce a Fourierexpansion method that is used throughout the paper. We will reproduce the well-known formula for the lasing linewidth in the low-pump limit and introduce the βfactor for lasers with low-quality cavities. These results will be used as a reference for comparison in the following sections, where we do take population fluctuations into account.
Neglecting population fluctuations, we replace population operators in Eq. (1b) by their mean valueŝ N e,g ≈ N e,g ,N ≈ N.
Thereby Eqs. (1a) and (1b) turn into a set of linear equationsȧ We expressâ(t) andv(t) and their corresponding Langevin forces through Fourier-component operatorŝ forα = {â,v,F a ,F v }, and obtain from Eqs. (10) linear algebraic equations for all α(ω) and find from them Coming back fromâ(ω) toâ(t) by an inverse Fouriertransformation, we calculate the mean number of photons in the cavity as where n(ω) is the spectral power density of the field in the lasing mode, or optical spectrum, which is related tô a(ω) as We will determine n(ω) and then find n by Eq. (13). In order to find n(ω) we must know the relevant diffusion coefficients [54,62]. After neglecting, as is usual, any thermal radiation in the lasing mode, since k B T ω, we take the diffusion coefficient 2D a + a = 0. When population fluctuations are neglected, the diffusion coefficient as shown in Appendix A. With these diffusion coefficients, we find the optical spectrum (16) This spectrum may either have one or two peaks. Two peaks occur when all emitters are collectively and strongly coupled to the lasing mode, under the condition and when P < P c , where P c is such that N (P c ) = −N c . The two peaks in n(ω) are then caused by and a signature of collective Rabi splitting [59]. Otherwise, n(ω) has a single peak, with full width at half maximum γ las , defined by n(γ las /2) = n(ω = 0)/2, with the value where the parameter r is given by For r 1, as obtained for pumping levels where N is close to N th , we expand Eq. (18) as a series in r and to first order in r obtain where γ c = 2κγ ⊥ /(2κ + γ ⊥ ). Examples of optical spectra calculated according to Eq. (16) are given in Fig. 2, showing both cases of single-and double-peaked spectra.
Our goal in the remainder of this section is to express the linewidth γ las in terms of familiar laser parameters. In order to do so we first determine the population inversion N . (Incidentally, the same procedure to calculate N will be used later again, when we also take population fluctuations into account.) From Eqs. (13) and (16) we find for the mean number of photons in the cavity By inserting this into the energy conservation law (4), we obtain a quadratic equation for the population inversion, where we introduced the parameters following Refs. [38] and [61], respectively. In the special caseβ c = 0 the two solutions of Eq. (21) for N coincide with the stationary population inversion found in semiclassical laser theory. The general solution of Eq. (21) withβ c = 0 is different, because the approach that led to Eq. (21) takes into account spontaneous emission into the lasing mode. Indeed, in the limit 2κ/γ ⊥ → 0, the coefficientβ c tends toβ, which was introduced in Ref. [61] as the ratio of the rate of spontaneous emission into the lasing mode to the rate of all other emission processes (i.e. background emission). Solving Eq. (21), we find the pump-dependent population inversion N (P ), given by Eq. (B1a) of Appendix B. By inserting N (P ) into Eq. (18) we obtain an explicit expression for the pump-dependent linewidth γ las (P ). Similarly, by inserting N (P ) into the energy conservation law (4) we obtain the pump-dependent photon number n(P ), given in Eq. (B1b) of Appendix B.
We can now express the linewidth (19) in terms of the laser output power W out = 2κ ω 0 n by using Eq. (20) to express the factor (1 − N/N th ) in terms of W out , giving where N sp = N e /N th is the so-called spontaneousemission factor [63]. Eq. (23) is the well-known result for the laser linewidth below threshold, which (apart from notations) coincides with, for example, results in Refs. [57,63]. It is generally accepted that the laser linewidth far above threshold is suppressed by a factor of two compared to Eq. (23) [57]. We notice, though, that recent work [64] challenges this result, based on a semiclassical analysis. In the next sections we will show that a fully quantum mechanical theory for the lineshape far above threshold agrees with adding the extra factor of 1/2 to Eq. (23), and that it can be ascribed to the effect of population fluctuations, in particular to relaxation oscillations induced by population fluctuations.

IV. LINEARIZATION OF EQUATIONS
Our aim is now to develop a theory for the optical spectrum of a laser, valid at any pump rate, without making the assumption that fluctuations in the lasing field and polarization are always small. Our linearization procedure is therefore different from the usual small-signal analysis as presented, for example, in Refs. [36,37,54,55]. The theory is still approximative and we shall seek to clearly identify the approximations made.
We begin with the linearization of Eq. (1b) by writing the population operators as the sum of their mean values N e,g and population fluctuations δN e,g , whereby the population fluctuations are defined. We consider a large number N 0 1 of emitters and suppose small fluctuations of populations δN 2 e,g 1/2 N e,g . We also suppose weak coupling, 2Ω 2 0 f /γ ⊥ γ 1, and, for superradiant lasers, low-Q cavities with 2κ ≤ γ ⊥ . The mean photon number for such a laser below and near the semiclassical threshold is of the order of unity or less, which we see in Fig. 5(a). So we do not assume that fluctuations of the lasing field and polarisation are small compared to their mean values.
We next insert Eq. (24) into Eqs. (1) and obtaiṅ We shall now show that instead of the conventional representation of the laser field in terms of amplitude and phase, it is convenient to represent the field and polarization by their quadratureŝ whereα stands forâ orv. In our stochastic approach, α x,p are represented by real-valued stochastic variables. The equations of motion for the quadratures follow from Eqs. (25), The Langevin forces in Eqs. (27a) and (27b) arê where α stands for a or v. We linearize Eqs. (27) in several steps. In the first step we consider only the nonlinear terms in Eqs. (27b) and (27c). In these terms, we neglect low-frequency fluctuations of the field and polarization, in a frequency range ∆ω around ω = 0. In more detail, we approximatê where the index µ stands for x or p, and the c-number V will be determined below. The "cut-off" frequency ∆ω is chosen such that â 2 µ n and v 2 µ V 2 , i.e.â µ andv µ are small perturbations relative to √ n and V . In Sec. VI we will calculate the optical spectra [see Eq. (49)] and find that the approximation (28) can be made in the high-excitation limit, where almost all energy of the lasing field resides in a narrow spectral peak n A (ω) of width γ las around the optical frequency ω 0 , with only a small part of the energy in the wide spectral background n S (ω) of width γ bg γ las . The peak n A (ω) and the background n S (ω) of the full optical spectrum are depicted in Fig. 3. The n A (ω) and n S (ω) correspond to (but do not exactly coincide with), respectively, the field W ε and the amplitude W δ spectra calculated, for example, in Ref. [35]. Figure 3 shows that ∆ω can be chosen from the interval γ las < ∆ω γ bg , and we shall see that the final results are independent of ∆ω. Later we will explain that the approximation (28) can safely be made in the nonlinear terms in Eq. (27b) also at low or moderate laser excitation.
We insert the approximation (28) into Eq. (27a) and find One may expect that in Eqs. (28), different c-numbers should be defined for the different quadratures, instead of only √ n and V . However, the different c-numbers reduce to √ n and V by the replacementsâ →âe iϕ andv →ve iϕ with real-valued phase ϕ, which is the constant phase of the lasing field. The solution of the initial MBE (1) does not depend on such a replacement, which adds only a constant phase multiplier to the Langevin forcesF a and F v . Therefore √ n and V in Eqs. (28) can be chosen to be the same for both quadratures in the general case. In the approximation (34),â S (ω) is extended to the frequency interval |ω| < ∆ω/2, which is the small dotted segment in the nS(ω) curve.
In the second linearization step, we insert the approximation (28) into only the nonlinear terms in Eqs. (27). Then, by neglecting the products of fluctuations in these terms, we obtain linear equations forâ x,p andv x,p . This step is different from the usual small-signal analysis, where linear equations should be written for the perturbationsâ x,p andv x,p .
In the third and final step of the linearization, the set of Eqs. (27) can be further simplified, if we extract the symmetric (S)α S and anti-symmetric (A)α A combinations of the quadratures and their high-frequency fluctuations (denoted by primes)  las . Direction fluctuations are much slower than length fluctuations, so after ∆t the difference angle θ 1. Also shown are the quadraturesâ x,p as well as the symmetriĉ a S and anti-symmetricâ A combinations that are defined by Eqs. (30) and used in the approximation (28). For θ 1, theâ S andâ A combinations determine, correspondingly, the amplitude and phase fluctuations of the field. The green (the red) dashed vectors denote the change of the field vector 1 due to only phase (amplitude) fluctuations.
After the three linearization steps, we obtain, finally, from Eqs. (27) as one of our main results the approximate linear equations for the S-and A-combinations of quadratureṡ and where α stands for a or v. The derivation of Eqs. (31) is given in Appendix C. We see that Eqs. (31) are split into two sets of equations: the two Eqs. (31a) and (31b) for the A-combinations of the quadratures and the three equations (31a), (31c) and (31d) for the S-combinations and the population fluctuations. Both sets can be solved independently, but their solutions are related through their common c-number variables N e and n. For clarity, let us summarize the approximations made in the derivation of the system of equations Eq. (31): we neglected low-frequency fluctuations in the approximation (28); in Eq. (31b) we neglected the small term a S with respect to the large term ∼ √ n; similarly, in Eq. (31c) we neglected a termâ A δN e since the termâ A N is much larger because δN 2 with respect to 2 √ n (Ω 0v S + κâ S ) and we approximated as described in detail in Appendix C. Now let us find and discuss the solutions of the equations Eqs. (31). We solve them by Fourier transformation, and we find â 2 S,A ≡ n S,A (N e ). In order to find the mean photon number n, we expresŝ insert this Eq. (35) into n = â +â and obtain Inserting Eq. (35) into [â,â + ] = 1 we find [â A ,â S ] = i/2, while for the mean photon number we obtain By inserting this Eq. (37) into the energy conservation law (4), we obtain an equation for N e . Though Eqs. (31a), (31c) for the A-combinations do not depend on population fluctuations explicitly, the Acombinations do depend on n S and, therefore, on δN e implicitly, through the energy conservation law (4), where n(N e ) is given by Eq. (37). Because of this, we will obtain different relations between n A , n S and n at low and at high excitation.
In order to justify the approximation (34), we note that the Fourier expansions ofâ S andv S by their definition in Eq. (28) have the same Fourier components asâ S andv S for frequencies |ω| > ∆ω, but instead have zero Fourier components for |ω| ≤ ∆ω. The Fourier components of a S andv S are smooth functions of ω, as confirmed by Fig. 3, so as an approximation we extend them to the small interval |ω| ≤ ∆ω, which is exactly the approximation (34).
The approximation (28) was made during the linearization of equations (31) in the high-excitation limit, when the spectrum of the lasing field has a narrow peak and a broad background, as in Fig. 3. Most of the energy is concentrated in the peak and the cut-off frequency ∆ω is about as small as the width γ las of the peak, soâ √ n andv V . These inequalities do not hold for low and moderate excitations. For such excitations, however, population fluctuations are so small that, as a zero-order approximation one may neglect them altogether, as in Sec. III. The linearized equations (31) are still approximate, but they are the next-order approximation at low excitation, where they describe population fluctuations incompletely, but do not neglect them. This is an argument for using the equations (31) for all excitation levels as the first approximation to incorporate small population fluctuations. This approach is analogous to the one in Ref. [38], where it was shown that in the rate-equation limit, the standard linear small-signal analysis, when carried out for large number of emitters and weak coupling, agrees very well with results of exact numerical nonlinear analysis also below the lasing threshold.
For these reasons, below we will use the linearized equations (31) at all levels of excitation.

V. SOLVING THE COMBINATIONS OF THE FIELD QUADRATURES
In this section we first solve the equations (31a) and (31c) for the A-combinations of the field and polarization quadratures, these equations being the ones that do not explicitly depend on population fluctuations. In doing so, well-known results for the laser linewidth in the high-pump limit will be reproduced.
We replaceâ A andv A in Eqs. (31a) and (31c) by their Fourier expansions (11), then solve the linear equations for the Fourier componentsâ A (ω) andv A (ω), and find Fromâ A (ω) we can calculate the corresponding optical spectrum n A (ω) by analogous to Eq. (14). The relevant diffusion coefficients are given by as calculated in Appendix A. With the help of Eq. (38), the diffusion coefficients (40), and the correlations of Langevin forces (2), we find the spectrum of the Acombinations of the photon field to be The photon numbers n A,S residing in the A-or Scombinations are in general given by the integrated spectra By carrying out this integration with n A (ω) given by Eq. (41), we find the number of photons in the Acombinations of the field quadratures as So here we find, perhaps not surprisingly, that n A can grow very large when the laser is pumped strongly and the population inversion N approaches its semiclassical stationary value N th . However, in Sec. VI we will see that in the same high-pump limit, the number of photons n S in the S-combinations of the quadratures is much smaller than n A in Eq. (43). Mathematically, this is the case because n S does not have a corresponding term ∼ 1/ (N th − N ). Physically, because n S is suppressed by relaxation oscillations introduced by population fluctuations. So almost all lasing photons reside in the Acombinations and n A ≈ n. Anticipating these results for the spectra and for the number of photons in the Acombinations, we now replace n A by n in Eq. (43), and use this to derive an expression for the linewidth as a function of the laser output power, as we did before in Eq. (23): in Eq. (43) we neglect the term 2κ/γ ⊥ , which above threshold is small compared to the large first term ∝ (N th − N ) −1 . The linewidth γ las is again defined as the full width at half maximum of the optical spectrum, i.e. n A (γ las /2) = n A (ω = 0)/2. This leads to exactly the same expressions (18) and (19) for the linewidth γ las as found previously for the spectrum (16), when expressed in terms of the average population inversion, N . However, the crucial difference is that the variation of N with the pump level P changes quantitatively as the laser threshold is passed, leading to different dependencies of the laser linewidth on power above and below the laser threshold.
To see this, we express (N th − N ) in terms of n and the laser output power W out . We insert the result into Eq. (19), and obtain for the laser linewidth in the highexcitation limit This looks a lot like Eq. (23) for the laser linewidth in the low-pump limit. It differs only by the prefactor of 1/2 in Eq. (44), and by a different expression for the spontaneous-emission factor, which for Eq. (44) reads The result (44) is the same as, for example, in Ref. [65]. Our approach gives a new interpretation to this remarkable result: the linewidths (23) and (44) at low and at high pump rates are different due the different roles of population fluctuations below and above threshold. As shown in the next section, above threshold the population fluctuations reduce the number of photons in S-and increase it in A-combinations, which govern only slow frequency fluctuations and thereby narrow the linewidth. We will also see that the total number of photons n(P ) in the field practically does not depend on population fluctuations.

VI. CALCULATION AND ANALYSIS OF OPTICAL SPECTRA
Equations (31a), (31b) for the symmetric combinationsâ S andv S of the field and polarization and Eq. (31d) for population fluctuations lead to algebraic equations for the Fourier-component operatorsâ S (ω), v S (ω), and δN e (ω). By solving them, we obtain the Fourier-componentâ S (ω) of the S-combinations of the field quadratures. Then, using relation (39) with indices "S" instead of "A", we find the spectrum of the S-combinations as where γ P ≡ γ (P + 1) and is the squared relaxation oscillation frequency at high pump in the rate-equation limit 2κ γ ⊥ [37]. For the derivation of Eq. (46) we used the diffusion coefficients 2D a S a S = 2D a A a A and 2D v S v S = 2D v A v A given by Eq. (40), and the diffusion coefficient We have found that correlations between the Langevin forces representing fluctuations of the carrier population and the polarization give only a very small contribution to the spectrum for a large number of emitters N 0 1 and at weak coupling Ω 0 /(2κ+γ ⊥ ) 1, which is the case considered here. For that reason we shall in the following neglect these correlations and put the corresponding diffusion coefficients to zero (i.e. 2D v A,S Ne = 0). The full spectrum of the lasing mode is n(ω) = n A (ω) + n S (ω) + n AS (ω).
where n A (ω) and n S (ω) are given by Eqs. (41) and (46). So we are left with calculating n AS (ω). By inserting the Fourier expansions (11) forâ S andâ A into the commutator [â A ,â S ] in Eq. (36), we find for n AS (ω) the relation We calculate n AS (ω) approximately by noting that the contribution of n AS (ω) is only important when the number of photons is small, in which case population fluctuations and their contribution to S-combinations are also small. So we calculate n AS (ω) while neglecting population fluctuations. By setting δN e = 0 (that is the case considered in Section III) in Eqs. (31), we see that the equations forâ S andv S are identical to the equations for a A andv A , so that n A (ω) = n S (ω) and therefore where n(ω)| δNe=0 is obtained in Sec. III and given by Eq. (16). By inserting Eq. (50) into Eq. (49) and applying there Eqs. (16) and (41), we obtain the lasing spectrum with n S (ω) still given by Eq. (46). The spectrum n(ω) depends on the population inversion N , which can be found from the energy conservation law (4) in the same way as in Sec. III. We express n [entering n S (ω) through Eqs. (47)] through N by the same energy conservation law (4). Then, the population inversion N is the only unknown variable in Eq. (4). Written in terms of the numbers of photons in the S-and Acombinations, the conservation law becomes Here, n A (N ) is given by Eq. (43), and n S (N ) is found by integrating Eqs. (42) with n S (ω) given by Eq. (46). We find the mean population inversion N by solving the integral equation (52) numerically. In our calculation examples we choose parameters close to typical ones for photonic crystal nanolasers with quantum-dot active media [38]: for the wavelength of the lasing transition we pick λ 0 = 1.55 µm, for the background refractive index n r = 3.3, the cavity mode volume V c = 10(λ 0 /n r ) 3 with N 0 = 100 emitters; a population relaxation rate γ = 10 9 s −1 ; a vacuum Rabi frequency Ω 0 = (d/n r )[ω 0 /(ε 0 V c )] 1/2 with a dipole moment of the lasing transition d = 10 −28 Cm so that Ω 0 = 34γ ; the average atom-lasing mode-coupling factor f = 1/2; finally, we choose the cavity quality factor Q = 1.2 · 10 4 so that 2κ = 100γ .
In the examples below we vary the dephasing rate γ ⊥ and the pump P while keeping all other parameters fixed. The value for γ ⊥ is varied between γ min ⊥ = 50 GHz (so that 2κ/γ min ⊥ = 2) to γ max ⊥ = 1.5 THz (with 2κ/γ max ⊥ = 0.07). This is a realistic region of γ ⊥ for quantum dots [66]. Within this range for γ ⊥ , the conventional beta-factor β varies from 0.98 to 0.6, while the beta-factorβ c varies from 15 to 1.4, so lasers with the chosen parameters have significant amounts of spontaneous emission into the lasing mode.

A. Photon numbers and population inversions
Calculations of mean values of photon numbers and population inversions are helpful for the identification of different lasing regimes and for understanding the role of population fluctuations. Our procedure to find the mean photon numbers and population inversions is as follows: first we calculate n S (N ) by inserting the spectrum n S (ω) from Eq. (46) into Eq. (42). Then we insert n S (N ) found from Eq. (42) and n A (N ) from Eq. (43) into the energy conservation law (52). By solving the latter equation, we can determine the population inversion N . By inserting this N back into n S,A (N ), we find the mean photon number n from Eq. (37).
The red curve in Fig. 5(a) shows the mean photon number n(P ) for a superradiant laser, while the corresponding black curve is for a non-superradiant laser. For low pump rates the superradiant laser exhibits "subradiance" and excitation trapping [22,24], i.e. the photon number is smaller than for a non-superradiant laser. At high pump rates, by contrast, the superradiant laser is seen to generate more photons than the non-superradiant laser. There are two reasons for this: first, excitation trapping becomes weaker as the average emitter population grows and is suppressed when population inversion is achieved; second, in our case the superradiant laser has a smaller polarization relaxation rate γ ⊥ than the non-superradiant laser.
Panel 5(b) shows the population inversion N (P ) for a superradiant laser (red curve). The blue dashed curve depicts N (P ) as found by neglecting population fluctuations, which should be a good approximation for low pump rates. The dashed green curve shows N (P ) in the approximation n ≈ n A , which should be valid for high pump rates. The exact red curve indeed approaches the blue (green) curves at low (high) pump rates.
By following the red curve in Fig. 5(b) and observing where it approaches its asymptotics at low and at high pump rates, we can roughy identify three regions: the LED region at small pump rates, where fluctuations of populations are negligible; the lasing region at high pump rates, when almost all photons reside in the Acombinations of quadratures; and the remaining inter- mediate region between the lasing and the LED regions. These regions are separated by vertical dotted lines in Fig. 5(b). Fig. 6(a) shows the linewidth γ las (P ), given by Eq. (18) with N determined from Eq. (52) for a superradiant laser (red curve). The LED, intermediate and lasing regions shown in Fig. 6 are the same as found in Fig. 5(b). The effects of the same approximations as in Fig. 5 are now shown for the linewidth: the blue dashed curve in Fig. 6(a), given by Eq. (23) for a laser below threshold [57,63], represents the approximation of neglecting population fluctuations, which is again shown to be valid for small pump rates, i.e. in the LED region. The green dashed curve in Fig. 6(a) is given by the laser linewidth Eq. (44) for a laser far above threshold [65] that was found by taking n ≈ n A . As before in Fig. (5), this approximation is shown to be accurate at large pump rates, i.e. in the lasing region. Fig. 6(b) shows the photon numbers (n S,A − 1/4) of the S-and A-combinations of the field quadratures for a superradiant laser. The reason to display (n S,A − 1/4) on the vertical logarithmic axis is that in the low-pump limit P → 0, when n S → n A , we have n S,A → 1/4, in accordance with Eq. (52). We see that population fluctuations (from now on abbreviated as PF) have different influences on n S and n A . Removing δN e from Eqs. (31a)-(31c), we find that these equations are the same for Aand for S-combinations, so that without PF, the n A and n S would be identical and follow the blue dashed curve in Fig. 6 (b). Preserving δN e in Eqs. (31a)-(31c) and plotting n S,A (P ) [the black and red curves in Fig. 6(a)], we see that photon numbers n S,A are hardly affected by PF in the LED region, but PF lead to n S > n A in the beginning of the intermediate region. From the end of the intermediate region and onwards, the PF suppress n S making n S n A in the lasing region, where we also see that indeed n A n, as we anticipated in Sec. V to derive the central result Eq. (44) for the linewidth. Thus, due to PF above threshold, n S is strongly suppressed and n A increased by a factor of two. Fig. 7(a,b) shows the analogous results for a conventional (non-superradiant) laser, to be contrasted with the case of Fig. 6(a,b). The linewidth γ las (P ) of the conventional laser departs more gradually from its asymptotics at low pump to the asymptotics for high pump. An even more conspicuous difference with superradiant lasers is that for the non-superradiant laser n S never exceeds n A . This is an indication that PF have a stronger effect on the superradiant lasers. But n S and n A are both clearly affected by PF for both types of lasers, and for both lasers n S 0 and n A n at high pump rates. By contrast, both for the superradiant laser in Fig. 6(b) and for the conventional laser in Fig. 7(b), we see that the total number of photons practically does not depend on PF.
In the next section we will see that the larger PF of superradiant lasers will make their spectra qualitatively different from those of non-superradiant lasers.

B. Optical spectra
The different shapes of optical spectra for low and for high pump rates reflect different physical effects. It is therefore convenient to consider the spectra for high and for low pump rates separately. Figure 8(a) shows optical spectra n(ω), given by Eq. (51), for a superradiant laser, with 2κ/γ min ⊥ = 2, for high pump rates P ≥ 2. Fig. 8(b) shows the same for a non-superradiant laser, with 2κ/γ ⊥ = 0.2. Only in Fig. 8(a) for the superradiant laser do we see sideband peaks in the spectra. These sideband peaks (or spikes) in n(ω) for Curves 2 to 5 in Fig. 8(a) have the same nature as relaxation oscillations in lasers with 2κ γ ⊥ [36] that are described by rate equations. Sideband peaks appear because the carrier population reacts with some delay to changes in the field and polarization. The delay causes oscillatory energy exchange between the field, polarization and population with a resonance at the relaxation oscillation frequency. For 2κ γ ⊥ , the relaxation oscillation frequency is given by ω ro , as defined in Eq. (47) [36]. Such resonances cause well-known sidebands in the intensity fluctuation spectra [36,41], see also Fig. 10(a,b) below. Analogous sideband peaks due to relaxation os- cillations are not resolved in optical spectra of the nonsuperradiant laser with 2κ/γ ⊥ = 0.2 1 in Fig. 8(b). We attribute this to the fact that such a laser has smaller population fluctuations than a superradiant laser, as we have seen above in the analysis of the mean photon numbers, comparing Figs. 6(b) with 7(b).
For comparison, the dashed versions of Curves 1, 2 and 3 show the corresponding spectra if population fluctuations are neglected. While Curves 1 with and without population fluctuations are practically identical (apart from a small structure in the center that we will discuss below), the dashed and solid Curves 2 and 3 in Fig. 8(a) are qualitatively different: no sideband peaks are observed when population fluctuations are neglected. By contrast, solid and dashed curves are quite close to each other in Fig. 8(b) for the conventional (non-superradiant) laser, where population fluctuations are smaller than in the superradiant laser.
The two-peak structure in n(ω) for Curve 1 in Fig. 8(a) has a different origin than the sideband peaks in curves 2-5 in this Figure. The two broad peaks in Curve 1 are due to collective Rabi splitting (CRS), which occurs when a large number of emitters exhibit Stark shifts in the lasing field. The parameters used for Curve 1 satisfy the conditions for CRS, in particular P < P c , see Sec. III after Eq. (17). For Curves 2 to 5, we have P > P c , so CRS is absent. We described CRS in more detail in Ref. [59], but without taking population fluctuations into account. Figs. 9(a,b) show optical spectra for low pump rates P ≤ 2. The two broad peaks in the curves in Fig. 9 (a) are due to CRS. At these lower pump rates than previously considered in Fig. 8, population fluctuations lead to small features in the center of the optical spectra, namely a peak in Fig. 9(a) and a dip in Fig. 9(b). These small Physically, the photons emitted at the frequencies of the two CRS peaks participate in a nonlinear scattering process: they induce population fluctuations, which couple to the polarization, and thereby have a back action on photon emission. In other words, photons may be re-absorbed and then emitted again by other emitters, which constitutes a nonlinear photon scattering process. Fig. 9 (a) shows that photons from both CRS peaks are re-absorbed and re-emitted most effectively near the center of the spectrum, leading to the small central peak in the optical spectra in Fig. 9(a). Fig. 9(b) shows the modification of the spectrum n(ω) due to the nonlinear photon scattering in a conventional (i.e. non-superradiant) laser, when conditions for CRS are not satisfied and CRS peaks are absent. Here, most emitters are near the center of the spectrum, and photons emitted in this spectral region are absorbed and reemitted by all other emitters away from the center. By nonlinear scattering the energy of the field is thus taken from the center of the spectrum, where we see a dip, and re-emitted far from the center. Indeed we see that away from the dip, the solid curves in Fig. 9(b) lie slightly above the dashed curves, which do not take into account the nonlinear photon scattering.
For increasing pump rates, the main lasing peak grows, and the tiny nonlinear structures in the center of the lasing spectrum disappear.

C. Spectra of intensity fluctuations
The radio-frequency spectrum of laser intensity fluctuations is important for applications of lasers in optical communications and can be measured by direct photodetection of the lasing field as fluctuations of the photocurrent [36]. Here we will derive an approximate expression for the intensity fluctuation spectra of our nanolasers, restricting ourselves to the high-pump limit, where almost all energy resides in the A-combinations of the quadratures.
Slow fluctuations of A-combinations, or equivalently phase fluctuations according to Fig. 4, can be neglected in calculations of photon number fluctuations. In this high-excitation limit, according to Eqs. (28a) and (30), In Eq. (53) the operatorâ S describes broadband amplitude fluctuations, and we takeâ S =â S as in Eq. (34). Then from Eq. (53) we obtain the approximationâ ≈ ( √ n +â S )e iπ/4 , which leads to the photon number (or intensity) fluctuation spectrum where n S (ω) is given by Eq. (46). In Fig. 10, photon number fluctuation spectra are shown for the same parameter values as in Fig. 8. In Figs. 10(a,b) one observes the wellknown relaxation oscillation peaks [37,67]. When γ ⊥ in- Similarly, the sideband peaks in the optical spectra of the superradiant laser shown in Fig. 8(a) disappear upon increasing γ ⊥ , finally arriving to the case of the conventional laser with optical spectra as shown in Fig. 8(b). This similarity confirms that sideband peaks in optical spectra of a superradiant laser are caused by strong population fluctuations leading to strong relaxation oscillations.

VII. DISCUSSION
In this section we discuss our approach to linearize the Maxwell-Bloch equations, followed by a discussion of the main results.
In the high-pump limit, Eqs. (1) can be linearized around their mean values, when the field and polarization are sums of their coherent parts, which is a solution of the semiclassical Equations (5), plus small fluctuations [36,54,55]. Such a linearization is equivalent to the small-signal analysis that is well-known in electrical engineering [68] and which has also been applied to laser rate equations before [37,69]. Rate equations can not be applied for superradiant lasers, where polarization is a dynamical variable. The standard small-signal analysis of Eqs. (1) cannot be applied in the low-pump limit, where the mean values of the lasing field and polarization vanish. In order to linearize the Maxwell-Bloch equations in the lowpump limit, we made a first approximation within our new method in Sec. III by neglecting population fluctuations with respect to the large mean-value population of the large number of emitters. This approximation gives satisfactory results when calculating mean values. An important example is the mean photon number, which is well described both below, near and above the semiclassical laser threshold. Also, the laser linewidth is accurately accounted for, at least below threshold. However, it is one of the main points of this paper that accounting for population fluctuations is necessary in order to correctly account for the linewidth of the laser above threshold, as well as the detailed structure of the laser spectrum below as well as above threshold.
We therefore made an improved analysis, where both the population fluctuations and the dynamics of the material polarization are taken into account. First we considered the high-pump limit in Sec. IV, but again using an approach that differs from and extends the conventional small-signal analysis. The first difference is that we separate the mean values and fluctuations only in the nonlinear terms in the laser equations (1). A further difference is that we use the symmetric (S) and the antisymmetric (A) combinations (30) of the quadratures. In Fig. 4 we showed that S-combinations are related to the amplitude and A-combinations to the phase of the lasing field. We arrived at two linear sets of equations: a set of two equations (31a) and (31c) for A-combinations and a set of three equations (31a), (31b) and (31d) for S-combinations of the quadratures coupled to population fluctuations. Our finding that population fluctuations are explicitly coupled to the S-but not to the A-combinations, is in agreement with the well-known fact [70] that variations in the carrier density change the gain, which subsequently changes the amplitude. If the alpha-parameter (linewidth enhancement factor) is zero, as in our case, the gain changes do not affect the phase.
Next we put forward the hypothesis that the linearized equations (31) as derived in the high-pump limit can be used in the intermediate-and in the low-pump limits as well, where they give satisfactory approximate predictions. First of all, if we neglect population fluctuations, then Eqs. (31) become identical to the equations (10) that we derived for the low-pump limit in Sec. III.
Having discussed our new approach in detail, in the remaining part of this section we will discuss our new results. We demonstrated that the well-known factor 1/2 difference between the laser linewidth Eq. (23) for the low-pump limit and Eq. (44) for the high-pump limit arises due to population fluctuations. If we would neglect population fluctuations, then we would obtain the linewidth (23) for arbitrary pump rates. Furthermore, we calculated the laser linewidth at arbitrary pump rates, taking population fluctuations into account, and demonstrated a smooth transition between the limits for low and high pump rates, as shown in Figs. 6(a) and 7(a).
We calculated the full optical spectrum, given by Eqs. (49) or (51). In the high-pump limit, the optical spectrum features a narrow peaked spectrum n A (ω), with a width determined by the phase fluctuations, on top of a broad background spectrum n S (ω), related to the amplitude fluctuations. Similar sharp "coherent" peaks and "incoherent" wings are seen in experiments, for example in Fig. 1 of Ref. [63] or, for superradiant lasers, in Fig. 4(b) of Ref. [20]. The usual analysis completely neglects this broad background emission in the optical spectra [35].
We found sidebands in the optical spectrum n(ω), shown in Fig. 8(a) for a superradiant laser, characterized by a high beta factor and a low-quality cavity (2κ/γ ⊥ ≥ 1). In this case, population fluctuations are especially important and the laser displays strong relaxation oscillations, which show up as sidebands in the optical spectra.
Another mechanism that gives rise to satellite peaks in optical spectra is well-known for semiconductor lasers. It arises because the refractive index depends on the level of excitation, as described by the so-called α-parameter in the rate equations [33]. However, this is not the mechanism leading to the sideband peaks in our model, where the frequency of the lasing field is exactly on resonance with the lasing transition and the α-parameter is zero.
For weak pumping, when the main lasing peak has not yet appeared, we predict small features in the center of the lasing spectra as displayed in Fig. 9. These small spectral peaks and dips only appear when taking population fluctuations into account. With our present theory we only make qualitative predictions about these small features, because we neglect the second-order correlations (33), although at low excitation these are of the same order of magnitude as the linear terms (34) in Eq. (31d). From calculations that do account for these correlations (33), not shown in this paper, we find that the small central peaks and dips are preserved, but their heights and widths are changed. Below we give an interpretation why the central peak [dip] appears for the lasers parameters of Fig 9(a) [Fig 9(b)].
Population fluctuations (PF) at low excitation have a relatively high spectral power density in the range γ γ ⊥ , 2κ near zero frequency, where the PF are mostly due to the pump and the population decay noise. Far from the center of the spectrum on the other hand, the PF are caused by the interaction with the field and polarization, and the spectral density of PF is relatively small. Since γ is much smaller than the width of the spectra in Fig. 9, the PF can produce either a peak or a dip of width ∼ γ in the center of the optical spectrum.
When collective Rabi splitting (CRS) is absent, as in Fig. 9(b), then the maxima in the linear and nonlinear parts of the polarization coincide and interfere destructively in the frequency region of width γ around the center of the optical spectrum. The linear part of the polarization (i.e. the part that is independent of population fluctuations) induces the birth of a photon. This leads to a PF, which decreases the population of the upper level. Such a PF, in turn, leads to a polarization fluctuation, which suppress the birth of the next photon. This is a reason for a dip in the center of the optical spectra in Fig. 9(b).
The situation is different when the maxima of the linear part of the polarization are not in the center of the optical spectrum, as it happens when collective Rabi splitting occurs. Then the PF in the frequency region near the center of the spectrum are less affected by the lasing field, compared to the case without CRS. The field now resides mostly in CRS peaks far from the center, where the spectral power density of the PF is relatively small. The pump, therefore, increases the PF near the center of the spectrum, leading to a central peak, as shown in Fig. 9(a).
The new sideband peaks and fine structures of spectra of superradiant laser have not been seen in the experiments of Refs. [20-22, 24, 25]. Based on our theory, we predict the kind of lasers and the parameter region, given after Eq. (52), in which to observe such features. The active medium should work in a three-level (effectively two-level) scheme. The conditions for a laser to be superradiant must be satisfied, i.e.β c 1 and 2κ ≥ γ ⊥ . For example, one can reduce γ ⊥ , so that 2κ ≥ γ ⊥ , by lowering the temperature of a nanolaser with q-dots as the active medium [66]. Or vice versa, by increasing the temperature one can go from the superradiant to the nonsuperradiant regime in the same laser. The pump rate must be in the range from γ to 20γ , corresponding to the intermediate region [see Fig. 5(b)], where population fluctuations have their maximal effect, and the sideband peaks in optical spectrum [see Fig. 8(a)] can be observed. Lasers with possibly larger relaxation oscillation peaks in the intensity fluctuation spectrum, as in Fig. 10(a), are good candidates for observing sideband peaks in the optical spectrum as predicted here.
In our theory we have neglected inhomogeneous broadening of the active medium. This implies that for our current theory to be directly appliccable, the actual inhomogeneous broadening must be much smaller than the distance between sideband peaks in the field spectra in Fig. 8(a). It was shown in Ref. [71] that collective effects synchronize emitters, even if the cavity linewidth is smaller than the inhomogeneous broadening. So we expect that some features in the spectra of SR lasers as discussed here will also show up in the presence of inhomogeneous broadening. More precise criteria for the observation of sideband peaks in the optical spectra of superradiant lasers will be derived from detailed investigations of the optical spectrum (51) in the future.

VIII. CONCLUSION
We presented a quantum theory for the spectra and fluctuations of a single-mode homogeneously broadened two-level laser. We developed a new approach, which does not employ the common approximation of eliminating the polarization adiabatically, meaning that the theory also applies to superradiant lasers with a low-quality cavity and a high beta factor. We linearised the equations, solved them and obtained analytical results.
We identify the LED region, where the laser linewidth described by Eq. (23) is wide as in Refs. [57,63], the lasing region with a narrow laser linewidth Eq. (44) in agreement with Ref. [65], and the intermediate region in between them.
Different from the rate-equation approaches of Refs. [33][34][35][36], we describe a smooth transition between the linewidth Eq. (23) in the LED region and the linewidth Eq. (44) in the lasing region, as in Figures 6(a) and 7(a). Furthermore, we calculated optical spectra, including phase and amplitude fluctuations, which respectively lead to a narrow peak and a wide background of the spectra at high excitation, as in Fig. 8.
For superradiant lasers in the intermediate region, we predict two sideband peaks in the optical spectra n(ω), as a consequence of strong relaxation oscillations induced by strong population fluctuations. The sideband structure in n(ω) dissapears for a non-superradiant laser with smaller relaxation oscillations and weaker population fluctuations.
In the LED region, we predict a structure (peak or dip) in the center of the lasing spectrum that is caused by the interference of the linear and nonlinear parts of the polarisation. The interference is constructive, leading to a small peak in the center of the spectrum, if the linear part of polarisation displays collective Rabi splitting [59]. Otherwise the interference is destructive and then we predict a corresponding small dip in the center of the spectrum.
We expect that our approach and results will be useful also for further theoretical and experimental studies of spectra, fluctuations and correlations in lasers for a wide range of parameters, at any pump values and when the active-medium polarisation cannot be adiabatically eliminated, as in superradiant lasers.
In this paper we used our method only to calculate mean values and lasing spectra. In the future, our approach could also be used to calculate higher-order correlations, for example the second-order correlation function g 2 of the lasing field. while the remaining diffusion coefficients forâ,â + ,v,v + all vanish. From Eqs. (A3) -(A5) we obtain the diffusion coefficients (40) of the main text.
Using the operator relationsN eNe =N e andN eNg = 0, we obtain 2D NeNe in Eq. (48) of the main text.

Appendix B: Stationary solutions
In this Appendix we derive expressions for the stationary population inversion and photon number both below and above threshold.
First we consider the situation below threshold. We solve Eq. (21) and find the stationary population inversion N (P ) = N th 2(P + 1) where M b = (P − 1)N 0 /N th − P − 1 −β c and Q b = M 2 b + 8Pβ c N 0 /N th . After that we determine the stationary photon number n(P ) from the energy conservation law Eq. (4) and find n(P ) = 1 4β Next we do the analogous analysis in the highexcitation limit. In this case, almost all photons are in the A-combination of quadratures, so in Eq. (43) we set n A = n and we neglect the constant term 2κ/γ ⊥ with respect to the term ∼ 1/(N th − N ) that gets large when N → N th at high pump. By inserting this photon number n into Eq. (4) we again arrive at a quadratic equation for N , Solving this equation and using the energy conservation law Eq. (4), we find solutions N and n that are again given by Eqs. (B1a), but now with the coefficients M a = (P − 1)N 0 /N th − P − 1, Q a = M 2 a + 2(P + 1)β c (N 0 /N th + 1), instead of M b and Q b , respectively.
We rewrite Eqs. (C2) in terms ofâ A,S ,v A,S andâ A,S , v A,S introduced in Eq. (30). In Eq. (C2c) we use expression (29) for V and replacev x +v p = 2v S ,â x +â p = 2â S . After this we obtain from Eqs. (C2) a A,S = −κâ A,S + Ω 0vA,S +F a A,S , Equation (C3d) depends onâ S andv S , so we must find them. According to Eq. (30), the Fourier-components of a S andv S are zero in the frequency region |ω| < ∆ω/2 and otherwise are identical to the Fourier-components of a S andv S for |ω| > ∆ω/2. Therefore, in order to find a S andv S we solve Eqs. (C3a) and (C3c) by Fouriertransform and set the Fourier componentsâ S (ω) and v S (ω) to zero for |ω| < ∆ω/2. After that the inverse Fourier transform givesâ S andv S . The "cut-off" frequency ∆ω is the same as in the approximation (28) and as shown in Fig. 3.