Soizic Terrien
email: soizic.terrien@univ-lemans.fr

Bernd Krauskopf

Neil G R Broderick

Venkata A Pammi

Rémy Braive

Isabelle Sagnes

Grégoire Beaudoin

Konstantinos Pantzas

Sylvain Barbay

MERGING AND DISCONNECTING RESONANCE TONGUES IN A PULSING EXCITABLE MICROLASER WITH DELAYED OPTICAL FEEDBACK PREPRINT

Excitability, encountered in numerous fields from biology to neurosciences and optics, is a general phenomenon characterized by an all-or-none response of a system to an external perturbation of a given strength. When subject to delayed feedback, excitable systems can sustain multistable pulsing regimes, which are either regular or irregular time sequences of pulses reappearing every delay time. Here, we investigate an excitable microlaser subject to delayed optical feedback and study the emergence of complex pulsing dynamics, including periodic, quasiperiodic and irregular pulsing regimes. This work is motivated by experimental observations showing these different types of pulsing dynamics. A suitable mathematical model, written as a system of delay differential equations, is investigated through an in-depth bifurcation analysis. We demonstrate that resonance tongues play a key role in the emergence of complex dynamics, including non-equidistant periodic pulsing solutions and chaotic pulsing. The structure of resonance tongues is shown to depend very sensitively on the pump parameter. Successive saddle transitions of bounding saddle-node bifurcations constitute

Introduction

Excitability refers to the spiked or pulsed response of a system at rest to an external perturbation when the perturbation amplitude exceeds the so-called excitable threshold, while no response occurs for smaller perturbations [START_REF] Izhikevich | Neural excitability, spiking and bursting[END_REF]. This general phenomenon has been described in a variety of scientific fields from biology [START_REF] Kca Wedgwood | Robust spike timing in an excitable cell with delayed feedback[END_REF] to neurosciences [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF] and optics [START_REF] Turconi | Control of excitable pulses in an injection-locked semiconductor laser[END_REF], and typically results from the interplay between different internal timescales in the system. After a (short) excitable pulse has been triggered, the system enters a refractory period, during which it is either impossible or much more difficult to trigger another excitable response [START_REF] Selmi | Relative refractory period in an excitable semiconductor laser[END_REF]; note that the refractory period is significantly larger than the duration of the pulse itself. In the presence of delayed feedback, an excitable system can regenerate its own output: after the first excited pulse, the output is reinjected after the feedback delay time τ which triggers the next pulse. As the process repeats, this results in a periodic pulsing regime whose period is directly related and close to the delay τ . This general mechanism for selfpulsations has been demonstrated in a variety of optical systems [START_REF] Garbin | Topological solitons as addressable phase bits in a driven laser[END_REF][START_REF] Romeira | Regenerative memory in time-delayed neuromorphic photonic resonators[END_REF][START_REF] Terrien | Asymmetric noise sensitivity of pulse trains in an excitable microlaser with delayed optical feedback[END_REF], as well as in an excitable biological cell [START_REF] Kca Wedgwood | Robust spike timing in an excitable cell with delayed feedback[END_REF].

From a more general point of view, the introduction of a delay to an excitable system can induce a wealth of complex dynamics beyond regular self-pulsing [START_REF] Garbin | Interactions and collisions of topological solitons in a semiconductor laser with optical injection and feedback[END_REF][START_REF] Terrien | Bifurcation analysis of the yamada model for a pulsing semiconductor laser with saturable absorber and delayed optical feedback[END_REF][START_REF] Thomas | Influence of time-delayed feedback on the dynamics of temporal localized structures in passively mode-locked semiconductor lasers[END_REF]. This includes a high degree of multistability [START_REF] Yanchuk | Delay and periodicity[END_REF] and an enhanced dynamical complexity, such as, for example, quasiperiodic [START_REF] Munsberg | Topological localized states in the time delayed adler model: Bifurcation analysis and interaction law[END_REF] and chaotic regimes. In particular, recent experimental and numerical investigations have demonstrated that an excitable microlaser with delayed feedback can sustain multistable periodic pulsing regimes [START_REF] Terrien | Pulse train interaction and control in a microcavity laser with delayed optical feedback[END_REF][START_REF] Terrien | Equalization of pulse timings in an excitable microlaser system with delay[END_REF]. Depending on the ratio between the internal timescales of the excitable microlaser and on the delay time of the feedback loop or external cavity, these include pulsing patterns with equidistant pulses of equal amplitude or with non-equidistant pulses of different amplitudes [START_REF] Terrien | Pulse-timing symmetry breaking in an excitable optical system with delay[END_REF]. The emergence of equidistant pulsing patterns is well understood [START_REF] Yanchuk | Delay and periodicity[END_REF][START_REF] Krauskopf | Bifurcation Study of a Semiconductor Laser with Saturable Absorber and Delayed Optical Feedback[END_REF][START_REF] Terrien | Equalization of pulse timings in an excitable microlaser system with delay[END_REF]. The emergence of non-equidistant pulsing regimes, on the other hand, has been suggested to originate in resonance phenomena associated with locked periodic orbits on stable tori [START_REF] Terrien | Pulse-timing symmetry breaking in an excitable optical system with delay[END_REF] -yet this was still to be investigated, which motivated the work presented here.

We adopt a dynamical systems point of view to investigate experimentally and numerically the emergence of multi-frequency dynamics in an optical realisation of an excitable system with delayed feedback. We consider an excitable microlaser subject to delayed optical feedback [START_REF] Selmi | Relative refractory period in an excitable semiconductor laser[END_REF][START_REF] Elsass | Control of cavity solitons and dynamical states in a monolithic vertical cavity laser with saturable absorber[END_REF], whose study is motivated by potential applications to neuromorphic photonic computing [START_REF] Va Pammi | Photonic computing with single and coupled spiking micropillar lasers[END_REF][START_REF] Shastri | Photonics for artificial intelligence and neuromorphic computing[END_REF]. Its feedback-induced dynamics is investigated both experimentally and with a suitable mathematical model -the Yamada equations with delayed feebback [START_REF] Yamada | A theoretical analysis of self-sustained pulsation phenomena in narrow-stripe semiconductor lasers[END_REF][START_REF] Krauskopf | Bifurcation Study of a Semiconductor Laser with Saturable Absorber and Delayed Optical Feedback[END_REF], which take the form of three coupled delay differential equations (DDEs) with two slow and one fast variables. Compared to ordinary differential equations, solving DDEs requires specific numerical methods due to their infinite-dimensional nature [START_REF] Roose | Continuation and bifurcation analysis of delay differential equations[END_REF]. We use the Matlab-based numerical continuation software DDE-Biftool [START_REF] Engelborghs | Ddebiftool: a matlab package for bifurcation analysis of delay differential equations[END_REF][START_REF] Engelborghs | Numerical bifurcation analysis of delay differential equations using dde-biftool[END_REF][START_REF] Sieber | Dde-biftool manual-bifurcation analysis of delay differential equations[END_REF] to perform an indepth bifurcation analysis in three parameters of practical importance: the feedback delay τ , the feedback strength κ and the pump parameter A. The results presented here unveil how very large, experimentally accessible locking regions emerge in the parameter space, leading to an increased and observable degree of multistability. We show that this phenomenon involves several transitions in the (τ, κ)-plane that change the structure of regions of locked dynamics. These occur in such a very small A-interval that the associated switch in observed pulsing may seem instantaneous from an experimental perspective. From the mathematical point of view, on the other hand, we are able to distinguish and identify these transitions. In particular, we find two generic cases of resonance tongues interacting locally at an extremum of the rotation number on a torus bifurcation curve. These lead to the disconnection and disappearance of resonance tongues in a parameter plane, respectively, as a third parameter is changed. As we also show, extrema of the rotation number emerge naturally, including the Yamada model with delayed feedback, when two torus bifurcation curves reconnect differently at a saddle transition.

The article is organised as follows. The experimental device is described in Section 1.1, and experimental observations are presented and discussed. Background on the mathematical model is provided in Section 1.2. In Section 2, multi-frequency dynamics, including non-equidistant periodic pulsing regimes, quasiperiodic regimes and chaos are investigated through time-domain simulations of the model. In Section 3, a bifurcation analysis demonstrates that resonance tongues in the (τ, κ)-plane play a key role in the emergence of the observed multi-frequency dynamics. The sensitivity of the structure of resonance tongues to the experimentally relevant pump parameter A is investigated in detail in Section 3.2, and the phenomenon of disconnecting and disappearing resonance tongues is the focus of Section 3.3. We draw some conclusions in Section 4.

Experimental device and observed pulsing regimes

In our optical realisation we use a semiconductor micropillar laser with integrated saturable absorber, consisting of a 5 µm diameter pillar laser with an original design [START_REF] Selmi | Relative refractory period in an excitable semiconductor laser[END_REF][START_REF] Elsass | Control of cavity solitons and dynamical states in a monolithic vertical cavity laser with saturable absorber[END_REF][START_REF] Barbay | Excitability in a semiconductor laser with saturable absorber[END_REF].

In particular, its microcavity includes both a gain and a saturable absorber section. This microlaser emits light perpendicularly to its surface at the cavity resonance wavelength of 980 nm and is optically pumped at around 800 nm. The pump intensity is set just below the self-pulsing threshold such that, in the absence of feedback, the microlaser is in the excitable regime [START_REF] Selmi | Relative refractory period in an excitable semiconductor laser[END_REF][START_REF] Selmi | Spike latency and response properties of an excitable micropillar laser[END_REF][START_REF] Erneux | Two distinct excitable responses for a laser with a saturable absorber[END_REF]. Short optical pulse perturbations (of 80 ps duration) can be sent by an external Ti:Sa mode-locked laser to trigger excitable responses that consist of optical pulses of approximately 200 ps duration.

