Efficient stochastic simulation of rate equations and photon statistics of nanolasers

Based on a rate equation model for single-mode two-level lasers, two algorithms for stochastically simulating the dynamics and steady-state behaviour of micro- and nanolasers are described in detail. Both methods lead to steady-state photon numbers and statistics characteristic of lasers, but one of the algorithms is shown to be significantly more efficient. This algorithm, known as Gillespie's First Reaction Method (FRM), gives up to a thousandfold reduction in computation time compared to earlier algorithms, while also circumventing numerical issues regarding time-increment size and ordering of events. The FRM is used to examine intra-cavity photon distributions, and it is found that the numerical results follow the analytics exactly. Finally, the FRM is applied to a set of slightly altered rate equations, and it is shown that both the analytical and numerical results exhibit features that are typically associated with the presence of strong inter-emitter correlations in nanolasers.


I. INTRODUCTION
Optical cavities on the micro-and nanometer scale can reduce the number of available modes for light emission and increase the coupling of spontaneously emitted light into the cavity mode(s) [1]. This can be useful for laser devices, as it allows for low power consumption and high modulation speeds [2][3][4][5]. These features, combined with the small footprint of the devices, make them promising candidates for on-chip optical interconnects [6], as well as for many other uses, like chemical and biochemical sensing and detection [7][8][9][10]. Therefore, micro-and nanolasers and their properties are currently a prominent subject of research.
One major area of interest in regards to nanolasers is the photon noise and photon statistics: With a relatively small number of intra-cavity photons and a large fraction of the spontaneously emitted photons ending up in the cavity mode(s), the associated quantum noise becomes increasingly important [11][12][13][14][15]. Recently, stochastic methods have been used to simulate nanoscale lasers with large values of the spontaneous emission β-factor [16][17][18], and it was found that the noise and statistics of nanolasers can be captured surprisingly well by including only shot noise; that is, the noise associated with the discreteness of the photons and emitters in the cavity.
In this work we will show that another stochastic approach, often used in the fields of chemistry and biochemistry [19][20][21][22], can be used to describe the laser dynamics, and can decrease the computation times by several orders of magnitude compared to previous approaches [16][17][18]. Additionally, certain numerical assumptions made in [18] can be avoided with the new method, leading to a robust and versatile stochastic approach with numerous applications. Here we use it to consider the intra-cavity * mwubs@fotonik.dtu.dk photon probability density functions as a lasing system transitions from thermal to coherent emission, and we will investigate the prospects of introducing effectively altered rates of emission in the rate equations as a way to qualitatively describe the effects of collective interemitter correlations.
The paper is structured as follows: In Section II we will introduce the basic laser model that we consider in most of the paper, and give a brief overview of the laser rate equations. In Section III we will review and discuss the stochastic simulation algorithm used in [18] and then introduce another, more efficient algorithm for stochastically simulating the laser rate equations. Using the new algorithm, in Section IV we will examine the intra-cavity photon distributions of the laser and in Section V we will examine the effects of introducing an effective asymmetry between the spontaneous and stimulated emission rates. Finally, we will summarize and discuss our results in Section VI.

II. LASER MODEL AND RATE EQUATIONS
We consider a model for a nanolaser consisting of n 0 distinct two-level emitters interacting with a single cavity-mode, as shown schematically in Fig. 1. We assume for simplicity that all emitters couple to the cavity mode with the same coupling strength, so it it is convenient to work with the total level populations n e and n g for the excited and ground levels, respectively. Since we consider two-level emitters, we have n e + n g = n 0 . The intra-cavity photon number, i.e. the population of the cavity mode, will be denoted by n p . The cavity mode has a decay rate per photon given by γ c , and the emitters have a rate of loss into non-lasing modes given by γ bg . The rate per emitter of emission into the lasing mode is denoted γ r , giving a total decay rate γ t = γ r + γ bg from the excited state. The dependence of γ r on material and device parameters is discussed in [18], and includes Pur-

