Stochastically sustained population oscillations in high-beta nanolasers

Non-linear dynamical systems involving small populations of individuals may sustain oscillations in the population densities arising from the discrete changes in population numbers due to random events. By applying these ideas to nanolasers operating with small numbers of emitting dipoles and photons at threshold, we show that such lasers should display photon and dipole population cycles above threshold, which should be observable as a periodic modulation in the second-order correlation function of the nanolaser output. Such a modulation was recently reported in a single-mode vertical-cavity surface-emitting semiconductor laser.


I. INTRODUCTION
Nanolasers, such as semiconductor microcavity lasers [1][2][3] or plasmonic lasers [4][5][6], display cavity volumes of the order of a fraction of a cubic wavelength and thus contains few emitters, while their lasing threshold is attained with very few photons in the cavity, of the order of β −1/2 [7] at any given time, where β is the fraction of spontaneous emission funnelled in the "useful" mode. In a nanocavity, β may be close to 1, so that the number of photons at threshold is very small, typically 1 to 10. For such small numbers, fluctuations are large (of the order of the average values), occur only in integral multiples of a base value (corresponding to a change of a single unit), are asymmetric (as excursions into negative values are impossible) and cause the nanolaser output to deviate from that predicted by the traditional laser rate equations [8,9]. The reason for this deviation is that these real-number differential equations have been developed for conventional macroscopic lasers, in which the discrete nature of the number of emitters and photons may be ignored (termed "the thermodynamic limit" [7]) and are therefore inadequate for describing the operation of nanolasers.
Theoretical methods based on full quantum electrodynamic treatments have been developed [10] to account for phenomena observed in single-atom or few-atom microlasers and micromasers, which have been studied experimentally since the 1980s [11][12][13]. However, in these systems the emitting atoms are relatively well isolated from their environment, have sharp optical spectra, and generally couple strongly to the cavity mode. Semicon- * Electronic address: isabelle.robert@lpn.cnrs.fr ductor (or plasmonic) nanolasers, on the other hand, operate at the opposite limit: The emitters (generally semiconductor quantum dots) are embedded in the cavity material and thus their radiating dipoles interact extensively with their environment, undergoing very fast dephasing processes. At the same time, under the strong pumping conditions required for lasing, the emitters are usually highly-excited, containing several electron-hole pairs (dipoles) and producing very broad emission spectra [14]. These conditions are contrary to the assumptions of the few-atom microlaser theories, as they preclude strong coupling between the emitters and the cavity mode, while the incoherent fluctuation processes cause decay of the atomic polarization and rapidly wash out any phase-related phenomena. In this limit, the full quantum equations reduce to the traditional laser equations that involve only photon and dipole numbers [7]. However, these are continuous-variable differential equations that describe the evolution of the mean photon and dipole populations, but ignore their discrete nature: In each individual realization of the lasing process only integer values of photon and dipole numbers may be obtained, a feature that may be ignored in large lasers but becomes important in nanolasers where these numbers are small. As described below, the traditional laser equations can be adapted to the discrete nature of the photon and dipole populations through discrete simulation methods in which small increments in time are made, and within each increment a discrete but random number of constituents is chosen, obeying the statistics of the physical process being modelled.
This procedure is quite similar to that used in modelling small predator-prey or epidemiological ecosystems [15,16], where it has been shown that "demographic stochasticity", that is population fluctuations due to random birth, death, predation or contagion events, produces cyclical population variations. These cycles are due to the discrete nature of the individuals in the population and cannot be accounted for through traditional deterministic models for the populations, based on realnumber continuous differential equations, such as the Lotka-Volterra equations, whose domain of validity corresponds to the cases in which discretization may be ignored, such as in infinite ecosystems, or when extensive ensemble-averaging is possible.
The traditional method of dealing with fluctuations in lasers consists of inserting randomly-fluctuating Langevin forces in the laser rate equations and linearizing about the steady-state photon and dipole populations [18]. However, this method is not suitable when small populations are involved, since in that case the magnitude of the fluctuations is of the same order as the mean values (and thus linearization is not an adequate approximation) while, at the same time, the small-population statistics deviate significantly from the Gaussian statistics of the Langevin forces. Thus, while the traditional treatment may be adjusted so that the variance of the Langevin forces reproduces the observed variance of the fluctuations, it cannot account adequately for the higherorder moments of their distribution. In a recent publication [9] we measured the higher-order intensity fluctuation statistics of a pulsed nanolaser (up to the fourth order) and, through the implementation of a discrete simulation method, we showed that the discrete nature of the photon number at threshold induces a jitter in the pulse timing, which accounts for the observed photon statistics to all orders.
In this paper, we examine a continuously pumped nanolaser and show that the "granularity" of the photon and dipole populations and the "discretization noise" it introduces lead to large self-sustained oscillations in these populations, often reaching extinction periodically for the case of photons, even above threshold. While this behavior is similar to that of small predator-prey ecosystems, it can also be understood within the traditional laser theory as arising from the laser's relaxation oscillations driven by the discrete jumps associated with the spontaneous and stimulated emission events. The photon population cycles in nanolasers should be observable as a periodic modulation of g (2) (τ ), the intensity autocorrelation function of the nanolaser.
Such a periodic modulation has recently been reported by Takemura et al. [17] in experiments involving a single-mode vertical-cavity surface-emitting laser (VC-SEL). While VCSELs are not nanolasers, the sample examined displayed β ≈ 10 −4 , corresponding to β −1/2 ≈ 100 photons in the cavity at threshold, so that a periodic modulation of g (2) (τ ), albeit of small amplitude, could be measured in the vicinity of the threshold. The discretization noise model accounts well for this modulation, thus validating its use in solid-state nanolasers.

