Open Access
Add to BibSonomy
Abstract
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.
© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. 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–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–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–15]. Recently, stochastic methods have been used to simulate nanoscale lasers with large values of the spontaneous emission $\beta$-factor [16–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–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–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 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 inter-emitter correlations.
The paper is structured as follows: In Section 2 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 3 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 4 we will examine the intra-cavity photon distributions of the laser and in Section 5 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 6.
2. 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 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 $\gamma _{c}$, and the emitters have a rate of loss into non-lasing modes given by $\gamma _{bg}$. The rate per emitter of emission into the lasing mode is denoted $\gamma _{r}$, giving a total decay rate $\gamma _{t} = \gamma _{r} + \gamma _{bg}$ from the excited state. The dependence of $\gamma _{r}$ on material and device parameters is discussed in [18], and includes Purcell enhancement. Finally, the emitters are pumped at a rate per emitter of $\gamma _{p}$.
In this work we will focus on so-called Class B lasers, where the dephasing rate $\gamma _{2}$ of the emitters is much larger than the photon loss rate $\gamma _{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]:
Table 1. Average rates of event types in the laser rate equations.
The rate Eqs. (1)–(2) can be solved analytically in steady state, yielding expressions for the mean photon number $\langle n_{p}\rangle$ and number $\langle n_{e}\rangle$ of excited emitters as functions of the pump rate, $\gamma _{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 $\beta \rightarrow 1$, these features become less distinct, and instead the photon statistics are used to characterise the transition to lasing [11–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 $\langle n_{p}\rangle$ and the variance $\langle\Delta n_{p}^{2}\rangle = \langle n_{p}^{2}\rangle - \langle n_{p}\rangle^{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) [26] as
3. 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} \geq 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.
3.1 Fixed time-increment stochastic algorithm
The basic idea behind the stochastic simulations is to assume that the laser rate Eqs. (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 1.
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 $\mu$ 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_{\mu }$ 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 $\tau$-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_{\mu }$ for all event types $\mu$ according to Table 1 using current particle numbers.
- 3. Determine the number of each type of event $\mu$ occuring during $dt$ by random draws from Poisson distributions with parameters $a_{\mu }dt$.
- 4. Update $n_{p}$ and $n_{e}$ accordingly, set $t \rightarrow 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 intra-cavity 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 $\langle n_{p}\rangle$, excited-emitter number $\langle n_{e}\rangle$, and the photon variance $\langle \Delta n_{p}^{2} \rangle$. In practice, the size of $dt$ can be chosen at a given pump rate as some fraction $dt_\textrm {{frac}}$ of the smallest reciprocal rate appearing in the rate Eqs. (1)–(2), computed using the analytical steady-state 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_\textrm {{frac}} / \gamma _{r}\langle n_{e}\rangle$. This ensures that most iterations in the low-pump limit involve at most one event of any type, and convergence is typically obtained for $dt_\textrm {{frac}}$ on the order of $10^{-2}$ or less. For all FTI simulations in this work we have used $dt_\textrm {{frac}} = 10^{-2}$. 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.
3.2 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 Eqs. (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 $\tau _{\mu }$ until this event occurs. Once all the tentative times $\tau _{\mu }$ have been determined, we can choose the shortest, $\tau _{\mu _{0}}$, as the actual time until the first event (or reaction) occurs, and the corresponding type of event, $\mu _{0}$, as the event that occurs. The tentative time $\tau _{\mu }$ can be efficiently generated by drawing a random number from an exponential distribution with parameter $a_{\mu }$, 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)
- Initialize: Set system parameters, set initial number of photons $n_{p}(t=0)$ and excited emitters $n_{e}(t=0)$.
- Calculate rates $a_{\mu }$ for all event types $\mu$ according to Table 1 using current particle numbers.
- For each $\mu$, generate a tentative time increment $\tau _{\mu }$ by a random draw from an exponential distribution with parameter $a_{\mu }$.
- Determine the event type $\mu _{0}$ for which $\tau _{\mu _0}~=~\min _{\mu }\{\tau _{\mu }\}$.
- Update $n_{p}$ and $n_{e}$ according to event type $\mu _0$, set $t \rightarrow t+\tau _{\mu _0}$.
- 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. Therefore, there is no need for additional checks of convergence in terms of such parameters. In this sense, the FRM is an exact method to stochastically simulate the laser dynamics. In addition, exponentially distributed random numbers $\tau _{\mu }$ with parameters $a_{\mu }$ can be quickly and inexpensively generated using uniformly distributed random numbers $r_{\mu }$ on the unit interval $(0,1)$ as $\tau _{\mu } = -\log (r_{\mu })/a_{\mu }$ [29,30]. 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. We note that the very small mean photon numbers obtained for the blue curves in the low-pump range of Fig. 2 lead to relatively large statistical uncertainties. Therefore, results below a certain pump rate have been discarded, like in the corresponding figure of [18].
Fig. 2. (a)–(c) Steady-state analytical [full lines] and numerical [markers] results for the mean photon number, relative intensity noise and photon auto-correlation versus pump rate. For all plots, the cavity decay rate is $\gamma _{c}~=~ 10^{11}\;s^{-1}$ and the total decay rate is $\gamma _{t}=10^{10}\;s^{-1}$. The $\beta$-factor and number of emitters, $(\beta , n_{0})$, are varied between $(10^{-3},2\cdot 10^{4})$ [blue], $(10^{-1},2\cdot 10^{2})$ [red] and $(1,2\cdot 10^{1})$ [black]. For $\beta =10^{-3}$ [blue], low-pump results have been discarded because of statistical uncertainty. (d) Computation time using the FTI [dashed lines] and FRM [full lines] algorithms for different lengths of simulated time at the pump rates indicated in (a). Computation times have been normalized to the lowest computation time, which is 0.01 $s$, obtained using MatLab on a PC with a 2.7 GHz quad-core processor and 16 GB DDR4 RAM. The parameters are the same as those used for the $\beta =0.1$ plots in (a)–(c).
Fig. 3. Mean photon number (a), photon auto-correlation (b) and Fano-factor (c) versus pump rate for a laser with $\beta =1$ [black curves] and $\beta =10^{-3}$ [blue curves]. (d)–(f) Normalised photon number distribution for the laser with $\beta =1$ at the pump values indicated with black crosses in (a)–(c). (g)–(i) Normalised photon number distribution for the laser with $\beta =10^{-3}$ at the pump values indicated with blue circles in (a)–(c). The parameters used to make these plots are the same as those used to produce the curves in Fig. 2. In (d,g) the orange curves indicate Bose-Einstein distributions with mean values determined by the mean photon numbers $\langle n_{p} \rangle$. In (f,i) the orange dashed curves indicate Poisson distributions with parameters given by the mean photon number $\langle n_{p} \rangle$, and the red curves indicate Gaussian distributions with means and variances given by the mean photon number $\langle n_{p} \rangle$ and photon number variance $\langle \Delta n_{p}^{2} \rangle$, respectively.
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.
4. 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 $\beta =1$ and $\beta =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 $\beta =1$ and $\beta = 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 at high pump rates, the photon PDF approaches a Gaussian distribution which is wider than a Poisson distribution with the same mean value. Mathematically, we can explain this by the photon variance approaching a value slightly different from the characteristic Poissonian $\langle \Delta n_{p}^{2} \rangle = \langle n_{p} \rangle$. Equivalently, the Fano factor $F = \langle\Delta n_{p}^{2} \rangle / \langle n_{p} \rangle$ 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
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.
5. 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 [26]. 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, $\gamma _{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,34–39]. 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 $\gamma _{r}$ in Eqs. (1)–(2) by two separate rates, $\gamma _{r}^\textrm {sp}$ and $\gamma _{r}^\textrm {st}$, 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 Eqs. (1)–(2) can be rewritten using this quantity. We can solve the rate equations including $\xi$ in the same way as we did for the traditional rate Eqs. (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 $\xi$-dependent analytical expressions for $\langle \Delta n_{p}^{2} \rangle$, 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_{\mu }$ to include $\xi$. 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 $\xi$, 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 auto-correlation 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 $\xi$ 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 $\xi$ from 1 leads to changes in one of the derived parameters $\beta = \gamma _{r}^\textrm {sp}/(\gamma _{r}^\textrm {sp} + \gamma _{bg})$ or $n_{th} = \gamma _{c}/\gamma _{r}^\textrm {st}$. 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, $\xi = \gamma _{r}^\textrm {sp}/\gamma _{r}^\textrm {st} = 1$.Fig. 4. Photon number (a,b) and photon auto-correlation function (c,d) as functions of the pump rate for different values of the ratio $\xi = \gamma _{r}^\textrm {sp}/\gamma _{r}^\textrm {st}$. All plots are produced using $n_{0}=20000$. In (a,c) $\gamma _{r}^\textrm {st}$ and hence $n_{th}= \gamma _{c} / \gamma _{r}^\textrm {st}=10000$ are kept fixed, while $\gamma _{r}^\textrm {sp}$ changes with $\xi$. In (b,d) $\gamma _{r}^\textrm {sp}$ and hence $\beta = 10^{-3}$ are kept fixed, while $\gamma _{r}^\textrm {st}$ and hence $n_{th}$ change with $\xi$.
From Fig. 4 we see that decreasing the rate of spontaneous 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 sub-thermal 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 increased rate of stimulated emission is also the reason for a larger mean photon number in the high-pump limit seen in Fig. 4(b). 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,38,39]. Figure 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 [38,39], meaning excitations in the system preferably reside in the emitters rather than the cavity mode. This effect, which is related to subradiance [40,41], 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 $\xi \langle1$, we have plotted in Fig. 5 the time-resolved photon number for two values of $\xi$. As shown in Figs. 5(a) and (b), we examine the case where the steady-state mean photon number is the same, namely $\langle n_p\rangle = 0.2$, while the photon distributions and statistics are either thermal or super-thermal. From (c,d) it is clear that the case where $\xi \langle1$, 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.
Fig. 5. (a) Photon auto-correlation function versus mean photon number for $\xi =1$ [full blue line] and $\xi =1/4$ [dashed red line]. Circles indicate simulated results with steady-state mean photon number $\langle n_p\rangle = 0.2$. (b) Simulated photon distribution for the two cases indicated by the circles in (a). (c,d) Time-resolved photon number for the two cases indicated in (a) during 10 $ns$ of simulated time after reaching steady state.
6. 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 intra-cavity 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 $\tau$-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. Indeed, optimizing the size of the time-increment may lead to computation times for the FTI algorithm which are closer to those for the FRM algorithm in the high-pump limit, at a low cost to 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,42], 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 $\tau _{\mu }$ and only updating certain rates $a_{\mu }$, 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 [43]. 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.
A. Derivation of the FTI and FRM algorithms
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. 3, 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 Eqs. (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_{\mu }$ for each type of event $\mu$, which depend only on the physical properties of the system (not particle numbers), such that
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 that occurs. There are several ways of doing this, but Gillespie’s first reaction method is based on the following idea: For each $\mu$ we can find a tentative time $\tau _{\mu }$ until an event of type $\mu$ would occur, assuming no other event occurs before, and then we choose the smallest of these tentative times, $\tau _{\mu _0}$, as the time until an event actually occurs. The corresponding event type, $\mu _{0}$, is chosen as the event type that happens. To generate the random $\tau _{\mu }$’s, we should consider the probabilities at time $t$ for an event of type $\mu$ to happen between times $t+\tau$ and $t+\tau +dt$. We may compute such a probability as the product of the probability that no event of type $\mu$ happens between $t$ and $t+\tau$, given by
and the probability $a_{\mu } dt$ that an event of type $\mu$ does occur during an interval of length $dt$ [19]. In other words, the probabilities that we are looking for are of the form To generate the random tentative times $\tau _{\mu }$, we can therefore draw random numbers from the distributions $P_{\mu }(\tau )$, which are essentially exponential distributions with parameters $a_{\mu }$.Funding
Danmarks Grundforskningsfond (DNRF147); Villum Fonden (8692).
Acknowledgments
The authors wish to thank I.E. Protsenko and A.V. Uskov for stimulating discussions regarding collective effects in nanolasers. J.M. wishes to thank G.L. Lippi for helpful discussions regarding stochastic simulations.
Disclosures
The authors declare no conflicts of interest.
References
1. E. M. Purcell, “Spontaneous emission probabilities at radio frequencies,” Proc. Am. Phys. Soc. 69, 681 (1946). [CrossRef]
2. H. Yokoyama and S. D. Brorson, “Rate equation analysis of microcavity lasers,” J. Appl. Phys. 66(10), 4801–4805 (1989). [CrossRef]
3. H. Altug, D. Englund, and J. Vučković, “Ultrafast photonic crystal nanocavity laser,” Nat. Phys. 2(7), 484–488 (2006). [CrossRef]
4. T. Suhr, N. Gregersen, K. Yvind, and J. Mørk, “Modulation response of nanoLEDs and nanolasers exploiting Purcell enhanced spontaneous emission,” Opt. Express 18(11), 11230–11241 (2010). [CrossRef]
5. C.-Y. A. Ni and S. L. Chuang, “Theory of high-speed nanolasers and nanoLEDs,” Opt. Express 20(15), 16450–16470 (2012). [CrossRef]
6. D. A. Miller, “Attojoule optoelectronics for low-energy information processing and communications,” J. Lightwave Technol. 35(3), 346–396 (2017). [CrossRef]
7. M. C. Gather and S. H. Yun, “Single-cell biological lasers,” Nat. Photonics 5(7), 406–410 (2011). [CrossRef]
8. G. Shambat, S.-R. Kothapalli, J. Provine, T. Sarmiento, J. Harris, S. S. Gambhir, and J. Vuckovic, “Single-cell photonic nanocavity probes,” Nano Lett. 13(11), 4999–5005 (2013). [CrossRef]
9. A. H. Fikouras, M. Schubert, M. Karl, J. D. Kumar, S. J. Powis, A. Di Falco, and M. C. Gather, “Non-obstructive intracellular nanolasers,” Nat. Commun. 9(1), 4817 (2018). [CrossRef]
10. R.-M. Ma and R. F. Oulton, “Applications of nanolasers,” Nat. Nanotechnol. 14(1), 12–22 (2019). [CrossRef]
11. C. Gies, J. Wiersig, M. Lorke, and F. Jahnke, “Semiconductor model for quantum-dot-based microcavity lasers,” Phys. Rev. A 75(1), 013803 (2007). [CrossRef]
12. W. W. Chow, F. Jahnke, and C. Gies, “Emission properties of nanolasers during the transition to lasing,” Light: Sci. Appl. 3(8), e201 (2014). [CrossRef]
13. S. Kreinberg, W. W. Chow, J. Wolters, C. Schneider, C. Gies, F. Jahnke, S. Höfling, M. Kamp, and S. Reitzenstein, “Emission from quantum-dot high-β microcavities: transition from spontaneous emission to lasing and the effects of superradiant emitter coupling,” Light: Sci. Appl. 6(8), e17030 (2017). [CrossRef]
14. A. Moelbjerg, P. Kaer, M. Lorke, B. Tromborg, and J. Mørk, “Dynamical properties of nanolasers based on few discrete emitters,” IEEE J. Quantum Electron. 49(11), 945–954 (2013). [CrossRef]
15. M. Marconi, J. Javaloyes, P. Hamel, F. Raineri, A. Levenson, and A. M. Yacomotti, “Far-from-equilibrium route to superthermal light in bimodal nanolasers,” Phys. Rev. X 8(1), 011013 (2018). [CrossRef]
16. A. Lebreton, I. Abram, N. Takemura, M. Kuwata-Gonokami, I. Robert-Philip, and A. Beveratos, “Stochastically sustained population oscillations in high-β nanolasers,” New J. Phys. 15(3), 033039 (2013). [CrossRef]
17. G. P. Puccioni and G. L. Lippi, “Stochastic simulator for modeling the transition to lasing,” Opt. Express 23(3), 2369–2374 (2015). [CrossRef]
18. J. Mørk and G. L. Lippi, “Rate equation description of quantum noise in nanolasers with few emitters,” Appl. Phys. Lett. 112(14), 141103 (2018). [CrossRef]
19. D. T. Gillespie, “A general method for numerically simulating the stochastic time evolution of coupled chemical reactions,” J. Comput. Phys. 22(4), 403–434 (1976). [CrossRef]
20. D. T. Gillespie, “Exact stochastic simulation of coupled chemical reactions,” J. Phys. Chem. 81(25), 2340–2361 (1977). [CrossRef]
21. M. A. Gibson and J. Bruck, “Efficient exact stochastic simulation of chemical systems with many species and many channels,” J. Phys. Chem. A 104(9), 1876–1889 (2000). [CrossRef]
22. J. Pahle, “Biochemical simulations: stochastic, approximate stochastic and hybrid approaches,” Briefings Bioinf. 10(1), 53–64 (2009). [CrossRef]
23. F. T. Arecchi, G. L. Lippi, J. R. Puccioni, and G. P. Tredicce, “Deterministic chaos in laser with injected signal,” Opt. Commun. 51(5), 308–314 (1984). [CrossRef]
24. G. Bjork and Y. Yamamoto, “Analysis of semiconductor microcavity lasers using rate equations,” IEEE J. Quantum Electron. 27(11), 2386–2396 (1991). [CrossRef]
25. P. R. Rice and H. J. Carmichael, “Photon statistics of a cavity-QED laser: A comment on the laser–phase-transition analogy,” Phys. Rev. A 50(5), 4318–4329 (1994). [CrossRef]
26. L. A. Coldren, S. W. Corzine, and M. L. Mashanovitch, Diode lasers and photonic integrated circuits (John Wiley & Sons, 2012).
27. R. Loudon, The Quantum Theory of Light (Oxford University Press, 2000).
28. D. T. Gillespie, “Approximate accelerated stochastic simulation of chemically reacting systems,” J. Chem. Phys. 115(4), 1716–1733 (2001). [CrossRef]
29. C. P. Robert and G. Casella, Monte Carlo Statistical Methods (Springer New York, 2004).
30. V. Krishnan, Probability and Random Processes (John Wiley & Sons, 2006).
31. A. V. Kozlovskii and A. N. Oraevsky, “Non-Poisson photon statistics in an n-level (n= 2–4) laser above the threshold,” Quantum Electron. 24(3), 255–260 (1994). [CrossRef]
32. A. V. Kozlovskii, “Super-poissonian fluctuations of number of photons in the radiation of a single-mode laser operating near active particles number threshold,” Opt. Spectrosc. 116(1), 115–121 (2014). [CrossRef]
33. N. J. van Druten, Y. Lien, C. Serrat, S. S. R. Oemrawsingh, M. P. van Exter, and J. P. Woerdman, “Laser with thresholdless intensity fluctuations,” Phys. Rev. A 62(5), 053808 (2000). [CrossRef]
34. E. C. André, I. E. Protsenko, A. V. Uskov, J. Mørk, and M. Wubs, “On collective Rabi splitting in nanolasers and nano-LEDs,” Opt. Lett. 44(6), 1415–1418 (2019). [CrossRef]
35. J. G. Bohnet, Z. Chen, J. M. Weiner, D. Meiser, M. J. Holland, and J. K. Thompson, “A steady-state superradiant laser with less than one intracavity photon,” Nature 484(7392), 78–81 (2012). [CrossRef]
36. D. Meiser, J. Ye, D. R. Carlson, and M. J. Holland, “Prospects for a millihertz-linewidth laser,” Phys. Rev. Lett. 102(16), 163601 (2009). [CrossRef]
37. D. Meiser and M. J. Holland, “Steady-state superradiance with alkaline-earth-metal atoms,” Phys. Rev. A 81(3), 033847 (2010). [CrossRef]
38. H. A. M. Leymann, A. Foerster, F. Jahnke, J. Wiersig, and C. Gies, “Sub-and superradiance in nanolasers,” Phys. Rev. Appl. 4(4), 044018 (2015). [CrossRef]
39. F. Jahnke, C. Gies, M. Aßmann, M. Bayer, H. A. M. Leymann, A. Foerster, J. Wiersig, C. Schneider, M. Kamp, and S. Höfling, “Giant photon bunching, superradiant pulse emission and excitation trapping in quantum-dot nanolasers,” Nat. Commun. 7(1), 11540 (2016). [CrossRef]
40. R. H. Dicke, “Coherence in spontaneous radiation processes,” Phys. Rev. 93(1), 99–110 (1954). [CrossRef]
41. M. Gross and S. Haroche, “Superradiance: An essay on the theory of collective spontaneous emission,” Phys. Rep. 93(5), 301–396 (1982). [CrossRef]
42. L. Chusseau, J. Arnaud, and F. Philippe, “Rate-equation approach to laser light statistics,” Opt. Spectrosc. 94(5), 746–754 (2003). [CrossRef]
43. I. Protsenko, E. C. André, A. Uskov, J. Mørk, and M. Wubs, “Collective effects in nanolasers explained by generalized rate equations,” – arXiv:1709.08200 (2017).