Event type
aµ aµ aµ Average rate Stimulated Em.
In this work we will focus on so-called Class B lasers, where the dephasing rate γ 2 of the emitters is much larger than the photon loss rate γ c , so that the emitter polarization can be adiabatically eliminated from the analysis [23]. For such devices, the laser dynamics may be described through rate equations for the photon number n p and the number of excited emitters n e [2,18,24,25]: Here F p and F e are stochastic, zero-mean Langevin forces taking into account the shot noise in the particle numbers due to the discrete nature of photons and excited emitters. The laser rate equations (1)-(2) give a simple and intuitively pleasing description of laser dynamics, since each term can be easily associated with a process occurring in the laser, as summarized in Table I. Note that the pumping term explicitly takes into account the effect of pump blocking [18]. The spontaneous emission β-factor is defined as the fraction of spontaneously emitted photons which end up in the cavity mode, i.e.
The rate equations (1)-(2) can be solved analytically in steady state, yielding expressions for the mean photon number n p and number n e of excited emitters as functions of the pump rate, γ p n 0 . These steady-state solutions can exhibit characteristic features of laser operation, namely emitter population clamping at large pump rates and a sudden jump in the photon number as a function of the pump rate. These features are characteristic of a system crossing the lasing threshold, but for certain parameter values, typically associated with β → 1, these features become less distinct, and instead the photon statistics are used to characterise the transition to lasing [11][12][13][14]. The steady-state photon statistics of a laser are typically measured by the first two statistical moments of the photon number, or the mean value n p and the variance ∆n 2 p = n 2 p − n p 2 . They are often summarised using the zero-delay second-order photon auto-correlation function g (2) (0) [26,27]. This can be expressed through the relative intensity noise (RIN) [27] as (4) Lasing is then characterised by a transition from thermal (LED) light, for which g (2) (0) = 2, to coherent (laser) light, for which g(2)(0) ≈ 1.

III. STOCHASTIC SIMULATION
In [18] it was shown that stochastic simulations based on the algorithm presented in [17] agreed very well with analytical results obtained from a small-signal analysis of the rate equations that is accurate even down to a few emitters (n 0 ≥ 10). In this section we will start by reviewing the specific method used for the stochastic laser dynamics simulations in [18], and then we will introduce a different method to produce the same results with significantly increased efficiency.

A. Fixed Time-Increment Stochastic Algorithm
The basic idea behind the stochastic simulations is to assume that the laser rate equations (1)-(2) describe a collection of discrete particles, namely photons and excited emitters, the numbers of which fluctuate due to several processes: loss of particles (cavity or background losses), exchange of particles (emission or absorption) or addition of new particles (pumping). Individual events occur randomly, but their average rates are given by the terms on the right-hand sides of eqs. (1)-(2), as summarized in Table I.
In Appendix A, we give a more detailed derivation of the simulation algorithm used in [18], but essentially it can be boiled down to the iteration over many time increments dt, where one asks the question "how many events of type µ happened during dt?" for each event type: stimulated emission, absorption, spontaneous emission, cavity loss, background loss and pumping. In general, this is not an easy question to answer, since each event that occurs will immediately change the current number of photons and/or excited emitters; and this in turn affects the rates a µ and the probability of subsequent events. However, in [18] a few simplifying assumptions were made, which make it possible to approximate the answer: It was assumed that there is some order in which event types happen, in the sense that all the events of a particular type occur together, before all events of the next type; it was assumed that one occurrence of an event does not affect the probability of any other events of the same type occurring; and it was assumed that the binomial distributions for the number of occurring events can be approximated by Poisson distributions. These assumptions are all reasonable when the value of dt is sufficiently small. The algorithm is similar to the so-called τ -leap method, described in detail in [22,28].
Explicitly, the algorithm can be written as follows: Fixed Time-Increment Method for Laser Rate Equations (FTI) 1. Initialize: Set system parameters, set initial number of photons n p (t = 0) and excited emitters n e (t = 0).

2.
Calculate rates a µ for all event types µ according to Table I using current particle numbers.
3. Determine the number of each type of event µ occuring during dt by random draws from Poisson distributions with parameters a µ dt.
4. Update n p and n e accordingly, set t → t + dt.

