Delay-induced dynamics and jitter reduction of passively mode-locked semiconductor lasers subject to optical feedback

We study a passively mode-locked semiconductor ring laser subject to optical feedback from an external mirror. Using a delay differential equation model for the mode-locked laser, we are able to systematically investigate the resonance effects of the inter-spike interval time of the laser and the roundtrip time of the light in the external cavity (delay time) for intermediate and long delay times. We observe synchronization plateaus following the ordering of the well-known Farey sequence. Our results show that in agreement with the experimental results a reduction of the timing jitter is possible if the delay time is chosen close to an integer multiple of the inter-spike interval time of the laser without external feedback. Outside the main resonant regimes the timing jitter is drastically increased by the feedback.


Introduction
Passively mode-locked (ML) semiconductor lasers are of broad interest as sources of ultrashort picosecond and sub-picosecond optical pulses with high repetition rates. These sources are among other applications needed for data communication, optical clocking, high-speed optical sampling, all-optical clock recovery and microscopy [1][2][3]. Compared to actively [4] and hybrid ML lasers they are simpler to fabricate and to handle. However, their major drawback is their relatively large timing jitter due to the absence of an external reference clock [5]. Since timing fluctuations degrade the performance of the laser, a great deal of efforts have been made to reduce them. An efficient and inexpensive method for timing jitter reduction based on the use of all-optical feedback from an external cavity will be discussed in this paper.
One basic timescale of the passively ML laser is the inter-spike interval time T ISI,0 , which is the inverse of the pulse repetition frequency. The optical feedback introduces an additional timescale into the system, namely the roundtrip time of the light in the external cavity τ . By measuring the radio frequency linewidth of a passively ML quantum well laser as done already in 1993 by Solgaard and Lau [6], it was shown that optical feedback can have a stabilizing or a destabilizing effect on the laser dynamics depending on the ratio of τ and T ISI,0 . Further, the authors observed a pulling of the repetition rate of the pulses by the variation of τ . Both effects have recently been observed in experiments with quantum dot lasers [7][8][9][10]. In addition, the dependence of the timing jitter on the pump current and the feedback strength has been studied for a fixed ratio of τ and T ISI,0 [8,11] and for variable τ [12].
In contrast to the wealth of experimental results, only a few theoretical works have been published on this subject. The bifurcation scenarios of a free-running quantum dot laser with optical feedback have been studied in [13][14][15]. In [16] the authors studied the amplitude and timing jitter of a quantum dot ML laser subject to optical feedback with τ = 5T ISI,0 (intermediate feedback delay time). In [17] the authors discuss the impact of a variation of the ratio τ/T ISI,0 on the dynamics of an ML quantum well laser for τ -values ranging from τ ≈ 0.3T ISI,0 to τ ≈ 3.1T ISI,0 (short feedback delay time). They identify different dynamical regimes of ML laser operation: main resonances are found when τ is an integer multiple of T ISI,0 , higher order resonances take place when τ/T ISI,0 is a rational number and nonresonant regimes are detected for the values of τ/T ISI,0 in between two resonant regimes. The dynamics in each regime will be discussed in detail in section 4 of this paper. For a short external cavity, the impact of different ratios of τ/T ISI,0 on the pulse characteristics of a quantum dot ML laser model has recently been studied in [18].
In the theoretical studies mentioned above, finite-difference traveling wave models have been used. Our study is based on a delay differential equations (DDE) mode-locking model, which has the advantage of a strongly reduced computational cost. On the one hand, this permits us to study the regime of longer external cavities with τ up to 70T ISI,0 . This is of experimental interest, because in experiments optical feedback is often realized by long optical fiber loops [12]. On the other hand, it enables us to continuously study the dependence of the timing jitter on feedback strength and delay time in intervals around the exact resonance condition and not only for the exact resonance discussed in [16].
The paper is organized as follows. In section 2, an extension of the DDE model that takes into account feedback from an external mirror is discussed. The dynamics of the laser without feedback is reviewed in section 3 and the impact of the optical feedback is studied in section 4. Finally, in section 5 the reduction of the timing jitter by optical feedback is discussed, and concluding remarks are made in section 6.