To realise an external cavity providing the delayed feedback, part of the signal emitted by the microlaser is trans-mitted through a beamsplitter (R/T=70/30) and reflected back by a distant mirror after passing through a 5 cm focal length lens. The resulting optical feedback delay time τ can be set to between ∼ 5 to 10 ns. The light reflected by the beamsplitter is analyzed using an avalanche photodiode, amplified with a high bandwidth RF-amplifier and recorded with an oscilloscope. The microlaser is temperature controlled close to room temperature thanks to a Peltier cooler.

In the presence of feedback, it has been shown that a first excitable pulse can regenerate itself when reinjected in the micropillar after the delay τ , provided that the feedback strength is sufficiently large [START_REF] Terrien | Equalization of pulse timings in an excitable microlaser system with delay[END_REF]. A periodic pulsing regime results with fundamental repetition frequency close to τ -1 ; this regime can coexist with harmonic pulsing solutions with several regularly timed pulses in the feedback cavity. It has been shown that individual pulse trains can be added or erased by single external optical perturbations [START_REF] Terrien | Pulse train interaction and control in a microcavity laser with delayed optical feedback[END_REF][START_REF] Terrien | Equalization of pulse timings in an excitable microlaser system with delay[END_REF]. Under certain experimental conditions, non-regular pulse trains can also be emitted following a seemingly pulse-timing symmetry-breaking phenomenon [START_REF] Terrien | Pulse-timing symmetry breaking in an excitable optical system with delay[END_REF].

Figure 1 illustrates the diversity of pulsing dynamics observed experimentally. The left column represents time series of the measured intensity I and the right column shows their pseudo-space representation [START_REF] Ft Arecchi | Two-dimensional representation of a delayed dynamical system[END_REF]: temporal traces are folded at (approximately) multiples of the delay time τ and stacked vertically. In this representation, the x-axis represents the delay line (realized by the external feedback cavity) and the y-axis the recorded number of roundtrips in the external feedback cavity. Figure 1(a) shows a periodic regime with two equidistant pulses per feedback roundtrip, triggered by an appropriate sequence of two successive external perturbations [START_REF] Terrien | Equalization of pulse timings in an excitable microlaser system with delay[END_REF]. The amplitude of the pulses is quite irregular, while the interpulse timing repeats consistently roundtrip after roundtrip. Such laser intensity fluctuations can be attributed to pump noise and also to detection noise since the microlaser has a very low output power [START_REF] Terrien | Asymmetric noise sensitivity of pulse trains in an excitable microlaser with delayed optical feedback[END_REF]. In Figure 1(b), the system is in the symmetry-broken regime [START_REF] Terrien | Pulse-timing symmetry breaking in an excitable optical system with delay[END_REF]: it sustains two pulses per roundtrip with unequal but well-defined interpulse timings. In the pseudo-space representation, it appears as two non-equidistant pulses in the external feedback cavity. Figure 1(c) illustrates a regime reminiscent of quasiperiodic dynamics. In particular, it displays a strong modulation of the pulse amplitude. Note that the pulse timing is also affected, as the pseudo-space representation shows: the group of pulses shows a negative shift (towards the left) when the pulse amplitude increases. This is explained by the strong amplitude-timing coupling in excitable systems: the response time to an external perturbation gets shorter when the perturbation amplitude is increased [START_REF] Terrien | Pulse train interaction and control in a microcavity laser with delayed optical feedback[END_REF]. Importantly, all the regimes shown in panels (a-c) are multistable and coexist with the off-state (non-lasing equilibrium) of the laser. Finally, Figure 1(d) illustrates a complex pulsing regime reminiscent of chaotic dynamics. It is recorded for a value of the pump parameter slightly above the first lasing threshold at which the off-state of the laser loses stability. As such, this regime is not triggered by external pulse perturbations, but rather by noise. The pseudo-space representation clearly shows multiple competing pulses in the feedback cavity. This complex regime cannot be observed for a long period of time because the microlaser heats up, which results in it switching off.

Overall, Figure 1 which possibly display slightly different internal parameters (including different recombination rates of carriers in the gain and absorber sections) due to fabrication.

Background on the Yamada rate equations with delay

The dynamics of the experimental system under consideration is investigated by using a straightforward extension of the original Yamada rate equation model [START_REF] Yamada | A theoretical analysis of self-sustained pulsation phenomena in narrow-stripe semiconductor lasers[END_REF], a well known system of three ordinary differential equations (ODEs) for single-mode, Q-switched lasers. It has been studied extensively, in particular, through a complete numerical bifurcation analysis which highlighted all its possible dynamics [START_REF] Dubbeldam | Self-pulsations of lasers with saturable absorber : dynamics and bifurcations[END_REF][START_REF] Otupiri | The yamada model for a self-pulsing laser: Bifurcation structure for nonidentical decay times of gain and absorber[END_REF]. This showed that this model exhibits an excitable regime for a large range of parameters, below the lasing threshold at which the non-lasing (i.e., zero-intensity) equilibrium becomes unstable. We consider here the Yamada model with an additional delayed optical feedback term [START_REF] Krauskopf | Bifurcation Study of a Semiconductor Laser with Saturable Absorber and Delayed Optical Feedback[END_REF]. This model has been shown to produce a wealth of pulsing regimes [START_REF] Krauskopf | Bifurcation Study of a Semiconductor Laser with Saturable Absorber and Delayed Optical Feedback[END_REF][START_REF] Terrien | Bifurcation analysis of the yamada model for a pulsing semiconductor laser with saturable absorber and delayed optical feedback[END_REF] and to describe accurately a range of dynamics observed experimentally in the micropillar laser with integrated saturable absorber and delayed optical feedback considered in this article [START_REF] Terrien | Asymmetric noise sensitivity of pulse trains in an excitable microlaser with delayed optical feedback[END_REF][START_REF] Terrien | Pulse train interaction and control in a microcavity laser with delayed optical feedback[END_REF][START_REF] Terrien | Equalization of pulse timings in an excitable microlaser system with delay[END_REF][START_REF] Terrien | Pulse-timing symmetry breaking in an excitable optical system with delay[END_REF]. This includes a variety of stable pulsing periodic regimes with different numbers of equidistant and non-equidistant pulses in the feedback cavity. Overall, the model is written as a system of three coupled DDEs for the dimensionless gain G, absorption Q and intensity I:

Ġ = γ G (A -G -GI); Q = γ G σ (B -Q -aσQI); İ = (G -Q -1)I + κI(t -τ ). (1)
Here, A is the pump parameter, B describes the linear absorption, a is the saturation parameter, γ G is the recombination rate of carriers in the gain section, and σ is the ratio between the recombination rate of carriers in the gain and absorber sections; furthermore, time is rescaled in [START_REF] Ft Arecchi | Two-dimensional representation of a delayed dynamical system[END_REF] to the photon lifetime in the cavity (which is on the order of 1-2 ps). Importantly, γ G is usually small in semiconductor lasers and σ is between 0.5 and 2 for the microlaser we consider. As such, system (1) is a slow-fast dynamical system with two slow variables G and Q, and one fast variable I. In the intensity equation, the delayed term describes the incoherent delayed optical feedback with feedback strength κ and delay time τ . The influence of the feedback parameters κ and τ has been investigated through an extensive bifurcation analysis [START_REF] Krauskopf | Bifurcation Study of a Semiconductor Laser with Saturable Absorber and Delayed Optical Feedback[END_REF][START_REF] Terrien | Bifurcation analysis of the yamada model for a pulsing semiconductor laser with saturable absorber and delayed optical feedback[END_REF]. This highlighted, in particular, an important and increasing level of multistability as the delay τ is increased, with a large number of coexisting stable periodic pulsing solutions.

The parameters describing material properties of the laser are fixed here throughout to B = 2, γ G = 0.01, a = 5.5, and σ = 1.8; these values are chosen to match the parameters considered in previous work [START_REF] Terrien | Pulse-timing symmetry breaking in an excitable optical system with delay[END_REF]. The pump parameter A, the feedback strength κ and the feedback delay τ , on the other hand, can be changed during the experiment and are considered as bifurcation parameters. Importantly, for all the parameters combinations considered in this article, the solitary laser for κ = 0 (i.e., the model without the feedback term) is in the excitable regime: the non-lasing equilibrium, corresponding to the laser off-state, is stable but the system can release a single intensity pulse when subject to an external perturbation with sufficiently large amplitude [START_REF] Dubbeldam | Self-pulsations of lasers with saturable absorber : dynamics and bifurcations[END_REF].

Multi-frequency dynamics in time-domain simulations

The experimentally observed regimes in Figure 1 can be identified in time-domain simulations of the Yamada model with feedback (1). Figure 2 illustrates such simulations with time series of the intensity I for fixed feedback delay τ = 335 and pump parameter A = 2.7, and the four values of the feedback strength κ = 0.007, κ = 0.0159, κ = 0.04 and κ = 0.1. We checked that, for the considered value of A, the solitary laser (for κ = 0) is indeed in the excitable regime. Moreover, in the presence of feedback, the nonlasing (i.e., zero-intensity) equilibrium solution of ( 1) is still stable for all the values of κ considered in Figure 2: hence, the intensity I remains zero in the absence of external perturbations. In the simulations, an initial external perturbation is accounted for by setting initial conditions (given for a DDE by a history segment over [-τ, 0]) with suitable non-zero intensity. Specifically, these are set to the (unstable) equilibrium solution of ( 1) with a non-zero intensity I, which corresponds to the continuous-wave regime of the laser [START_REF] Terrien | Bifurcation analysis of the yamada model for a pulsing semiconductor laser with saturable absorber and delayed optical feedback[END_REF].

