Self-generation of optical frequency comb in single section Quantum Dot Fabry-Perot lasers: a theoretical study

Optical Frequency Comb (OFC) generated by semiconductor lasers are currently widely used in the extremely timely field of high capacity optical interconnects and high precision spectroscopy. Very recently, several experimental evidences of spontaneous OFC generation have been reported in single section Quantum Dot (QD) lasers. Here we provide a physical understanding of these self-organization phenomena by simulating the multi-mode dynamics of a single section Fabry-Perot (FP) QD laser using a Time-Domain Traveling-Wave (TDTW) model that properly accounts for coherent radiation-matter interaction in the semiconductor active medium and includes the carrier grating generated by the optical standing wave pattern in the laser cavity. We show that the latter is the fundamental physical effect at the origin of the multi-mode spectrum appearing just above threshold. A self-mode-locking regime associated with the emission of OFC is achieved for higher bias currents and ascribed to nonlinear phase sensitive effects as Four Wave Mixing (FWM). Our results are in very good agreement with the experimental ones.


Introduction
An Optical Frequency Comb (OFC) is a light emission characterized by equally spaced narrow optical lines with same intensity, low phase noise and low mode partition noise [1]; an OFC laser is therefore different from a multi-wavelength laser, because, in the latter, the lasing lines may be not phase locked and the coupling between them may lead to important mode partition noise.
In semiconductor lasers, OFC have been obtained with mode-locking techniques (either active mode-locking via RF electrical modulation of the laser gain/losses or passive massive mode-locking with reversely biased saturable absorber in laser cavity) or gain-switching. A breakthrough has been represented by the experimental evidences of self-mode-locking in single section lasers characterized by the absence of any active or passive modulation; self-generated OFC has been experimentally reported in different types of Fabry-Perot lasers based on Quantum Well (QW) [2], Quantum Dash (QDash) [3], QD [4] active materials and in Quantum Cascade Lasers (QCLs) [5]. Self-generated OFCs in 1.3 µm and 1.5 µm wavelength bands have been extensively investigated in QD and QDash FP laser under CW injection [6,7,8,4,9,10,11,12,13]; pulse trains have been reported directly at the laser output [6,7,8,4,9] or after compensation of the group velocity dispersion trough, e.g., a dispersive optical fiber of suitable length [10,11,12,13]. Despite the observation of pulses is an unambiguous fingerprint of the OFC formation, we evidence that OFC lasers generate lasing lines that are in general phase locked, even when there is no regular phase relation between the lines [5]. For many applications, as those discussed below, the phase locking condition (leading to narrow RF beat-note and narrow optical linewidths) is the sole requirement and regular output pulses are not strictly needed.
These very simple OFC sources have attracted an impressive interest in the rapidly growing field of high-data rates DWDM optical interconnection where the OFC laser diode feeds the silicon photonics optical modulators to realize a compact and low cost transmitter [14,15,16]. Another interesting application is in the field of photonic generation of sub-THz signals via the mixing of the comb lines, spaced of hundreds of GHz, on a fast photodetector [17]. QCL in the mid-infrared range are used for high precision and high speed spectroscopy based on dual comb techniques [18,5].
The understanding of the fundamental physical effects at the origin of self-generation of OFC in FP semiconductor lasers is still matter of present research. Multi-wavelength lasing close to threshold is generally ascribed to fast carrier gratings with sub-wavelength spatial variation which is generated by the optical field standing wave pattern in the FP cavity (see e.g. [19]). This longitudinal modulation of the carriers (spatial hole burning, SHB) destabilizes the single longitudinal mode emission and favors the onset of several longitudinal cavity modes with different spatial configuration. SHB is rather week in QW lasers because the carrier diffusion tends washing out the carrier gratings while the cross-coupling between longitudinal cavity modes is responsible for the high mode partition noise [20,21].
In low dimensional active materials as QDs and QDashes, the carriers are trapped in the quantum box and cannot diffuse, therefore the SHB is much stronger than in the QW counterpart. From this viewpoint, QD and QDash laser diodes are similar to QCLs, where carrier diffusion time is longer than the typical carriers life time [22]. Despite SHB alone can explain the multi-wavelength lasing, it is not sufficient to enable the phase-locking among the lasing spectral lines [23,24,25].
Recent theoretical analysis carried on for the QCL case (see e.g. [26]) have attributed the phaselocking and the consequent OFC generation to nonlinear phase sensitive effects, like FWM, which are inherent to the saturable character of radiation-matter interaction in SCLs and impose phase-matching conditions for sufficiently high intra cavity field (ie: high CW bias current) and when the dispersion is limited. For standard QW laser diodes, the rigorous description of the role played by the FWM effect in the phase locking of three cavity longitudinal modes has been presented in [27]. In [27], the Authors also identified the operating conditions where the FWM is able to compensate the non uniform frequency spacing among the cold cavity modes, thus leading to self-phase-locking.
In our previous work [28], we have studied the formation of OFC in FP single section QD lasers using a Time Domain Traveling Wave model. In that work, starting from the assumption that the carrier gratings originated by the standing wave pattern could lead to a compression of the gain at the wavelength of the lasing mode [29,30], we have introduced an empirical self-gain compression factor ( -parameter) of the QD material gain and we have demonstrated that it led to a multi-wavelength emission and, in some cases, to spontaneous OFCs generation.
Here, upgrading that model [31,28], we present a TDTW model that eliminates the need of an empirical -parameter and includes in a rigorous way the effect of SHB. With rigorous numerical simulations and analysis of the results, we demonstrate that the SHB is, in QD systems, the fundamental physical effect at the origin of the multi-wavelength lasing spectra always observed experimentally at any current just above threshold. Moreover, thanks to our model, we can identify in this work intervals of the bias current where, together with multi-wavelength emission, we have also an actual OFC thanks to FWM which is responsible for mutual injection locking of the cavity modes.
The paper is organized as follows: in section 2, we describe the TDTW model used for simulating the multi-mode dynamics of the QD laser and, in section 3, we present the results of numerical simulations for a standard QD single section FP laser. In section 4 we compare our observations on comb self-generation with the experimental evidences reported in the existing literature and we provide a physical interpretation of the self-phase locking. Finally we draw our conclusions in Section 5.