Delay model for a mode-locked laser
We consider the DDE model for a passively ML ring cavity laser proposed in [19][20][21]. This model generalizes the model of Haus [22,23] to the case of large gain and loss per cavity roundtrip, which is typical of semiconductor ML lasers. The model has been extended to take into account a quantum dot material model [24,25], hybrid mode-locking [26,27] and external optical injection [28]. Further, for a quantum dot model with ground-and excited state lasing, a multi-section modification of the DDE model was proposed, which permits us to take into account more symmetric cavity designs and tapered waveguides [29][30][31] (see also [32], where a multi-section model for a linear-cavity ML laser was developed).
In this paper we extend the DDE model to study delayed optical feedback from a passive external cavity. A sketch of the model setup is shown in figure 1(a), while figure 1(b) depicts a sketch of the corresponding configuration for a linear-cavity Fabry-Perot ML two-section laser that is often considered in experimental setups. The passively ML ring laser consists of an actively pumped gain section (GAIN) and a saturable absorber section (SA). Spectral filtering is taken into account within a lumped element approach [33] by an infinitely thin Lorentzian filter (red bar in figure 1(a)). At the interface between the gain section and the absorber section, a part of the light is coupled out into an external ring cavity. The model assumes unidirectional lasing and the same group velocity of the light v in all sections. The carrier dynamics in the gain section is described by the saturable gain G(t), which is proportional to the inversion in the gain section. The carrier dynamics in the slow saturable absorber is described by the saturable loss Q(t), which is anti-proportional to the inversion in the absorber section (see the appendix for the exact definitions). The term 'saturable' is used in this context to describe that G and Q are strongly varied during one roundtrip of the pulse in the cavity due to light-matter interactions. Thus, at the arrival of the pulse the saturable gain G is depleted due to induced emission and the saturable absorption Q is decreased due to induced absorption, i.e. the absorber is bleached. , an external cavity (EC) and a Lorentzian filter (red bar). The roundtrip time of the light in the laser is T , and τ is the delay time per roundtrip in the external cavity. The strength of the feedback is denoted by K , and κ is the attenuation factor modeling accumulative nonresonant and out-coupling loss per cavity roundtrip. (b) Sketch of a linear-cavity Fabry-Perot ML two-section laser considered in experiments that we aim to describe with the DDE model.
The final set of three coupled delay differential equations for the slowly varying field amplitude E, the dimensionless saturable gain G and the dimensionless saturable loss Q reads where describes amplification and losses of the field E during one roundtrip in the ring cavity. The linewidth enhancement factors (α-factors) in the gain and absorber sections are denoted by α g and α q , respectively. Out-coupling and internal losses at the interfaces between the different sections are lumped together with the intensity transmission coefficient of the interface between the gain section and the external cavity into the attenuation factor κ (see the appendix for details). The frequency spacing of the cavity modes is 1/T , where T ≡ v/L is the cold-cavity roundtrip time of the light given by the ratio of the group velocity v and the cavity length L. The finite width of the gain spectrum is taken into account by a bandwidth-limiting element, which is described by a Lorentzian-shaped filter function f (t) ≡ γ e (−γ +i )t of full-width at half-maximum γ . It is centered at the frequency ≡ max − 0 taking into account a possible detuning of the frequency of the gain maximum max from the optical frequency of the closest cavity mode 0 . In the following, we assume zero detuning = 0. The sum in the second line in equation (1a) sums up the contribution of the field after l = 1, . . . , ∞ roundtrips in the external cavity. During each roundtrip in the external cavity, the field is delayed by the delay time τ . This is modeled by the terms R(t − T − lτ ) in equation (1a). The phase C ≡ 0 τ describes the phase shift of the light due to the finite delay time. Since the optical frequency 0 is large (THz) a tiny variation of τ causes a variation of C over its full range [0, 2π ], while the variation of τ in the delayed terms R(t − T − lτ ) and E(t − T − lτ ) has only a negligible effect on the solution for E(t). This is why we consider C and τ as independent parameters. The feedback strength is given by K l ≡ r (1 − r ec ) l/2 (1 − r ) l/2−1 , with the intensity reflection coefficients of the interface between the gain section and the external cavity r and of the interface in the external cavity r ec (see figure A.1 in the appendix). In the following, we assume the feedback to be small and neglect all but one roundtrip in the external cavity, i.e. we keep only the term with l = 1 in the sum of equation (1a) and define K ≡ K 1 . Spontaneous emission noise of strength D is modeled by a complex Gaussian white noise term ξ(t), i.e.
In the carrier equations, equations (1b) and (1c), G is pumped with the pump rate J g , which is proportional to the injection current in the gain section and J q is the unsaturated absorption parameter, which models carrier losses due to a reverse bias that is applied to the absorber section (see the appendix for the exact definition). The carrier decay rates in the gain and absorber sections are denoted by γ g and γ q , respectively. For semiconductor materials typically the inequality γ g γ q holds [34]. The last terms in equations (1b) and (1c) model the light-matter interactions. They describe how G and Q are depleted by the pulse that travels in the cavity. The factor r s is proportional to the ratio of the linear gain coefficients of the absorber and gain sections (see the appendix for the exact definition).