II. DISCRETIZATION NOISE
We consider an optical cavity with one only "useful" (or reference) mode, while all other modes are considered as "leaks". The cavity contains an ensemble of N excited emitters (dipoles) each emitting photons at a spontaneous rate γ . We assume that the emitters interact extensively with their material environment, so that they undergo very rapid dephasing processes (as is the case, for example, in solid-state room-temperature lasers) in a time short compared with the characteristic time of the coupling with the electromagnetic field.
Under these assumptions, the quantum mechanical Maxwell-Bloch equations (or, equivalently, the Jaynes-Cummings model), which describe the dipole-field interaction, may be simplified by adiabatically eliminating the emitter polarization, so that the emission kinetics in such a system is adequately described by the traditional laser rate equations [7]: where P (t) is the incident pumping rate; s(t) is the number of photons in the lasing (reference) mode; 1/Γ c and 1/γ are respectively the lifetime of the reference mode (also termed cavity lifetime) and the dipole lifetime; and β is the fraction of spontaneous emission funnelled into the reference mode. These are continuous-variable differential equations that describe the evolution of the mean photon and dipole populations, but ignore the "granularity", that is the discrete nature of these populations. In any particular realization of the lasing process, the photon and dipole populations can assume only integer values, with discrete transitions between the different population configurations. The physical processes behind these transitions are the excitation of a dipole through "pumping", the emission of a photon by a dipole inside or outside the reference mode, or the escape of a photon from the reference mode to the outside. The discrete nature of the photon and dipole numbers can be taken into account rigorously (and under the same assumptions of rapid loss of phase memory) through the master equation approach, which considers the probability of each realization of the discrete photon and dipole populations, as well as the transitions between the different population configurations [7]. It can therefore account for the full statistics of the emission process. Within the master equation approach the discrete transitions between different population configurations can be regarded as probabilistic partitions of the photon and dipole numbers produced by the processes of pumping, emission, and escape, governed by the statistical parameters of the master equation.
When the master equation is applied to the calculation of the mean photon and dipole numbers ( s and N ), it yields the traditional laser rate equations (1) and (2), under the assumption that the fluctuations in the dipole and photon numbers are uncorrelated ( N s = N s ) [7,8], which is a good approximation when the number of photons and dipoles is large (so that their relative fluctuations are small) or, more precisely, in the "thermodynamic limit", that is when 1/β → ∞.
Calculation of the laser photon statistics by use of the master equation requires calculating at each time-step a matrix of (N max × s max ) interdependent probabilities (where the subscript "max" indicates a cut-off value chosen large enough so as to encompass all population values with a significant probability), a formidable task when the dimensionality of the matrix is larger than a few tens. Alternatively, the master equation may be solved through stochastic simulation methods [8,15], whereby the system evolves from one discrete state to another through a jump event in each time-step, thus describing a trajectory in state-space. The average of a large number of trajectories converges to the solution of the master equation. This greatly simplifies the calculation, as each trajectory is calculated independently of the others, and convergence of the average is reasonably rapid. Beyond the mathematical convenience it offers in the calculation of the master equation, the calculation of individual trajectories permits a direct visualization of the time-evolution of the system parameters and thus provides physical insight into the ongoing processes.
The laser rate equations provide a guideline in the calculation of the stochastic state-space trajectories, since these equations deal with the mean values that may be calculated through the master equation. Thus, at each step, the mean of the stochastic evolution of a given initial state must follow the laser rate equations.
To this end, the rate equations (1) and (2) can be solved numerically by combining an iterative procedure over a time-step dt, small compared with 1/Γ c and 1/γ , and a quantum jump approach to ensure that N and s are integers. The use of a small time-step ensures that this iterative procedure follows well the fine structure of the time-evolution of the system, while at the same time it preserves the correlations between the fluctuations in the photon and dipole populations. Thus, re-writing the laser rate equations, the discrete evolution of N and s follows at each iteration step: (4) This presentation of the laser rate equations permits us to interpret each increment or decrement of N and s as resulting from a random drawing from an appropriate distribution, such that its average is the corresponding expression in these equations. Thus, denoting by P D[µ] a random number out of a Poisson distribution with mean value µ, and by BD[n, p] a random number drawn from a binomial distribution of n trials, each of which yields success with probability p, we can obtain the number of additional dipoles excited by the pumping process during time-step dt as P D[P (t)dt], while N d (t), the number of dipoles having decayed between t and t + dt by spontaneous or stimulated emission, is given by N d (t)=BD[N (t), γ (1 + βs(t))dt]. Each photon emitted by these dipoles has a probability (1 + s(t))β/(1 + βs(t)) to be funneled in the cavity mode, and thus the number of emitted photons entering the cavity mode is BD[N d (t), (1 + s(t))β/(1 + βs(t))]. Finally, the number of photons that escape to the outside due to cavity losses, is drawn from the binomial distribution BD[s(t), Γ c dt]. This procedure accounts for the evolution of the photon and dipole populations as integer numbers. It thus takes account of the "discretization noise" without any need to introduce it in the continuous-variable laser rate equations externally. Traditionally, noise (including quantum fluctuations) is re-introduced in the laser equations (1) and (2) by linearizing around a (large) mean value of the populations and adding a Langevin driving force with Gaussian statistics [18]. However, this approach is valid only in the case of large populations [15], where the relative fluctuations are small and obey Gaussian statistics thanks to the central limit theorem.