The time series in Figure 2(a1-d1) are shown after a few dozens or hundreds of roundtrips to ensure that any transient dynamics has died down. The displayed dynamical regimes thus correspond to stable pulsing regimes of (1). Panels 2(a2-d2) show the corresponding pseudo-space representation of these time series [START_REF] Ft Arecchi | Two-dimensional representation of a delayed dynamical system[END_REF], as explained above. These simulation results illustrate the diversity of stable pulsing regimes observed over a small range of the feedback strength κ.

Figure 2(a) for κ = 0.007 shows a periodic pulsing regime with four equidistant pulses in the external feedback cavity, that is, over the span of one delay time τ that constitutes the feedback loop. Note that the amplitude also repeats exactly here because (1) does not feature (pump or other) noise. The positive slope of the intensity pulse trains in the pseudo-space representation in panel (a2) shows that the period of pulsing is slightly larger than τ /4, which is due to the latency time of the system to the re-injected perturbation [START_REF] Yanchuk | Delay and periodicity[END_REF][START_REF] Terrien | Bifurcation analysis of the yamada model for a pulsing semiconductor laser with saturable absorber and delayed optical feedback[END_REF][START_REF] Krauskopf | Bifurcation Study of a Semiconductor Laser with Saturable Absorber and Delayed Optical Feedback[END_REF].

Figure 2(b-c) for increasing values of κ (where all the other parameters are fixed) show examples of dynamics on a torus, which may be quasiperiodic or locked to an attracting periodic orbit. Locked periodic solutions are found inside resonance regions or resonance tongues, which are regions of a parameter plane that are bounded by curves of saddle-node bifurcations of periodic orbits. Resonance tongues emerge from resonance points, which are points along a torus (or Neimark-Sacker) bifurcation curve where the rotation number is rational [START_REF] Ya Kuznetsov | Elements of applied bifurcation theory[END_REF][START_REF] Keane | Delayed feedback versus seasonal forcing: resonance phenomena in an el niño southern oscillation model[END_REF]. Quasiperiodic solutions, on the other hand, are found in the parameter plane along curves 'in between' infinitely many and generally very narrow resonance tongues. Figure 2(b) for κ = 0.0159 is past a torus bifurcation and shows a locked periodic orbit on a torus. Note that the period of this periodic regime is not directly related to an integer submultiple of the feedback delay time τ . Rather, a single period of the periodic solution displays several intensity pulses with different amplitudes, which is characteristic of a locked periodic orbit on a stable invariant torus; see panel (b1). Figure 2(c) for κ = 0.04 shows a quasiperiodic regime (or a periodic regime with a very large period); in particular, the time series does now not repeat exactly but displays roughly five pulses in the feedback loop, with a deeply modulated pulse amplitude. The pseudo-space representation in panel (c2) highlights the strong amplitude-time coupling of the pulses [START_REF] Terrien | Pulse train interaction and control in a microcavity laser with delayed optical feedback[END_REF]. The fact that a stable quasiperiodic regime is observed in panels (c) for a nearby value, yet larger value of κ support the interpretation that panels (b) indeed show a locked periodic solution on a stable torus.

Finally, Figure 2(d) for the yet larger value of κ = 0.1 shows an example of a chaotic regime. Here no clear structure is observed in neither the repetition rate nor the pulse amplitudes, which is visually clear especially in panel (d2); see also panel (d1) and the enlargement of the intensity time series over two roundtrips in the inset; note, in particular, that the modulation of the amplitude is no longer periodic.

Transitions between the different stable pulsing regimes in Figure 2 as the feedback strength κ is changed are explored further in Figure 3. Panels (a) and (b) show one-parameter bifurcation diagrams for increasing and decreasing values of κ, respectively, where observed dynamics is represented by the pulses amplitudes observed over one single roundtrip of duration τ . These diagrams are obtained from sweeped simulation as follows. For the smallest values of κ, the initial condition/history is set to the (unstable) lasing equilibrium solution; for each subsequent value of κ, the previously calculated solution (i.e., for a slightly smaller or larger value of κ in Figure 3(a) and (b), respectively) is considered as initial history. In all cases, the simulation is run over several hundreds of roundtrips so that any transient phenomenon are disregarded, before the pulses amplitudes over one single roundtrip are recorded. In each case, we show the attractor in (G, Q, I) space and in a Poincaré section represented by the (G, Q)-plane. It is important to note here that, due to the infinite-dimensional nature of DDEs, these are both projections: onto the threedimensional physical space of the variables G, Q and I for the trajectory and, similarly, onto the (G, Q)-plane for the Poincaré map, which is defined by the fixed value I eq of the intensity at the unstable equilibrium point [START_REF] Krauskopf | Bifurcation Study of a Semiconductor Laser with Saturable Absorber and Delayed Optical Feedback[END_REF]; specifically, I eq = 1.383 in panel (c1), I eq = 1.412 in panel (c2), I eq = 1.444 in panel (c3), I eq = 1.572 in panel (c4), and I eq = 1.796 in panels (c5) and (d). We also show the corresponding radio-frequency (RF) spectrum (the power spectral density) of the real-valued intensity time series.

As κ is varied, several transitions are observed between the pulsing regimes displayed in Figure 2. For the smallest values of κ in Figure 3(a), a single value is observed for the pulse amplitudes. This corresponds to a periodic regime with a fixed number of equidistant pulses of equal amplitude per roundtrip [START_REF] Terrien | Pulse-timing symmetry breaking in an excitable optical system with delay[END_REF]. This periodic regime is illustrated in panel (c1) for κ = 0.002, where one observes a periodic orbit in (G, Q, I) space and two corresponding points in the (G, Q)-plane representing the Poincaré section; the RF spectrum is discrete, featuring a main peak with large harmonics owing to the strongly pulse-like nature of the oscillation. Figure 3(a) shows that increasing κ leads to a tran- sition from this periodic to a regime with an increasingly stronger modulation of the pulse amplitude. This is clear evidence of a torus bifurcation with subsequent quasiperiodic or locked dynamics. Panels (c2) for κ = 0.013 show the corresponding stable torus in projection onto (G, Q, I)space; its intersection with the Poincaré section is a densely filling closed curve in the (G, Q)-plane, and its RF spectrum displays a large number of peaks at incommensurate frequencies; compare with Figure 2(b). Hence, the dynamics on the torus is quasiperiodic (or of very high period).

(c1) (c2) (c3) (c4) (c5) (b) (d) I G Q (c1) G 0.6 1.8 Q 0 0.4 ν 0 0.1 0 I(ν) 7e4 I G Q (c2) G 1.15 1.5 Q 0.1 0.4 ν 0 0.1 0 I(ν) 4e4 I G Q (c3) G 1.12 1.6 Q 0.1 0.4 ν 0 0.1 0 I(ν) 5e4 I G Q (c4) G 0.6 1.8 Q 0 0.7 ν 0 0.1 0 I(ν) 3e4 I G Q (c5) G 0.6 1.8 Q 0 0.7 ν 0 0.1 0
I(ν) 4e4 I G Q G 0.6 1.8 Q 0 0.7 ν 0 0.1 0 I(ν) 12e4 (d)
When κ is increased in Figure 3(a), different (usually very narrow) resonance tongues are briefly crossed, where the dynamics is locked to a periodic orbit on the stable torus; see also Figure 2(c). Figure 3(c3) shows such a locked periodic regime for κ = 0.025, which lies in a larger resonance tongue corresponding to a wider range of κ. The trajectory in (G, Q, I)-space shows a closed trajectory, which does not cover the entire stable torus but rather winds six times around it. The Poincaré section in the (G, Q)-plane shows the corresponding six intersection points, which confirms that the system is in a 1:6 locking region. When κ is increased further, a transition is observed in Figure 3(a) from quasiperiodic (or high-period) dynamics on a torus to a chaotic regime; this is detected by a sudden growth of the maximal recorded amplitude. The chaotic regime is illustrated further in Figure 3(c4-c5); compare with the chaotic intensity time series in Figure 2(d). Here, neither the (projected) trajectory nor the Poincaré section display a clear structure; moreover, the RF spectrum is characteristic of a chaotic regime: although some clear peaks are visible, there are multiple frequency bands in which no clear structure can be seen. This chaotic regime is observed over a large range of κ in Figure 3(a) before the system jumps back to a stable periodic pulsing regime as κ is increased further.

Starting from this periodic solution, Figure 3 As can be seen from its RF spectrum, this periodic regime features a different number of equidistant pulses in the feedback cavity compared to the periodic regime in panel (c1). It is important to note that the non-lasing equilibrium solution of ( 1) is stable as well over the entire range of κ considered in Figure 3(a,b). Which regime is observed in practice in the presence of the overall multistability depends on the considered initial condition/history [START_REF] Terrien | Equalization of pulse timings in an excitable microlaser system with delay[END_REF].

3 Bifurcation analysis in the (τ, κ)-plane

The complex dynamics highlighted in Figures 2 and3 for the continuation of families of equilibrium and periodic solutions of delay-differential equations, as well as their bifurcations. The results are presented in the form of bifurcation diagrams consisting of bifurcation sets in the (τ, κ)-plane, formed by curves of codimension-one bifurcations that bound regions of different dynamics.

Bifurcation diagram for A = 2.7