QD laser TDTW model in presence of carrier grating
We consider a single section Quantum-Dots-in-a-Well (DWELL) InAs/GaAs laser emitting from the ground state (GS) around 1258 nm. The length of the FP cavity (L) is a few hundred microns. We sketch in Fig. (1) the typical device configuration (a) and the electron dynamics (b) as considered in our model, whereas the main material and device parameters are summarized in Table 1.
The SHB is introduced in the model adopting a multiple scale approach where a slowly varying envelope approximation (SVEA) is adopted for describing the spatio-temporal evolution of the optical electric field; this approximation allows discretizing the cavity length with spatial steps of few microns including, nevertheless, the effects of the fast standing wave pattern on carrier dynamics [23,19,32,33].
In the SVEA, the optical electric field of the fundamental TE mode is written as: where φ(x, y) is the transverse mode profile that we assume independent on z, ω 0 is our reference angular frequency and k 0 is the propagation constant of the cold cavity TE mode at ω 0 . The terms E(z, t) ± represent the slowly varying spatio-temporal variation of the forward and backward components of the electric field that obey the following equations ( [19]): In the FP configuration the forward and backward field envelops satisfy the boundary conditions: where r R,F are the reflectivities at the rear (z = 0) and front (z = L) facets, respectively. The variable P (z, t) is the complex macroscopic polarization slowly varying in time (having extracted the component exp(jω 0 t)), but still including the fast spatial variation in the z-direction. The symbol < · > denotes in (1) a spatial average over many wavelengths along the z-direction; v g is the group velocity of optical mode in the waveguide, Γ xy is the transverse optical confinement factor in the active region, α wg are the waveguide losses per unit length and S ± sp is the noise source associated with the process of spontaneous emission. The relation between the macroscopic polarization P and the microscopic polarizations p m (z, t) due to the GS and Excited State (ES) levels of each QD of the ensemble is: where N D is the number of QDs per unit area, h QD is the thickness of a single QD layer, µ m is the degeneracy of the m-th level (with m = GS, ES) and d m is the dipole matrix element associated to the optical transition from level m-th. The variables p m (z, t), with m = GS, ES, satisfy the following two-level like dynamical equations: wherehω m is the energy associated to the GS or ES transition and 1/Γ is the dephasing time which in our model is set to fit the gain dispersion of the QD ensemble. The variables ρ m (z, t) represent the electron occupation probability in the m-th level and keep the fast spatial variation of the carrier density due to the SHB; we assume that the hole dynamics will follow the electron dynamics (excitonic approximation). We adopt in the following the scalings: where w is the ridge area and η L is the number of QD layers, such that the quantity |E ± (z, t)| 2 represents the total power across the device transverse section.
The dynamics of the carriers in the GS, ES and in the wetting layer (WL) is described via a set of rate equations that read [31]: where the stimulated emission rate per unit area has the following expression: where η is the effective refractive index. In Eqs. (4)-(5), ρ W L is the occupation probability of the WL which is represented as a single energy level with degeneracy per unit area D W L , η L is the number of QD layers, Λ = I η i /(e D W L η L w L) is the carrier injection rate in the WL, η i is the internal quantum efficiency, and I the pump current. Referring to Fig. 2, the characteristic times τ GS C,e represent the capture (C) and the escape time (e) between GS and ES. Similarly, the capture/escape processes between WL and ES is quantified by the time constants τ ES C,e . The carrier lifetime due to spontaneous and non-radiative recombinations in the states is quantified by τ m sp,nr with m = GS, ES, W L. Since we study here laser emission from the GS only, we will assume that the photon emission process occurs only in the GS.
To couple the rate equations (3)-(5) with the slowly varying traveling wave equation of the electric field in Eq. (1), we expand in Fourier series the spatial variation at the wavelength scale of both the polarizations and the occupation probabilities: where p ± n,GS (z, t), ρ n,m (z, t) with m = GS, ES, W L and n = 0, 1, 2... are slowly varying functions of z. Inserting these expression in Eqs. (1)-(5) and neglecting the terms with spatial frequency higher than 2k 0 , we finally get the following system of differential equations: where we used the notation: In the dynamical equations (6a)-(6h), we do not include the effect of higher order terms such as the coherence grating (p 1,GS ). These terms would leave qualitatively unchanged the system multi-mode behavior [24].
We observe that, if we neglect the carrier grating at the sub-wavelength scale by forcing (ρ + m = 0 with m = GS, ES, W L), the model reduces to the TDTW model presented in [31,28].