III. POPULATION OSCILLATIONS IN NANOLASERS
To gain an insight into the effects of discretization noise, we have calculated state-space trajectories at different pumping powers, by solving Equations (3) and (4) numerically through the above procedure, using the values of β = 0.25, 1/γ = 400 ps and 1/Γ c = 10 ps (corresponding to a quality factor of 13000 at λ = 1.5µm), which are typical for quantum dot nanolasers [1-3]. Figure 1 (left) displays the time-dependence of the number of intra-cavity photons and excited dipoles at a pumping power of 4.5P th , where P th = Γ c /β is the nominal threshold pumping power. At this pumping power, the mean number of dipoles has reached its pinning value of Γ c /βγ ≈ 160, yet the photon and dipole populations display well-defined oscillations of the same frequency, as can be seen in the corresponding power spectra (Fig.  1 right). These oscillations are synchronized (but in quadrature), in a manner analogous to the population cycles driven by demographic stochasticity observed in finite predator-prey ecosystems [15].
The synchronization of the photon and dipole population cycles as a function of pumping power can be studied through the covariance of these two quantities, as can be seen in Fig. 2 (left) (black continuous line), where the covariance curve for a low-β laser (red continuous line) is also given for comparison. Below threshold, the photon and dipole populations fluctuate randomly and independently of each other, as the presence of less than one photon (on average) in the cavity produces no change in the dipole emission dynamics and spontaneous emission outside the cavity mode destroys any correlation between these two populations. At higher powers, as stimulated emission grows, a larger fraction of the emission goes into the reference mode, so that any decrease in the dipole population is matched by a synchronous increase of the intra-cavity photon population, thus producing a strong negative correlation. The difference between highand low-β lasers is essentially the speed with which the rescaled covariance reaches its limiting value of -1. At a given pumping power (normalized with respect to threshold), a high-β laser has a smaller number of intra-cavity photons, so that its stimulated emission rate is smaller implying that a smaller fraction of the emission goes into the reference mode. As a consequence, the rescaled covariance drops more gradually to its limiting value of -1.
The power dependance of the covariance (Fig. 2 left) and of the central frequency of the population oscillations (Fig. 2 right) are quite similar to those calculated through the traditional linearized ocillation relaxation theory [19] ( Fig. 2 dotted lines), indicating that the photon and dipole population cycles in a nanolaser can be understood also as arising from relaxation oscillations driven by the noise due to the granularity of photon emission. This mechanism is similar to that of the origin of the excess linewidth in semiconductor lasers [19] and of the excess noise in class-B (i.e. Γ c > γ ) lasers [18], however with an important difference. In traditional low-β lasers, the (spontaneous and stimulated) emission events are distributed essentially uniformly in time, corresponding to a white noise spectrum that gets filtered by the transfer function of the system to produce a damped oscillation. On the other hand, in high-β lasers, because of the small mean value of the intra-cavity photons, emission events are more frequent at the top of the photon population cycle, and thus constitute a periodic driving force that produces a resonant amplification sustaining the photon and dipole population cycles.
The photon population cycles should be directly observable in the second-order photon correlation function g (2) (τ ), which is measurable experimentally through a Hanbury Brown and Twiss setup. As can be seen on Fig.  FIG. 3: Second-order photon correlation function g (2) (τ ) of a nanolaser (β = 0.25) undergoing photon population cycles, for P = 0.5P th (red line) and P = 4.5P th (black line). The dashed horizontal line corresponds to g (2) (τ ) = 1, the value in continuous-wave lasers.
3 (red line), below threshold g (2) (τ ) displays the usual behavior for chaotic light. That is, when the two detected photons are in temporal coincidence it has the value g (2) (0) = 2, while when the two photons are separated by a time interval τ the value of g (2) (τ ) drops to 1 exponentially e −Γcτ , where Γ c is the cavity energy decay rate. This is the familiar Siegert relation, where g (1) (τ ) is the first-order correlation function and corresponds to the Fourier tranform of the optical spectrum. On the other hand, above threshold, the presence of photon population cycles in nanolasers causes the value of g (2) (0) to be larger than 1, in contrast to traditional lasers which display g (2) (τ ) = 1 for all time intervals τ . At finite time-intervals, g (2) (τ ) tends towards the value of 1, in a damped oscillation fashion dropping periodically below 1 for a few cycles (see Fig. 3 black line). This feature arises from the presence of population cycles in nanolasers (even under constant pumping) and is a deviation from the Siegert relation (Eq. (5)), which does not pemit g (2) (τ ) values smaller than 1.