Figure 4 shows the bifurcation diagram in the (τ ,κ)-plane for the value A = 2.7 of the pump parameter considered in Figures 2 and3; its panel (a) shows the range of τ from 0 to 550, and panel (b) is an enlargement that focuses on the resonance structures we consider here. A main feature of Figure 4 is a single curve H of Hopf bifurcation of the lasing (i.e., non-zero intensity) equilibrium of (1): when this curve is crossed, for instance increasing τ at constant κ, a small-amplitude periodic solution emerges. More precisely, the bifurcating periodic solution is locally stable when the Hopf bifurcation is supercritical, indicated by bold parts of the curve H, and it is locally unstable along thin parts of H where the Hopf bifurcation is subcritical. The criticality of H changes at codimension-two degenerate Hopf bifurcation points DH, indicated by blue dots in Figure 4. We also find codimension-two Hopf-Hopf bifurcation points HH in Figure 4; they arise here due to self-intersections of the curve H and also lead to changes of criticality of the Hopf bifurcation. When τ is increased for a (sufficiently large) fixed value of κ, several stable periodic solutions emerge successively from the supercritical parts of the curve H, and their periods are close to submultiples of the delay τ [START_REF] Krauskopf | Bifurcation Study of a Semiconductor Laser with Saturable Absorber and Delayed Optical Feedback[END_REF][START_REF] Terrien | Bifurcation analysis of the yamada model for a pulsing semiconductor laser with saturable absorber and delayed optical feedback[END_REF]. Past the Hopf bifurcations where they emerge, these periodic solutions correspond to pulsing patterns with different numbers of equidistant pulses in the feedback loop [START_REF] Terrien | Bifurcation analysis of the yamada model for a pulsing semiconductor laser with saturable absorber and delayed optical feedback[END_REF]; these numbers are shown in boxes in Figure 4(a) for the respective supercritical parts of the curve H. The increasing level of multistability with the delay τ is typical for delay systems [START_REF] Yanchuk | Delay and periodicity[END_REF], and has been discussed in the literature for this particular system [START_REF] Ruschel | The limits of sustained self-excitation and stable periodic pulse trains in the yamada model with delayed optical feedback[END_REF].

The codimension-two points DH and HH in Figure 4 also give rise to additional bifurcation curves in the (τ ,κ)-plane. From each degenerate Hopf bifurcation point DH, where the Lyapunov coefficient (third-order normal form coefficient) vanishes, emerges a curve S of saddle-node bifurcation of periodic orbits [START_REF] Krauskopf | Bifurcation Study of a Semiconductor Laser with Saturable Absorber and Delayed Optical Feedback[END_REF][START_REF] Terrien | Bifurcation analysis of the yamada model for a pulsing semiconductor laser with saturable absorber and delayed optical feedback[END_REF]; for sake of keeping this exposition focused, Figure 4 only shows the curve S that emerges from a point DH for low κ near τ = 300, which we require for the discussion in Section 3.2. Moreover, at each Hopf-Hopf bifurcation point HH, where the linearization of system (1) at the bifurcating lasing equilibrium has two pairs of complex conjugate eigenvalues on the imaginary axis, two curves of torus bifurcation T can typically emerge [START_REF] Ya Kuznetsov | Elements of applied bifurcation theory[END_REF]. For sake of clarity, only two pairs of curves T are shown in Figure 4: those emerging from the two points HH for low values of κ around τ = 335. Indeed, many Hopf-Hopf bifurcation points are encountered as the delay τ is increased, resulting in a very complex bifurcation diagram with many more curves of torus bifurcation; this has been discussed in the literature [START_REF] Krauskopf | Bifurcation Study of a Semiconductor Laser with Saturable Absorber and Delayed Optical Feedback[END_REF][START_REF] Terrien | Bifurcation analysis of the yamada model for a pulsing semiconductor laser with saturable absorber and delayed optical feedback[END_REF][START_REF] Ruschel | The limits of sustained self-excitation and stable periodic pulse trains in the yamada model with delayed optical feedback[END_REF] and is beyond the scope of this article.

Connecting resonance tongues

We now focus on the two shown pairs of torus bifurcation curves T. They are shown in Figure 4(b) in the enlarged parameter area of interest, where they are clearly seen to emerge from the two respective Hopf-Hopf bifurcation points; see also the further enlargement in the inset. From each torus bifurcation curve T bifurcates a smooth invariant torus on which the (multi-frequency) dynamics is either quasiperiodic or locked. Resonance tongues in the (τ ,κ)plane emerge from p:q resonance points, some of which are highlighted on the curves T by black dots: at these points, the pair of critical complex conjugate Floquet multipliers on the unit circle e ±2πα has the rational rotation number α = p q , where p and q are relatively prime [START_REF] Ya Kuznetsov | Elements of applied bifurcation theory[END_REF]. A pair of curves S of saddle-node bifurcation of periodic orbits emerge from each p:q resonance point; they form the boundaries of a p : q resonance tongue in which the dynamics is p:q locked on the torus, that is, periodic rather than quasiperiodic; see Figure 3(c3) for an example of such locked dynamics. Figure 4(b) also shows curves S bounding the resonance tongues of some lower-order resonances; we remark that the pairs of curves S bounding narrow resonance tongues are sometimes so close to each other that they appear as single curves in the parameter plane. We find that these pairs of curves S connect different resonance points on the pair of torus bifurcation curves T emerging from the respective Hopf-Hopf bifurcation point HH; specifically a given p:q resonance point on one of the curve T connects to a p:(p + q) resonance point of the other curve T. To the best of our knowledge, such connecting resonance tongues have been observed only once before, namely in a constructed system of DDEs with two state-dependent delays [START_REF] Calleja | Resonance phenomena in a scalar delay differential equation with two state-dependent delays[END_REF]. The Yamada model with a single feedback term and with a constant (non-state-dependent) delay, hence, constitutes the first physically relevant and yet also simplest mathematical model featuring connecting resonance tongues near a Hopf-Hopf bifurcation point. Our results and those in Ref. [[4]] suggest that this phenomenon is generic and related to the occurence of chaotic dynamics near Hopf-Hopf bifurcation. The argument is as follows. Close to the resonance points, the dynamics is either p:q or p:(p + q) locked on the smooth invariant torus that emerges from the respective torus bifurcation. The fact that the bounding saddle-node bifurcation curves S connect means that the p:q locked periodic solutions smoothly evolves into p:(p + q) locked solutions. This is perfectly possible in a phase space of dimension at least three, but it cannot happen while these periodic orbits lie on a smooth invariant two-dimensional torus throughout (there cannot be two periodic orbits of different winding number on one and the same two-dimensional torus). Therefore, the smooth invariant torus bifurcating from the respective torus bifurcation T necessarily breaks up at some point when going through the connecting resonance tongue. Such torus break-up can typically lead to the creation of chaotic dynamics through a quasiperiodic route to chaos [START_REF] Dg Aronson | Bifurcations from an invariant circle for two-parameter families of maps of the plane: a computer-assisted study[END_REF][START_REF] Green | Bistability and torus break-up in a semiconductor laser with phase-conjugate feedback[END_REF].

The simulations shown in Figure 3 and, in particular, the appearance of chaotic dynamics can, therefore, be interpreted as follows. For a fixed value of τ and for small κ, a stable periodic solution (corresponding to a fixed number of equidistant pulses in the feedback loop) is observed. When κ is increased, it looses stability through a torus bifurcation, which leads to the emergence of a stable invariant torus with dynamics that is quasiperiodic or p:q locked with very high q; see Figure 3(c2). As κ is increased, the rotation number on the torus evolves and different resonance tongues are crossed, only the larger of which with low q can be identified in Figure 3(a); see Figure 3(c3) for an example of such locked dynamics. As κ is increased further, the torus breaks up and chaotic dynamics is observed over a large range of κ; as is illustrated in Figure 3(c4,c5). Indeed, this overall transition provides an explanation for experimentally observed chaotic pulsing regimes as in Figure 1(d), which are reported here in this microlaser system for the first time. Curves S near this minimum that bound resonance tongues that connect two p:q resonance points are shown in lighter magenta. The inset is an enlargement near the 1:3 point.

Influence of decreasing the pump parameter A

We now investigate how the structure of resonance tongues highlighted in Figure 4(b) evolves when the pump parameter A is decreased down towards a value A = 2. The motivation for this is two-fold. First, in the experiment, the pump parameter is related to the amount of energy provided to the microlaser through optical pumping and, as such, is the main control parameter in the experiment. Second, a recent study has shown that for A = 2 multistable periodic regimes are found, which include non-equidistant pulses in the feedback cavity [START_REF] Terrien | Pulse-timing symmetry breaking in an excitable optical system with delay[END_REF]. Each of these periodic non-equidistant pulsing regimes exists in a very large region of the (τ, κ)-plane of feedback parameters, which is bounded by saddle-node bifurcation curves S of periodic solutions. It has been suggested that these regions are in fact very large resonance tongues, and we now investigate the mechanism for the emergence of such unusally large locking regions. For clarity, we focus on resonance tongues associated with one pair of torus bifurcation curves only, namely the one that emerges from the Hopf-Hopf point HH with the smallest value of τ in Figure 4(b).

Rotation number near Hopf-Hopf point

Figure 5 shows the bifurcation diagram of (1) for A = 2.42 in the relevant region of the (τ, κ)-plane, featuring the Hopf bifurcation curve H with the Hopf-Hopf bifurcation point HH, the associated two torus bifurcation curves T, and resonance tongues arising from selected p:q resonance points; the latter are bounded by pairs of curves S of saddlenode bifurcations of periodic orbits. The lower curve T in Figure 5 has two fold points (with respect to τ ) close to the Hopf-Hopf point HH; moreover, the rotation number α first decreases along this curve away from HH, reaches a minimum of below 1 11 near (but not exactly at) the leftmost sharp fold, and subsequently increases again. This is in contrast to the properties of the corresponding lower torus bifurcation curves T for A = 2.7 (see the leftmost part of Figure 4(b)), which lacks the leftmost sharp fold and along which α increases monotonically from the point HH. The non-monotonicity of α along this curve T in Figure 5 is associated with an intriguing phenomenon: the lighter colored saddle-node curves S connect p:q resonance points on the decreasing part to the left of the minimum with p:q resonance points on the increasing part to the right of the minimum. Beyond HH, on the other hand, we find connections by pairs of curves S between p:q on the upper curve T and p:(p+q) resonance points on the other curve T; this situation is as discussed above for A = 2.7. Another interesting difference with the case A = 2.7 lies in the shape of the 1:4 resonance tongue in Figure 5. Notice that the 1:4 locking region in the (τ, κ)-plane is significantly larger than for the other resonances, and it is bounded by a single curve S; as the inset in Figure 5 shows, this curve passes very close to the 1:3 resonance point, but does not connect to this point.