Results of dynamical simulations
In this Section we report numerical results on the multi-wavelength QD laser dynamics obtained by integration of the PDEs (Eqs. (6a)-(6h)) using an optimized finite difference algorithm described in details in the Appendix. Since the temporal variations of the carrier densities occur in the time scale from hundreds of femtoseconds to few picoseconds (which are slower than the polarization dephasing time 1/Γ set to 15 fs) we can solve the QD active material frequency dependent optical response as an infinite impulse response filter as described in the Appendix. We discuss first the key role of the fast carrier grating in generating the multi-wavelength lasing regime at all currents above threshold and then we focus on the conditions for having the phase-locking of the lasing lines. In the following analysis, we start from the total field at the laser output facet E out (t) = 1 − r 2 R E(L, t) that allows to calculate the total power emitted by the QD laser as P (t) = |E out (t)| 2 . By numerically filtering E out (t) with Hanning windows of spectral width W h = 0.8F SR centered at the relevant peaks of the optical spectrum, we are able to isolate the temporal evolution of the field of each q-th lasing line E q (t) and then retrieve the power P q (t) = |E q (t)| 2 and the instantaneous phase φ q (t) = E q (t) .
In order to characterize the multi-mode dynamics in the different simulated regimes, we introduce convenient parameters that will measure the intensity noise and the amount of phase locking. We will use the notation < · > to indicate the temporal average, the symbols µ q and σ q to indicate mean and standard deviation in time for each line, and M to indicate the average over all the N optical lines within the −10 dB optical bandwidth. These parameters are: • The average output power µ P q =< P q (t) >, the intensity noise δP q (t) = P q (t) − µ P q and the r.m.s of the intensity noise σ δPq = < (δP q ) 2 > of the q-th line; we will use the average value over the N lines M σ δP = 1 N Σ N q=1 σ δPq as a measure of the average intensity noise.
• The average optical frequency of the q-th line µ f q =< 1 2π d dt φ q (t) > and the associated phase fluctuations Φ q (t) = φ q (t) − 2πµ f q t.
• The differential phase fluctuation between two adjacent lines ∆Φ q (t) = Φ q+1 (t)−Φ q (t) and its av- is a measurement of the differential phase noise, whereas M σ∆Φ = 1 N Σ N q=1 σ ∆Φq is the average differential phase noise.
• The frequency separation between couples of adjacent lines ∆µ f q = µ f q+1 −µ f q ; M ∆µf represents the average value for the modes within the −10 dB optical bandwidth. Our model does not include the waveguide dispersion, but, because of the QD active material dispersion, the longitudinal mode separation under current injection might be different from the cold FSR and it might change among couple of modes and with current.
Respect to parameters defined above, an ideal OFC would have for any q-th line • σ δPq =0 indicating that the mode partition noise is completely suppressed • σ ∆Φq = 0 as consequence of the phase-locking, i.e. no phase fluctuation and therefore extremely narrow optical linewidth • the same value of ∆µ f q for any q-th line because, by definition, OFC requires equally spaced optical lines.
In order to better compare our results with the experimental evidences so far reported in the literature (see e.g. [34]), we will also provide in this Section a calculation of more conventional (and experimentally more easily accessible) quantifiers of the degree of coherence in the system. These are: • The Relative Intensity Noise (RIN) spectrum calculated for the total output power and the power of each line as, respectively: where F{·} denotes the Fourier transform operator.
• The integrated RIN in the electrical bandwidth B ranging from 10 MHz to 21 GHz calculated for the total output power and for each line as, respectively: the quantity M iRIN represents the average value of iRIN q , calculated for each q line in the −10 dB optical bandwidth.

