1. 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, meaningful only for large enough averages [9]) rather than those intrinsic to the physics of the process near threshold.

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 [10, 11] –, 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.

2. 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.

The chain of physical processes which are mimicked by the SS is the following (cf. Fig. 1):

 

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. Absorption is neglected as we consider an ideal system (e.g., perfect four-level laser). Its inclusion is expected to mainly raise the laser threshold.

Download Full Size | PDF

As detailed in Table 1, the probability distribution for the pump is taken to be a Poissonian (ℙ) [12], while all other processes are simulated through a binomial distribution (𝔹) [12], 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”.

Table 1. 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 (D_L); thus no explicit process is associated to this channel (we assume, for simplicity, all relaxations to be radiative). Each event (column 3) is computed from the probability distribution (column 4) whose first argument is defined in column 5.

View Table

The following recurrence relations define the simulator’s rules:

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 1). 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 1) 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).

3. Results

 

Fig. 2 Left panel: Photon number ℘ vs. injection current i computed from the steady-states of the rate equations (eq. (25) in [2]) with the following parameter values: γ = 1 × 10⁻¹⁰s⁻¹, ξ = 0.1, τ_sp = 1×10⁻⁹s, τ_nr = 1×10⁻¹⁰s. Right panel: Average photon number (<S>) as a function of pump (P) obtained from the SS. <S> is computed as the average of the temporal average photon number $\overline{S}=\frac{1}{N}{\sum }_{\ell =1}^N{S}_{\ell }$ over ten series of events (issued from different values input into the random number generator). The corresponding error bars represent the standard deviation of the individual averages S̄ and give a measure of the variability of the photon number, due to the stochastic nature of the conversion process. Parameter values: γ_‖ = 2.5 × 10⁹s⁻¹, Γ_c = 1 × 10¹¹s⁻¹, Γ_o = 5 × 10¹³s⁻¹. The time step used for evaluating the probabilities (last column of Table 1) is $\mathrm{\Delta }\tau =\frac{1}{10\times {\mathrm{\Gamma }}_o}=2\mathit{fs}$, thus τ_r = r Δτ (r already defined). Notice that the curve with β = 10⁻⁷ is rescaled by a factor 10⁻² both in the horizontal and vertical axes for graphical purposes.

Download Full Size | PDF

Figure 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⁻⁷). 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 [13]. Thus, the SS can be used to model class B lasers [14] of any size. An eplanation for the bump appearing in the laser response for β = 0.1 will require additional work.

Figure 3 offers other predictive aspects of the SS. The left panel shows the average photon number, as a function of 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 ≈ .1 for the parameters of this simulation. In addition, we remark that far below threshold $\frac{<R_L>}{<R_o>}\approx \beta$, 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 $\frac{<R_L>}{<S>}\ll 1$. It is important to remark that while the stimulated component (blue curve) steadily grows, both spontaneous components remain clamped above threshold. The residual growth – which then decays – of both spontaneous emission components in the “S”-portion of the stimulated emission remains for the moment unexplained.

 

Fig. 3 Static and dynamical predictions for β = 0.01. Left panel: Average photon number for stimulated photons (blue curve), on-axis (brown) and off-axis (magenta) spontaneous photons as a function of pump. Center panel: dynamical evolution of <S> for three different values of pump, marked by the corresponding colors on the steady-state response (right panel).

Download Full Size | PDF

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 [15]) on a desktop PC AMD Phenom X6 1090T, and defining an efficiency $\eta =\frac{\text{physical}\hspace{0.17em}\text{simulation}\hspace{0.17em}\text{time}}{\text{computer}\hspace{0.17em}\text{running}\hspace{0.17em}\text{time}}$, we find $0.1\frac{\mathit{ns}}{s}\lesssim \eta \lesssim 0.25\frac{\mathit{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 fs accuracy, the equivalent of the trace displayed on an oscilloscope (typically 50ns long for good resolution) in as little as 20s of computation!

4. 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.