Dynamics of a solitary mode-locked laser
In this section the dynamics of the laser without feedback (K = 0) and without noise (D = 0) is discussed. The parameters for the simulations are given in table 1, unless stated otherwise. In a passively ML laser, the interplay of saturable gain G and saturable loss Q leads to a fixed phase relation, i.e. phase locking, of many cavity modes resulting in narrow pulses with large intensity. The number of locked modes can be roughly estimated by the ratio of the bandwidth of the spectral filter γ and the frequency spacing of the cavity modes 1/T . This means that the number of modes participating in the ML process usually increases with γ , which results in a shorter pulse width τ p ∝ γ −1 [20]. If γ is large compared to the decay times of gain and loss media γ g and γ q , respectively, the dynamics of the electrical field amplitude E is much faster than the dynamics of the carriers (equations (1b) and (1c)), due to the presence of the small factor γ −1 1 in the left-hand side of equation (1a). This means that the ML laser acts in this regime as a typical slow-fast system [35]. In the ML regime a regular pulse train is obtained. The width τ p of the pulses is short compared to the inter-spike interval time , which is close to the cold-cavity roundtrip time T . Figure 2(a) shows time traces of the absolute value of the amplitude |E| (red solid line), the saturable gain G (blue dashed line), the total loss Q t ≡ Q + |ln κ| (green dash-dotted line), which is the sum of saturable and nonresonant losses per cavity roundtrip, and the net gain G ≡ G − Q t (thin black dash-dot-dotted line). The latter is the difference of the overall gain and total loss per cavity roundtrip. It is illustrative to calculate the net gain in the simplest case of a real-valued field equation α g = α q = 0 for a continuous wave (cw) solution of the model equations. Inserting the ansatz (E(t), G(t), Q(t)) = (E s , G s , Q s ) with |E s | = 0 into equations (1a)-(1c), we see from equation (1a) that above the linear laser threshold G s ≡ G s − (Q s + ln κ) = 0 holds. Thus G = 0 corresponds to the gain-clamping above threshold observed in single-mode one-section laser models with only one carrier type [35]. A bifurcation analysis of the rotating wave solutions of equations (1a)-(1c) (cw emission) with complex E, i.e. nonzero α-factors, and K = 0 was performed in [36].
In phase space the system evolves on a stable limit cycle as one can see in figure 2(c) in a projection onto the (G, |E|)-plane. From the power spectrum in figure 2(d) we see that the pulse repetition rate ν ISI,0 ≡ 1/T ISI,0 is close to 40 GHz. Figures 2(a) and (c) illustrate that the dynamics can be divided into a slow evolution close to the slow manifold E ≡ 0 in between the pulses (slow stage) and a fast and short excursion through the phase space during the pulse (fast stage). In each of the two stages, approximate analytical solutions for the dynamical system (equations (1a)-(1c)) can be found [20]. During the slow stage G and Q recover after the passage of the previous pulse. Throughout the recovery process the total losses predominate the gain, thus G is negative. This can be seen in figure 2(a) where Q t > G holds in the slow stage. The pulses have stable background, i.e. the pulse train remains stable with respect to perturbations as long as the perturbations are too weak to open a net gain window during the slow stage. This means that New's stability criterion [37] is fulfilled [20]. The leading edge of a pulse depletes G and Q on the fast timescales of the light-matter interaction given by vg g and vg q for the gain and absorber sections, respectively, where g g is the linear gain coefficient for the gain section and g q is that for the absorber section (see the appendix for details). For r s > 1 the depletion of Q is faster than the depletion of G, which leads to a short net gain window with G > 0 triggering the pulse (gray shaded region in figure 2(b)). The condition r s > 1 means that the saturation energy of the absorber section is smaller than that of the gain section (see the appendix for details); thus the absorber bleaches first opening the net gain window. The faster decay of Q can also be seen from equations (1b) and (1c). In the fast stage the linear terms J g − γ g G and J q − γ q Q can be neglected [20]. Since r s > 1 and Q > G > 0 hold at the beginning of the fast stage, the saturation term in equation (1c) that is proportional to |E| 2 is larger than the corresponding term in equation (1b). During the net gain window the center part of the pulse is amplified.