Disappearance of Hopf-Hopf point

The Hopf-Hopf bifurcation point HH in Figure 5 arises because the curve H has a little loop. Comparison with Figure 4(b) shows that this loop becomes smaller when the pump parameter A is decreased. Figure 6 shows that the loop and the point HH disappear when A is decreased further; moreover, it explains what this means for the associated curves T of torus bifurcation. The rotation number α along the torus bifurcation curves is represented by a color scale in the larger view of panels (a1) to (c1), while panels (a2) to (c2) are further enlargements. The observed configuration in Figure 6(a) for A = 2.39 is topologically as that for A = 2.42 in Figure 5. However, the point HH in Figure 6(a) is now at a lower value of κ and the loop of the curve H is even smaller; similarly, the part of the lower torus curve T along which the rotation number α decreases is significantly smaller than for A = 2.42.

Growth and merging of resonance regions

We now show how the resonance tongues identified thus far grow into large and experimentally relevant regions of different observable locked periodic dynamics. To this end, Figure 7 shows how the bifurcation diagram of (1) in the relevant larger region of the (τ, κ)-plane changes when the pump parameter A is decreased from A = 2.29 in panel (a) to A = 2.23 in panel (d). The focus here is on curves S of saddle-node bifurcation of periodic orbits that bound the (shaded) regions with stable 1:4 dynamics, that is, where one observes periodic pulsing with four (non-equidistant) pulses in the feedback cavity.

Figure 7(a) for A = 2.29 shows that the 1:4 resonance tongue bounded by a single curve S, which was identified in Figure 5, has grown; note that this curve S in Figure 7(a) still connects to and from a 1:4 resonance point on the torus bifurcation curve T associated with the leftmost point HH that just disappeared. Close by, there is another, large 1 : 4 resonance region that is bounded by another curve S and exists for larger values of τ ; we remark that this resonance region does not connect to any resonance point and extends to very large values of τ (as we checked by computing its bounding curve S beyond the shown range). When the pump parameter is changed to A = 2.28 as shown in Figure 7(b), these two resonance regions have merged into a single and very large 1 : 4 locking region emerging from the left-most 1:4 resonance point, which happens to be near the right-most fold on the associated curve T. This change occurs in a saddle transition of the bounding curves S of saddle-node bifurcations of periodic orbits, where the two nearby fold points of the two curves S meet and the curves connect differently. other with decreasing A, and in Figure 7(c) for A = 2.25 they now connect differently; also after a saddle-transition, now of the curves T. This means, in particular, that the two 1 : 4 resonance points now lie on the same curve T. As A is decreased further, these two points meet and then disappear, as is shown in Figure 7(d). Simultaneousy, the two resonance tongues meet, and their bounding curves S connect differently. This leads to a very large resonance region (shaded) of 1:4 locked pulsing, bounded by a curve S, which is not connected anymore to any of the torus bifurcation curves found in the considered parameter range. Inside this large four-pulse region, one finds a smaller region, bounded by a second curve S, where the locked solution does not exist. Notice further that pairs of cusp points have occured in Figure 7(c) and (d) on the upper curve S and in panel (d) on the lower intermediate curve S; these cusp points appear locally in codimension-three swallowtail bifurcations and lead to small regions with additional 1 : 4 locked solutions that are not discussed further in this article.

The changes to the curves T and the bounding curves S of the 1:4 resonance regions occur over a small parameter range of A and are difficult to see on the scale of Figure 7. This is why they are illustrated further in Figure 8 with enlargements near the two 1:4 resonance points involved. Figure 8(a) for A = 2.29 (an enlargement of Figure 7(a)) clearly shows the respective two resonance tongues emerging from the 1:4 resonance points near the folds of the two different torus bifurcation curves T. Panel (b) for A = 2.2605 is just before the saddle transition of these curves T: the two folds are now extremely close together and there are two 1:5 resonance points practically on them. When the curves T are connected differently, just after the saddle transition for A = 2.25 as shown in Figure 8(c) (an enlargement of Figure 7(c)), the two 1:4 resonance points are now clearly on the same curve T. This means, in particular, that there must be an extremum of the rotation number. Different resonance points are indicated along the curves T to give information on the rotation number along them, and this shows that the point m is a minimum near 1 5 . Similarly, there is a maximum M of the rotation number, also near 1 5 , on the other torus bifurcation curve T in panel (c). The value of the rotation number at the minimum m increases as A is decreased; in the process, the two 1 : 4 resonance points move towards each other until they collide and disappear when the minimum moves through 1 4 , as discussed previously. The result is the situation for A = 2.23 in Figure 8(d) (an enlargement of Figure 7(d)), where the bounding curves S of the associated resonance tongues are now connected differently: they lie either side of the point m on T and no longer interact with this curve of torus bifurcation. Note that, throughout in Figure 8, the lower parts of curves T are very close to the curve S that emerges from the point DH near the minimum of the Hopf bifurcation on the left.

Finally, as A is decreased further, the region without locking inside the large 1:4 pulsing region in Figure 7 shrinks and eventually disappears from the shown range of τ through a minimax transition of the bounding curve S (this is not illustrated here). Overall, this leads to a single, very large locking region, which is not connected to any resonance point along a torus bifurcation curve and which extends to very large values of the delay time τ . Similar transitions of other resonance tongues exist as well (but are not shown in Figures 7 and8); this agrees with the observation of associated regions with stable periodic solutions with different numbers of non-equidistant pulses [START_REF] Terrien | Pulse-timing symmetry breaking in an excitable optical system with delay[END_REF]. From a practical point of view, this means that, for these lower values of the pump parameter A, locked periodic solutions corresponding to such non-equidistant pulses in the external feedback cavity can be observed over very large ranges in the (τ, κ)-plane. This agrees with the behaviour observed recently in both the model and an actual experiment [START_REF] Terrien | Pulse-timing symmetry breaking in an excitable optical system with delay[END_REF], which confirms that the non-equidistant pulsing periodic regimes originate in a resonance phenomenon.

While the resulting regions of stable pulsing solutions are large in the (τ, κ)-plane, we stress that the transitions of resonance tongues that generate them happen in a very narrow range of A. Hence, the system is highly sensitive to small changes of the pump parameter. This is of practical importance since the pump parameter A is a main control parameter in the actual experiment: even a small change in A may result in the disappearance of the locked solution (corresponding to non-equidistant pulses in the feedback cavity) and, as a consequence, the appearance of qualitatively different long-term dynamics [START_REF] Terrien | Equalization of pulse timings in an excitable microlaser system with delay[END_REF]. This Figure 9: Sketches of the two scenarios for disconnecting (a) and disappearing (b) p:q resonance tongues near an extremum of the rotation number on a curve T of torus bifurcation. Shown is the case for a minumum m where, from top to bottom, m < p q , m = p q and m > p q ; column (b) also shows a connecting r : s resonance tongue with m < r s .

observed change with A appears to be sudden from an experimental and practical perspective; mathematically, on the other hand, it is explained by the sequence of transitions of resonance tongues presented above.

Disconnecting and disappearing resonance tongues near extrema of the rotation number

While Figures 7 and8 illustrate the case of 1 : 4 resonance, disconnecting resonance tongues are a generic phenomenon that occurs for any pair of p:q resonance points near an extremum of the rotation number on a torus bifurcation curve T. We now discuss in more detail the disappearance of pairs of resonance points and the consequences for the associated resonance tongues. To our knowledge, this generic phenomenon of codimension three (meaning that it happens at a specific point in a three-dimensional parameter space, as the (τ, κ, A)-space of (1) considered here) has not been reported in previous literature. In fact, there are actually two different cases, as is sketched in Figure 9 for the case that the extremum on T is a minimum m in a two-parameter plane. In both columns (a) and (b) of Figure 9 a pair of p:q resonance points with m < p q collides when m increases and reaches the value m = p q as a third parameter is changed, and subsequently does not exist any longer on the curve T when m > p q . As a function of the third parameter (A for (1)), this corresponds to a fold point of p:q resonance points on the curve T. Column (a) shows the situation encountered in Section 3.2.3, where the resonance tongues 'point away' from the minimum m; as a result, the pairs of bounding curves S then meet at m when m = p q and form two curves on either side of the curve T for m > p q . Since the bounding curves S are then no longer connected with T at a point of resonance, we refer to this generic case as that of disconnecting resonance tongues. In contrast, the p:q resonance tongues in Figure 9(b) connect near the minimum m, meaning that the curves S bound a single small resonance region. This p:q resonance region shrinks to a point when m = p q and has disappeared for m > p q ; left are then only r : s resonance tongue with m < r s , of which one is shown in the sketch. We refer to this generic case as that of disappearing resonance tongues. Figure 9 illustrates the codimension-three transitions of disconnecting and disappearing resonance tongues locally for a single pair of resonance points. As the value of m increases smoothly and monotonically with a third parameter, the same local scenario takes place at any point on the torus bifurcation T where the rotation number is rational. Since this set is dense in T, any change of the parameter leads to infinitely many disconnecting or disappearing resonance tongues. In other words, Figure 9 represents a countably infinite sequence of such transitions of resonance tongues as the third parameter is changed and the minimum m increases. We remark that the case of a maximum M on T is completely analogous: here the resonance tongues merge or disappear at rational points when M decreases with a third parameter. With this interpretation, we claim that columns (a) and (b) of Figure 9 represent the unfolding of the two cases of disconnecting and of disappearing resonance tongues.