IV. OBSERVATION OF POPULATION OSCILLATIONS
In a recent publication [17], observation of a fast periodic modulation in the photon correlation g (2) (τ ) was reported for a single-mode vertical-cavity surface-emitting semiconductor laser displaying a spontaneous emission coupling factor of β ≈ 10 −4 (as can be seen on Fig. 2 of Ref. [17]). The modulation was clearly observable at a pumping current of 2.47 mA, corresponding to a power of 1 % above threshold and ∼ 150 photons in the lasing mode. However, the modulation was weaker at a pumping current of 2.72 mA (11 % above threshold and ∼ 1500 photons in the lasing mode) and was below instrumental sensitivity at currents of 3.00 mA (23 % above threshold and ∼ 2300 photons in the lasing mode). It was attributed by the authors to relaxation oscillations excited by (generic) noise in the laser.
We argue that in that experiment the oscillations are driven by discretization noise, as our model can fit the experimental second order autocorrelation function g (2) (τ ), with no adjustable parameters other than the cavity lifetime 1/Γ c = 3.4 ps and the spontaneous emission lifetime 1/γ = 4 ns extracted from [17]. Figure 4 shows the experimental second order autocorrelation function g (2) (τ ) (black points) for a pump current of I = 2.6 mA (corresponding to a pumping power 6.5 % above threshold and ∼ 1000 photons in the cavity) and the fit using the discretization noise model is represented as a continuous red line.