Role of SHB in the multi-mode dynamics
In the following, we report simulation results for a 250 µm long FP laser with cold cavity FSR of 178 GHz.
To study the effect of the carrier grating, in Fig. 2 we have prepared the laser with an injected current of 200 mA above the threshold current and, for t < 0, we have forced ρ + m = 0 with m = GS, ES, W L thus neglecting Eqs. (6d), (6f) and (6h). As expected, the system has a stable Continuous Wave (CW) solution consisting in only one lasing mode (q = 0) that is the one closest to gain peak [19,28]. At t ≥ 0 we "switch-on" the fast carrier grating in the system. We plot the temporal evolution of P (t) in Fig. 2a, P q (t) in Fig. 2b, and a zoom of the transient of some exemplar optical lines in Fig. 2c. The standing wave pattern induced by the lasing mode q = 0 generates a carrier grating at the wavelength scale which has two main effects. First, it compresses the gain of the mode q = 0, triggering a significant reduction of its modal amplitude (as shown in Fig. 2b,c) and a consequent increase of the carrier density (ρ 0,m , with m = GS, ES, W L) above threshold (SHB). Second, this spatial modulation of the carriers allows the side lines gain (q = 0) to overcome the cavity losses in a cascading mechanism that brings lines with larger distance from the central one to be later excited (Fig. 2b,c). In agreement with the experimental findings (see e.g. [9]), we observe that SHB makes the laser operating in multiwavelength regime just above the lasing threshold. Figure 3 shows the power versus current curve for the simulated QD laser from threshold (∆I = I − I th = 0 mA) to ∆I = 400 mA and, in the insets, two representative optical spectra corresponding to ∆I = 40 mA and ∆I = 200 mA. In this current range we see an increment of the −10 dB optical bandwidth from ∼14 nm (just above threshold) to ∼18 nm.
After having identified in the SHB the mechanism responsible for the multi-wavelength lasing, we analyse in the rest of this Section two different dynamical regimes encountered by increasing the current above threshold: a multi-mode regime with phase-unlocked modes (indicated in the following as unlocked regime) and a multi-wavelength regime with phase-locked modes (indicated in the following as locked regime). The latter corresponds to the actual OFC generation. These regimes will be mainly characterized by mapping the temporal evolution of the quantities P q (t), ∆Φ q (t) and estimating the OFC indicators defined in the previous paragraph.