Proving that these unfoldings are generic is a challenge that is left for future work, especially since it would require the simpler setting of an ODE rather than a DDE as studied here. However, the correctness of these two unfoldings is supported by numerical evidence in the (τ, κ)-plane of (1). Specifically, Figure 10 highlights the local situations before and after the respective transition for the disconnecting resonance tongues in column (a) and for To complete this section, we present how extrema arise on a torus bifurcation curve in a two-parameter plane. As row (a) of Figure 11 shows, the generic mechanism is that of a saddle-transition of a pair of torus bifurcation curves, which necessarily needs to respect and agree with the respective changes of the rotation number along these curves. Note that generically the rotation number is not constant, and Figure 11(a) shows the case that it increases as indicated by the arrows on the curves labeled T. As a third parameter is changed, the saddle transition occurs: here the two curves meet as shown at a singular point (which is a saddle point of the surface of torus bifurcations in the three-dimensional parameter space). The locus of torus bifurcations connects differently past the saddle transition, as is shown. On the level of just the curves T, this is a standard way of reconnecting the respective branches of torus bifurcation. However, we can conclude more here: since the rotation number is unique and common to all branches at the singular point, a minimum m on one curve T and a maximum M on the other curve T are necessarily created in the process. The sketch in Figure 11(a) represents the generic case of the saddle transition of a locus of torus bifurcations. As we discussed in Section 3.2.3, such a saddle transition occurs in the (τ, κ)-plane of (1), namely between the two situations shown in Figure 8(b) and (c). This is illustrated more clearly in Figure 11(b) with enlargements near the saddle transition that also shows the rotation num-ber along the respective curves T. Notice that the situation in panel (b1) is very close to the saddle-transitions, where the curves T connect.

Conclusion

We investigated the emergence of complex multi-frequency dynamics in an excitable system subject a delayed feedback loop -specifically a microlaser with optical feedback from an external mirror. The Yamada model with delayed feedback has been shown to reproduce accurately a range of phenomena observed experimentally in an excitable microlaser subject to delayed optical feedback. This includes quasiperiodic dynamics, locked periodic dynamics on a torus and chaotic dynamics. A bifurcation analysis in the physically relevant variables, the feedback delay τ , the feedback strength κ and the pump parameter A, unveiled the key role played by resonance tongues in the (τ, κ)plane, which arise from torus bifurcation curves that, in turn, emerge from Hopf-Hopf bifurcation points on the main curve of Hopf bifurcations. In particular, the emergence of chaotic dynamics in the DDE results from the mechanism of torus break-up. This is associated with the observation of connecting resonance tongues, which connect p:q resonance points with p:(p + q)resonance points.

Our results confirm that non-equidistant pulsing periodic solutions originate in resonance phenomena. Unexpectedly large locking regions in which such pulsing regimes are observed emerge in the (τ, κ)-plane from initially much smaller resonance tongues in an intriguing sequence of transitions when the pump parameter A is changed. In particular, this involves several re-arrangements of pairs of bifurcation curves in saddle transitions, as well as pairs of p:q resonance points coming together at points on a torus bifurcation curve where the rotation number has an extremum. We presented (conjectural) unfoldings of the two generic cases of disconneting and of disappearing resonance tongues, and also showed how extrema on torus bifurcation curves arise. In the Yamada model with delayed feedback the respective transitions all occur in quick succession, showing that the system is very sensitive to small changes in the pump parameter A. Because A is a main control parameter this is of practical interest in an actual experiment -the bifurcation analysis presented here explains the consistency and genericity of what may otherwise be interpreted as a non-obvious sudden jump from stable single-pulse behavior to pulsing with several equidistant and/or non-equidistant pulses.

From a more general point of view, the DDE model considered here has only two main ingredients: excitability and feedback. As such, our results are expected to be of more general practical interest, beyond the particular laser device considered in this paper. Examples with these same ingredients are other optical systems with different types of feedback, as well as biological systems with delayed feedback, which show considerable similarities, as confirmed by a recent study with an excitable cell [START_REF] Kca Wedgwood | Robust spike timing in an excitable cell with delayed feedback[END_REF].

Figure 1 :

 1 Figure 1: Experimental intensity time series (left, a1-d1) and their pseudo-space representation (right, a2-d2) for: (a) equally spaced two pulses per feedback roundtrip (τ = 5.47 ns), (b) symmetry-broken two pulses per roundtrip (τ = 8.22 ns), (c) modulated quasiperiodic regime (τ = 4.78 ns), and (d) complex dynamics above the lasing threshold (τ = 4.77 ns).

Figure 2 :

 2 Figure 2: Time-domain simulations of (1) for A = 2.7 and τ = 335, showing intensity time series (left, a1-d1) and their pseudo-space representation (right, a2-d2); here κ = 0.007 in (a), κ = 0.0159 in (b), κ = 0.04 in (c), and κ = 0.1 in (d). The inset in panel (d1) shows an enlargement of the intensity time series over two roundtrips in the external feedback cavity.

  Figure 3 also illustrates individual dynamics at selected values of κ that are indicated by vertical lines in panels (a) and (b).

Figure 3 :

 3 Figure 3: One-parameter bifurcation diagrams of (1) for A = 2.7 and τ = 335 where κ is increased (a) and decreased (b). Also shown are representative trajectories in (G, Q, I) space (left), their intersection with a Poincaré section represented by the (G, Q)-plane (middle), and corresponding RF spectra (right), for κ = 0.002 (c1), κ = 0.013 (c2), κ = 0.025 (c3), κ = 0.07 (c4) and κ = 0.14 (c5,d).

  (b) shows the one-parameter bifurcation diagram of (1) for decreasing κ. The comparison between the bifurcation diagrams in panels (a) and (b) highlights a region of multistability, where the chaotic regime displayed in Figure 3(c5) coexists with the stable periodic regime shown in panels (d).

  is now investigated further by means of numerical bifurcation analysis, where the feedback delay τ and the feedback strength κ are considered as the main two bifurcation parameters. This is performed with the Matlab-based continuation software DDE-Biftool[START_REF] Engelborghs | Numerical bifurcation analysis of delay differential equations using dde-biftool[END_REF][START_REF] Engelborghs | Ddebiftool: a matlab package for bifurcation analysis of delay differential equations[END_REF][START_REF] Sieber | Dde-biftool manual-bifurcation analysis of delay differential equations[END_REF], which allows

Figure 4 :

 4 Figure 4: Two-parameter bifurcation diagram of (1) in the (τ ,κ)-plane for A = 2.7 (a), showing a curve H of Hopf bifurcation (blue), curves T of torus bifurcation (green), a curve S of saddle-node bifurcations of periodic orbits (red), as well as codimension-two points DH of degenerate Hopf bifurcation points (blue dots) and HH of Hopf-Hopf bifurcation (green dots). The enlargement in panel (b) additionally shows some of the (lowest-order) resonance points along the torus bifurcation curves and the curves S of saddle-node bifurcations of periodic orbits (magenta) bounding the corresponding resonance tongues. The inset shows a further enlargement of the framed area in panel (b), highlighting the emergence of the torus bifurcation curve from a Hopf-Hopf bifurcation point.

4 HHFigure 5 :

 45 Figure5: Bifurcation diagram in the (τ, κ)-plane for A = 2.42 near the leftmost Hopf-Hopf bifurcation point HH (green dot), showing the self-intersecting Hopf bifurcation curve H (blue), the emerging pair of torus bifurcation curves T (green), and curves S (magenta) of saddle-node bifurcations of periodic orbits that bound resonance tongues arising from selected resonance points (black dots); the open circle on the lower curve T indicates the minimum of the rotation number. Curves S near this minimum that bound resonance tongues that connect two p:q resonance points are shown in lighter magenta. The inset is an enlargement near the 1:3 point.

Figure 6 :

 6 Figure 6: Bifurcation diagram in the (τ ,κ)-plane illustrating the disappearance of the Hopf-Hopf bifurcation point HH from Figure 5 for A = 2.39 (a), A = 2.38 (b) and A = 2.35 (c). Panels (a1)-(c1) show a larger view of the Hopf bifurcation curve H and the torus bifurcation curves T, where the rotation number α along T is indicated according to the color bar and low-order resonance points are marked. Panels (a2)-(c2) are further enlargements.