Dynamics of a mode-locked laser subject to optical feedback
In the previous section we have seen that an important timescale of the solitary laser is the inter-spike interval time T ISI,0 . By including external feedback a new timescale is introduced, namely the delay time τ . The facet between the ring laser and the external cavity splits the pulse stream into two parts: one part propagating in the ring laser and another part propagating through the external cavity. The latter is, after one roundtrip in the delay line, coupled back into the internal cavity multiplied by the feedback strength K < 1 and a phase shift e −iC . The delay time τ and phase C determine the relative temporal position and the phase relation of the pulses coupled back from the external cavity to the pulses in the laser. We first describe the resonances of the system for the case C = 0, before discussing the influence of the feedback phase C on the dynamics. Throughout this study we concentrate on the case α g = α q = 0, which shows the clearest resonance structure. In combination with the assumption = 0 this yields a real field equation. For simplicity, we use constant initial conditions Re (E(t)) = 0.4, Im (E(t)) = 0, G(t) = 4 and Q(t) = 1 for t ∈ [−τ, 0], unless stated otherwise. Here Re and Im denote the real and imaginary parts, respectively.
Synchronization effects can be observed if an integer multiple of τ matches an integer multiple of T ISI,0 , i.e.
where the relation p < q holds in the regimes of intermediate and long delays discussed here.
The resonances described by equation (3) can be further sub-classified into the main and higher order resonances. In the main resonances ( p = 1), the pulses traveling in the external cavity match perfectly with the pulses in the laser cavity, and only one pulse per laser cavity roundtrip is present. In the higher order resonances ( p > 1), p pulses may travel in the cavity, i.e. p pulses are found in an interval of length T ISI,0 : one main pulse stemming from the solitary laser and p − 1 secondary pulses with smaller amplitudes that are induced by the feedback. Note that in [17] the main and higher order resonances are referred to as the integer and fractional resonant cases, respectively. Here the notion of the main and higher order resonances is preferred to emphasize the ordering of the corresponding fractions p/q in a Farey tree, which will be discussed in section 4.1.
The maximum number of pulses in the cavity is limited by the cold cavity roundtrip time T and the width of the pulses, which is determined by γ . The longer the cavity and the smaller the pulse width the larger the number of pulses that can fit into the cavity. For the set of parameter values given in table 1, we find up to five pulses in the cavity, i.e. p 5. At the exact resonances, the main pulses and the feedback-induced pulses are well separated by low values of the intensity. Furthermore, the pulses are almost evenly spaced over the interspike interval. That is, if we have p pulses and the main pulse is located at t = 0, the p − 1 feedback-induced pulses are located close to lT ISI,0 / p for l = 1, . . . , p − 1. . In simulations with spontaneous emission noise (D = 0), we found that in the main resonances the stability of the pulse train with respect to noise is enhanced (see also section 5). At the higher order resonances, the peak intensities of the dominant pulses decrease in the presence of delay, because the total energy of the pulse stream is now divided between the p pulses in the cavity. Figure 4(a) depicts a one-parameter bifurcation diagram of the local maxima of the time trace of |E(t)| for exactly resonant optical feedback τ = 7T ISI,0 and K varies over its full range from 0 to 1. For each K -value the dynamic equations (1a)-(1c) have been integrated over 5 × 10 4 time units to avoid transient effects and then the local maxima of the time trace of |E| were collected for 100 time units. The bifurcation diagram has been obtained by up-sweeping (orange dots) and down-sweeping (black dots) K taking the last τ -interval of the previous run as the initial conditions. (For the first runs of each sweep direction the constant initial conditions mentioned above were chosen.) At the critical feedback strength K crit = 0.55 the system becomes bistable (gray shaded region). For K K crit and up-sweeping K the system remains in the fundamental mode-locked (FML) regime with period T ISI,0 . The upper inset shows a time trace of |E(t)| in this regime for K = 0.8. However, for the constant initial conditions and the highest possible K -value of K = 1, the system exhibits harmonic modelocking (HML) with period T ISI,0 /2. For down-sweeping K , the HML branch of solutions remains stable until K = K crit , where the system abruptly jumps back to the FML branch. The lower inset in figure 4(a) depicts a time trace of |E(t)| for the HML solution and K = 0.8. This bistability could be introduced into the systems by a saddle-node bifurcation of limit cycles  (fold bifurcation) close to K = K crit , which creates a pair of new solutions: a stable node branch, which corresponds to the HML solution in the simulations, and an unstable saddle branch. Figure 4(b) depicts the dependence of K crit on the feedback phase C as observed by downsweeping K starting with constant initial conditions from K = 1. As seen from equation (1a) the system is 2π-symmetric with respect to C. The lowest value of K crit is observed for C = 0 and K crit increases strongly with C. In a symmetrical interval around C = π, the system remains in the FML regime and no bistability is observed (blue shaded region). For a shorter external delay time of τ = 3T ISI,0 (not shown), we find an additional interval symmetric around C = π, in which a secondary Andronov-Hopf (Neimark-Sacker) bifurcation introduces a new frequency into the system that is noncommensurate with 1/T ISI,0 . This leads to a quasiperiodic motion in  phase space and thus to unstable pulse trains. For higher feedback strength the FML solution is stabilized again. Thus for short, exactly resonant feedback the dynamics of the ML laser is similar to the dynamics of a single-frequency semiconductor laser subject to external optical feedback described by the Lang-Kobayashi equations [13,38,39] as already observed by Avrutin et al in [17] using a finite difference traveling wave modeling approach.
To present a more complete picture of the feedback-induced dynamics away from the exact resonances, figure 5 shows the number of pulses in the cavity (color code) as a function of the delay time τ (intermediate values) and feedback strength K for C = 0. For each pair of (τ ,K )values the system was integrated over 3 × 10 4 time units and then the number of maxima with different peak heights is counted for 25 time units.
For up-sweeping τ (figure 5(a)) for each K and the lowest value of τ the same constant initial conditions mentioned in the beginning of this section were chosen. Then, τ was increased stepwise from its lowest value taking an interval of length τ from the end of the time series of the previous run as the initial condition in each step. To catch bistability, the same procedure was performed by down-sweeping τ from its maximum to its lowest value ( figure 5(b)). For the case of intermediate delay times presented in figure 5, bistability is observed only for high K 0.5. (This becomes obvious by comparing the right boundaries of the white 'cones' with the blue regions in figures 5(a) and (b).) However, for larger delay times τ ∼ 70T ISI,0 bistability was also found for small K by varying τ (see section 4.2). The contour plot shown in figure 5 is almost periodic in τ with the period T ISI,0 corresponding to the distance between the vertical, white, dashed lines within the large cone-like blue areas around the main resonances τ = 7T ISI,0 , 8T ISI,0 and 9T ISI,0 . Within these blue areas we find complete synchronization, i.e. only one pulse travels in the cavity. In the following we call these intervals the main resonant regimes. For increasing K the resonant regimes broaden and the laser synchronizes also for larger deviations from the exact main resonances. The increase of their width is not linear in K ; instead a sudden broadening is observed for K ≈ 0.15, which we interpret as a transition from small to intermediate values of K . Away from the main resonant regimes, this transition is marked by the appearance of quasiperiodic pulse forms (white region) to be discussed below. For small K we find that the periodic orbit is deformed by the feedback in a nonlinear way (see figure 3 (middle column)). These deformations lead to τ -intervals around the exact higher order resonances (equation (3) with 2 p 5), where we find p pulses in the cavity that are well separated by low values of the intensity. They appear as red and green regions in figure 5 and in the following they are denoted as higher order resonant regimes. For very small values of K and p = 5, the smallest pulse cannot be detected numerically. This is why the brown regions indicating the presence of five pulses in the cavity do not reach down to the lowest K -value (K = 0.08 for τ/T ISI,0 = 37/5 and 38/5, as well as for τ/T ISI,0 = 42/5 and 43/5).
In between the resonant regimes the mismatch of the two timescales τ and T ISI,0 leads to overlapping pulses that are not well separated by low values of the intensity. We find broad (main) pulses that are not unimodal but have one to three side peaks. In the following, we call these τ -intervals between two adjacent resonant regimes nonresonant regimes. In the nonresonant regime in between the main resonant regime around τ/T ISI,0 = 7 and the fifthorder resonant regime around τ/T ISI,0 = 36/5, a broadened pulse with one to three sidepeaks is observed leading to the regions with two (red), three (light-green) and four maxima (darkgreen), respectively. (The same holds for τ values in between the main resonant regimes around τ/T ISI,0 = 8 and 9 and the adjacent fifth-order resonant regimes, respectively.) Also between adjacent higher order resonances, nonresonant regimes are observed, which are for example responsible for the light-green ( p = 3) and dark-green ( p = 4) colored regions between the second-order resonance regimes around τ/T ISI,0 = 15/2 and 17/2 (red) and the adjacent fifthorder regimes (brown), respectively. For K = 0.12, figure 6 shows time traces of |E|, G, Q t and G in the nonresonant regimes to the left (see figure 6(a)) and to the right (see figure 6(b)) of the second-order resonant regime around τ/T ISI,0 = 15/2, respectively. We observe broadened main pulses with shoulders at the leading (see figure 6(a)) and the trailing edge (see figure 6(b)) of the main pulse, respectively. The pulses in the nonresonant regimes have relatively low peak intensities compared to the pulses in the resonant regimes (compare figures 6 and 3 (left panel)). In addition the net gain G (dash-dot-dotted line) shows a broad shoulder with only slightly negative values between the small pulse and the leading edge of the next larger pulse. This is why the pulses in the nonresonant regimes are more easily destabilized by noise than the pulses in the resonant regimes (see section 5).
The resonance structure observed for small feedback strengths K 0.15 is robust under changes of the feedback phase C. However, in the regions of the (τ ,K )-plane where bistability is observed (compare figures 5(a) and (b)) the trajectory can switch from one attractor to the other if C is changed. As a consequence, the regions in the (τ ,K )-plane, where quasiperiodic motion is observed, depend on C.
Examples of the quasiperiodic behavior that is observed away from the main resonant regimes for intermediate and high values of K (white regions in figure 5) are depicted in figure 7 for K = 0.25. Figures 7(a) and (b) depict the dynamics at the second-and third-order resonance with τ/T ISI,0 = 15/2 and 22/3, respectively. We find two and three competing pulse trains, respectively, leading to a periodic modulation of the peak amplitude of the pulses. Like in the periodic regime for small K , at the exact resonances the inter-pulse temporal distance is given by approximately T ISI,0 /2 and T ISI,0 /3 resulting in dominant frequencies of 2ν ISI,0 and 3ν ISI,0 , respectively (see figure 7 (right panel)). The lowest frequency in the power spectra is the competition frequency ν c of neighboring pulse trains, which is small (∼ 1-3GHz) compared to the pulse repetition frequency ν ISI,0 ≡ 1/T ISI,0 . Thus, we plot only the maxima of the competing pulse trains in figures 7(a) and (b) (left column). For the second-order resonance (see figure 7(a)) two competing pulse trains appear that are labeled A and B. Their competition frequency ν c = 2.38 GHz corresponds to the temporal distance between the largest maximum of pulse train A and the largest maximum of pulse train B. For the third-order resonance (figure 7(b) (left)), we find three competing pulse trains A, B and C with a competition frequency of ν c = 2.55 GHz. The two (three) competing pulses can also be seen in the phase space projection onto the (G, |E|)-plane in figure 7(a) (middle) that shows projections of the trajectories for one competition period T c ≡ 1/ν c . The black dashed parts of the trajectory are the projection of one T ISI,0 -interval of the limit cycle. For high K > 0.73 we find HML in the second-order resonance, which results in a pulse train with pulse repetition rate 2ν ISI,0 .
In figure 8 the transition to quasiperiodic motion is shown in a one-parameter bifurcation diagram depicting the local maxima of the time trace of |E| for fixed τ/T ISI,0 = 15/2 and K varying from 0 to 1. The diagram has been obtained by up-sweeping (orange diamonds) and down-sweeping (black dots) K as described for figure 4(a). For small K we retrieve the periodic pulse train with p = 2 pulses in the cavity (see figure 3(b) and left inset of figure 8). At a feedback strength of K = 0.19 a new noncommensurate frequency is introduced by a secondary Andronov-Hopf (Neimark-Sacker) bifurcation leading to a quasiperiodic motion on a torus in phase space (compare with figure 7(a)). Eventually, at K = 0.865 a transition to a harmonic ML state with period T ISI,0 /2 takes place. In the right blowup of figure 8, a time trace of the amplitude in the harmonic ML state is shown for K = 0.9. For decreasing K this state persists down to K = 0.73.
For nonzero α-factors and α g < α q the resonance regimes discussed in this section subsist, but in certain main resonances we find bifurcation cascades leading to chaotic pulse trains as already mentioned by Avrutin and Russell [17]. For sufficiently large α g > α q , the laser without feedback already exhibits unstable pulse trains. A detailed discussion of the additional effects introduced by finite α-factors is beyond the scope of this paper.