FIG. 4:
Second order autocorrelation function g (2) (τ ) as function of time. Black dots: experimental data from [17]. Red continuous line: simulation through the discretization noise model.

V. POPULATION OSCILLATIONS AND NOISE
For small-β lasers (in which the number of photons at threshold is large), or when the nanolaser response is averaged over times much longer than the relaxation oscillation period (and thus corresponds to an ensembleaveraged response), the population cycles arising from granularity appear as fluctuations around the (mean) steady-state populations. These fluctuations are often characterized by the Fano factor [18], where σ 2 s is the variance in the photon number of mean values. By analogy to phase transitions, it is often expected that fluctuations should be maximum at threshold  [7,20], as can be seen in Fig. 5 for a low-β laser. However, for class-B lasers (Γ c > γ ) with large β, the Fano factor maximum occurs at powers significantly higher than the stimulated emission threshold (see Fig. 5 left), as shown by van Druten et al. [18]. While this discrepancy between the laser threshold and the fluctuation maximum has been interpreted as an indication that such lasers produce partially chaotic light [21] and become coherent only at pumping powers much higher than that of the Fano maximum, our results indicate that the increase of the Fano factor not due to the presence of a chaotic component in the laser output but arises from the large amplitude of the photon population cycles that are excited and sustained by discretization noise.
An alternative designation of the magnitude of the fluctuations, useful in engineering, is the normalized rootmean-square noise (nRMS), which is related to g (2) (0) by nRM S = σ 2 s s = g (2) (0) − 1 + 1 s As can be seen on Fig. 5 (right), for traditional low-β lasers and above threshold, the nRMS noise drops rapidly to very small values. As the signal-to-noise ratio is essentially the inverse of nRMS, this feature makes the output of low-β lasers an ideal carrier for high-fidelity information transmission. For high-β nanolasers, on the other hand, the nRMS noise remains quite high even at pumping powers much above threshold, both because of the presence of population cycles (which cause a large value for g (2) (0)) and the relatively small value ofs. As the signal-to-noise ratio is of the order of 1 even high above threshold, this implies that such lasers are ill-adapted to the transmission of information, unless particular noisereducing protocols are implemented taking into account the population cycles.

VI. CONCLUSION
In conclusion, the output of a continuously-pumped high-β nanolaser cannot obey the traditional (continuous, real-number) laser equations, because of the small numbers of intra-cavity photons and dipoles involved, which bring forth the discrete nature of these small populations. Such nanolasers should display population cycles both for the intra-cavity photons and the excited dipoles because of the intrinsic stochasticity of the photon and dipole granularity, analogous to what is observed in small ecosystems. These cycles can be understood also as resulting from relaxation oscillations driven synchronously by the sudden and discrete jumps in photon and dipole numbers, due to the elementary emission processes. The photon population cycles are observable as a periodic modulation of the second-order correlation function of the laser, as was recently reported by Takemura et al. [17]. The presence of population cycles in high-β nanolasers may appear as excess noise which may hinder information processing and transmission with such lasers.