Stochastic Simulator for modeling the transition to lasing

A Stochastic Simulator (SS) is proposed, based on a semiclassical description of the radiation-matter interaction, to obtain an efficient description of the lasing transition for devices ranging from the nanolaser to the traditional"macroscopic"laser. Steady-state predictions obtained with the SS agree both with more traditional laser modeling and with the description of phase transitions in small-sized systems, and provide additional information on fluctuations. Dynamical information can easily be obtained, with good computing time efficiency, which convincingly highlights the role of fluctuations at threshold.


I. INTRODUCTION
The progressive miniaturization of lasers, begun with the introduction of the VCSEL [1], has opened a wealth of questions related to the operation of very small devices. Problems with the definition of threshold were recognized very early on [2] and still persist [3], and difficulties with a quantitative interpretation of correlation functions have been more recently pointed out [4]. Failure of traditional continuous modeling based on rate equations has been widely recognized and is discussed in detail in [4] and various alternatives exist to predict the behaviour of small lasers. In essence, the difficulty arises from the strongly discrete and stochastic nature of the problem and the limited number of spontaneous emission channels in the cavity.
We can distinguish two classes of stochastic models: 1. those which describe the evolution of the ensemble, either predicting the static features or their statistical dynamics, and 2. those which predict the detailed dynamics. To the first class belong Monte-Carlo simulations, which have been recently applied to lasers with success [5], and Fokker-Planck models which have a long history in laser physics [6]. The Master Equation clearly belongs to the second class, since it allows for detailed predictions of the dynamical evolution of the system, but its use becomes unwieldy when the number of elements (modes of the e.m. field) exceeds a few units. A practical reduction to a random walk on a lattice has been proposed [7]. Recently developed quantum models [8] are to be ascribed to this second group.
The aim of our Stochastic Simulator (SS) is to provide a tool for rapid computation predicting the laser dynamics and its emergence from the spontaneous emission, without the mathematical assumptions normally associated to the derivation of a differential description. Rather, we set up an ensemble of very simple (semiclassical) rules reproducing the physics of the interaction between radiation and matter, and let the SS predict the sequence of events which will occur, on the basis of the randomness of each occurrence. The advantages of this choice are that: a. it provides extreme computational efficiency, thus offering the possibility of quick dynamical predictions and of performing averages over wide samples; and b. letting the operator observe emerging dynamical aspects which otherwise would be either washed out (in stationary approaches) or reproduce the statistical characteristics of the modelling choice (rate equations with Langevin noise) rather than those intrinsic to the physics of the process.
We remark that the present proposal for a SS aims at modeling only the basic ingredients necessary for a description of laser action. This choice is dictated by the interest in identifying those elements which are essentially intrinsic to the stochastic nonlinear process. Refinements are possible, and planned, to include specific effects which will allow for a closer description of particular kinds of lasers.
A discretization of the rate equations, fulfilling a similar role, has previously been proposed [4], but its scope was strongly reduced by the neglection of the spontaneous emission -inherent to the rate equations [9,10] -, which precludes the determination of thresholds or of the mixtures of coherent and incoherent photons. The absence of the latter entirely removes the interplay spontaneous/stimulated emission, thereby preventing one from correctly observing the rise of the stimulated component of the laser field intensity.

II. STOCHASTIC SIMULATOR
The Stochastic Simulator mimics a sequence of events as a succession of possible occurrences taking place at discrete times, with a time discretization (∆τ , cf. caption of Fig. 2) small compared to the fastest time scale. Indeed, since each process is represented by a probabilistic distribution determining whether and when the corresponding process takes place, the discretization must be sufficiently fine to allow for a correct description of even the fastest process (the escape of off-axis spontaneous photons -cf. below). This ensures that all processes will be described with sufficient accuracy, while most of the processes will have zero outcome at any given step, at least below and around threshold.
FIG. 1: Scheme of principle of the intracavity processes included in the SS: each dipole is pumped by an external source and may decay by emitting a photon in one of the three following channels -a. the off-axis spontaneous "mode", grouping all modes other than the lasing one (green photons); b. the on-axis spontaneous mode (blue photons); c. the (on-axis) stimulated mode (red photons). All on-axis photons are recycled by the cavity and are transmitted by the coupling mirror; the off-axis (green) ones exit the interaction volume laterally. All processes (including the mirror transmission) are described stochastically. In the schematic representation, we summarize the ensemble of processes leading to the population of the upper state -with possible population inversion -as a single upwards magenta arrow. The large bracket denotes the interaction processes occurring between each individual dipole and the radiation.
The chain of physical processes which are mimicked by the SS is the following (cf. Fig. 1): a. excitation of a dipole by a pumping mechanism; b. disexcitation of a dipole through one of the following channels: i. spontaneous relaxation into an off-axis cavity mode; ii. spontaneous relaxation into the on-axis cavity mode; iii. stimulated emission (into the on-axis mode).
c. escape of an off-axis spontaneous photon; d. transmission of an on-axis spontaneous photon through the cavity mirror; e. transmission of a stimulated photon through the cavity mirror.
As detailed in Table I, the probability distribution for the pump is taken to be a Poissonian (P), while all other processes are simulated through a binomial distribution (B), since they correspond to yes/no events [4]. All variables -dipole number (N ), stimulated photon number (S), on-axis spontaneous photon number (R L ), and off-axis spontaneous photon number (R o ) -need to be defined at all times, independently of their actual value (e.g., S = 0 far below threshold), since the SS can only accumulate population numbers in preexisting "categories". I: Synoptic of the processes handled by the SS. The off-axis spontaneous emission results from the difference between the population relaxation events (N d ) and the on-axis spontaneous emissions (DL); thus no explicit process is associated to this channel (we assume, for simplicity, all relaxations to be radiative). The last two columns explicitely define the procedures for computing each individual event probability.