Farey tree
The pulse trains found in the ML laser subject to feedback can be further classified by comparing them with the locking dynamics of externally driven dynamical systems. A generic example of such a system is the circle map with an external forcing [40], where the forcing introduces an additional frequency into the system and leads to dynamic evolution on a torus in phase space. We now briefly recall the concept of Farey tree ordering in these systems before we show that this concept can also be applied to the ML laser with delay. In general, the frequency of the external signal and the resulting resonance frequency of the system itself are noncommensurate, which results in a quasiperiodic motion on the surface of the torus. Locking takes place if the ratio of the frequency of the external signal and the resonance frequency of the system is a rational number p/q with p, q ∈ N. The fractions p/q are called winding numbers and can be ordered in a Farey tree by applying the Farey-sum operation [41]. This means that the denominators and the numerators of two neighboring fractions p/q and p /q are added up separately yielding a new fraction ( p + p )/(q + q ) in the next higher level of the tree (see figure 9(a)). With increasing injection strength the intervals around the exact p/q-resonances, at which locking takes place, increase, forming the well-known Arnold tongues. For fixed injection strength, one finds a hierarchy of the sizes of the locking plateaus according to their ordering in the Farey tree, i.e. the plateau sizes decrease from lower to upper levels, resulting in a self-similar pattern. This generic type of frequency locking has been found in a wealth of nonlinear dynamical systems, for instance, modulated external cavity semiconductor lasers [42,43].
In the case of an ML laser the τ -intervals of the resonant regimes form plateaus with p well-separated pulses in the cavity (compare with figures 3(b)-(d) (left column)), where p defines the level of the plateau in the Farey tree. In the numerical simulations we consider two subsequent pulses to be well separated if between the pulses a time point t min exists with |E(t min )| < 1 × 10 −3 . Figure 9 . It can be seen that plateaus similar to the locking plateaus in the Farey tree form. The width of the plateaus shrinks with increasing the number of pulses p, i.e, with increasing the tree level, as one would expect from the study of the generic circle map. A similar scenario has also been observed for a pair of delay-coupled oscillators [44].
For a higher feedback strength of K = 0.25 the width of the plateaus increases, which is similar to the broadening of the Arnold tongues for increasing control parameters in the case of externally driven systems. As already discussed in figures 5 and 7 quasiperiodic motion is observed away from the main resonances for the higher feedback strength (K = 0.25). Nevertheless, we still find p = 2 and 3 well-separated-pulses in the cavity in the secondand third-order resonant regimes, respectively. They are indicated by open black diamonds in figure 9(b).

