Bifurcation analysis of the Yamada model for a pulsing semiconductor laser with saturable absorber and delayed optical feedback
Bifurcation analysis of the Yamada model for a pulsing semiconductor laser with saturable absorber and delayed optical feedback
Soizic Terrien
Abstract
Semiconductor lasers exhibit a wealth of dynamics, from emission of a constant beam of light, to periodic oscillations and excitability. Selfpulsing regimes, where the laser periodically releases a short pulse of light, are particularly interesting for many applications, from material science to telecommunications. Selfpulsing regimes need to produce pulses very regularly and, as such, they are also known to be particularly sensitive to perturbations, such as noise or light injection.
We investigate the effect of delayed optical feedback on the dynamics of a selfpulsing semiconductor laser with saturable absorber (SLSA). More precisely, we consider the Yamada model with delay – a system of three delaydifferential equations (DDEs) for two slow and one fast variable – which has been shown to reproduce accurately selfpulsing features as observed in SLSA experimentally. This model is also of broader interest because it is quite closely related to mathematical models of other selfpulsing systems, such as excitable spiking neurons.
We perform a numerical bifurcation analysis of the Yamada model with delay, where we consider both the feedback delay, the feedback strength and the strength of pumping as bifurcation parameters. We find a rapidly increasing complexity of the system dynamics when the feedback delay is increased from zero. In particular, there are new feedbackinduced dynamics: stable quasiperiodic oscillations on tori, as well as a large degree of multistability, with up to five pulselike stable periodic solutions with different amplitudes and repetition rates. An attractor map in the plane of perturbations on the gain and intensity reveals a Cantor setlike, intermingled structure of the different basins of attraction. This suggests that, in practice, the multistable laser is extremely sensitive to small perturbations.
Key words
Semiconductor lasers with saturable absorber, Delaydifferential equations, Bifurcation analysis, Multistability.
1 Introduction
Semiconductor lasers are very efficient sources of light, widely used for many different applications, such as telecommunications [37], cryptography [33], optical data storage [36] and eye surgery [55]. This popularity is easily explained by the features of semiconductor lasers: they are small, cheap, and efficient [2]. Semiconductor lasers have thus attracted considerable attention in the last decades, and their dynamics have been studied widely, both from theoretical and experimental points of view; see for example [20, 24, 33] as entry points to the extensive literature on laser dynamics. Apart from stable emission of a continuous beam of light, semiconductor lasers have been shown to produce a wealth of other dynamics. This includes excitability [7, 13, 51, 27], which corresponds to a all or none, wellcalibrated response to a perturbation, periodic oscillations [1, 6, 11], quasiperiodic modulation of the light amplitude [39, 57], as well as chaotic dynamics [35, 48]. Pulselike periodic solutions, where the laser periodically produces a very short, highamplitude pulse of light, are particularly interesting from an application point of view. However, trains of pulses may be very sensitive to perturbations, such as noise [12], optical feedback [29], or injection of the light produced by another laser [56]. A small amount of internal or external noise can induce timingjitter of the pulses [17, 38, 52]: pulses are then produced with a nonconstant repetition rate, which is detrimental to many applications.
We focus here on a micropillar semiconductor laser with saturable absorber (SLSA) subject to delayed optical feedback. Figure 1 is a sketch of the device: the laser consists of an optical resonant cavity of a few microns bounded by two parallel, highreflectivity mirror stacks. The top facet where the light is emitted is of such small diameter (3 ) that lasing occurs only in the fundamental mode. The microcavity embeds a gain medium, consisting of two quantum wells [8]. The energy provided to the system through optical pumping is stored in that gain section, which can then act as a light amplifier. The key feature of the device is that it includes a saturable absorber section in the cavity, consisting of a single quantum well [8]. In such an absorber, the optical transmission increases when the light intensity is larger than a given threshold. Without feedback, this device has been shown, both experimentally and theoretically, to exhibit excitability [4, 7, 45] and selfpulsations [8, 11] arising through an homoclinic bifurcation [6].
Different physical mechanisms can induce selfpulsations in lasers. Modelocking relies on a phase synchronisation of the different longitudinal modes of the laser cavity, leading to the periodic emission of extremely short pulses of light with high repetition rate [42]. On the other hand, a laser as discussed here is essentially a single longitudinal mode system, so that modelocking cannot take place. For this type of device, selfpulsations occur through a passive Qswitching mechanism [11], which results from variable losses in the cavity and leads to the emission of a train of highintensity light pulses. In the case of passive Qswitching, variable losses are induced by the saturable absorber. The physical mechanism can be explained schematically as follows. At the first stage, the intracavity intensity is low, and the losses due to the saturable absorber are consequently high. The energy provided by pumping the gain medium is then almost entirely stored in the gain section. When the gain becomes sufficiently high to overcome the losses, the intracavity intensity starts to grow, and eventually saturates the absorber, which results in a sudden drop of the cavity losses. As the energy stored in the gain section is high, the intensity can grow quickly, until all the stored energy is released as a pulse of light. In the process, the intracavity intensity drops back to a low value, and the same process can repeat. However, the next pulse can only be triggered after a given amount of time, called the refractory period, during which the gain section recovers; see, for example, [45]. Because the electric field in the cavity evolves on a shorter timescale than the carrier population in gain and absorber media, one obtains a train of short light pulses separated by long periods with almostzero intensity.
The Yamada model [58], a system of three ordinary differential equations (ODEs) for the gain , the absorption and the intensity , is a well known model for Qswitched singlemode lasers. This model has been extensively studied; in particular, a complete bifurcation analysis has been performed in [6], where all the possible different dynamics have been determined. It has been shown recently to describe accurately, both qualitatively and quantitatively, a range of dynamics of the micropillar laser sketched in figure 1 in the case without feedback: this includes selfpulsations, excitability and temporal summation [4, 45, 46].
We consider in this study the effect of delayed optical feedback on the SLSA dynamics. Since feedback can either occur naturally in laser systems through unintentional reflections, or be introduced in an attempt of control in many physical systems [41], the effects of feedback on a laser’s dynamics have attracted considerable attention in the last decades; see, for example, [20] as an entry point to the extensive literature. Different kinds of feedback can be considered, including filtered or phaseconjugate optical feedback [23]. However, we consider here the simplest case of conventional optical feedback from a regular mirror, where a fraction of the output light is reinjected into the laser cavity after a given delay ; see the sketch in Figure 1. More precisely, a beam splitter separates the output light into two beams: the first one constitutes the output of the SLSA, and the other part is reinjected into the cavity, after a delay time that is directly related to the length of the feedback loop. An optical isolator prevents unwanted backward reflections. In Figure 1, the feedback loop is sketched with an optical fiber for sake of clarity; however, the delay loop is realised in the experiment through propagation in free space [32]. Our focus is the effect of delayed optical feedback on the excitable regime, which has been shown to occur, in a SLSA, just below the laser threshold at which the laser starts to emit light and then produces selfpulsations [6]. The motivation for adding a feedback loop into the excitable SLSA relies on the following simple idea: when an excitable pulse in reinjected after a delay , it will trigger a new pulse if the reinjection strength is large enough. As this process repeats, a regular pulse train will be produced, where the repetition rate is directly related to the delay time , and the pulse amplitude and shape are wellcalibrated. By adding a delayed optical feedback loop, we aim here at achieving a better control of the selfpulsing solutions: this includes the reduction of timingjitter, as well as the control of amplitude and repetition rate.
We consider here the Yamada model with an additional delayed term to describe the delayed optical feedback, as already introduced in [29]. In dimensionless form, it can be written as the system of three delaydifferential equations for the gain , the absorption and the intensity :
(1) 
Here, is the pump parameter of the gain, describes the nonsaturable losses, is the saturation parameter, and is the recombination rate of carriers in the gain and absorber medium (which are assumed here to be the same). In the feedback term of the intensity equation, is the feedback strength, and is the feedback delay. For physical reasons, all these parameters are positive. The time variables are rescaled to the cavity photon life time, which has been estimated to be of the order of a few picosecond for the device we consider here [4].
Many models for semiconductor lasers with delay include the complex electric field and not only the intensity as discussed here. One speaks of LangKobayashitype rate equation models [17, 31, 38]. We consider here a situation in which the SLSA is either off, in an excitable regime, or produces pulses. In between such very short light pulses, the electric field inside the laser is thus almost zero for long periods of time, its only contribution being the spontaneous emission noise (which is not considered in (1)). Such a configuration is conceptually very different from optical feedback introduced into a laser that operates in or near a continuouswave regime. More precisely, the response of the excitable SLSA to a sufficiently large external perturbation is a highamplitude pulse of approximately 200 duration [45]. When an additional feedback loop is considered, the feedback delay is in practice considerably larger than the pulse duration (typically between 4 and 60 ). When the excitable pulse is reinjected, the intracavity field has thus come back to an almostzero value and has low coherence. Consequently, its phase is actually not properly defined. Indeed, it has been shown numerically that the feedback phase is not relevant in the limit of low to medium pumping [52]. In this specific configuration, one can thus reasonably consider that the electric field at the instant does not interfere with the reinjected, delayed field . Conceptually, this configuration is very similar to the case, without feedback, when the excitable SLSA is subject to an external pulse perturbation. It has been shown experimentally that then the system dynamics only depends on the amplitude of the perturbation: the response of the device to a perturbation is qualitatively and quantitatively the same whether this perturbation is to the gain or to the intensity [46]. In particular, the feedback is incoherent and linear in the intensity [40, 50]. The Yamada model with an additional feedback term in the intensity equation thus emerges as a natural first minimal model for the device considered here. Because it is mathematically quite simple, it is possible to study analytically the steadystate solutions [29], as well as to perform an indepth numerical bifurcation analysis, highlighting a wide range of the different possible, complex dynamics.
Although they are derived from the same physical principles, the Yamada model with feedback and LangKobayashitype models for semiconductor lasers subject to feedback rely on different assumptions related to the specifics of the considered device; as such they considerably differ from each other [31, 38]. In fact, equations (1) can be related directly to the LangKobayashi equations: considering only the intensity dynamics and taking into account the assumption of noninterference between the delayed and instantaneous fields results in a feedback term that is formally different but gives the same results. A more detailed comparison between these Yamadatype and LangKobayashitype models for pulsing lasers with (self)feedback is beyond the scope of this article and will be discussed elsewhere. In particular, we remark that any of these laser models with feedback may show the phenomenon of practically independent pulses in the external feedback loop, providing the delay is sufficiently large (larger than considered here). Such solutions are also referred to in the literature as temporal cavity solitons or localised structure; see for example [34, 43].
The Yamada model should be of more general interest, beyond the specific device sketched in figure 1. In fact, it can be seen as an optical analogue to excitable spiking neurons: such systems exhibits related behaviours, including excitability, and are modelled accurately with very similar equations that describe the dynamics of a fast voltage, see for example [16]. In particular, these models do not have a phase variable either. More generally, the Yamada model can be seen as a paradigmatic model for coupled pulsing systems, where the configuration considered here corresponds to the particular case of selfcoupling.
In this paper, we adopt a dynamical system point of view and perform a bifurcation analysis of the Yamada model with delayed optical feedback. By considering both the feedback strength and the feedback delay as bifurcation parameters, we aim at studying extensively the effect of the delayed feedback on the laser’s dynamics. In a first bifurcation analysis of this model [29], bifurcation diagrams were computed in the plane, for two different values of the pump parameter , corresponding to two qualitatively different behaviour of the system without feedback. In the first case, the laser without feedback is in the excitable regime, while it is in the nonexcitable offstate in the second case [6]. The study in [29] has shown interesting new feedbackinduced dynamics, including bistability of two different pulsing solutions. However, it focused on very low values of the feedback delay of . We extend here the bifurcation analysis for larger values of the feedback parameters. This is physically relevant, as such large values of the delay are usually considered experimentally [32]. Moreover, after considering the bifurcation diagrams in the plane, we also consider bifurcation diagrams in the plane of pump current and feedback strength. Again, this is of particular interest from an experimental point of view: the value of the delay being directly related to the length of feedback loop (see figure 1), it is fixed during a given experiment. Conversely, the feedback strength can be tuned quite easily, and the pump current is the main control parameter.
Compared to ordinary differential equations (ODEs), solving delaydifferential equations (DDEs) requires to specify as an initial condition not only the state at time , but its values on a whole time interval . Consequently, a DDE has an infinitedimensional phase space (see, for example, the review [44]), and one has to use specific numerical methods for continuation and bifurcation analysis. We use here the numerical continuation software DDEBiftool [9, 10, 49], which has been used to investigate many different delay systems [14, 21, 54]. As explained above, the intensity , on the one hand, and the gain and absorption , on the other hand, evolve on different timescales. In equations (1), this is reflected by a small value of parameter , which is characteristics of a slowfast system; see for example [19]. The analysis performed in this paper is thus also an example of a bifurcation analysis of a slowfast system with a single delay, which can be of broader interest than only for laser applications.
New important information in this study includes the precise loci of the bifurcations of periodic solutions; such bifurcations can now be continued numerically for DDEs [49], thus providing more complete bifurcation diagrams. For both the case of the SLSA in an excitable regime () and in a nonexcitable offstate (), we perform a bifurcation analysis in the plane. Both cases lead to qualitatively different bifurcation diagrams. Both bifurcation diagrams have in common a selfrepeating structure, meaning that they become more and more complex as the delay is increased, which is characteristic feature of delay systems [59]. In particular, new feedbackinduced dynamics is observed in both configurations for sufficiently large values of . This includes multistability between up to five pulselike periodic solutions, with different pulse amplitudes and repetition rates, as well as stable quasiperiodic oscillations on tori. In the case of multistability, we show that the periods of the different periodic solutions (i.e. the repetition rate of the pulses) are directly related to the value of the delay time . Considering three fixed values of , we then perform a bifurcation analysis in the plane of the pump parameter and feedback strength. It shows that the increasing complexity and wealth of dynamics observed in the plane can also be accessed in the plane. Finally, we investigate the consequences of feedbackinduced multistability. In particular, we show that the basins of attraction associated with the different stable pulsing solutions have an intricate, Cantor setlike structure in the plane of perturbations on the gain and the intensity . This highlights the extreme sensitivity of the multistable laser to small perturbations.
The paper is organized as follows. The main results of the first bifurcation study performed in [29] are briefly summarized in section 2 as a starting point. In section 3, an extended bifurcation analysis of the model in performed in the plane of feedback delay and feedback strength, for two different values of the pump parameter corresponding to two qualitatively different configurations. Bifurcation diagrams in the plane are studied in section 4, for three fixed values of . Finally, the consequences of multistability are investigated in section 5. We summarize and point out directions for future work in section 6.
2 Background on the Yamada model with delayed feedback
The bifurcation analysis of the Yamada model in [6] and [29], demonstrated, for , the existence of an equilibrium:
(2) 
Regardless of the other parameters values, this equilibrium has zerointensity, and thus corresponds physically to the nonlasing or off solution. In presence of feedback, two other equilibria and appear for a specific value of the feedback strength through a saddlenode bifurcation [29]. Its locus is given by:
(3) 
If , then is the only equilibrium of system (1). Conversely, if , both , and are solutions of the system. is a saddle, is a node, and they are given by:
(4) 
(5) 
where
(6) 
The equilibria and display constant, nonzero intensity, and thus correspond physically to continuouswave emission. As demonstrated in [29], undergoes a transcritical bifurcation when the feedback strength is further increased to the value :
(7) 
At this point, exchanges stability with . Meanwhile, the value of given by (6) becomes negative: enters the halfplane with , and this solution is thus mathematically stable, but not physically relevant anymore. More precisely, as long as , the equilibrium is attracting and is of saddletype, with positive intensity. Conversely, is of saddle type for . This value , which directly depends on the value of the pump parameter , thus corresponds to the laser’s threshold. Other bifurcations of equilibria, as well as bifurcations of periodic solutions, are found in system (1). However, they cannot be calculated analytically, and one has to use dedicated continuation software for DDEs [9, 10, 49].
3 Bifurcation study in the ()plane
We now perform a bifurcation study of equations (1) in the plane of feedback parameters, for two different values of the pump parameter , where for the laser is excitable or has a unique nonlasing equilibrium, respectively. To investigate the influence of the optical feedback, we consider the fixed values , and , as in [29].
3.1 Case of the SLSA in the excitable regime
We first consider the case of a fixed pump parameter . In this configuration, the laser without feedback is excitable [6], which corresponds to the phase portrait 2 in figure 3. Figure 2(b) summarizes the main results of the bifurcation analysis performed in [29] for this configuration. The whole bifurcation diagram has been recomputed and includes new information, with extra curves representing the bifurcations of periodic solutions.
Each curve of this twoparameters bifurcation diagram corresponds to the locus of a specific bifurcation, and the different curves thus divide the ()plane into regions with qualitatively different dynamics. An overview of these different dynamics is provided in figure 3, which represents twodimensional projections, onto the plane, of the phase portraits associated with the different regions. Such a projection is relevant because, for small values of the delay, the dynamics effectively takes place in a twodimensional attracting surface [29]. It is worth noting that the numbering of the different regions and phase portraits is consistent with that in [6] and [29], the only phase portrait not reproduced here is that in region 1, where the equilibrium is a global attractor (which is not actually found in the bifurcation diagram in figure 2).
Overall, we find the following bifurcation curves and points of system (1) (see, for example, [30] for definitions and properties of these bifurcations):