Process
Physical event Event Math expression P

Excitation Pump Absorption
Np On-axis spont. em. Ssp The following recurrence relations define the simulator's rules: where N P represents the pumping process, N d the spontaneous relaxation processes which reduce the population inversion N , E S the stimulated emission processes which also reduce N , L S the leakage of stimulated photons through the output coupler, S sp accounts for the seed starting the first stimulated emission process, D L represents the fraction of spontaneous relaxation processes which enter the on-axis mode (and therefore superpose to the stimulated emission), L L the losses for the on-axis fraction of the spontaneous photons through the output coupler, and L o the losses for the off-axis fraction of the spontaneous photons exiting (laterally) the cavity volume. The definitions of the quantities appearing on the r.h.s. of eqs. (1-4) are given in Table I, where the parameters appearing in the last column are the usual ones appearing in laser models [4]. Specifically, P represents the pump (average of the Poissonian), γ is the population relaxation rate, Γ c the cavity losses (on-axis mode), Γ o the cavity losses for the off-axis "mode" (i.e., the average lifetime for spontaneous photons emitted in modes other than the on-axis one), β represents the fraction of spontaneous emission coupled into the on-axis mode [2]. At variance with the simple discretization introduced in [4], the SS handles the different time scales over which the dynamical evolution of the processes take place by introducing independent, explicit probabilities (cf. last column in Table I). This is achieved by writing individual probability distribution for each process where the time τ r (r = i, j, k, m) represents the time accumulated from the previous event. In other words, time is reset to zero, independently for each kind of process r, every time an event has taken place.
Finally, since the SS accounts for the initiation of stimulated emission randomly starting from a spontaneous event, we introduce the seed S sp (cf. Table I) which activates (producing S sp = 1) only in the absence of stimulated photons (S = 0), in the presence of spontaneous photons in the on-axis mode (R L = 0), and with probability determined by a random number α and a set level K. The value of the constant K tunes the probability of transforming a spontaneous photon into a stimulated one, thus initiating the amplification process.
With the values of the laser parameters defined above, the SS is started with the desired value of the pump P and initial values for the variables (N , S, R L , and R o ). The choice of the initial values is not crucial as long as they do not differ by more than a couple of orders of magnitude from the actual averages belonging to the set of parameters chosen, since the SS is run with a transient τ t to let the variables relax. The first determination of reasonable starting values has to be empirically obtained with sufficiently long runs (monitoring the averages).  Table I) is ∆τ = 1 10×Γo = 2f s, thus τr = r ∆τ (r already defined). Notice that the curve with β = 10 −7 is rescaled by a factor 10 −2 both in the horizontal and vertical axes for graphical purposes. Fig. 2 shows the predicted laser output as a function of the pump (right panel) which gives a good qualitative agreement with the equivalent curves computed from the rate equations (cf. caption for details). The SS is capable of not only producing a meaningful laser response, but also predicting the deviations (error bars) for nanolasers (0.01 β ≤ 1), for microlasers (0.0001 β 0.01) and even for macroscopic lasers (e.g., β = 10 −7 ). The progressive sharpening of the transition confirms the well-known trend with system size and agrees with the prediction that in the thermodynamic limit the laser threshold becomes a true phase transition [11]. Thus, the SS can be used to model class B lasers [12] of any size.  pump, for the three different kinds of photons: spontaneous off-axis, spontaneous on-axis and stimulated. According to "conventional wisdom", which defines threshold as the pump value for which the average number of stimulated photons equals the average sum of all the spontaneous photons, one can read off the graph the threshold value as P ≈ 0.1 for the parameters of this simulation. In addition, we remark that far below threshold <R L > <Ro> ≈ β, while this ratio increases by a factor up to 5 with growing < S >. This shows that the cavity amplification acts also on the (on-axis) spontaneous photons, although <R L > <S> 1. It is important to remark that while the stimulated component (blue curve) steadily grows, both spontaneous components remain clamped above threshold. The residual growthwhich then decays -of both spontaneous emission components in the "S"-portion of the stimulated emission remains for the moment unexplained.

III. RESULTS
The center panel shows the time-dependent stimulated emission output by the laser for different values of pump (averaged over a ∆t = 0.1ns time window, to simulate the finite bandwidth of a typical detector). This convincingly shows that in the threshold region (cf. corresponding color-coded points in the right panel) the fluctuations are as large as the average (black and red curves), thus confirming the observations of [4].
As a final point, the running time of the SS is remarkably short. Running the SS (programmed in C using GSL routines [13]) on a desktop PC AMD Phenom X6 1090T, and defining an efficiency η = physical simulation time computer running time , we find 0.1 ns s η 0.25 ns s depending on the value of β and on the average values of the variables (the better efficiencies being attained for large β-values and low variable averages). Thus, we can predict, with f s accuracy, the equivalent of the trace displayed on an oscilloscope (typically 50ns long for good resolution) in as little as 20s of computation!

IV. CONCLUSIONS
In conclusion, the SS offers the possibility of obtaining very fast predictions on the dynamics of stimulated and spontaneous photons (both on-and off-axis) in lasers whose size ranges from the nanoscale, all the way to the conventional macroscopic devices. Thanks to the computational speed, average predictions are easily obtained and agree with theoretical considerations about phase transitions in small-size systems. As such, we can consider the SS as a valuable tool for providing complementary information to the one usually obtained from traditional modeling of lasers.
One of the authors (GLL) acknowledges helpful discussions with A. Beveratos.