Long delay
In this section we study a feedback loop that is ten times longer than that discussed in the previous sections. This regime of long delays is of experimental interest, because in experimental setups optical feedback is often provided by a long optical fiber loop [6,12]. Figure 10 figure 5). In the gray shaded parts of the main resonant regimes in figure 10(a), the peak values of the time trace of |E| are higher and the pulse width is shorter than for the solitary laser. Simulations with white noise (D = 0) reveal that in these regions the stability of the pulses with respect to fluctuations is increased. In contrast to the case of intermediate delay times discussed in the previous subsection, bistability of periodic orbits is now detected also for small feedback strengths.
In addition to bistability between FML and HML solutions (compare with figure 8) bistability is now also observed between solutions with one pulse in the cavity and solutions with two, three or four pulses in the cavity for τ -values ranging from 67.25T ISI,0 to 67.47T ISI,0 (shaded blue in figure 10(a)). In figure 10(b) exemplary time traces of |E| are plotted for two periodic orbits that are bistable for τ = 67.45T ISI,0 (vertical gray dashed line in figure 10(a)). By up-sweeping τ , a pulse train with only one pulse traveling in the cavity ( p = 1) is observed (orange dashed line), while by down-sweeping τ a pulse train with two pulses in the cavity ( p = 2) is detected (black solid line). Figure 10(c) depicts phase space projections of the trajectories onto the (G, |E|)-plane.