Unlocked regime
As an example, we report the results of an unlocked regime observed at 40 mA above threshold. Figure  4 shows P q (t) and ∆Φ q (t) for a few representative lines of the optical spectrum; we plot next to the right vertical axis also the temporal averages µ P q and µ ∆Φq (circular symbols) and the associated intensity noise σ δPq and differential phase noise σ ∆Φq (bars amplitudes). All the parameters have been calculated for an observation time window of 1 µs. The same statistical moments are also plotted in Fig. 4c for the lasing lines within the −10 dB optical bandwidth and, together with the average frequency separation between couples of adjacent lines (∆µ f q plotted in Fig. 4d), they represent convenient indicators of the self-generated OFCs. The evidence that differential phase fluctuations are uncorrelated in time, that differential phase standard deviation bars are quite large (> 1 rad/s) for all lines and that optical frequency separation is quite dependent on the mode number allows us to conclude that we are dealing with an unlocked regime.

Locked regime: self-generated OFC
For higher bias current, for example 200 mA above threshold, the system dynamics changes radically as shown in Fig. 5. Compared to the case in Fig. 4, we see a consistent reduction of both the intensity noise and the differential phase fluctuations and a significant uniformity of the optical line spacing; all indicating that the phase-locking with OFC generation has been achieved. Since the differential phases are practically constant with time, but not the same for all the lines, it is not possible to observe any optical pulse at the laser output (as it happens in passive mode-locking) despite that the modes are  phase-locked.
We compare the two regimes (locked and unlocked at 200 mA and 40 mA above threshold respectively) in terms of degree of coherence of the system by reporting quantifiers which can be more easily measured: the RIN spectra (Fig. 6a,b), the beat-note RF linewidth (Fig. 6c), and the RF spectrum in an extended frequency range (Fig. 6d). In the unlocked regime we observe high RIN (especially in the low frequency part of the spectrum between 1 MHz and 1 GHz) and very broad RF line at the beat-note. On the contrary, the locked regime is characterized by a reduced RIN (it decrease down to around −170 dBc Hz −1 for total power RIN and to around −140 dBc Hz −1 for RIN of mode q = 0) and extremely narrow beat-note. Finally, the RF spectrum in a broader frequency range reveals in both cases the presence of more than 10 lines in the first decade.
The temporal evolution of the total power in a locked regime is dominated by regular oscillations with ∼5 ps period (corresponding to the inverse of the beat-note) superimposed to sub picosecond oscillations associated with the high frequency peak around 2 THz in the RF spectrum of Fig. 6d. We associate this high frequency peak with the Rabi resonance [35] and we consider it as a signature of the coherent interaction between the intracavity electric field and the ensemble of "artificial atoms" represented by the QDs in the active material . Experimental and theoretical evidences of Rabi oscillations at room temperature in QD materials have been found in QD based semiconductor optical amplifier in condition of high intensity pulse injection [36,37] but, at the best of our knowledge, no signature of this phenomenon was reported in QD laser diodes. We however stress here that the observation of the Rabi resonance is in this context only the consequence of the excitation of the QD medium with the coherent broad band electric field self-generated by the FWM process, as discussed in the next Section. Therefore, contrary to the prediction of the Risken-Nummedal analysis for a two level system laser [38,35], the excitation of the Rabi resonance is not the key physical mechanism that triggers the phase-locking of the optical lines. Indeed, if the Risken-Nummendal instability was responsible for the phase-locking, we would have found only locked lasing lines separated by ∼2 THz [39].
In order to describe the transition from unlocked to locked regime, we plot in Fig. 7 the average OFC indicators (M σ δP , M σ∆Φ , ...) versus the bias current. The beat-note spectra for the same values of the electrical pumping are reported in the image of Fig. 8. The typical scenario observed by increasing ∆I from 0 mA (laser threshold) to 300 mA consists initially in a multi-wavelength unlocked regime (denoted with letter A in the figure) associated to large intensity noise (Fig. 7a and c), large differential phase fluctuations (Fig. 7b) and very unevenly spaced optical lines (Fig. 7d). The RF beat-note is also Figure 4: Unlocked multi-mode dynamics for a pump current 40 mA. Temporal evolution of the modal power (P q (t)) (a) and of phase difference between adjacent modes (∆Φ q (t)) (b) for selected modes.
The average values (µ P q and µ ∆Φq ) and standard deviations (σ δPq and σ ∆Φq ) are reported next to the right vertical axis. In panel (c) these statistical moments are reported for all the modes within −10 dB bandwidth in the optical spectrum. The temporal average of the frequency separations between couples of adjacent modes (∆µ f q ) (left vertical axis) and its difference with respect to the overall average over the selected modes (M ∆µ f ) (right vertical axis) are represented in panel (d).
rather broad (Fig. 8). This regime is then followed by clear phase-locking regime (denoted by the letter C in the figure) where both intensity noise and differential phase noise is drastically reduced. We also observe lower integrated RIN, constant optical line spacing and narrow RF beat-note linewidth. In the region between ∆I = ∼60 mA and ∆I = ∼100 mA (denoted by the letter B in the figure) the two configurations coexist for the same current parameter and the final steady state depends from the system history (initial conditions). After ∆I = 100 mA we observe self-generated OFC apart from few values of current where the system is unlocked. These results well agree with the experimental evidences reported in [34] for a single Section QDash laser where it is shown that an increase of the pump current above a certain value brings the laser to move from an unlocked regime (measured as high integrated RIN and large RF line at the beat-note) to a locked regime characterized by a significant drop of the integrated RIN and reduction of the RF linewidth. (c) (d) Figure 5: Self-generated OFC regime for a pump current of 200 mA. Temporal evolution of the modal power (P q (t)) (a) and of phase difference between adjacent modes (∆Φ q (t)) (b) for selected modes.
The average values (µ P q and µ ∆Φq ) and standard deviations (σ δPq and σ ∆Φq ) are reported next to the right vertical axis. In panel (c) these statistical moments are reported for all the modes within −10 dB bandwidth in the optical spectrum. The temporal average of the frequency separations between couples of adjacent modes (∆µ f q ) (left vertical axis) and its difference with respect to the overall average over the selected modes (M ∆µ f ) (right vertical axis) are represented in panel (d).