: a curve of transcritical bifurcation, located at , where the equilibria and exchange stability. This corresponds, for example, to the transition between phase portraits 6 and 7 in figure 3.

: a curve of saddlenode bifurcation, located at , where the equilibria and bifurcate. This corresponds, for example, to the transition between phase portrait 1 where is a global attractor, and phase portrait 4.

: curves of Hopf bifurcation, where the equilibrium changes stability and a periodic solution of small amplitude is created or disappears. The periodic solution is attracting when the bifurcation is supercritical (bold parts of curve H in figure 2(a)), and it is of saddle type when the Hopf bifurcation is subcritical. This corresponds, for example, to the transition between phase portraits 7 and 8 in figure 3.

: a curve of homoclinic loop to saddle equilibrium . Crossing this curve, a periodic orbit of very large period is either created or disappears. This is illustrated in figure 3 by the transition between phase portraits 3 and 5.

: curves of saddlenode bifurcations of periodic orbits, where two periodic orbits collide and disappear (or are created). An example is the transition between the phase portraits 8 and 12.

: curves of torus bifurcations of periodic orbits, where a periodic orbits has two complex conjugate Floquet multipliers crossing the unit circle, and an invariant torus bifurcates. This corresponds, for example, to the transition between phase portraits 14 and 15.