Delay-induced frequency pulling
In general, the inter-spike interval time T ISI differs for finite k and τ from the inter-spike interval time T ISI,0 of the laser without feedback, except for τ -values fulfilling exactly the main resonance condition ( p = 1 in equation (3)). This dependence of T ISI on τ results in a delayinduced variation of the pulse repetition frequency ν ISI (τ, K ) ≡ 1/T ISI (τ, K ), which is known as frequency pulling and has been observed experimentally in quantum well lasers [6] and more recently in quantum dot lasers [7-10, 12, 45]. Figure 11 depicts the deviations ν ISI (τ, K ) ≡ ν ISI (τ, K ) − ν ISI,0 of ν ISI from the pulse repetition frequency in a laser without feedback ν ISI,0 ≡ 1/T ISI,0 as a function of the delay time τ for small feedback strength, K = 0.03 (red crosses), intermediate, K = 0.25 (white diamonds), and strong, K = 0.5 (orange circles), respectively. We only discuss the frequency pulling in the main resonant regime where one pulse travels in the cavity (blue regions in figure 5 and the two largest plateaus in figure 9(b)).
In the main resonant regimes, we have synchronization if the inter-spike interval time T ISI adapts to the external delay time τ according to qT ISI (q, τ ) = τ leading to a frequency  figure 11(a) that the main resonant regimes broaden and the frequency deviations ν ISI increase with K (see also figure 5). For small values of K the frequency detuning is very small and a nearly horizontal line is found (red) that reaches zero detuning (horizontal gray dashed line) when τ approaches the borders of the main resonant regime. Instead, for intermediate and large K (white diamonds and orange circles, respectively) ν ISI increases nearly linearly with the distance from the exact resonance, but remains below the detuning expected for synchronization (dash-dotted black lines in figure 11). From figure 11(b) it can be seen that the width of the main resonant regimes is increased for longer delay (note the different scales of the τ -axis in figures 11(a) and (b)). For small values of K (red crosses) we find a similar behavior of ν ISI as for intermediate delay. An increase of K (white diamonds), however, leads to a frequency entrainment that stays close to its value for synchronization (see equation (4)) over nearly the whole interval length of T ISI,0 (compare the slope of the white diamonds and the black dash-dotted line). Further, jumps of ν ISI take place close to the exact resonance conditions. In figure 11(b) such a jump is observed at τ = 67.85T ISI,0 when the system switches from the main resonant regime with q = 67 to the main resonant regime with q = 68 at τ = 67.9T ISI,0 . In between two subsequent main resonant regimes a short interval with more complex dynamics is found. The sawtooth-like shape of ν ISI in terms of τ corresponds well to the experimental findings for long delay lines [6,9,12]. For large values of K (orange circles) we find a different behavior. Now the system remains in one main resonant regime for τ -intervals larger than T ISI,0 . In figure 11(b) the laser stays in the main resonance with q = 67 from τ = 67T ISI,0 up to τ = 68.3T ISI,0 . Figure 11 has been obtained by stepwise increasing the τ -values (up-sweep). This explains the asymmetry of ν ISI with respect to zero detuning (horizontal gray dashed line) for long delay ( figure 11(b)). For down-sweeping τ we obtain qualitatively the same resonance behavior, but in this case the larger absolute values of ν ISI are observed for positive detunings. This means that for intermediate and large K -values, we have bistability between periodic orbits from the neighboring main resonant regimes (here q = 67 and 68). These orbits have the simple shape of the one depicted in figure 3(a) (left) (only one pulse in the cavity) and a frequency difference of ν = 1/τ , which is the vertical distance of neighboring black dash-dotted lines in figure 11(b). When the length of the delay line is increased even further, these bistabilities become more pronounced and appear at lower values of the feedback strength.

Reduction of timing jitter by optical feedback
In this section we study the parameter dependence of the timing jitter for intermediate delay times near a main resonance. The dynamics of the ML laser with delay given by equations (1a)-(1c) is simulated with a noise strength of D = 0.05. To characterize the timing jitter we use the procedure proposed in [16] to determine the root-mean-square (rms) timing jitter. The specific value of the noise strength has been chosen to obtain values of the timing jitter that are comparable with those stated in [16]. However, simulations with larger noise strengths reveal that the reduction of the jitter in the main resonant regimes is robust against changes of the noise strength.
The This means that the timing fluctuations can be modeled by a nonstationary random process. The reason is that in contrast with active and hybrid mode-locking there is no external reference clock in the modeled passively ML laser. Thus the white noise applied to the field equation yields in good approximation a discrete Lévy-Wiener process, i.e. a random walk, of the timing fluctuations [16].
For the simulations pulse trains containing N = 3 × 10 4 pulses are used. By applying a discrete Fourier transform to each realization of the timing fluctuations { t} j and subsequently averaging over M noise realization, we obtain the timing spectral density where we have introduced the time span T ≡ N T ISI . Now a phase that expresses the position of the pulses in the pulse train with respect to the clock can be introduced in such a way that a temporal deviation of the pulse center from the ideal clock by T ISI results in a phase shift of 2π [46]. The timing spectral density is proportional to the phase noise spectrum The rms timing jitter is obtained by integrating the phase noise spectrum L(ν) over the frequency range from ν low to ν high where ν low and ν high are the minimal and the maximal frequency offset from the repetition frequency of the ideal clock 1/ T ISI (carrier frequency). The factor of 2 in equation (6)  from ν low = 1.5 MHz to ν high = 5 GHz. The results for the rms timing jitter σ rms are plotted as a color code in figure 13(b) as a function of the delay time τ and feedback strength K for τ -values in an interval of ±0.15T ISI,0 around the main resonance with τ = 7T ISI,0 . Figure 13(c) shows the timing jitter normalized to its value at K = 0. Here black corresponds to an increase and lighter colors to a decrease of σ rms with respect to σ rms at K = 0, thus underlining the regions where a jitter reduction by the external feedback is possible. Figure 13(a) shows several one-dimensional sections of figure 13(b) for different fixed values of the feedback strength K . The results of section 4 suggest that the stability of the pulses with respect to noise is increased in the main resonances. Indeed, in the density plot of figure 13(b) one sees that the lowest values of the timing jitter σ rms < 100 fs are found in a cone around the exact resonant condition (black dashed line). Thus, the cone-shaped main resonant regime (blue area in figure 5) is qualitatively reproduced. For τ -values away from the exact resonance and intermediate feedback strengths 0.08 < K < 0.16, we find a dramatic increase of the jitter (see dark-red and black regions in figure 13 for τ > 7.1T ISI,0 and for τ < 6.9T ISI,0 , respectively). Comparing this result with those of figure 5 it becomes clear that in this regime the transition to quasiperiodic motion takes place, leading to less stable dynamics. For higher K the main resonant regimes are already relatively broad and the jitter decreases nearly over the full length of the studied τ -interval.
For longer delay times the τ -intervals in which the jitter is reduced broaden. This makes it easier to access these regions experimentally. Our findings are in good agreement with recent experimental results [12] showing that a reduction of the timing jitter is possible for a large range of τ -values. In our terminology this corresponds to broad main resonant regimes, which we found for long delay (see figure 11(b) and K = 0.25 (white diamonds)). In between these intervals a sudden increase of the timing jitter above the value of the solitary laser was observed experimentally [12]. This is in agreement with the results of our numerical simulations indicating the existence of complex dynamics between two neighboring main resonant regimes.