Discussion
The numerical evidences reported in the previous Section regarding the phase-locking of the lasing lines might be interpreted by extending to several modes the results derived in [27] for the study of self-pulsing DBR semiconductor lasers (based on quantum well active medium) where three cavity longitudinal modes are locked. Whereas in that case the three lasing lines were initially unevenly separated because of the DBR dispersion, in the case we are simulating this non-uniform separation derives from the dispersion of the QD active medium. The nonlinear interaction between the electric field and the gain medium (Four Wave Mixing) generates sub-bands close to each lasing line. These extra lines act as internally self-generated optical injection terms, pulling the instantaneous angular frequency of the main optical lines. For sufficiently high bias currents, when the intracavity field is strong enough, this optical injection process converges to the only possible equilibrium configuration, where all the lines are all equally spaced in frequency and exhibit a constant phase difference. This configuration corresponds to the OFC self-generation shown in Fig. 5. On the contrary, in the low current regime, the power of the generated side band lines is not strong enough, and the FWM lines will not be able to force the locking: as a result, the system remains unlocked. In this condition, the high RIN of each single line in the low frequency region (reported for example as the red curve of Fig. 6a) is the consequence of the beating between the side-bands and the main lasing line.
This behavior is clearly observable in Fig. 9, where the optical line q = 0 are presented for injected currents ranging from 40 mA to 200 mA above threshold. For low values of the injected current (Fig. 9a) the FWM generated side lines around the main lasing peak are clearly visible, indicating that the modes power is not high enough for the FWM to induce a clean locking of the optical lines. Increasing the current (and thus the power) the side lines amplitudes reduce (Fig. 9b,c) and completely disappear in the locked condition achieved (Fig. 9d). A similar self-mode-locking mechanism has been also introduced in [40,26] to interpret the comb formation in broad band in quantum cascade lasers .
Finally, on a more fundamental level, we may identify the global indicators of the degree of coherence in the considered multi-wavelength QD laser, like M σ δP and M σ∆Φ plotted in Fig. 7, as "order" parameters for this complex physical system. Hence, their discontinuous transition across the boundary between the unlocked ("disordered") and locked ("ordered") regimes ("phases") when bias current is varied is a manifestation of the Landau symmetry principle [41], according to which a phase transition involving a spontaneous symmetry breaking is associated to an order parameter which is non-analytical at the critical (transition) point. In particular, in first order transitions the discontinuity affects the   order parameter itself, as seems to occur in our case. This analogy with first order phase transition might also help explaining the coexistence of the locked and unlocked "phases" for a given interval of pump current (region B) that would play the role of an inverse temperature ( [42,35]) and might indicate that the transition from unlocked to locked regimes occurs through a nucleation mechanism.