5.
If the maximal simulated time has been reached, end. Otherwise, go to step 2.
Solving the laser rate equations with this algorithm leads to equidistant time series for the number of intracavity photons and excited emitters, (n p (t i ), n e (t i )), and statistical convergence is ensured once reduction of the time increment dt no longer affects the steady-state mean photon number n p , excited-emitter number n e , and the photon variance ∆n 2 p . In practice, the size of dt can be chosen at a given pump rate as some fraction dt frac of the smallest reciprocal rate appearing in the rate equations (1)-(2), computed using the analytical steadystate values for the particle numbers. For instance, if the mean steady-state rate of spontaneous emission is the largest among the steady-state rates in eqs. (1)-(2), then dt = dt frac /γ r n e . Convergence is typically obtained for dt frac on the order of 10 −2 or less. The mean photon and emitter numbers obtained using the FTI show the behaviour expected for lasers, and the steady-state photon statistics, as computed using eqs. (4), also exhibit the characteristic transition from chaotic to coherent light. This is shown in Fig. 3 of [18].
While the algorithm is stable and very easy to implement, it also has a few downsides. First, the computationally expensive act of drawing six Poisson-distributed random numbers at each iteration leads to long computation times, since many iterations are needed to obtain reliable photon statistics. Second, the assumption that all events of one type happen independently of all other event types introduces some ambiguity in the ordering of the events which happen during each time-increment: For instance, if there are n e (t) excited emitters at time t and one is de-excited due to spontaneous emission into the lasing mode, then the number of excited emitters available for stimulated emission should in principle be n e (t) − 1. However, if the stimulated emission event occurred first, then the situation is opposite. This ordering issue may only play a small role in the outcome of the simulations as long as the time increment dt is chosen sufficiently small, but this again leads to the question of how small is sufficient. All these challenges can be ameliorated by use of another stochastic simulation algorithm, which will be described in more detail below.

B. Gillespie's First Reaction Method for the Laser Rate Equations
To obtain a faster algorithm, which has the added benefit of avoiding ordering issues as well as the need to define a time-step size, we can change the fundamental question asked at each iteration to "How long time until the next event occurs?" and "Which type of event will occur?". In this way, the size of the time increment at each iteration will differ, and only one event will happen during each time increment. Changing the basic point of view in this way corresponds to applying the algorithm known as Gillespie's First Reaction Method (FRM) [19], which is well-known in the chemistry and biochemistry communities [22]. It is used in numerical calculations regarding chemical reactions involving several species of molecules with finite populations, assuming that each type of reaction that can occur is a stochastic Markov process. The rate equations (1)-(2) describe an analogous type of system, where the photons and excited emitters are analogous to the particles, and the events of emission, absorption, loss and pumping are analogous to the reactions that change the particle populations.
Changing the operational question of the simulation in this way means that instead of generating six integers corresponding to the change in particle numbers due to the six different event types, we are interested in generating two numbers: One corresponding to the time until the next event, and one corresponding to the type of event. As shown in [19], one way of doing this involves assuming for each of the different event types that it happens before any of the others, and then generating a tentative time τ µ until this event occurs. Once all the tentative times τ µ have been determined, we can choose the shortest, τ µ0 , as the actual time until the first event (or reaction) occurs, and the corresponding type of event, µ 0 , as the event that occurs. The tentative time τ µ can be efficiently generated by drawing a random number from an exponential distribution with parameter a µ , which is proved in Appendix A and in [19].
Specifically, the algorithm can be stated as follows: Gillespie's First Reaction Method for Laser Rate Equations (FRM) 1. Initialize: Set system parameters, set initial number of photons n p (t = 0) and excited emitters n e (t = 0).
2. Calculate rates a µ for all event types µ according to Table I using current particle numbers. 3. For each µ, generate a tentative time increment τ µ by a random draw from an exponential distribution with parameter a µ .
5. Update n p and n e according to event type µ 0 , set t → t + τ µ0 .
6. If the maximal simulated time has been reached, end. Otherwise go to step 2.
An important property of the FRM is that there are no numerical parameters, like the time increment dt in the FTI algorithm, and therefore there is no need for checks of convergence. In this sense, the FRM is an exact method to stochastically simulate the laser dynamics. In addition, exponentially distributed random numbers can be quickly and inexpensively generated using uniformly distributed random numbers on the unit interval (0, 1) [29]. Therefore, step 3 in the FRM involves six draws from a uniform distribution, instead of the six draws from six different Poisson distributions needed in the FTI algorithm. This drastically reduces the computation time of the FRM compared to the FTI, while the end results for the steady-state mean particle numbers and photon statistics are the same for both algorithms. This is shown in Fig. 2 by reproducing the results of Fig. 3 in [18] through the use of the FRM, and comparing computation times versus simulated times for the two algorithms.
It is clear from Fig. 2(d) that the computation time increases approximately linearly with the simulated time for both the FTI and FRM algorithms, as expected. However, as illustrated, the computation time at a fixed simulated time also increases for increasing pump rates. This is because the time delays between events become smaller at higher pump rates, as the larger particle numbers lead to larger rates of the different events of emission, absorption and loss. Smaller time increments imply more iterations needed to obtain the same simulated time, and hence larger computation times. In all cases, it is clear that the FRM algorithm leads to significantly shorter computation times at all pump rates: Indeed, improvements of several orders of magnitude can be obtained by using the FRM algorithm instead of the FTI algorithm.
We see that the computation time for the FRM algorithm in the low-pump case is approximately constant for the shortest simulated times. This seeming lower bound on the computation time is related to the time for initialization rather than having to do with the specific algorithm.
Since the FRM significantly reduces computation times, new features of the laser dynamics can easily be studied, which would otherwise be extremely time consuming. One such feature is the evolution of the photon number distributions within the laser cavity as the system transitions from chaotic, LED-like behaviour to coherent, laser-like behaviour. This will be examined in the next section.