Figure 6 (

 6 Figure 4(b) shows that this loop becomes smaller when the pump parameter A is decreased. Figure6shows that the loop and the point HH disappear when A is decreased further; moreover, it explains what this means for the associated curves T of torus bifurcation. The rotation number α along the torus bifurcation curves is represented by a color scale in the larger view of panels (a1) to (c1), while panels (a2) to (c2) are further enlargements. The observed configuration in Figure6(a) for A = 2.39 is topologically as that for A = 2.42 in Figure5. However, the point HH in Figure6(a) is now at a lower value of κ and the loop of the curve H is even smaller; similarly, the part of the lower torus curve T along which the rotation number α decreases is significantly smaller than for A = 2.42. Figure 6(b) for A = 2.38 is just past the transition of codimension three where the loop disappears. It suggests that this type of degenerate Hopf-Hopf bifurcation occurs when the point HH and the leftmost fold of the lower curve T coincide at a cusp point of the curve H. As a result, for lower values of A, as can be observed clearly in Figure 6(c) for A = 2.35, the two torus bifurcation curves have merged into a single smooth curve T, which is now not connected to the Hopf bifurcation curve. Notice that the rotation number α changes smoothly along this single torus bifurcation curve and has a minimum close to the leftmost fold point of T.

Figure 7 :

 7 Figure 7: Bifurcation diagram of (1) in the (τ, κ)-plane for A = 2.29 (a), A = 2.28 (b), A = 2.25 (c), and A = 2.23 (d), where the region of stable four-pulse solutions is shaded; compare with Figure 5.

Figure 8 :

 8 Figure 8: Disconnecting 1:4 resonance tongues of (1) in the relevant region of the (τ, κ)-plane, showing the relevant bifurcation curves for A = 2.29 (a), A = 2.2605 (b), A = 2.25 (c), and A = 2.23 (d); compare with Figure 7.

Figure 10 :

 10 Figure 10: Disconnecting (a) and disappearing (b) resonance tongues in the (τ, κ)-plane of (1) near a minimum m of the rotation number on the torus bifurcation T; here A = 2.25 with m = 0.234 (a1), A = 2.23 with m = 0.257 (a2), A = 2.42 with m = 0.032 (b1), and A = 2.46 with m = 0.104 (b2). Compare panels (a) with Figure 8(c)-(d) and panels (b) with Figure 5.

Figure 11 :

 11 Figure 11: Saddle transition of torus bifurcation curves generating a minimum m and a maximum M of the rotation number. Panels (a) are sketches, while panels (b) show the (τ, κ)-plane of (1) just before (b1) and after (b2) the saddle transition, for A = 2.25 and A = 2.23, respectively. Compare with Figure 8(b)-(c).

Acknowledgments

VAP, RB, IS, GB, KP and SB acknowledge partial support from the French Network of Nanotechology facilities Renatech.


Soizic Terrien
email: s.terrien@auckland.ac.nz

Venkata A Pammi

Bernd Krauskopf

Neil G R Broderick

Sylvain Barbay

Pulse-timing symmetry breaking in an excitable optical system with delay

Keywords:

Excitable systems with delayed feedback are important in areas from biology to neuroscience and optics. They sustain multistable pulsing regimes with different number of equidistant pulses in the feedback loop. Experimentally and theoretically, we report on the pulse-timing symmetry breaking of these regimes in an optical system. A bifurcation analysis unveils that this originates in a resonance phenomenon and that symmetry-broken states are stable in large regions of the parameter space. These results have impact in photonics for e.g. optical computing and versatile sources of optical pulses.

Introduction Time periodic regular pulsing regimes can emerge in many dissipative physical systems with delayed feedback [1,[START_REF] Erneux | Applied Delay Differential Equations[END_REF]. This phenomenon is encountered in various fields, from neurosciences [START_REF] Longtin | [END_REF][START_REF] Campbell | Understanding Complex Systems[END_REF] to optics and opto-electronics [START_REF] Ikeda | [END_REF][6][7][8], ecology [9] and chemistry [10,11]. Typically, these systems are multistable, with several coexisting regular periodic regimes [12,13]. Multistability has been shown to be of particular interest for all-optical processing capabilities, e.g. associative memories [14][15][16].

Here we consider an excitable optical system with delayed feedback, namely a micropillar laser with integrated saturable absorber and delayed optical feedback. In the excitable regime and for a sufficiently large delay time, the system regenerates its own output at regular time intervals [17][18][19][20]: if a short duration perturbation with sufficiently large amplitude is sent as an input, the system emits a light pulse which is re-injected by the feedback loop after a time delay τ . If the losses in the feedback loop are sufficiently low, the re-injected pulse is regenerated in the excitable medium. As the process repeats, this results in a periodic pulsing regime with period T slightly larger than τ due to the finite response time of the excitable medium.

When several perturbations are sent sequentially, the timing structure of the regenerated pulses persist in the short term. This can lead to interesting applications such as optical buffer memories [17][18][19][20]. However, in the long term, the system must converge to one of the coexisting stable pulsing regimes. These regimes generally have period T n close to sub-multiples of τ and, as such, consist of roughly n equidistant pulses in the feedback loop [16].

In this Letter we report on a different stable asymptotic behavior, where the long-term dynamics consists of nonequidistant pulses in the feedback loop. This is observed experimentally and in a model, namely a system of three coupled delay-differential equations (DDEs). A bifurcation analysis unveils that the observed non-equidistant pulsing dynamics results from destabilising bifurcations of the equidistant pulsing solutions; these occur when the delay τ is increased, provided that the recombination rate of carriers in the gain medium is faster than the one of the saturable absorber. We show that the emergence of stable pulsing patterns with n non-equidistant pulses is generic in the system and arises from symmetry breaking due to locked dynamics on invariant tori. Two stable non-equidistant pulses in the feedback loop emerge from a period-doubling bifurcation of the equidistant two-pulse solution, which is close to a 1 : 2 strong resonance. Because of the amplitude-timing coupling in the excitable system [21,22], the relative timings of the pulses are very strongly affected. This is responsible for the observed immediate symmetry breaking of pulse timings. Pulsing patterns with n equidistant pulses per feedback loop with n ≥ 3 destabilize in torus bifurcations that are close to 1 : n resonance. This results in large resonance tongues, i.e. regions of the parameter space where the dynamics on the torus is locked. Due to the amplitude-timing coupling, these locked 1 : n periodic orbits correspond, here, to higher-order symmetry-broken, that is, non-equidistant pulsing regimes. We show that non-equidistant and equidistant pulsing regimes can coexist, thus leading to a much increased level of multistability of pulsing patterns.

Model equations We consider the Yamada equations with incoherent delayed feedback [23,[START_REF] Krauskopf | Bifurcation study of a semiconductor laser with saturable absorber and delayed optical feedback[END_REF], a model of semiconductor laser written in the form of three DDEs for the dimensionless gain G, absorption Q and intracavity intensity I:

Ġ = γ G (A -G -GI); Q = γ Q (B -Q -sQI); İ = (G -Q -1)I + κI(t -τ ). (1)
The time variables are rescaled with respect to the cavity photon lifetime [START_REF] Barbay | [END_REF]. Here, A is the pump parameter, B is the non-saturable absorption, s is the scaled saturation parameter, γ G and γ Q are the recombination rates of the carriers in the gain and absorber media, respectively, and κ and τ are the feedback strength and delay, respectively. Unless stated otherwise, we consider the following parameters values: A = 2, B = 2, γ G = 0.01, (a) -1.2 γ Q = 0.055, s = 10, and κ = 0.2. The delay time τ is considered as a bifurcation parameter. This model has been shown to produce rich and complex dynamics [START_REF] Krauskopf | Bifurcation study of a semiconductor laser with saturable absorber and delayed optical feedback[END_REF]26], and to describe accurately the dynamics of an excitable micropillar laser with integrated saturable absorber subject to delayed optical feedback [16,20].

ℑ(µ) 1.2 -1.2 ℜ(µ) 1.2 (c 1 ) -1.2 ℜ(µ) 1.2 (c 2 ) -1.2 ℜ(µ) 1.2 (c 3 ) -1.2 ℜ(µ) 1.2 (c 4 ) P T 3 T 4 T 5 S 4 S 5 S 3 ✡ ✡ 2N 3N 4N ✁ ✁ 5N ✁ ✁ 2 ❅ ❅
Bifurcation analysis Figure 1(a) shows the oneparameter bifurcation diagram of system (1) in the delay time τ , where solutions are represented by their maximum value I p of intensity I. When τ is increased from zero, successive Hopf bifurcations (H) are encountered, from which several branches of coexisting periodic solutions emerge. Far from the Hopf bifurcations, these solutions correspond to the periodic emission of short light pulses, in between which the intensity is practically zero, with periods close to submultiples of the delay [26]. Each of these solutions thus corresponds to a fixed number of equidistant pulses in the feedback loop, as indicated by the numbering in Figure 1(a). The fundamental solution with one pulse per feedback loop appears at τ = 51.7 and is stable for any larger value of τ . On the other hand, all the n-pulses solutions with n ≥ 2 emerge unstably from a Hopf bifurcation, subsequently stabilize in a torus bifurcation (T) when τ is increased, and finally destabilize through a second bifurcation. All these solutions coexist with the zero-intensity equilibrium solution (i.e. the non-lasing solution), which is stable over the entire range of τ in Figure 1(a).

Figure 1(b) present enlargements of panel (a) near the second (destabilizing) bifurcations of the pulsing regimes with two to five equidistant pulses (points P, T 3 , T 4 and T 5 , respectively), and panels (c) show the Floquet multipliers of the pulsing solutions at these points. The loss of stability of the two-pulse solution occurs at point P through a period-doubling bifurcation, with one Floquet multiplier crossing the unit circle at -1 [panel (c1)]. This bifurcation is close to a 1 : 2 resonance point, where two Floquet multipliers are equal to -1 [START_REF] Kuznetsov | Elements of applied bifurcation theory[END_REF]; this has been checked by varying slightly the value of the feedback strength κ (not shown). The three-, four-and five-pulse solutions destabilize at points T 3 , T 4 and T 5 , respectively, at torus bifurcations, where two complex conjugate Floquet multipliers cross the unit circle. Panels (c2)-(c4) show that these critical Floquet multipliers are extremely close to (but not quite equal to) e ±i2π/n , which means that the torus bifurcations that destabilize the equidistant solutions are very close to 1 : n resonance points with n =3, 4 and 5, respectively. As a result, the bifurcating stable multifrequency dynamics on the invariant torus is almost immediately 1 : n locked; this happens at saddle-node bifurcations S of periodic solutions.

As Figure 1(b1) shows, these two bifurcation mechanisms lead to the emergence of additional periodic solutions. A branch of stable (period-doubled) periodic solutions emerges from the period-doubling bifurcation of the two-pulse solution. This solution of period close to τ corresponds to two non-equidistant pulses per feedback loop, and thus appears as a pulse-timing symmetry broken state. This is illustrated in Figure 2(a), which shows the evolution of the amplitudes and relative timings of pulses with respect to τ , along the branches of the two equidistant and non-equidistant solutions. After the period doubling bifurcation at τ = 472, one observes both a splitting of the pulse amplitudes in panel (a1) (as expected), but also a strong splitting or symmetry breaking of the relative pulse timings [panel (a2)], which is due to the strong time-amplitude coupling [16,[START_REF] Yanchuk | [END_REF].

We use numerical simulations to further assess how the regime with two non-equidistant pulses is accessed. plied by depleting G slightly (not shown here), the phaseshifted, symmetric version of this solution is obtained, with the first (green) and second (gray) pulses converging to the low-amplitude and high-amplitude state, respectively. Although this leads seemingly to the same long-term dynamics, both of these different states occur, one being a phase-shifted version of the other. We also point out that the convergence is very slow and occurs over several thousands of delay times, showing that the stable non-equidistant solutions are only weakly attracting.

The bifurcation mechanism leading to the emergence of non-equidistant pulsing regimes with more than three pulses is slightly different. As shown in Figure 1(b1), a pair of (stable/unstable) periodic solutions emerges from a saddle-node bifurcation, for example at τ = 663 for n = 3. This bifurcation forms the boundary of the 1 : n resonance tongue associated with the destabilizing torus bifurcation of the n-pulse solution. The emerging periodic solutions have a period close to τ , compared to the period close to τ /n of the n-pulse solution undergoing the torus bifurcation. Here, the stable 1 : n locked periodic solution corresponds to a pulsing regime with n nonequidistant pulses of different amplitude in the feedback loop. As such, the 1 : n resonance tongues are identified here as the stability regions of non-equidistant pulsing solutions. Their emergence leads to a rapidly increasing level of multistability.

Figure 3(a) shows the intensity profiles of the coexisting stable periodic solutions for τ = 1000. Here the non equidistant two-, three-, four-and five-pulse solutions [panels (a2)-(a5)] coexist with the stable one-pulse solution [panel (a1)], but also with the stable solutions with 5, 6 and 7 equidistant pulses in the feedback loop [panels (a5)-(a7)]. Overall, when τ is increased, more and more of the equidistant pulsing regimes become unstable, while more and more stable periodic solutions with non-equidistant pulses in the feedback loop appear. Typically, for sufficiently large τ , solutions with lower numbers of (at least two) non-equidistant pulses coexist with solutions with larger numbers of equidistant pulses. In Figure 3, all the periodic solutions with 1 to 7 pulses per feedback loop coexist, but the ones with two to five pulses already underwent the resonance tongue transition and, thus, correspond to non-equidistant pulsing patterns.

Figure 3(b) presents the regions of stability in the (τ, κ)-plane of feedback parameters of the different pulsing regimes with one up to eight (equidistant and non equidistant) pulses per feedback loop. Here the regions of stability of the non-equidistant pulsing solutions are resonance tongues bounded by saddle-node bifurcations. The respective stability regions of both types of solutions extend over large areas of the (τ, κ) parameter plane. Moreover, they show a large degree of overlap, which is why we show them in individual panels (b1)-(b8) for one up to eight pulses per feedback loop. This represents the high degree of multistability between all the different solutions represented in panels (a); indeed, the long-term convergence to one or the other pulsing solution depends on the chosen initial conditions. Interestingly, for n ≥ 3 both the solutions with n equidistant pulses and with n non-equidistant pulses may coexist and be stable for the same parameter values. As shown in Figure 1(b1), this results from the fact that the 1:n resonance tongues are entered (at points S) slightly before the n-pulse solutions destabilize at torus bifurcation points T. Hence, in these ranges of τ , depending on the initial condition, one observes a pattern with n either equidistant or nonequidistant pulses. In particular, in Figure 3(a5) both the solutions with 5 equidistant and 5 non-equidistant pulses coexist for the considered parameters; see also panel (b5).

Experimental realization We compare the results of the bifurcation analysis with experimental measurements of an excitable micropillar laser. It consists of two gain and one saturable absorber (SA) quantum wells [29,30], is optically pumped at 800 nm and emits light around 980 nm. Part of the output light is reinjected into the microlaser after a delay τ , through free-space propagation and after reflection on a mirror. The micropillar laser is perturbed by short optical perturbations of 80 ps duration from a mode-locked Ti:Sa laser emitting around 800 nm. In the absence of feedback, the micropillar laser is in the excitable regime [31], where the steady state intensity I is zero, but a single high-amplitude, short pulse of light can be emitted in response to an external perturbation of sufficient amplitude [START_REF] Barbay | [END_REF][START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF]. In the presence of feedback, an excitable pulse is regenerated when it is reinjected by the delay loop after the delay τ , thus resulting in the regular emission of light pulses at a period close to τ [19,20].

In system (1) the ratio of the recombinations rates γ G and γ Q of the gain and SA media, respectively, play a crucial role in the pulsing dynamics [21,[START_REF] Dubbeldam | [END_REF]34]. In particular, the pulse-timing symmetry breaking is observed only for a faster gain recombination, that is, for γ G > γ Q . Experimentally, this parameter regime can be accessed by selecting a suitable micropillar laser (from many on the same chip) by taking advantage of the spread of physical parameters in the course of the nanofabrication process.

Figure 4 shows the evolution of relative pulse timings of experimental pulse trains over several hundreds of roundtrips in the feedback loop, in the same representation as in Figure 2(b2). The microlaser is started (after suitable external perturbations) with two [panels (a)] and three [panels (b)] almost equidistant pulses in the feedback loop. We observe that the pulse-timing information is preserved in the short term [16,20]. On the other hand, in the long term, the system slowly converges towards non-equidistant pulsing patterns with well-defined and different interpulse relative timings. These interpulse timings then stay very stable over a large number of roundtrip. It was not possible to monitor the amplitude difference in the final state due to the limited signal to noise ratio -the emitted pulse energy is on the order of only 100 fJ. On the other hand, in agreement with Figure 2, even a small difference in amplitude is associated with a large interpulse interval difference in the non equidistant pulsing regime. Overall, the experimental observations show excellent agreement with the dynamics predicted by the bifurcation analysis of the model.

Discussion and conclusions

We demonstrated that an optical excitable system with delayed feedback can sustain stable pulsing patterns with different numbers of non-equidistant pulses in the feedback loop. These arise from stable solutions with n equidistant pulses via torus bifurcations and associated 1 : n resonances, which manifest themselves as a swift breaking of the timingsymmetry due to the strong amplitude-timing coupling of the excitable system. We find stable non-equidistant pulsing in large resonance tongues in the parameter space, bounded by saddle-node bifurcations of periodic solutions. As the delay is increased, there is a high and increasing degree of multistability between both symmetric and symmetry-broken pulsing patterns in the feedback loop. Which long-term behavior is observed depends on the initial condition. We demonstrated that non-equidistant pulsing can be observed reliably in an experiment with an excitable micropillar laser.

Our results are reminiscent of the pulsing dynamics of models describing delay-coupled neuron by either two limit cycle oscillators coupled through a time dependent synaptic response [35,36] or a pulsing oscillator with delayed self-coupling [37,38]. In our case however, oscillations do not pre-exist and originate from the delayed feedback itself and their period is intimately linked to the delay time. This further illustrates that our results are expected to be generic and to extend beyond optics. Moreover, a mathematical connection between temporal dissipative solitons in spatially extended systems and pulsing regimes in delay systems has been recently suggested [START_REF] Yanchuk | [END_REF]. This raises open questions on possible connections between non-equidistant pulsing regimes and soliton molecules, which are bound states of pulses [START_REF] Grelu | Lecture Notes in Physics[END_REF][START_REF] Krupa | [END_REF].

Beyond their fundamental interest for the nonlinear dynamics of delay systems, our results may contribute to the realization of non-conventional pulsing and reconfigurable optical sources [41], and to optical computing schemes [42][START_REF] Peng | Conference on Lasers and Electro-Optics[END_REF][START_REF] Pammi | [END_REF][45] that rely on the large phase space available in delay systems [46].

Figure 1 .

 1 Figure 1. (a) Bifurcation diagram of(1), showing the pulse intensity Ip with respect to τ , with the number of pulses per feedback loop along each stable periodic solution branch. (b1) Enlargement of the framed area in (a), with further enlargments around point T3 (b2) and S3 (b3). Stable equidistant (E) and non-equidistant (N) pulse solutions are represented in dark and light blue, respectively, and unstable E and N solutions in dark and light orange, respectively. The dots indicate Hopf (H), torus (T), period doubling (P), saddle-node (S) and homoclinic (L) bifurcations. (c1)-(c4) Floquet multipliers at points P, T3, T4 and T5, with critical multipliers highlighted in red.

Figure 2 (Figure 2 .

 22 Figure 2. (a) Maximum Ip of pulse intensity (a1) and relative interpule timings tp (a2) along the branches of two equidistant and non-equidistant pulses, with respect to τ . Stable and unstable solutions are represented in blue and red, respectively. (b) Simulation of (1) for τ = 1000 (grey lines in panel (a)) with initial condition very near the (unstable) two-pulse solution, showing the long-term evolution of Ip (b1) and of tp (b2). The subpanels in (b1) show the intensity time series during the two first and two last roundtrips through the feedback loop; the dots and arrows indicate the amplitudes and relative timings as represented in (b1) and (b2), respectively.

Figure 3 .

 3 Figure 3. (a) Intensity profiles of stable periodic solutions of (1), for τ = 1000. (b) Regions of stability, in the (τ, κ)-plane of feedback parameters, of the families of equidistant (E) and non-equidistant (N) periodic solutions with 1 to 8 pulses per feedback loop. The number of pulses is indicated in the colored regions, and the star indicates the parameter point (τ, κ) = (1000, 0.2) of the time series in panels (a).

Figure 4 .

 4 Figure 4. Evolution over the roundtrip number of the relative interpulse timing tp of experimental pulse trains following two (a) and three (b) external perturbations, for a feedback delay of τ = 8.2ns: just after the external perturbation (panels 1), during the convergence towards non-equidistant pulsing patterns (panels (2)) and in the long-term (panels 3)).

Acknowledgments S.B. and V.A.P. acknowledge support from the CNRS Renatech Network of Technology for the nanofabrication of the samples.