Conclusions
We simulated the multi-mode dynamics of a single section Fabry-Perot QD laser using a TDTW approach that correctly describes the material gain and dispersion and takes into account for the carrier grating generated by the optical standing wave pattern formed in the laser cavity (SHB) and that is not washed out by diffusion as it happens in 1D or 2D semiconductor active regions. The SHB is at the origin of the single mode instability that leads to irregular multi-wavelength emission close to the lasing threshold (unlocked regime). For higher currents, degenerate and non degenerated FWM assure mode proliferation and phase-locking typical of OFC emission (locked regime). We characterized the degree of coherence in the system using both a spectrally resolved analysis through the estimation of the optical lines amplitude and the phase temporal fluctuations and more conventional, experimentally accessible, global quantifiers as the integrated RIN and the RF beat-note linewidth. Plotting these quantities vs. the pump current we were able to identify a value of bias current above which a quite sharp transition between the two regimes occurs and we observe low RIN and narrow inter-modal linewidth configurations associated with OFC self-generation. Our results are in good agreement with recently reported experimental evidences of self-mode-locking in QDash and QD single section lasers and add a physical insight into the understanding of these self-organization phenomena with undoubted impact on its applications in the field of e.g. high capacity optical communication.

Appendix
The coupled differential equations (6a)-(6h) can be numerically solved using a finite differences scheme, discretizing both in time (with step ∆t) and space (with step ∆z = v g ∆t). In order to significantly reduce the computational time required by the numerical integration, it is possible to observe that the dynamics of the occupation probabilities ρ 0 GS and ρ + GS is slower than the characteristic dephasing time 1/Γ of the corresponding inter-band transition [31]. In this adiabatic limit, we obtain for the microscopic polarization the expression p ± k,j = j d GS hΓ (2ρ 0 GS k,j − 1)I k,j + 2ρ ± GS k,j I ∓ k,j where the terms I ± = e −Γ∆t I ± k,j−1 + Γ∆t 2 e −Γ∆t E ± k,j−1 + Γ∆t 2 E ± k,j represent the numerical approximation of the convolution integral between the field components and the material Lorentzian optical susceptibility and the subscripts k and j are the space and time steps counters, respectively. The approximation E ± k,j of the forward and backward components of the electric fields E ± (k∆z, j∆t) in Eq. (6a) can be expressed as obtaining therefore with G 0;k,j = g 0 (2ρ 0;k,j − 1) and G ± k,j = 2g 0 (ρ ± k,j ); for sake fo brevity the terms related to the spontaneous emission, calculated as in [31], are not reported in Eqs. (7) and (8).
Introducing the spatial and temporal dependent coefficients A 0;k,j = ∆z (Γ∆tG 0;k,j − α wg ) /4, A + k,j = ∆zΓ∆tG + k,j /4, A − k,j = A + * k,j and C ± (k∓1,k),(j,j−1) = ∆z 2 − α wg 2 E ± k±1,j−1 + G 0;k,j e −Γ∆t I ± k,j−1 + which has to be solved in each longitudinal slice to calculate the new approximation of the fields E ± k,j . It has to be observed that the coupling between the forward and backward components of the field depends only on the coefficients A + k,j , related to the spatial carrier grating.

Acknowledgments
This work has been supported by the Fondazione CRT under the action "La ricerca dei Talenti". The Authors would like to acknowledge Prof. I. Montrosset for the fruitful discussions and helpful suggestions.