IV. PHOTON DISTRIBUTIONS
Once the laser settles into steady-state operation, the mean values for observables no longer change, but the particle numbers still fluctuate in time. By examining the simulated time series of the intra-cavity photon number in steady state, we can obtain the discrete photon probability density function (PDF) by evaluating how much of the total simulation time is spent with 1 photon in the cavity, how much is spent with 2 photons in the cavity, etc. In Fig. 3 the photon PDF obtained from the stochastic simulations is shown for three different values of the pump rate using parameters corresponding to the cases of β = 1 and β = 10 −3 from Fig. 2. In both cases there is a clear transition from thermal statistics giving a Bose-Einstein distribution (d,g), through an intermediate distribution near threshold (e,h), to near-Poissonian statistics (f,i).
As shown in Fig. 3(b), the photon auto-correlation g (2) (0) seems to indicate that Poissonian statistics are obtained at large pump rates for both β = 1 and β = 10 −3 , the simulated values being g (2) (0) = 1.0015 and g (2) (0) = 1.0001, respectively. However, we see in (f,i) that the photon PDF approaches a Gaussian distribution at high pump, rather than a Poisson distribution. Mathematically, we can explain this by the photon variance approaching a value slightly different from the characteristic Poissonian ∆n 2 p = n p . Equivalently, the Fano factor F = ∆n 2 p / n p does not reach a value of 1 at high pump rates. From the analytical expression in [18] for the photon variance we can show that at large values of the mean photon number we have F = ∆n 2 p n p ≈ 1 + c, with c = 1 + n th /n 0 2(n 0 /n th − 1) , for n p ≫ 1. Here n th ≡ γ c /γ r is the semi-classical threshold value of the emitter population inversion. The result in eq. (5) is illustrated in Fig. 3(c), where we have plotted the Fano factor versus the pump rate, with the constant c = 0.75 for both values of β. Physically speaking, this deviation from exact Poissonian statistics has its roots in the laser level-scheme: Similar, non-Poissonian photon variances for lasers above threshold have been reported in earlier theoretical and experimental work [30][31][32], where two-, three-and four-level lasers are considered. In particular it is argued in [25,30,31] that the two-level emitter scheme we consider here should generally lead to super-Poissonian statistics above threshold due to depletion of the emitter ground states. These effects are what we see in our analytical and numerical results.
With the FRM algorithm it is possible to obtain photon distributions which match the analytical predictions almost perfectly, while the computation time is kept relatively short -less than 2 hours on a commercially available laptop. This demonstrates once more the extreme efficiency and exactness of the FRM applied to the laser rate equations. In the next section we will use this to show how small alterations in the laser rate equations can lead to results which exhibit features that are typically obtained from much more intricate laser models.