: degenerate Hopf bifurcation points, where the Hopf bifurcation changes criticality and a curve SL emerges.

: neutral saddle points along the homoclinic loop curve L, where the periodic orbit associated to the homoclinic loop changes stability, and a curve SL emerges.

: FoldHopf bifurcation points, lying at the tangential intersection of the curve S and the curve H, where the system linearized around the equilibrium has a zero eigenvalue and two complex conjugate, purely imaginary eigenvalues.

: a BogdanovTakens point BT, lying along the curve S, where the system linearized around the equilibrium has a zero eigenvalue of multiplicity two. Generically, a curve H of Hopf bifurcation and a curve L of homoclinic loop emerge from this point. The periodic orbit arising from the Hopf bifurcation undergoes the homoclinic bifurcation.

: HopfHopf bifurcation points, where two curves of Hopf bifurcations intersect, or a single Hopf bifurcation curve selfintersects. The system has two pairs of complexconjugate, purely imaginary eigenvalues. Curves TR of torus bifurcation can emerge from these points.
Note that the saddlenode bifurcation of equilibria involving and does not occur in figure 2: for as considered in this section, the value defined by equation (3) is negative, and , and consequently all exist in the physically relevant halfplane with .
Figure 2(a) shows the precise locations of the curves of saddlenode bifurcations of limit cycles, which have now be continued numerically with DDEBiftool [9, 49]. As theory predicts, they emerge from degenerate Hopf points DH and from codimensiontwo points N, where the Hopf bifurcations change criticality and the periodic orbit involved in the homoclinic loop changes criticality, respectively. Panel (b) represents an enlargement of the bifurcation diagram for low values of . In this region, the curves H and L are found to cross each other twice, near and . The criticality of L changes twice at codimensiontwo points N located at and . The curve SL emerging from the first point N approaches the limit , while the curve SL emerging from the second point N connects to the degenerated Hopf bifurcation point DH located near . Between its two boundaries, this curve SL is just below the curve L and stays very close to it. Together with curves H and L, the curve SL bounds a very small region 13 of the bifurcation diagram, whose phase portrait is represented in figure 3. In this phase portrait, both the fact that the two periodic orbits are almost superimposed, and really close to the saddle equilibrium , clearly highlights the proximity of both the saddlenode bifurcation of limit cycles and the homoclinic bifurcation. Note that phase portrait 13 highlights an additional small intensity peak along both periodic orbits, when the intensity is close to zero. This can be explained by the detuning between the delay and the period : approaching the homoclinic bifurcation curve L, the period becomes much larger than the delay, and the second peak is thus due to a secondary excitation by the feedback. However, this feedback reexcitation occurs during the refractory period, when the gain has not entirely recovered: the next highamplitude pulse cannot be triggered, and the returning pulse only results in a small increase and relaxation of the intensity.
Even for relatively small values of the delay, the bifurcation diagram 2(a) highlights new feedbackinduced dynamics. This can be seen clearly when one considers the different dynamics encountered along a crosssection of the diagram, for a fixed value of the feedback strength . For small values of the delay (), in region 9, two equilibria exist and are physically relevant: is attracting, and is unstable. When the delay is increased, the Hopf bifurcation curve H is crossed to enter region 7: the bifurcating equilibrium becomes unstable and a smallamplitude periodic solution is created; see figure 3 for the phase portrait. Increasing the delay further, the curve H is crossed a second time and undergoes a second (subcritical) Hopf bifurcation: entering region 8, a new periodic orbit of saddle type is created, while becomes stable again. The same process repeats when regions 11 and 12 are entered: the Hopf bifurcation curve is crossed two additional times, leading each time to the birth of a new periodic orbit, which is alternately stable and of saddle type. Different stable periodic solutions with different amplitudes thus coexist in regions 10, 11 and 12. Physically, this means that the laser can display different qualitative and quantitative behaviors for the same parameter values: depending only on the initial conditions, the output is either a light beam with constant intensity, or pulselike selfoscillations with different amplitudes and periods.
Extending the range of the delay
Although the values of considered in figure 2(a) can seem large compared to the internal time scale of the SLSA, these values remain small from an experimental point of view. Indeed, optical paths of up to several meters are commonly used for the feedback loop in the experiment [32]. Figure 4(a) presents an extended bifurcation analysis for larger values of . It is a wellknown property of delay systems, that periodic solutions reappear infinitely many times for different values of the delay [59]. The consequence here is that the curve H starts to selfintersect for sufficiently large values of the delay. Each of these intersection points is a HopfHopf bifurcation point HH, where the system has two pairs of complexconjugate, purely imaginary eigenvalues. From each of these codimensiontwo points, two curves TR of torus bifurcation of periodic orbits can emerge. When a periodic orbit undergoes a torus bifurcation, two complexconjugate Floquet multipliers cross the unit circle, and a (stable or saddletype) invariant torus bifurcates from the periodic solution [44]. These additional bifurcation curves make the bifurcation diagram in figure 4(a) increasingly complicated as the delay is increased. In practice, it becomes impossible to map out the dynamics in each of the many different regions of the ()plane. Considering, as before, a crosssection of the bifurcation diagram for a fixed value of highlights new feedbackinduced dynamics. Figure 4(b) represents an enlargement of the bifurcation diagram in panel (a) for , and figure 5 presents some of the additional phase portraits encountered along this section while increasing .
Starting from region 12 where we stopped in section 3.1, increasing leads to the crossing of the Hopf bifurcation curve H to enter region 14, where we find the coexistence of three stable periodic solutions. When the torus bifurcation curve TR, which separates regions 14 and 15, is crossed one of the unstable periodic orbits existing in region 14 undergoes a torus bifurcation: in region 15 it thus has two additional Floquet multipliers outside the unit circle. Numbers near unstable periodic orbits in figure 5 indicate the number of Floquet multipliers outside the unit circle, when this number is larger than one. It is not possible with DDEBiftool to determine the criticality of the torus bifurcation; this thus has to be investigated through timedomain simulations, it will be discussed in section 3.1.3. The transition from region 15 to 16 is of the same nature as the transition between regions 11 and 12 previously described: a new unstable periodic orbit is thus observed in region 16. Up to this point, increasing the delay for the fixed value has led to coexistence of more and more stable periodic orbits. However, when the curve SL of saddlenode bifurcation of periodic orbits is crossed to enter region 17, the stable and the unstable periodic orbits with the largest amplitudes collide and disappear. In region 17, the system thus goes back to multistability between only two periodic orbits and equilibrium . Transition between regions 17 and 18, through the crossing of a torus bifurcation curve, is very similar to the transition from region 14 to region 15. The curve H is crossed to enter region 19, and the phase portrait thus again shows the coexistence of three different stable periodic orbits.
Multistability: features of periodic orbits
Figure 6 represents, for fixed , the evolution of the period of periodic orbits with respect to in panel (a) and the time series of periodic orbits at in panels (b), (c) and (d). For the parameter set considered here and for all the branches of periodic solutions, is never smaller than at the Hopf bifurcation point, and is always at the saddlenode bifurcation constituting the other end of the branch. Physically, these values are related to the recovery and saturation times of the gain and absorber media; that is to say, it is determined by the value of in equations (1).
In figure 6(a), the period evolves practically linearly with the delay , with different slopes depending on the considered periodic solution branch. Along the first pair of periodic orbit branches – one attracting which is born at , and one of saddle type which is born at – the period stays really close to the value of the delay. In other words, the period evolves with a slope of 1 with respect to the delay along these branches. The period of the two following periodic orbits – one attracting and one of saddletype, which are born at and at , respectively – stays close to half the delay along their branches. In the same way, the next pair of periodic orbits have, along their corresponding branches, a period close to a third of the delay, and so on. This corresponds to one, two or three pulses per roundtrip of the external feedback loop.
For a fixed value of the delay of , we now compare the time series corresponding to the different periodic orbits. For the considered value of , figure 6(a) shows the coexistence of three pairs of periodic orbits (b), (c) and (d), each pair consisting of an attracting and a saddletype periodic solution. Time series of the intensity are shown in panels (b1), (c1) and (d1) for the stable periodic orbits, and in panels (b2), (c2) and (d2) for the saddle periodic orbits. Beyond the differences in terms of the repetition rate of the pulses already discussed above, this highlights significant differences in the temporal shape of the laser output. Panels (b1) and (b2) show highamplitude, narrow pulses, with low repetition rates of about . Compared to point (b), point (c) in figure 6(a) is closer to the respective Hopf point where the periodic solution branches are born. Compared to times series (b1) and (b2), the pulselike oscillations in time series (c1) and (c2) shows wider pulses, with smaller amplitude and higher repetition rate of about . In the same way, at point (d) we are closer to the Hopf bifurcation points where the considered periodic orbits are born than for point (c). In particular, the periodic orbit of saddletype is born at the Hopf point at . The corresponding time series for in panel (d2) thus shows, as expected, smallamplitude and almost sinusoidal oscillations. The stable periodic orbit is born at ; at point (d), we are thus further away from the Hopf bifurcation than in the case of the saddletype periodic orbit. The corresponding time series in panel (d1) is therefore more pulseshaped, but the pulses are wider, and their amplitude smaller than for times series (b1) and (c1).
Bifurcating tori
So far, we only considered multistability between periodic orbits and equilibria. We now investigate transitions through torus bifurcations curves TR, which appear for sufficiently large , for example, between regions 22 and 23 and at the transition into region 28. In both cases, the periodic orbit undergoing the torus bifurcation is stable before the bifurcation (i.e. for smaller values of the delay, along the cross section defined by ), and thus becomes unstable in the torus bifurcation. Because invariant tori cannot be continued for delay systems, we run timedomain simulations with dde23 [47], starting from the bifurcating periodic orbit as initial condition.
Figure 7 represents, in panels (a) and (b), the stable tori in regions 23 and 28, respectively. In both cases, the time series in panels (a1) and (b1) suggest either quasiperiodic oscillations or highperiod locked periodic solutions on a torus. Projections of these solutions onto the physical space in panels (a2) and (b2) highlight the two associated invariant tori. The spectra of both time series, in panels (a3) and (b3), are clearly non harmonic, and both of them are composed of discrete peaks at frequencies equal to linear combinations of two base frequencies. Altogether, figure 7 confirms the quasiperiodic feature of the oscillations as taking place on an invariant torus.
Both in regions 23 and 28, the system thus displays multistability between a stable invariant torus and different periodic orbits (three in region 23, and four in region 28). The stable invariant torus in region 23 disappears when the same torus bifurcation curve TR is crossed again, which occurs at for ; see figure 4(b). At this point, the bifurcating periodic orbit regains stability. For the considered values of , we do not observe a coexistence of the two stable invariant tori. However, the increasing complexity of the bifurcation diagram for large values of suggests that the coexistence of several stable invariant tori and stable periodic orbits is to be expected. The dynamics of the Yamada model without delay has been shown to take place in a globally attracting twodimensional surface [6]. The existence of stable tori clearly shows that the dynamics of the SLSA with feedback is no longer twodimensional for sufficiently large .
3.2 Case of the SLSA in a nonexcitable off state
From a practical point of view, the influence of the pump parameter is particularly interesting, because it can easily be changed experimentally. We now consider the influence of the delayed optical feedback on the dynamics of the Yamada model for a fixed value of . In this case, the SLSA without feedback is in a nonexcitable regime: the nonlasing equilibrium is the unique steadystate [6]. The influence of the delayed optical feedback on this phase portrait has been investigated in [29], but for small values of only. For the same configuration, figure 8 (a) represents an extended bifurcation diagram in the plane. The curve S, where equilibria and bifurcate, and the curve T, where and exchange stability, are located at and , respectively; see section 2. Contrary to the excitable case in section 3, multiple disjoint curves H of Hopf bifurcation are found, which either involve or . Moreover, codimensiontwo bifurcation points occur along the curve S: two FoldHopf points FH, and a BogdanovTakens point BT; see section 3 for their description. A curve H and a curve L of homoclinic loop (involving the periodic orbit that appears at the Hopf bifurcation) emerge from BT, and they both end at the left limit of the bifurcation diagram. As before, curves TR emerge from HopfHopf points HH where curves H intersect, and SL curves emerge from degenerate Hopf bifurcation points DH where a Hopf bifurcation change criticality. Note that now, the different curves SL tend to an asymptotical value of when is increased: continuing one of the branch SL up to clearly shows that it remains located between and . These qualitative differences with the case of the SLSA in the excitable regime in figure 4 has multiple consequences for the system’s dynamics. Since curve S in figure 8 is in the physicallyrelevant halfspace , the equilibrium is the only stable solutions for . Moreover, remains stable for , that is to say, in most regions of the bifurcation diagram. Here, multistability thus involves not only pulsing regimes and the continuous wave solution , but also the offstate .
As before, the complexity of the bifurcation diagram in figure 8(a) makes it impossible in practice to map out the dynamics in all the different regions. However, the enlargement of the bifurcation diagram for the crosssection near , in panel (b), highlights some of the new dynamics. Along this crosssection, no new dynamics are found for , and the different regions correspond to phase portraits already observed for the case of the SLSA in the excitable regime; see figures 3 and 5. For the selected regions where new dynamics are found for , figure 9 represents the additional phase portraits. Namely, a curve H of Hopf bifurcation involving the equilibrium is crossed to enter region 29, and a smallamplitude periodic orbit is created around , which can be seen in the phase portraits of regions 31 and 32. For the range of parameters considered in figure 8(a), this Hopf bifurcation is subcritical, and the corresponding periodic orbit stays unstable when varying and . Additionally, no curve SL is crossed along the crosssection , as they all tend to an asymptotical value close to . On the other hand, more and more curves H are crossed when increasing , and more and more periodic orbits are thus created. As a consequence, the phase portraits in regions 31 and 32, in figure 9, show the coexistence of the stable nonlasing equilibrium with four and five stable periodic orbits, as well as five and six unstable periodic orbits, respectively. In contrast, no more than four coexisting stable periodic orbits were found, for the same range of , for the case of the SLSA in the excitable regime. As discussed in section 3.1.2, the period of the different periodic orbits increases linearly with . Here, this means that more and more stable periodic solutions are observed when is increased, which correspond to narrow, highamplitude light pulses. Moreover, a stable torus is born while entering region 31, which is computed by timedomain simulation and represented in black in figure 9(a). This torus being very thin, the corresponding pulsing solution only shows a very weak amplitude modulation. From a practical point of view, it can thus be easily mixed up with a periodic orbit.
Although the introduction of a delayed optical feedback leads to qualitatively different behaviours for the SLSA in an excitable and a non excitable regime, several common features can be highlighted. In both cases, the bifurcation diagrams quickly becomes intricate when the feedback delay in increased beyond value of studied in [29], with many bifurcations occurring in small parameters regions. While the case of the SLSA in the non excitable regime corresponds to very simple dynamics, the introduction of delayed feedback leads to a complicated bifurcation diagram, even for relatively small delay times compared to the experimental configuration [32]. In both cases considered so far, multistability considerably increases with the value of , so that coexistence of several pulsing solutions is always observed when the delay is large enough. Finally, the systematic presence of torus bifurcations demonstrates that, in contrast to the system without delay [6], the dynamics is not twodimensional anymore from intermediate values of . This also suggests that coexistence of several stable tori is to be expected for large values of the delay.
4 Bifurcation study in the plane
To study the influence of the feedback, we considered so far bifurcation diagrams in the plane. Firstly to understand the influence of both feedback parameters, and secondly because the DDE for is a regular perturbation of the ODE for , whose dynamics has been thoroughly studied in [6]. The comparison of the bifurcation diagrams in the plane for two values of the pump parameter A (in figures 4 and 8) has shown that the bifurcation scenario strongly depends on , not only quantitatively but also qualitatively. We focus in this section on bifurcation diagrams of system (1), in the plane of feedback strength and pump parameter . These are particularly relevant from a practical point of view: the pump A is the main control parameter of the laser, and can easily be tuned with an attenuator. The delay time , on the other hand, cannot be changed easily in an experiment.
Figure 10 shows the three bifurcation diagrams in the plane for in panel (a), in panel (b) and in panel (c). For small value of , in panel (a), the bifurcation scenario remains relatively simple, and the dynamics in all the different regions can easily be mapped. For an intermediate value of , in panel (b), both the bifurcation diagram, and the sequence of bifurcations encountered when increasing A for a fixed value of , become more intricate. Nevertheless, it remains possible to identify the dynamics in the different regions, especially for reasonably low values of , considered in practice. On the other hand, when the delay is increased further, in panel (c), the bifurcation diagram has considerable complexity: many different bifurcation curves and regions with different dynamics are found in small parameter ranges, so that it is hard to see the details of the bifurcation scenario. Figure 10 provides an additional and complementary representation of the increasing complexity of the dynamics when becomes larger, by illustrating how crosssections of the bifurcation diagrams in figures 4 and 8 for fixed , evolve with the pump parameter A.
More precisely, for in figure 10 (a), most of the phase portraits have already been observed for the system without feedback [6]. A new feature, however, is the occurence of two cusp bifurcations of saddlenode bifurcations curves SL: one cusp is found close to , and the second one is highlighted in the inset of panel (a). Here, bistability occurs between equilibria and a single periodic orbit; however is not sufficiently large to observe multistability between several periodic orbits. We conclude that, although the sequence of bifurcations encountered when increasing depends on , some constant features are observed. Namely, the SLSA is necessarily in the offstate for small (region 1), and on the continuouswave solution for large (region 9). In between these two regions, selfpulsations can occur in regions 5, 6, 7 and 8. However, all of them except for 7 show bistability between a pulsing solution and at least one equilibrium. From a practical point of view, this means that region 7 is the only one where selfpulsations are necessarily observed.
For the larger value of in panel (b), the bifurcation diagram becomes more complex: the accumulation of both cusp points and curves H, SL, and TR around makes it impossible in practice to distinguish between all the different regions. For this value of , new dynamics is found compared to the case without feedback, including multistability between two different periodic orbits in regions 10 and 11. As a consequence, for large values , the bifurcation scenario observed when increasing or decreasing A becomes more complicated: it includes jumps between different periodic solutions, hysteresis phenomena, and many different dynamics in a tiny parameter range. On the other hand, for small values , which are most often considered in experiments, no new dynamics is observed compared to the case without feedback. Moreover, the bifurcation sequence observed when increasing or decreasing A for fixed values of is similar as that for .
For the delay of in panel (c), the bifurcation diagram clearly illustrates the increasing complexity of the system dynamics. Namely, an additional curve H is found, along which both points DH and HH occur. The emergence, from these points, of additional curves SL and TR, respectively, leads to a particularly intricate diagram. In particular, the bottom inset shows an enlargement near the point , where more and more bifurcation curves and codimensiontwo bifurcation points accumulate when is increased. In the same way, the top inset is an enlargement near the BogdanovTakens point BT: it clearly shows that, for this value of , the bifurcation diagram becomes intricate even for small values of , with many bifurcation curves in a small parameters range. Some of the numerous new regions found in figure 10 (c) show multistability between different pulsing solutions. It is possible to observe up to three stable periodic solutions, corresponding to light pulses of different amplitudes and repetition rates. On the other hand, in spite of the existence of many curves TR, no stable torus is found for this value of .
All together, the bifurcation diagrams in the plane and in the plane show that the SLSA’s dynamics becomes more and more complex as the delay is increased. Even for relatively small values of , the feedback triggers complex bifurcation scenario, with a wealth of new interesting dynamics, including multistability of pulselike periodic solutions with different amplitude and repetition rate. It is of practical importance that these new dynamics, which were first highlighted in the plane, are also accessible in the plane for small values of . Indeed, in an experiment, only values of can usually be considered [32]. This means that the new regimes described above can be observed in an experiment when is changed.
5 Consequences of multistability
We focus here on the consequences of the observed high level of multistability. Specifically, we consider the example of phase portrait , in figure 11(a) for . Region is adjacent to region 15 in the bifurcation diagram in figure 4, for , below the curve T of transcritical bifurcation. The stable solutions in this region include the equilibrium and three periodic orbits with different amplitudes and periods. As the three periodic orbits are very close to each other where their intensity is small, one can expect the system to be very sensitive to small perturbations. To investigate the features of multistable dynamics in more detail, we focus on the basins of attraction associated with the different attractors. In particular, knowing the basins’ structure highlights low or high sensitivity of the system to the initial conditions.
In the experiment, when the laser is initially off, pulses are obtained by providing an external perturbation either on the gain or on the intensity [45]. From a practical point of view, investigating the basins structure corresponds to determining the attractor on which the laser settles when it is perturbed from the off state. This question is addressed by timedomain simulations, in which we consider the equilibrium as initial condition and focus on the influence of small initial perturbations on the gain and on the intensity , so that:
(8) 
For the parameters considered here, the system is close to the threshold at which becomes unstable, and a small but nonzero perturbation is enough to leave the basin of the offstate . Figure 11(b) represents in the plane the stable solution to which the system eventually settles in timedomain simulations, each solution being represented by the color of the corresponding attractor in the phase portrait in figure 11(a). As we deal here with a delay system, the basins associated with the different attractors live in infinitely many dimensions, and the attractor map in figure 11(b) is the trace of the basins structure in the plane. The blue region is the basin of the equilibrium solution : if the intensity perturbation is smaller than a threshold which depends on the value of , the laser relaxes back to its offstate without producing a pulse. The red region is the basin of the periodic orbit with the largest amplitude. Hence, a sufficiently large intensity perturbation always results in pulsations of the laser with a pulse repetition rate close to the delay time . Again, the exact value of from which this basin is reached depends on the value of . The same solution is observed for a sufficiently large, positive perturbation , as long as is not too small. In between these two regions, the structure of the basins is much more intricate. The attractor map in figure 11(b) suggests that it is quite complicated in practice to predict both the transient and the longterm behaviour of the SLSA. The computation of the map in figure 11 is particularly timeconsuming. Not only is there a need for a fine mesh on both and to highlight the basins structure, but long transients are observed before the system settles to an attractor. In order to determine the attracting solution, simulations are thus run for a duration . As a matter of fact, the computation of the map in figure 11(b) takes approximately a week with the method of steps for DDEs (implemented in C), on a desktop machine with three processors in parallel.
Knowing the trace of the basins of attraction in the plane is not only interesting from an application point of view, but also from the perspective of the theoretical analysis of the model. Indeed, the complicated structure of the basin boundaries, highlighted in figure 11, gives insight into the structure of the stable manifolds of unstable solutions and, thus, a more global overview of the system dynamics. This is of particular interest because we deal here with a system with delay, with an infinitedimensional phase space, meaning that it is not possible to compute the invariant manifolds with methods for ODEs [25, 26].
The attractor map in figure 11(b) clearly highlights the stripe structure of the basins, and suggests that it selfrepeats at different scales when is varied. To study the basins structure in more details, three crosssections for , and are computed with a much finer mesh size of 0.0001 along . These three sections, represented in panels (c), (d) and (e) of figure 11, clearly have a similar structure, which repeats at different scales for the three different values of . Moreover, this highlights the selfrepeating structure of the basins boundaries along , with the same pattern repeating along the sections. In panel (c) for example, the pattern found for , composed of a large red stripe and the upper green and orange stripes, repeats many times at smaller scales when decreases. The basins structure is thus particularly intricate for small , where the pattern repeats at very small scales. In this pattern, the largest green stripe is bounded by two orange stripes, one at the top and the other at the bottom. Each of these orange stripes is in turn the lower and the upper limit, respectively, of other green stripes. Enlarging the crosssection in this region reveals that the same structure repeats again and again at increasingly smaller scale, leading to thiner and thiner stripes in very small region of the plane. This clearly shows that the boundaries of the basins associated with the three attracting periodic solutions accumulate on each other in a Cantor structure. It is interesting to note that this Cantor structure does not involve the basin of the nonlasing equilibrium, which has a single boundary curve with the other basins in the plane. Note however, that the three basins boundaries of the periodic orbits accumulate on the basin boundary of the offsolution. This selfrepeating property along , which is qualitatively independent of , implies that the basins’ boundaries are intermingled, which is a property that can be found in other complex systems without delay, including the Lorenz system [22, 5]. The most obvious consequence is that the prediction of the asymptotic state of the laser may become practically impossible. Moreover, it also implies that trajectories can take a very long time before approaching a given attractor.
The three similar crosssections in figure 11 also allow us to determine the scale at which the pattern repeats along the whole range of . In figure 11(b), the boundaries between the different stripes are fitted by thirdorder polynomials, given by:
(9) 
From the precise representation of the basins along the three crosssections with fixed , we estimate the set to obtain the best fit for each stripe boundary. We find that the stripe structure is represented by only two oneparameter families of thirdorder polynomials (9), a family having fixed values of , and , and a single changing parameter . More precisely, the family with represents the uppermost stripes (i.e., the green and orange stripes with the largest values of ), and the family with represents all of the other stripes. For each family, two polynomial boundary curves, corresponding to two different values of the coefficient , are plotted in white in figure 11(b).
The parameters and can then be use to extend the accurate crosssection to the whole range of . In this way, we obtain a refined version of the attractor map without running CPUconsuming simulations on a finer mesh for the whole range of and . The resulting much finer map of the basins structure is shown in figure 12. Clearly, the SLSA is less sensitive to perturbations for negative values of . Conversely, for the range of in which the different basins are intricate is smaller; on the other hand in this range, the system is extremely sensitive to small perturbations and the behaviour of the SLSA thus easily becomes practically unpredictable.
Associated with the Cantorset structure of the basins is a wealth of transient dynamics exhibited by the system. Figure 13 presents three time series observed for and , and in panels (a), (b) and (c), respectively. Each of the three values of leads to a given stable periodic solution of the phase portrait found in figure 11(a). The comparison between the three transients highlights the different time scales on which they take place. Moreover, as shown in the inset of figure 13(a), complicated dynamics can occur during the transients, including multipulses, where several pulses occur before the intensity drops back to zero. Both the structure of the basins of attraction and the complex transient dynamics observed here in a deterministic configuration suggest that even a small amount of noise, as necessarily present in an experiment, may induce complicated, even practically unpredictable dynamics. Moreover, in an actual experiment, the delay time can be up to [32], and an even more complex configuration, with more stable and unstable periodic orbits, is to be expected [59]. One can then expect intermingled basins boundaries, involving many more basins, each of which is associated with one of the many attracting periodic orbits.
6 Discussion
We investigated the effect of delayed optical feedback on the dynamics of the Yamada model for selfpulsations in a semiconductor laser with saturable absorber (SLSA). The bifurcation analysis in the plane of the delay time and feedback strength shows that the bifurcation diagram quickly becomes very complex when the delay is increased. The case of the SLSA in the nonexcitable regime is particularly interesting: the zerointensity equilibrium being the only solution when , one might expect that the feedback only has small consequences on the laser’s dynamics. However, for nonzero feedback, the laser shows a wealth of new dynamics from intermediate values of . This includes multistability of several pulselike periodic solutions with different amplitudes and periods, as well as stable quasiperiodic oscillations on tori. In particular, the level of multistability increases significantly with the delay time : more and more attracting selfpulsing solutions coexist, where repetition rates of the pulses are submultiples of . We showed that these feedbackinduced dynamics can be found in the plane of the main control parameters, which means that they are accessible experimentally. These results show that, despite its simplicity, the Yamada model with feedback is able to produce a wealth of complex dynamics. Some of the scenarios observed in the Yamada model with feedback appear to correspond well to dynamics observed in an actual experiment, which at present time have been reported only in [32, 28]; an actual comparison between the model dynamics and experimental observations is the subject of ongoing research.
We especially focused on features of multistability, and showed evidence of the Cantorset like structure of the basins of the different pulsing solutions, by considering their trace in the plane of perturbations of gain and intensity . It follows that the system is very sensitive to perturbations, especially when the intensity is small. A small perturbation can thus be sufficient to induce the system to jump between different basins. This suggests that for large values of the delay, as often considered in practice, the basins of the many coexisting stable periodic solutions might have intermingled boundaries. In particular, this means in practice that a small amount of noise may trigger unpredictable dynamics.
Highsensitivity to noise might explain complex phenomena observed in an experiment, including multipulses and jumps between different stable solutions [32]. An ongoing study suggests that the varying life times of pulse trains observed in the experiment are explained by noiseinduced transitions from a stable pulsing solution to the offstate equilibrium. The Yamada model with feedback and additional noise terms accurately reproduces this phenomenon, showing good agreement between the experimental and numerical escape times [3]. In future work, we will focus on a more direct comparison between the experimental observations and the results from the theoretical analysis of the model. In particular, ongoing experiments investigate some of the complex dynamics that we discussed above. It aims at showing experimental evidence of multistable dynamics, and associated characteristics of the attractors’ basins. Although this last point is particularly challenging to demonstrate in an experiment, it can nevertheless be directly related to the observed sensitivity of the SLSA to noise and small perturbations. Ongoing experimental investigations also focus on quasiperiodic dynamics: because, in some cases, they may look similar to noisy periodically pulsing solution, such dynamics is quite difficult to identify experimentally [3].
Overall, our results suggest that the Yamada model accurately describes pulsing lasers with feedback as long as there is no other internal dynamics than pulsing or excitability. In future work, it would be interesting to compare the dynamics of the Yamada model with feedback with other types of pulsing lasers subject to feedback, in particular with passively modelocked lasers. Although the underlying physical mechanism responsible for selfpulsations in passively mode locked lasers is different from Qswitching, which implies the use of different mathematical models, the dynamics of both kinds of lasers display interesting similarities when subject to feedback. As an example, the increasing level of multistability observed when the delay becomes larger was also highlighted for a passively modelocked ring cavity laser with optical feedback: more and more pulsing solutions coexist, which correspond to different numbers of pulses in a cavity roundtrip [38, 18]. A recent study focused on multistable dynamics between several pulsing solutions and a nonlasing equilibrium in a passively mode locked laser subject to feedback [34]. It showed that, when the delay becomes larger than the lasers internal timescales, pulses can be triggered which almost behave as if they are independent from each other: arbitrary sequences of pulses can thus be addressed, where pulses are not necessarily regularly spaced. Such dynamics is often referred to, in the literature, as temporal cavity solitons, or localized structures. Despite the model we consider here being quite different, a preliminary study suggests that both the Yamada model subject to feedback with large delay and an actual experiment display, phenomenologically, this type of dynamics [3]. Moreover, some aspects of the dynamics of the Yamada model are very reminiscent of phenomena recently observed in neuromimetic excitable, slowfast photonic systems subject to feedback [43]. This type of device has also been shown to produce multiple localized structures. Overall, it suggests that the Yamada model might also be a minimal model suitable for the investigation of such phenomena.
An indepth comparison between these different systems would be valuable to better understand the underlying mechanisms of such common dynamics, as well as to provide insight regarding the domain in which the Yamada model is valid. In particular, conversely to modelocked lasers which may feature other intrinsic dynamics [17, 53], the device we consider here only displays excitability or pulsing in the absence of feedback. In this respect, the Yamada model is very similar to certain models of pulsing neurons [15, 45].
The extended bifurcation analysis of a pulsing SLSA with selffeedback presented here might, thus, be of wider interest to the applied dynamical systems community. Especially, it can be seen as a first step in the investigation of more general delaycoupled pulsing systems. The next step would be the case of two delaycoupled pulsing SLSA, where the system studied here appear as a special case where the two lasers are identical. While coupled lasers have been extensively studied both experimentally and theoretically, it would be particularly interesting to focus on the effect of delayed coupling when the lasers are operating strictly in the pulsing regime. Again, such a configuration is closely related to that of connected excitable neurons.
Acknowledgement
The authors thank Sylvain Barbay and Kathy Luedge for stimulating discussions, as well as the referees for helpful comments and suggestions. The research of S. T. has been funded by a Postdoctoral Fellowship of The DoddWalls Center for Photonic and Quantum Technologies.
Appendix A Parameters value
Fixed parameters: ; ; .
figure  A  figure  A  
2  6.5  5 panel 18  6.5  0.38  80.5  
3 panel 3  6.5  0.16  2.6  5 panel 19  6.5  0.38  90.5 
3 panel 5  6.5  0.25  38  6 panel (a)  6.5  0.38  
3 panel 6  6.5  0.25  25  6 panels (b), (c) and (d)  6.5  0.38  100 
3 panel 7  6.5  0.38  21  7 panels (a)  6.5  0.38  117 
3 panel 8  6.5  0.38  33  7 panels (b)  6.5  0.38  140 
3 panel 10  6.5  0.25  52  8  5.9  
3 panel 11  6.5  0.38  43  9 panel (a)  5.9  0.56  102 
3 panel 12  6.5  0.38  55  9 panel (b)  5.9  0.55  117 
3 panel 13  6.5  0.0152  26.8  10 panel (a)  25  
4  6.5  10 panel (b)  40  
5 panel 14  6.5  0.38  63  10 panel (c)  58  
5 panel 15  6.5  0.38  68  11  6.5  0.29  80 
5 panel 16  6.5  0.38  75  12  6.5  0.29  80 
5 panel 17  6.5  0.38  79  13  6.5  0.29  80 

