Monte Carlo modeling of the dual-mode regime in quantum-well and quantum-dot semiconductor lasers

: Monte Carlo markovian models of a dual-mode semiconductor laser with quantum well (QW) or quantum dot (QD) active regions are proposed. Accounting for carriers and photons as particles that may exchange energy in the course of time allows an ab initio description of laser dynamics such as the mode competition and intrinsic laser noise. We used these models to evaluate the stability of the dual-mode regime when laser characteristics are varied: mode gains and losses, non-radiative recombination rates, intraband relaxation time, capture time in QD, transfer of excitation between QD via the wetting layer... As a major result, a possible steady-state dual-mode regime is predicted for specially designed QD semiconductor lasers thereby acting as a CW microwave or terahertz-beating source whereas it does not occur for QW lasers.


Introduction
CW dual-mode semiconductor lasers are key components for low cost microwave or terahertz beating generation [1]. Regarding the optical beat source, several ways were explored. Besides the use of independent lasers that is set out of the focus of this paper, difference frequency generation was demonstrated using devices with independent gain regions [2,3], coupled-cavity microlasers [4], or edge-emitting lasers [5]. Such devices all require a precise control of injection currents in the different gain regions to achieve a balanced dual-mode operation. This is specially true if the two modes inject each other in the same optical path, leading to instabilities and decreasing seriously the tuning range. Devices with gain medium shared by the two modes [6,7] are however much more attractive for compactness and simplified pumping scheme. It is well known that even if the laser cavity allows two simultaneous modes with similar quality factors and gains, the gain medium non-linearity sometimes forbids their simultaneous emission as a result of a bistable behavior [8]. A stable dual-mode operation is thus only allowed when the coupling factor between the modes is not strong enough.
An illustrative stability analysis was proposed in [9] to account for the degeneracy of the two polarization modes of VCSELs, a question analogous to the dual-mode operation invoking two longitudinal modes. A fairly comprehensive system behavior was reported, including the polarization switching of a VCSEL based on the gain saturation in a manner consistent with what is observed experimentally.
We have recently addressed analytically the question of dual-mode stability of semiconductor lasers [10]. This has been proposed within the framework of rate-equation models. Bulk and quantum-well lasers with same gain and losses on the two modes were shown bistable. The physical origin is the short intraband relaxation time that strongly couples the carrier populations. To overcome this limitation, we proposed a quantum-dot active region to solve the mode competition because the difference between homogenous broadening due to temperature and inhomogeneous broadening due to growth process dispersion helps decoupling the modes [11][12][13]. This concept has been under experimental investigation in a dedicated device combining a photonic crystal membrane and a vertical Fabry-Perot cavity [14].
In this paper, we largely supplement our previous analytical work [10] by a numerical one using specially designed Monte Carlo markovian models. Considering QD lasers for example, electrons and photons are treated as particles exchanged in the course of time between two families of QDs and a common barrier. Such a model is completely different from the rate-equations model used in [10]. It does not make any assumption on mode self-or crosssaturations and thereby on their coupling strength. As a result, answer on stability will arise only from the carrier-carrier modeling through the chosen energy levels and from the carrierphoton interactions. Mode competition is thus taken into account ab initio, and simultaneously the model yields steady-state solutions and noise of both modes. To our knowledge, very few works have addressed the Monte Carlo simulation of laser with similar microscopic description of the physical system. The closest technique uses cellular automata and was shown to retrieve the observed transverse pattern formation, linking them to oscillatory behaviors at the laser output [15]. Aside from this latter work, the usual technique mostly involves Maxwell-Bloch equations to describe laser dynamics [16,17]. The laser description proposed here, which uses Poisson point processes for both carrier and photons, is new and was only introduced by us [18,19].  The dual-mode QW semiconductor laser model is derived from our initial Monte Carlo semiconductor laser simulator [18]. Let us recall as depicted in Fig. 1 that the energy diagram involves two bands with equally spaced energy levels accepting only one electron per level to figure out the Pauli exclusion principle. As such, electron-lattice processes are accounted for Table 1. Summary of notations and rates Notations δ energy level separation in bands n k occupancy at energy level k 1, 2 labels for the high (mode#1) and low (mode#2) energy mode m, m 1 , m 2 total, mode#1 and mode#2 photon count inside the cavity L 1 , L 2 labels of the lasing levels in conduction band ℓ 1 , ℓ 2 labels of the lasing levels in valence band α 1 , α 2 cavity losses for mode 1, 2 λ 1 , λ 2 fraction of pumping taken by mode 1, 2 (λ 1 + λ 2 = 1) Rates pn k downward thermal relaxation from level k to level k − 1 q δ pn k upward thermal excitation from level k to level k + 1

Laser model
photon escape from cavity, mode 1, 2 only with the closest neighbors. Thermal de-excitation is obtained by coupling each level with its closest downward neighbor, while thermoionic excitation is obtained with its closest upward neighbor with a rate depending on the lattice temperature. This chain of processes yields the correct Fermi-Dirac distribution in each band when no optical events occur. Electron-photon processes are accounted for using two pairs of levels (one in conduction band and one in valence band) at which lasing occurs. It figures out the two possible modes as selected by the dual-mode optical cavity. At that point, events are all considered as poissonian elementary processes, including pumping if an optical pump is considered. As a major result, the whole system behaves like a Markov process and may evolve dynamically as driven by the rates attached to each process, some of them being dependent on another process (Cox process). Finally, it can be treated either analytically by solving the corresponding Langevin equation for each energy level or it can be solved using a true Monte Carlo code [20] that exhibits a built-in account of every fine dynamics of either carrier-carrier or carrier-photon interaction, such as for example mode competition and spectral hole burning. Both are detailed in what follows.

Analytical solution for the steady-state regime
We are here interested in finding a steady-state dual-mode regime analytical solution that is consistent with the model drawn in §2.1. Establishing an analytical model based on the energy scheme defined above requires to write one rate equation per level to account for all the elementary processes. If B is the band size, we end up with 2B coupled equations indexed as given in the right part of Fig. 1. Fortunately these equations are highly coupled and chained, so that their resolution becomes possible analytically. Details about notations and rates are given in Tab. 1 where p and q = exp(−1/k B T ) are constants that parametrize thermal bath coupling and temperature respectively.
It should be noticed that this model is completely different from the rate-equations previously used in dual-mode stability analysis [10]. In the latter, carrier interaction between modes is not included, thereby prohibiting any conclusion about mode stability without the adjunction of self-and cross-saturation coefficients for optical gains. These effects are here included ab initio.
When the steady-state regime is reached, both the stability of photonic population for each mode and the stability of electronic population in each band are ensured. Combining relations (4) and (5) of the Appendix A yields a new relation between lasing levels occupancies for the same mode i = 1, 2.
For the sake of simplicity we assume here that band levels are symmetric with respect to the gap. We have therefore n L i = 1 − n ℓ i , and (1) thus reads On the other hand, as seen in Appendix A, occupancies of consecutive levels situated between the two lasing levels in conduction band are linked by a homography, which may be iterated from level L 1 to level L 2 to obtain an analytic relation between the occupancies of these levels Finally plugging (2) in (3) provides us with a relation that only involves variables λ i and laser parameters, which is the desired necessary condition for stationarity. Let us consider a dual-mode QW laser with B = 800 energy levels in each band and a energylevel spacing of δ = 1 meV. Lattice temperature is assumed to be 300 K. Eqs. (2) and (3) have been solved numerically using Mathematica, seeking for a perfect balance between the two modes (λ 1 = λ 2 = 1/2). In a first attempt, no solution appears while keeping exactly the same parameters for the two modes, i.e. either the mode gains g i or the mode losses α i . The stable dual-mode regime of this structure thus requires a slight adjustment of one cavity with respect to the other. In the second attempt, we choose to play with cavity losses and keep the same mode gain g 1 = g 2 = 1 ns −1 . Total losses were prescribed at a constant value α 1 + α 2 = 1.2 ns −1 , thereby ensuring that the same total photon count m is always obtained at the same pump intensity. Then Eqs. (2) and (3) were solved simultaneously for α 1 at various pumping values J .
The solid line of Fig. 2(a) plots the required cavity losses α 1 satisfying a perfectly balanced dual-mode operation, while the dotted lines figure out the limits where one mode becomes two times greater than the other. The α 1 -value clearly increases with m. Accordingly, any cavity adjustment is power dependent, a feature that seriously precludes any practical application. Another weakness is the very small range in α 1 , given in Fig. 2(a) by the spacing between the two dashed lines, that really allows a reasonable co-existence of the two modes. As a consequence, no care adjustment of cavity losses can be robust enough against pump fluctuations or any other external laser parameter variation such as temperature or optical feedback. This is illustrated in Fig. 2(b) with λ 2 calculated as a function of τ, the intraband relaxation time. In the framework of our model we have τ = p(1 − q δ ), and it has been varied here in a wide range of a few tenth to a few hundred of fs. The shorter τ, the steeper the transition from one mode to the other occurring at λ 2 = 1/2. As a major result, at τ ≈ 100 fs, a value that is consistent with semiconductor measurements, a nearly complete mode switch is observed with a variation of α 1 of only a few percent. These calculations agree with the well-known high difficulty to observe a stable dual-mode regime experimentally with QW semiconductor lasers.

Numerical results
In a second step, we built a Monte Carlo simulator following the model description of §2.1.
The program was derived from [18], following the prescription for a true Monte Carlo code described in [20]. We selected simulation parameters identical to those of the analytical analysis: B = 800, the level count per band, δ = 1 meV, the energy separation, T = 300 K, the lattice temperature, g 1 = g 2 = 1 ns −1 , the optical gain of both modes, α 1 and α 2 such as α 1 + α 2 = 1.2 ns −1 , the total optical losses including photon escape events in each cavity mode, 1/p = 20 fs, the thermalization rate on each energy level, and a 4 meV energy splitting between the two allowed optical modes corresponding to a 1 THz frequency beating. Monte Carlo runs were then conducted using as starting point an analytical solution of stable dual-mode operation previously determined, namely m = 350 and α 1 = 0.584. We then depart from these values by changing both α 1 and the lattice temperature. Numerous attempts were collected and are summarized in Fig. 3. In this first set of calculations, only the simultaneous existence of the two modes was sought, illustrated with a color code. These Monte Carlo simulations are in perfect agreement with the previous analytical analysis. Moreover, this example emphasizes a new dependence of α 1 with the lattice temperature, thus enforcing our previous remark of the extreme weakness of the dual-mode operation with respect to any of the laser parameters.
Monte Carlo simulations can also be exploited more extensively in order to extract features unavailable with analytical calculations, such as carrier dynamics or laser noise. This is first illustrated in Fig. 4(a) with the occupancy in the conduction band that reveals a complex carrier dynamics with strongly marked spectral hole. Surely, each optical mode depletes slightly the population, which otherwise would obey Fermi-Dirac statistics. As a result, the coupling between the two modes is conspicuous in that the more intense the depletion created by the mode#1, the less population remains for mode#2. This coupling is here brought by the Pauli exclusion principle. For the sake of completeness, it can be shown (without any assumption about bands symmetry or equidistant levels) that the difference between left and right slopes at level L i is equal to λ i J τ ln(q).
On the optical side, our Monte Carlo program allows to evaluate the normalized photodetection spectral density of the emitted photon flux [18]. The Fig. 4(b) gives these calculations for the total emitted power and for both modes which are supposed ideally filtered. Notice that a quiet pump was considered here, which gives rise to the strong subpoissonian laser statistics of the total emission. This quasi null spectral density at null frequency increases, as expected, up to a poissonian noise at very high frequency. Nevertheless the situation is completely different when considering each modes separately. They are both very noisy at low frequency with exactly the same values. By means of separate calculations we have seen that they differentiate when the two modes become unbalanced, for example by increasing the pumping. All these spectral density features are the signature of mode-partition noise. The Fig. 5 uses the point of view of photon histogram within the cavity to illustrate the sensitivity of the dual-mode regime to laser parameters and conditions. Plots are given as a function of α 1 and m and emphasize again that a very small change in the dual-mode power balance appears quickly when a parameter departs from its nominal value: Here, only 2 10 −3 in α 1 or 15% in pumping power are sufficient to discriminate the two modes that otherwise exhibit the same shape and position at α 1 = 0.584 and m = 350. Moreover, in perfect agreement with previous results on spectral density, the width and skewness of the histograms are also dependent on the laser parameters, indicating changes in noise characteristics for the laser modes.

Monte-Carlo model
Accounting for QD semiconductor laser requires another model be built in view of Monte Carlo simulation. The chosen one is described in Fig. 6. It differs from our previous QW laser model mainly because it is no more necessary to account for Pauli exclusion principle and thermalization within bands. Actually, QD lasers can be better viewed as set of quasi-atoms that behave just like an Ehrenfest urn model. This is well-known to be particularly easy to model and simulate using Markov chains. In practice, we assume only two distinct QD families centered on the two wavelengths of interest. They are supposed to operate on discrete energy levels. All events are governed by elementary Poisson processes whose corresponding rates are calculated from the known physical characteristics of the device (for instance, p i accounts for the carrier capture by the QD family i while q i p i accounts for the barrier thermoionic refilling).
On the other hand, the InGaAs or InP barrier and corresponding pumping mechanisms are also modeled with discrete levels involving very high recombination rates that figure the ultra-rapid dynamics of carriers in these bands. The same kind of rates applies to optical cavities to govern photon count inside both modes and photon escape. The model design is intimately based on coupled rate equations [10,21], but it has implicit assumptions that have to be justified. Firstly, accounting for the barrier model only by two energy levels is rough but can be justified by the very high internal relaxation that is taken in practice. Its validity is also very easy to prove during the course of the simulation, by verifying that relaxation events in barrier are predominant in count and that there is not carrier accumulation. Secondly, introducing only two QD families is a bit simplistic but the optical cavity is expected to filter two such categories among the whole set of QDs of the active layer. Within this framework, the other QDs are equivalently accounted for by an increase of the nonradiative recombination in InP, with great benefit for CPU time. Getting more inside model details, one notices that a fixed number of elementary cells is accounted for instead of a fixed number of QDs. This simplification was mandatory for our Monte Carlo model elaboration, and it is solved during calculations by only using reasonable pumping levels that do not saturate the proportion of excited QDs.
Other parameters were all dictated by the known data for InAs QD [11,[22][23][24]. Optical gain constants g 1 and g 2 were taken at 1 ns −1 like for QW lasers, the cavity losses α 1 & α 2 correspond to a 5 ps photon lifetime, the thermal relaxation in barrier corresponds to a 100 fs intraband relaxation time in InP, typical values for the capture time in QDs was chosen at 3 ps and thermoionic refilling obtained in accordance with Boltzmann law. Finally, a 1 THz splitting was chosen between modes, and in each run we simulated 10 6 elementary cells as depicted in Fig. 6.

Numerical results
First simulations were conducted in order to evaluate the steady-state regime when pumping is varied. Namely, we computed I(V ) curves with the above-defined parameters that are identical for both modes, just except for their gap energies that only affect the thermoionic refiling of QD#2 family with respect to QD#1. Because we select a non-radiative recombination rate p NR corresponding to a characteristic time of 7 ps, the I(V ) curve must present a threshold just like a true laser. Results are given in Fig. 7. The most important conclusion here is that both modes are lasing whatever J and that their outputs are yet superimposed. The latter remark is justified by the same parameters used for gain and losses of the two modes. Furthermore, all plots show the same threshold occurring at J th = 0.12 ns −1 . The first major conclusion with this QD laser model is thus a complete decoupling between optical modes that makes the dual-mode operation much more stable than for QW laser. A close-up view of what occurs just above threshold and far above threshold is given in Fig. 8. In two cases we have plotted the photon histograms in the cavity for both modes and for the total emission, as well as the corresponding normalized spectral densities. Histograms emphasize the simultaneous emission of both modes with similar bell-shaped curves centered at the same m-value, that is with the same intensity. The total emission also exhibits a bell-shaped curve, but centered at twice the value of each mode as expected from the sum of intensities, and with a slightly enhanced width corresponding to a more noisy emission. This proves that modes are mostly uncoupled and that their intrinsic noises add up in the total emission.
The comparison between the two pumping conditions, 1.33J th and 2.5J th , reveals the reduction of the ratio between the histogram width to its mean value when pumping increases. In other words, the Fano factor is reduced when J increases, in complete agreement with the well-known reduction of optical noise when departing from threshold. The plots of the normalized spectral density S (Ω) confirm the more noisy characteristic of the dual-mode laser close to threshold with much higher values in Fig. 8(b) than in Fig. 8(d). Moreover, it is evident for both pumping conditions that S (Ω) for the total emission is always twice that of one mode. It confirms the decoupling of the two modes by the canceling of mode-partition noise that was observed in QW laser.
An extensive parametric study has been conducted seeking for the CW stability conditions of such dual-mode microlasers as a function of external parameters (pumping, temperature) as well as intrinsic characteristics (gain dispersion, non-radiative recombination in QD and barrier, excited state exchange via the wetting layer or thermoionic refilling of the barrier). Whatever the conditions, the numerical Monte Carlo calculations show that the QD laser always operates dual-mode. This is illustrated in the following by two examples.
In Fig. 9(a) we consider the influence of mode gains g 1 and g 2 on the occurrence of a stable dual-mode regime. The intensity plot of each mode is given as a function of g 1 while keeping the product g 1 g 2 constant. When g 1 increases, mode#1 intensity increases and mode#2 intensity decreases. It ends up with the extinction of mode#2 at g 1 ≈ 2.5 ns −1 because the gain difference with mode#1 becomes too big. Nevertheless, the stable dual-mode emission is observed over this entire wide range of variation without any sudden mode extinction, the balance between In a second time, we consider in Fig. 9(b) the evolution of the power balance between the two modes as a function of τ, the carrier capture time in QD. We introduce an initial asymmetry between the two mode powers for illustration convenience and calculate then the steady-state photon values when τ is increased from 25 fs to 5 ps. In spite of this extreme range, the initially introduced asymmetry remains identical whatever the τ-value, up to the extinction of the two modes when τ becomes too high. The influence of the confinement energy was also explored numerically by varying the energy confinement from ≈ 50 meV to ≈ 0.5 eV. Again, the dual mode regime was obtained for all values and moreover the mode intensities and balance are almost the same on this whole range of variation.
As a final conclusion, Monte Carlo simulations conducted in the framework of our model ensure the robustness of the dual-mode regime of QD laser with respect to any laser parameter variation.

Conclusion
We have shown from Monte Carlo calculations that QW and QD lasers behave very differently considering their intrinsic dynamics. This was illustrated considering the mode competition that may preclude the possible stable dual-mode regime allowed by doubly resonant cavities. Similar structures and parameters were evaluated with these two kind of active materials. QW lasers were shown bistable, that is single-mode, for almost every set of material and cavity parameters considered. This is particularly the case for devices with balanced gain and loss parameters for the two modes. Stable two-mode regime was however demonstrated with unbalanced gain and loss parameters. Nevertheless, it occurs in a so tiny range of applied pumping and temperature that it should probably not be observed experimentally.
To the contrary, when considering QD lasers a stable dual-mode regime was obtained in every case. It proves that coupled-cavity laser including QD may operate CW in dual-mode far more easily than devices including QW. To our eyes, this stability comes from the carrier capture process in QDs that is independent for each dot. The coupling between two QDs requires that a captured carrier in the first QD get excited to the barrier and then captured by the second QD. This process is quite indirect even with small energy barriers. The corresponding process that couples two lasing levels in QW is far much direct and involves only a cascade of relaxation that are hardly linked by the Pauli exclusion principle. This is commonly expressed in other words by the observation of a nearly perfect overlap between inhomogeneous and homogeneous gain broadening in QW, whereas they remain quite different for QD, at least up to room temperature. Together with an analytical model [10], these numerical Monte Carlo support the design guidelines of a new stable dual-mode semiconductor laser operating as a THz beating source [14].