V. BREAKING THE SYMMETRY BETWEEN SPONTANEOUS AND STIMULATED EMISSION RATES
In the traditional rate-equation description of the laser dynamics, the rate per emitter of spontaneous emission into the cavity mode is the same as the rate per emitter per photon of stimulated emission into the cavity mode, consistent with the Einstein relations [27]. This can be seen in eqs. (1)-(2), where the terms related to spontaneous and stimulated emission and absorption are all proportional to the rate of radiative decay, γ r . However, it is possible to imagine that certain physical effects, which are not accounted for in the traditional derivation of the rate equations, could give rise to an effective difference between the rates of spontaneous and stimulated emission. For instance, collective effects due to correlations building up between emitters can affect the dynamics and steady-state behaviour of lasers, as has been shown in recent theoretical and experimental works [13,[33][34][35][36][37][38]. Some of these phenomena could potentially be captured by an asymmetry between the rates of spontaneous and stimulated emission, which is what we will investigate next.
To examine the effects of introducing a difference between the spontaneous and stimulated radiative events, we can replace the common rate γ r in eqs. (1)-(2) by two separate rates, γ sp r and γ st r , for the spontaneous and stimulated events, respectively. The ratio quantifies the difference between the rates of spontaneous and stimulated emission (and absorption), and the rate equations (1)-(2) can be rewritten using this quantity. We can solve the rate equations including ξ in the same way as we did for the traditional rate equations (1)-(2), both analytically and numerically. If we perform the same small-signal analysis applied in [18] with the new radiative rates, we can also obtain corresponding ξ-dependent analytical expressions for ∆n 2 p , RIN and g (2) (0). Incorporating the asymmetry in the stochastic simulations using the FRM is simple, since we only need to change the rates a µ to include ξ.
In Fig. 4, we have plotted the steady-state mean intra-cavity photon number and photon auto-correlation as functions of the pump rate for three different values of ξ, illustrating the effects of increasing or decreasing the spontaneous emission rate relative to the stimulated emission rate. We see several effects: For small pump rates, the reduction of the rate of spontaneous emission leads to a lower mean number of photons in the cavity as well as an increased value of the photon autocorrelation function to super-thermal values. Conversely, increasing the spontaneous emission rate leads to a larger mean intra-cavity photon number and sub-thermal photon statistics in the low-pump limit. In both cases, the results of the stochastic simulations follow the analytical results very well. In fact, we find that there is a simple relationship between the value of ξ and the value of the photon auto-correlation function at low pump rates, namely These effects occur in part because the alteration of the ratio ξ from 1 leads to changes in one of the derived parameters β = γ sp r /(γ sp r + γ bg ) or n th = γ c /γ st r . Note that the well-known result g (2) (0) = 2 at low pump rates is restored in the case where the spontaneous and stimulated emission rates are equal, ξ = γ st r /γ sp r = 1. From Fig. 4 we see that decreasing the rate of sponta- neous emission or increasing the rate of stimulated emission lead to the same qualitative effect, namely a reduced mean photon number and super-thermal photon statistics. Correspondingly, increasing the rate of spontaneous emission or decreasing the rate of stimulated emission both lead to a larger mean photon number and subthermal photon statistics. We can understand these effects qualitatively as follows: Less spontaneous emission will lead to fewer photons in the cavity in the low-pump regime where spontaneous emission dominates. This means that there will be more excited emitters which can be prompted to emit through stimulated emission by the few photons in the cavity, which in turn leads to a higher probability of several photons being present at the same time, hence a larger auto-correlation. Equivalently, an increased rate of stimulated events (stimulated emission, absorption) means that any emitted photons are more likely to be reabsorbed, giving fewer photons on average, while each photon has a larger probability to set off a stimulated emission from any excited emitter, resulting in a larger g (2) (0). The converse situation with a larger rate of spontaneous emission / lower rate of stimulated emission and absorption can be understood in a similar way.
From recent experiments and theories regarding nanolasers with collective correlations between the emitters, it is known that collective effects in small lasers can result in a reduced mean number of photons in the cavity and super-thermal photon statistics in the low-pump limit [13,37,38]. Fig. 4 shows that these phenomena also appear from the laser rate equations if one breaks the symmetry between the spontaneous and stimulated emission rates dictated by the Einstein relations. The effects are typically described theoretically by the emitter cross-correlations causing so-called excitation trapping at small pump rates [37,38], meaning excitations in the system preferably reside in the emitters rather than the cavity mode. This effect, which is related to subradiance [39,40], in turn gives rise to bursts of photon emission, similarly to what we described above in conjunction with the reduced rate of spontaneous emission. To see that such bursts of photons occur in our simulations for ξ < 1, we have plotted in Fig. 5 the timeresolved photon number for two values of ξ. As shown in Fig. 5(a,b), we examine the case where the mean photon number is the same, while the photon distributions and statistics are either thermal or super-thermal. From (c,d) it is clear that the case where ξ < 1, i.e. the rate of spontaneous emission is reduced, the photons often come in larger bursts, giving rise to the super-thermal statistics that we observe. This corresponds well with the theoretical descriptions given of the effects of emitter-emitter correlations in nanolasers, meaning an effectively reduced rate of spontaneous emission may provide a simple and intuitive way of further studying and understanding collective effects in nanolasers.

VI. SUMMARY AND DISCUSSION
In summary, we have presented Gillespie's First Reaction Method (FRM) in detail, and applied it to stochastic simulations of the traditional laser rate equations. We have shown that the FRM reduces computation times by around three orders of magnitude compared to earlier algorithms [18], while it also avoids potential numerical issues related to the size of time-increments and the ordering of events. Thus, the FRM is an efficient and exact algorithm for stochastic simulations of the laser rate equations. We have used the FRM to examine the intracavity photon distribution as the laser transitions from chaotic to coherent emission, and we have shown that the numerical results for the photon statistics follow the analytical results exactly. We have also applied the FRM to altered laser rate equations, in which the symmetry between the spontaneous and stimulated emission rates was broken. We found that by reducing the rate of spontaneous emission, or increasing the rate of stimulated emission, both the numerical and analytical results show features that are also found in theories and experiments regarding lasers with strong emitter-emitter correlations. This suggests that there may be a way to use a slightly altered version of the well-known laser rate equations to describe collective effects in micro-and nanolasers.
The FRM described here increases the efficiency of the stochastic simulations tremendously compared to methods like the FTI, while remaining exact and conceptually simple. One of the reasons why the exact FRM algorithm is more efficient than the approximate FTI algorithm for the laser rate equations is that there are relatively few event types and particles involved. Another reason is the conservative choice of time-increment size used in the FTI to obtain accurate results: In the FTI, and in the similar τ -leap method [22,28], faster computation times can always be obtained by choosing the time increments to be larger, but this naturally reduces the accuracy of the simulation. In the low-pump limit, where the photon statistics are very sensitive to the particle numbers, high accuracy is needed, implying small time-increments and hence long computation times. In the high-pump limit, where the particle numbers are large, it is possible that a less conservative choice of time-increment sizes could reduce the computation times of the FTI algorithm without significantly affecting the accuracy. We leave an investigation of this to future work.
Another exact and conceptually simple algorithm for stochastic simulation is Gillespie's Direct Method [19,20,41], which is similar to the FRM, but requires only two random draws per iteration. This could potentially increase the efficiency of the simulation further, but due to the relatively low number of event types in the rate equations, the difference in computation time is expected to be minor. We performed a few simple simulation tests using Gillespie's Direct Method, and for the simulated times needed to obtain reliable statistics, the preliminary results showed practically no difference in computation time compared to the FRM. Both the FRM and the Direct Method have been expanded upon, leading to even more efficient algorithms [21,22]. However, in most cases the increase in efficiency happens at the cost of the conceptual simplicity. For instance, the Next Reaction Method introduced in [21] reduces the computation time by storing and reusing the tentative times τ µ and only updating certain rates a µ , but it requires the introduction of additional concepts like dependency graphs and indexed priority queues. Since the rate equations for nanolasers deals with relatively few particles and event types, this added complexity would likely lead to fairly small improvements in efficiency. However, such algorithms could be implemented in stochastic simulations of the laser rate equations in future work.
Applying the FRM to stochastic simulations of the laser rate equations leads to short computation times, which makes it possible to experiment more easily with changes in system parameters and even alterations of the laser model. Using the FRM, one could potentially study other characteristics of the laser as well, like the emission spectrum and linewidth, and the response to different types of pumping. Additionally, if the method can be suitably modified, it can perhaps be used to further examine the effects of inter-emitter coupling in nanolasers through simulations using an extended set of rate equations [42]. In conclusion, the FRM is efficient, robust and versatile, and we hope that this detailed description of the algorithm and how to apply it to the laser rate equations will lead to more new research and insights in the future. In this appendix, we will give a slightly more detailed description of why the FTI and FRM algorithms have the specific forms described in the main text. We will describe in some detail the probability theoretical considerations that go into the derivation of the algorithms, though the full mathematical description is beyond the scope of this paper; see instead e.g. [19].
As mentioned in Sec. III, our aim is to simulate the interaction of a collection of discrete photons and emitters in a single-mode optical cavity, whose dynamics are governed by the set of rate equations (1)- (2). To perform these stochastic simulations, we make a fundamental assumption common to the stochastic formulation of chemical kinetics [19]: There exist a set of constants c µ for each type of event µ, which depend only on the physical properties of the system (not particle numbers), such that c µ dt = Average probability, to first order in dt, that one particular combination of photons and/or emitters will undergo (A1) an event of type µ during the next time increment dt.
For instance, the probability that a particular photon in the cavity mode will be lost through the cavity mirrors during a time increment dt, averaged over all photons in the cavity mode, is γ c dt + o(dt), where o(dt) are extra terms satisfying o(dt)/dt → 0 as dt → 0; In other words, c c = γ c . Likewise, the probability that a particular photon in the cavity mode will cause a specific excited emitter to emit in a stimulated emission event during dt is γ r dt + o(dt), so c st = γ r . Additionally, if at time t the numbers of photons and excited emitters are n p (t) and n e (t), respectively, we can define the functions h µ as h µ = Number of distinct combinations of photons and/or emitters that can undergo an event (A2) of type µ, given that the current particle numbers are (n p (t), n e (t)).
For instance, the number of photons which could potentially leak out of the cavity is the current photon number, i.e. h c = n p (t); and since each of the n p (t) photons can prompt a stimulated emission event from each of the n e (t) excited emitters, we take h st = n e (t)n p (t). Combining the two quantities defined in eqs. (A1) and (A2), we can define the rates, or propensity functions, a µ as a µ dt = h µ c µ dt = Probability, to first order in dt, that an event of type µ will (A3) occur during the next time increment dt.
In [18], the stochastic simulations are carried out by choosing a fixed time increment dt and using the quantities defined above to determine how many of each type of event µ that happen at every time step: It is suggested that the total number of events of type µ happening during the time increment dt follows a binomial distribution, Binom(h µ , c µ dt), based on the idea that each individual event which happens can be seen as a success in a Bernoulli trial, i.e. a yes/no-experiment, whose probability of success is c µ dt (Of course this assumes that dt is small enough that c µ dt ≤ 1). As an example, the number of photons lost through the cavity mirrors between t and t + dt is taken to be a random integer drawn from the binomial distribution Binom(n p (t), γ c dt), since the average probability of losing one photon is γ c dt and there are n p (t) "tries". In order to speed up the simulations, a further assumption is made in [18], namely that the time increment dt may be chosen sufficiently small that the binomial distributions Binom(h µ , c µ dt) may be replaced by Poisson distributions Poiss(h µ c µ dt) = Poiss(a µ dt). This removes the upper bound on the integers that may be drawn, introducing the risk that e.g. the number drawn for how many photons are lost through the cavity mirrors is greater than the number of photons currently in the cavity mode, but by choosing a sufficiently small time increment dt, this risk is minimal.
While the replacement of binomial distributions by Poisson distributions significantly reduces computation times for the algorithm of [18], the simulation must still run for several hours to obtain satisfactory photon statistics for all pump values. Implementing a more efficient algorithm is therefore highly desirable, and to do this we can change the fundamental viewpoint in the simulation: Instead of determining how many events of each type occur during a fixed time increment, we can try to determine the length of time until some event occurs and then determine which type of event occurs. There are several ways of doing this, but Gillespie's First Reaction Method is based on the following idea: For each µ we can find a tentative time τ µ until an event of type µ would occur, assuming no other event occurs before, and then we choose the smallest of these tentative times, τ µ0 , as the time until an event actually occurs. The corresponding event type, µ 0 , is chosen as the event type that happens. To generate the random τ µ 's, we should consider the probabilities at time t for an event of type µ to happen between times t + τ and t + τ + dt. We may compute such a probability as the product of the probability that no event of type µ happens between t and t + τ , given by P 0 (τ ) = exp(−a µ τ ), and the probability a µ dt that an event of type µ does occur during an interval of length dt [19]. In other words, the probabilities that we are looking for are of the form P µ (τ ) dt = a µ exp(−a µ τ ) dt.
To generate the random tentative times τ µ , we can there-fore draw random numbers from the distributions P µ (τ ), which are essentially exponential distributions with parameters a µ .