Footnotes
 The DoddWalls Centre for Photonic and Quantum Technologies, Department of Mathematics, The University of Auckland, New Zealand. corresponding author: s.terrien@auckland.ac.nz
 footnotemark:
 The DoddWalls Centre for Photonic and Quantum Technologies, Department of Physics, The University of Auckland, New Zealand.
References
 N. B. Abraham, L. A. Lugiato, and L. M. Narducci. Overview of instabilities in laser systems. JOSA B, 2(1):7–14, 1985.
 G. P. Agrawal and N. K. Dutta. Semiconductor Lasers. Electrical Engineering. Van Nostrand Reinhold, 1993.
 S. Barbay. personnal communication.
 S. Barbay, R. Kuszelewicz, and A. M. Yacomotti. Excitability in a semiconductor laser with saturable absorber. Optics letters, 36(23):4476–4478, 2011.
 E. J. Doedel, B. Krauskopf, and H. M. Osinga. Global invariant manifolds in the transition to preturbulence in the lorenz system. Indagationes Mathematicae, 22(3–4):111–240, 2011.
 J. L. A. Dubbeldam and B. Krauskopf. Selfpulsations of lasers with saturable absorber: dynamics and bifurcations. Optics communications, 159(4):325–338, 1999.
 J. L. A. Dubbeldam, B. Krauskopf, and D. Lenstra. Excitability and coherence resonance in lasers with saturable absorber. Phys. Rev. E, 60(6):6580–6588, 1999.
 T. Elsass, K. Gauthron, G. Beaudoin, I. Sagnes, R. Kuszelewicz, and S. Barbay. Control of cavity solitons and dynamical states in a monolithic vertical cavity laser with saturable absorber. The European Physical Journal D, 59(1):91–96, 2010.
 K. Engelborghs, T. Luzyanina, and D. Roose. Numerical bifurcation analysis of delay differential equations using DDEBIFTOOL. ACM Transactions on Mathematical Software (TOMS), 28(1):1–21, 2002.
 K. Engelborghs, T. Luzyanina, and G. Samaey. DDEBIFTOOL v. 2.00: a Matlab package for bifurcation analysis of delay differential equations. Technical report, Department of Computer Science, KU Leuven, Leuven, Belgium, 2001.
 T. Erneux. Qswitching bifurcation in a laser with a saturable absorber. J. Opt. Soc. Am. B, 5(5):1063–109, 1988.
 M. Georgiou and T. Erneux. Pulsating laser oscillations depend on extremelysmallamplitude noise. Physical Review A, 45(9), 1992.
 M. Giudici, C. Green, G. Giacomelli, U. Nespolo, and J. R. Tredicce. Andronov bifurcation and excitability in semiconductor lasers with optical feedback. Phys. Rev. E, 55(6), 1997.
 K. Green and B. Krauskopf. Global bifurcations and bistability at the locking boundaries of a semiconductor laser with phaseconjugate feedback. Physical Review E, 66(1), 2002.
 E. M. Izhikevich. Which model to use for cortical spiking neurons? IEEE transactions on neural networks, 15(5):1063–1070, 2004.
 Eugene M Izhikevich. Dynamical systems in neuroscience. MIT press, 2007.
 L. Jaurigue, A. Pimenov, D. Rachinskii, E. Schöll, K. Lüdge, and A. G. Vladimirov. Timing jitter of passivelymodelocked semiconductor lasers subject to optical feedback: A semianalytic approach. Phys. Rev. A, 92(5):053807, 2015.
 L. Jaurigue, E. Schöll, and K. Lüdge. Passively modelocked lasers subject to optical feedback: The role of amplitudephase coupling. In Numerical Simulation of Optoelectronic Devices, 2014, 2014.
 C. K.R.T. Jones. Geometric singular perturbation theory. In Dynamical systems, pages 44–118. Springer, 1995.
 D. M. Kane and K. A. Shore, editors. Unlocking dynamical diversity: optical feedback effects on semiconductor lasers. John Wiley & Sons, 2005.
 A. Keane, B. Krauskopf, and C. Postlethwaite. Delayed feedback versus seasonal forcing: Resonance phenomena in an el nino southern oscillation model. SIAM Journal on Applied Dynamical Systems, 14(3):1229–1257, 2015.
 Judy Kennedy. How indecomposable continua arise in dynamical systemsa. Annals of the New York Academy of Sciences, 704(1):180–201, 1993.
 B. Krauskopf, G. R. Gray, and D. Lenstra. Semiconductor laser with phaseconjugate feedback: Dynamics and bifurcations. Physical Review E, 58(6), 1998.
 B. Krauskopf and D. Lenstra, editors. Fundamental Issues of Nonlinear Laser Dynamics. American Institute of Physics, 2000.
 B. Krauskopf, H. M. Osinga, E. J. Doedel, M. E. Henderson, J. Guckenheimer, A. Vladimirsky, M. Dellnitz, and O. Junge. A survey of methods for computing (un) stable manifolds of vector fields. International Journal of Bifurcation and Chaos, 15(3):763–791, 2005.
 B. Krauskopf, H. M. Osinga, and J. GalanVioque, editors. Numerical continuation methods for dynamical systems. Springer, 2007.
 B. Krauskopf, K. Schneider, J. Sieber, S. Wieczorek, and M. Wolfrum. Excitability and selfpulsations near homoclinic bifurcations in semiconductor laser systems. Optics Communications, 215(46):367–379, 2003.
 B. Krauskopf, S. Terrien, N. G. R. Broderick, and S. Barbay. Quasiperiodic dynamics in a micropillar laser with saturable absorber and delayed optical feedback. In Australian Conference on Optical Fibre Technology. Optical Society of America, 2016.
 B. Krauskopf and J. J. Walker. Bifurcation Study of a Semiconductor Laser with Saturable Absorber and Delayed Optical Feedback, pages 161–181. WileyVCH Verlag GmbH & Co. KGaA, 2012.
 Y. A. Kuznetsov. Elements of applied bifurcation theory, volume 112. Springer Science, 2013.
 R. Lang and K. Kobayashi. External optical feedback effects on semiconductor injection laser properties. IEEE Journal of Quantum Electronics, 16(3):347–355, 1980.
 F. Lelievre. Etude de micropilier laser neuromimétique soumis à une rétroinjection avec délai (investigation of a neuromimetic micropillar laser with delayed feedback). Master’s thesis, Laboratoire de Photonique et Nanostructure, CNRS, Marcoussis, France, 2015.
 K. Lüdge, editor. Nonlinear Laser Dynamics: From Quantum Dots to Cryptography. John Wiley & Sons, 2012.
 M. Marconi, J. Javaloyes, S. Balle, and M. Giudici. How lasing localized structures evolve out of passive mode locking. Physical Review Letters, 112(22):223901, 2014.
 J. Mork, B. Tromborg, and J. Mark. Chaos in semiconductor lasers with optical feedback: theory and experiment. IEEE Journal of Quantum Electronics, 28(1):93–108, 1992.
 T. W. Mossberg. Timedomain frequencyselective optical data storage. Optics Letters, 7:77–79, 1982.
 J. Ohtsubo and P. Davis. Chaotic optical communication. In Deborah M. Kane and K. Alan Shore, editors, Unlocking dynamical diversity: optical feedback effects on semiconductor lasers. Wiley, 2005.
 C. Otto, K. Lüdge, A. G. Vladimirov, M. Wolfrum, and E. Schöll. Delayinduced dynamics and jitter reduction of passively modelocked semiconductor lasers subject to optical feedback. New Journal of Physics, 14(11), 2012.
 D. Pieroux and T. Erneux. Bridges of periodic solutions and tori in semiconductor lasers subject to delay. Physical Review Letters, 87(19), 2001.
 D. Pieroux, T. Erneux, and K. Otsuka. Minimal model of a classb laser with delayed feedback: Cascading branching of periodic solutions and perioddoubling bifurcation. Physical Review A, 50(2), 1994.
 K. Pyragas. Continuous control of chaos by selfcontrolling feedback. Physics letters A, 170(6):421–428, 1992.
 E. U. Rafailov, M. A. Cataluna, and W. Sibbett. Modelocked quantumdot lasers. Nat Photon, 1(7):395–401, 07 2007.
 B. Romeira, R. Avó, J. M. L. Figueiredo, S. Barland, and J. Javaloyes. Regenerative memory in timedelayed neuromorphic photonic resonators. Scientific reports, 6, 2016.
 D. Roose and R. Szalai. Continuation and bifurcation analysis of delay differential equations. In Numerical continuation methods for dynamical systems, pages 359–399. Springer, 2007.
 F. Selmi, R. Braive, G. Beaudoin, I. Sagnes, R. Kuszelewicz, and S. Barbay. Relative refractory period in an excitable semiconductor laser. Physical review letters, 112(18), 2014.
 F. Selmi, R. Braive, G. Beaudoin, I. Sagnes, R. Kuszelewicz, and S. barbay. Temporal summation in a neuromimetic micropillar laser. Opics Letters, 40(23):5690–5693, 2015.
 L. F. Shampine and S. Thompson. Solving ddes in matlab. Applied Numerical Mathematics, 37(4):441–458, 2001.
 K. A. Shore. Nonlinear dynamics and chaos in semiconductor laser devices. Solidstate electronics, 30(1):59–65, 1987.
 J. Sieber, K. Engelborghs, T. Luzyanina, G. Samaey, and D. Roose. DDEBIFTOOL v. 3.1 manual—bifurcation analysis of delay differential equations. Technical report, http://arxiv.org/abs/ 1406.7144, 2015.
 J. M. S. Solorio, D. W. Sukow, D. R. Hicks, and A Gavrielides. Bifurcations in a semiconductor laser subject to delayed incoherent feedback. Optics communications, 214(1):327–334, 2002.
 J. R. Tredicce. Excitability in laser systems: The experimental side. AIP Conference Proceedings, 548(1):238–259, 2000.
 G. H. M. van Tartwijk and M. San Miguel. Optical feedback on selfpulsating semiconductor lasers. Quantum Electronics, IEEE Journal of, 32(7):1191–1202, 1996.
 A. G. Vladimirov and D. Turaev. Model for passive mode locking in semiconductor lasers. Physical Review A, 72(3), 2005.
 A. G. Vladimirov, D. Turaev, and G. Kozyreff. Delay differential equations for modelocked semiconductor lasers. Optics Letters, 29(11):1221–1223, 2004.
 C. E. Webb and J. D. C. Jones. Handbook of Laser Technology and Applications: Laser design and laser systems. CRC Press, 2004.
 S. Wieczorek, B. Krauskopf, and D. Lenstra. A unifying view of bifurcations in a semiconductor laser subject to optical injection. Optics communications, 172(1):279–295, 1999.
 H. G. Winful, Y. C. Chen, and J. M. Liu. Frequency locking, quasiperiodicity, and chaos in modulated selfpulsing semiconductor lasers. Applied physics letters, 48(10):616–618, 1986.
 M. Yamada. A theoretical analysis of selfsustained pulsation phenomena in narrowstripe semiconductor lasers. Quantum Electronics, IEEE Journal of, 29(5):1330–1336, 1993.
 S. Yanchuk and P. Perlikowski. Delay and periodicity. Physical Review E, 79(4):046221, 2009.