Conclusion
In conclusion, we have extended the DDE model for a passively ML semiconductor laser to study the influence of external optical feedback on the laser dynamics. We find resonances between the inter-spike interval time T ISI,0 and the delay time τ that follow the ordering of the Farey sequence. If τ is close to an integer multiple of T ISI,0 only one pulse travels in the cavity (main resonant regime), while p pulses are found (higher order resonant regimes) if the ratio of τ and T ISI,0 is close to a rational number ( pτ ≈ qT ISI,0 ). In between the resonant regimes the mismatch of τ and T ISI,0 leads to pulse trains of broad and less stable pulses (nonresonant regime).
For small feedback strengths, the periodic orbit of the pulse train is deformed in a nonlinear way, but remains stable. For intermediate feedback strengths the main resonant regimes broaden, while quasiperiodic motion is observed outside of this regime. For high feedback strengths HML is found in the second-order resonant regimes.
By increasing the delay time we find, on the one hand, a broadening of the width of the main resonant regimes and, on the other hand, for long delay, bistability between pulse trains with different numbers of pulses in the cavity is observed.
Further, we have shown that the pulse repetition frequency is entrained by the external feedback in the main resonant regimes. For long delay and intermediate feedback strengths, a sawtooth-like shape of the frequency detuning as a function of the delay time τ is found as also observed in experiments. Moreover, we find a reduction of the rms timing jitter in the main resonant regimes, which is in good agreement with experimental results.
Sketch of the ring cavity laser with external optical feedback as used for the derivation of the DDE model. The laser consists of five sections: a saturable absorber section (SA), a gain section (GAIN), an external cavity (EC), a passive section (P) and a Lorentzian filter (red bar). The field amplitudes in the saturable absorber, in the gain section and in the passive section, are denoted by E q , E g and E p , respectively. The facets between the sections are located at z i for i = 1, . . . , 4. κ 1 and κ 2 are attenuation factors modeling nonresonant and out-coupling losses located at z 1 and z 2 , respectively. t, r and t ec , r ec are intensity transmission and reflection coefficients at the interface between the gain section and the external cavity and at the interface in the external cavity, respectively.
Rewriting equations (A.1a) and (A.1b) in the coordinates (t, ζ ), we obtain the following set of ordinary differential equations: ∂ ζ A r (t, ζ ) = 1 2 (1 − iα r )N r (t, ζ )A r (t, ζ ), (A.8a) ∂ t N g (t, ζ ) = J g (t) − γ g N g (t, ζ ) − N g (t, ζ )|A g (t, ζ )| 2 , (A.8b) ∂ t N q (t, ζ ) = −J q (t) − γ q N q (t, ζ ) −r s N q (t, ζ )|A q (t, ζ )| 2 , (A.8c) where we have introduced the rescaled amplitudes of the electric field A r ≡ vg g g E r , rescaled carrier densities N r ≡ vg r r [n r − n tr r ] that are proportional to the carrier inversion in the corresponding section, rescaled pump rates J g ≡ vg g g ( j g − γ g n tr g ) and J q ≡ vg q γ q q n tr q . Further,r s ≡ g q q /g g g = E sat,g /E sat,q , where E sat,g ≡ (v g g g ) −1 and E sat,q ≡ (v q g q ) −1 define the saturation energies of the gain and absorber sections, respectively, 0 is the vacuum permittivity and bg is the background permittivity. (The saturation energy is the energy of a short pulse that decreases the inversion to 1/e of its initial value.) The evolution of the field in the gain and absorber sections is obtained by integrating equation (A.8a) over these sections where dimensionless carrier densities To derive differential equations for G and Q equations (A.8b) and (A.8c) are integrated with respect to ζ over the gain and absorber sections, respectively. The last terms on the right-hand sides of the resulting equations can be rewritten in terms of the field intensities at the facets ζ 2 ,ζ 3 ζ 1 ,ζ 2 N q,g (t, ζ )|A q,g (t, ζ )| 2 dζ = |A q,g (t, ζ 2,3 )| 2 − |A q,g (t, ζ 1,2 )| 2 .  where we have introduced the averaged pump parameters and A(t) ≡ A q (t, ζ 1 ).
Next an integral equation for the evolution of the field amplitude during one round-trip is derived. Rewriting equation (A.5) in the time domain, we obtain from the right-hand side a convolution product ( ) of the filtering function in the time domain f and E p E q (t , z 1 + L) = √ κ 1 ( f E p )(t , z 4 ). By transforming equation (A. 19) to a co-moving frame A(t) = 1 √ κ 2 E(t) e i t , the final set of DDEs is obtained: where the rescaled ratio of saturation intensities in the gain and absorber sections r s ≡r s /κ 2 has been introduced [21]. In all simulations κ 2 is set to unity.