Stability diagrams for continuous wide-range control of two mutually delay-coupled semiconductor lasers

The dynamics of two mutually delay-coupled semiconductor lasers has been frequently studied experimentally, numerically, and analytically either for weak or strong detuning between the lasers. Here, we present a systematic numerical investigation spanning all detuning ranges. We report high-resolution stability diagrams for wide ranges of the main control parameters of the laser, as described by the Lang–Kobayashi model. In particular, we detail the parameter influence on dynamical performance and map the distribution of chaotic pulsations and self-generated periodic spiking with arbitrary periodicity. Special attention is given to the unfolding of regular pulse packages for both symmetric and non-symmetric configurations with respect to detuning. The influence of the delay –time on the self-organization of periodic and chaotic laser phases as a function of the coupling and detuning is also described in detail.

Apart from the mentioned practical applications, delay-coupled lasers are theoretically challenging because of the underlying infinite-dimensional equations of the paradigmatic Lang-Kobayashi phenomenological model normally used to describe them. For instance, the authors of [20][21][22][23] provided an in-depth analysis of compound laser modes (CLMs), the analogue of the external cavity modes (ECMs) in the context of external cavity-delayed semiconductor lasers. The case of two coupled external cavity lasers was considered in [24] and [25]. Rogister and Blondel [26] presented a detailed study of the control parameter space spanned by the laser coupling and detuning, showing how periodic and chaotic phases self-organize. The theoretical analysis provided by all of these works was done considering an adaptation for coupled lasers of the Lang-Kobayashi phenomenological model, whose results have shown good agreement with experiments [14,18]. (For firstprinciples derivation of this model see, e.g., [27].) The vast literature dealing with mutually delay-coupled semiconductor lasers provides much insight about the different dynamical behaviors of the system. However, for distinct reasons, the investigation of such systems has been frequently divided into two non-overlapping classes, referred to as the weak and the strong laser detuning limits. Furthermore, several publications deal with unstable phenomena, namely, phenomena that cannot be readily observed in laboratory experiments. Works dealing with observable stable phenomena frequently consider only fixed points, equilibria, and tame oscillations observed immediately after Hopf bifurcations. Little is known about the location, shape, and extension of stable laser phases corresponding to complicated laser self-pulsing of arbitrary periodicity or for phases of chaotic laser spiking.
The aim of this paper is to help bridge this lack of information by presenting a systematic numerical classification of all sorts of oscillations, of arbitrary periodicity or chaotic, which are supported by two mutually delay-coupled semiconductor lasers. Using the standard Lang-Kobayashi model, we classify laser self-pulsing for continuous coupling strengths and for detuning spanning from the weak to the strong laser detuning limits. More specifically, we report high-resolution stability diagrams for quite wide ranges of the main laser control parameters.
Two absolutely identical lasers with a coupling delay τ behave exactly like a single external cavity-delayed laser with delay time τ [28], a situation that we already considered in detail elsewhere [29] (see also the thesis work by Otto [30]). However, this somewhat idealized situation is very difficult to achieve in practice. For instance, variations in temperature can alter the lasers' optical angular frequencies [26], leading to some degree of detuning between the lasers. So, it is important to ask how detunings of this or of any other origin affect the dynamical performance of the coupled system. Are small detunings effective to destabilize laser operation? How robust are the several laser phases to changes in parameters? What do the main alterations that control parameter changes imply?
These are the main questions addressed in this paper. To this end, we report results of numerical simulations for a wide range of control parameter space of two mutually delay-coupled lasers, contrasting their performance with and without detuning. In addition, given the key importance of the delay parameter in this kind of system, we extend work done by Rogister and Blondel and investigate how the variation of the delay time τ affects the lasers under the simultaneous variation of detuning and coupling.

Model of two mutually delay-coupled laser diodes
As mentioned, pairs of single-mode semiconductor lasers mutually delay-coupled by their optical field were discussed previously, e.g., by Erzgräber et al [20][21][22][23], Li et al [24], Hicke et al [25], Rogister and Blondel [26], and, most recently, by Yang et al [3]. Such a configuration is represented schematically in figure 1. It can be modeled by a set of rate equations adapted from the Lang-Kobayashi model [31], providing the temporal evolution of the normalized slowly varying complex electric fields E t ( ) 1,2 and the normalized excess carrier numbers N t ( ) 1,2 for both lasers 1 and 2, respectively. In dimensionless units, the model is governed by the equations [14,26,27]  In these equations we use the optical angular frequency of laser 1 as a reference so that ω ω = 1 while the detuning frequency is given by δ ω ω ≡ − 2 1 . Time is normalized to the cavity photon lifetime ∼ ps ( 1 ), and T is the ratio of the carrier lifetime ∼ ns ( 1 )to the photon lifetime. The delay time τ represents the flying time of the light between the lasers (see figure 1), η controls the strength of the symmetric coupling, δ controls the magnitude of the detuning between the lasers, P is the normalized pump current above threshold, and α is the linewidth enhancement factor.
Note that equations (1)-(4) are written in the reference frame of the unchanged laser, i.e., laser 1. Accordingly, the optical fields of the lasers are represented by where ω 1 is the optical frequency of laser 1 operated solitarily at threshold. Apart from the difference in their solitary optical frequencies, both lasers are considered identical [20]. Unless stated otherwise, we use the following set of realistic laser parameters [29,32]: Here, we analyze the output of laser 1 (| | E 1 2 ), but the output behavior of laser 2 is found to be qualitatively identical to the one of laser 1, even in the asymmetric configuration. Therefore, the conclusions presented remain valid for both lasers. Before proceeding, we remark that the initial history used to start numerical integrations is an important piece of information to define the dynamics. Surprisingly, such information is traditionally omitted in the literature, so that comparisons are difficult to be made. The details about the initial configuration used in this paper are given in the appendix.

Dynamical performance of two mutually delay-coupled laser diodes
As already mentioned, for δ = 0 and when starting from the same parameter values and initial configuration, both lasers present the same output and behave exactly as a single laser with feedback from an external cavity of flying time τ. But when the optical angular frequency of the lasers are different, namely for δ ≠ 0, this symmetry is broken, and different dynamical behaviors arise.
We start by illustrating the typical temporal evolutions of | | E 1 2 for representative parameters of the coupled laser model. In  Figure 2(a) shows a periodic oscillation of the laser intensity with one peak per period, obtained for η = 0.1285. Increasing the coupling to η = 0.1342, the system enters a quasiperiodic regime, illustrated in figure 2(b). For η = 0.1600, one finds the so-called Regular Pulse Package (RPP) (figure 2(c)), a classical pulsation described most recently in [29]. In the RPP regime, the laser output is characterized by a sequence of packages composed by decreasing amplitude pulses with frequency related to the external cavity flying time. This type of regular pulsation has been extensively observed and described in semiconductor lasers with delayed feedback (see [29] and several references therein). An increase of the coupling to η = 0.1994 leads to a non-periodic regime, depicted in figure 2(d). This non-periodic solution reveals a sort of 'quasi-RPP' profile, very similar to the classical RPP but with small amplitude differences between the packages. Note, for example, the small differences between the large pulses of each individual package, magnified in the inset, which show an irregular distribution of amplitudes.
The appearance of RPPs in mutually delay-coupled lasers in a totally symmetric configuration was expected by the aforementioned equivalence between this system and the delayed feedback laser. The question we want to answer at this point is what happens if we somehow break the system symmetry. To this end, we consider a small difference between the optical angular frequencies of the lasers, setting δ ω ω = − = − 10 2 1 4 . For the value ω ω = 1 that we use (see equation (5)), the above choice of δ leads to a difference between ω 2 and ω 1 which is less than 1% of the value of ω 1 . Taking the relaxation frequency of laser oscillations as a reference, given here by 2 , one sees that a detuning of 10 −4 is clearly very small. But, as discussed below, it is not without measurable consequences.
The temporal evolution of laser 1 intensity shown in figures 2(e) and (f) were calculated with the same values of η used in figures 2(c) and (d) (0.1600 and 0.1994, respectively) but for δ = − 10 4 , representing the asymmetric detuning configuration. Differently from what was obtained in the symmetric case, here we see extremely irregular pulse packages in temporal evolution. Comparing the RPP of figure 2(c) with the non-periodic outputs of figures 2(d) (δ = 0) and (e) (δ = − 10 4 ), we see that in the symmetric case, the increasing of η induces a small destabilizations of the RPP attractor, leading to phases that, despite being non-periodic, are not far from RPPs (figure 2(d)), while the asymmetry produces a much more intense effect, totally destroying the RPP pattern (figure 2(e)).
Next, figures 3(a)-(c) display bifurcation diagrams of the local maxima of the laser amplitude | | E 1 2 (which, in this case, is identical of | | E 2 2 ) as a function of the bifurcation parameters τ, P, and η, for a perfectly symmetric configuration (δ = 0). The vertical lines in figure   laser [29]. In figure 3(a), there is a continuous sequence of transitions of laser regimes as the delay is increased: constant output (CW, orange), one-peak oscillations (P1, blue), quasiperiodic oscillations (QP, red), and an irregular alternation of regular pulse packages (RPP, green) and non-periodic oscillations (NP, black). The same regimes appearing, however, in a different sequence are observed when the pump current is varied, namely, CW, one-peak oscillations, quasiperiodicity, and the alternation of RPPs (green) and non-periodic oscillations (NP, black). Figure 3(c) shows the bifurcation diagram when increasing the coupling parameter η. In this case, the system presents a mix of behaviors seen in the variation of τ and P. For low values of η, the same transitions seen for the variation of τ take place. This kind of transitions hold until η ≈ 0.2, when the sequence CW, one-peak oscillation, and quasiperiodicity disappears, with only the irregular alternation between RPPs and non-periodic solutions remaining, similar to what is seen in figure 3(b), for large values of P.
Comparing the bifurcation diagrams of figures 3(d)-(f), calculated with δ = − 10 4 , with the ones shown in figures 3(a)-(c), calculated with δ = 0, we can clearly see that a small detuning is sufficient to generate a great alteration in the system dynamics. The other parameters and initial configuration were held the same.
Comparing the diagrams of figures 3(a) and (d), where the delay time is the bifurcation parameter, we note that the typical sequence of transitions, seen when increasing τ in the symmetric configuration, are maintained in the asymmetric case, but while the symmetric case presents only one-peak oscillations between CW and quasiperiodic phases, shown in blue in figure 3(a), in the asymmetric configuration, after the CW phase, an extremely narrow P1 window takes place (not labeled in the bifurcation diagram), followed by a phase with an additional branch created by continuous deformation of the output waveform [29,33,34], and two-peak oscillations (P2) take place in this region, as shown in light-blue in figure 3(d). Still, the small increasing in δ also generated a considerable enlargement and increase in the pulse amplitudes in the non-periodic phases, as can be seen by comparing the regions in black (NP) between figures 3(a) and (d).
Another interesting feature that we can observe comparing figures 3(a) and (d) is that, in the same range of the delay time τ, the diagram calculated using δ ≠ 0 presents more sequences of the transitions described above compared to the case with δ = 0, i. e., the transitions CW, P1 oscillations (P2 oscillations), quasiperiodic oscillations, and aperiodic oscillations (+RPP), happens in a shorter range of the delay when the system symmetry is broken. For the case of a single laser with delayed feedback from an external cavity, which is analogous to the symmetric configuration we are considering here, this sequence of transitions was explained before [29]: the increasing of τ leads to the creation and destabilization of ECMs, a solution with constant amplitude and inversion, and a linear variation of the phase with time, which strongly influences the system dynamics. In the case of mutually delay-coupled lasers, we have solutions which are analogous to what ECMs represent for single lasers with delayed feedback, the CLMs. However, as pointed by Erzgräber and Krauskopf [22], the relation between the creation and destabilization of such solutions and the lasers dynamics are much more complicated here than in the case of a single laser with delayed feedback, and it is very difficult to understand the overall picture of how the CLMs depend on the system parameters. Be that as it may, Erzgräber and Krauskopf gave a relation between the CLMs and parameters α, η, ω, τ, and the detuning δ. We believe that the increase in the number of transitions over the observed range of τ in the asymmetric configuration is a result of the influence of the detuning δ in the way that the delay time alters the number and/or the stability of CLMs, but a deeper theoretical analysis is needed to clarify this point. This analysis is beyond the scope of this paper and will be performed elsewhere.
The impact of increasing the pump current P is also strongly affected by the small difference in angular optical frequencies, as can be seen by comparing figure 3(e) with figure 3(b). The first big difference between these figures is the two-peak solution, indicated in light-blue in the asymmetric case, and which is absent in the symmetric case. Another interesting difference is that, after the quasiperiodic phases, shown in red, the irregular alternation between aperiodic and periodic (RPP) solutions, shown in black and green, respectively, in figure 3(b), does not take place in the asymmetric case, but only a completely chaotic phase arises , showing solutions with higher peaks. As for the variation of the feedback rate η, the impact of having δ ≠ 0 is even stronger. Comparison of figures 3(f) and (c) shows that, although small, the asymmetry was enough to completely destroy all periodic solutions, leaving behind a continuous chaotic phase with much higher output amplitudes.
Comparison of the bifurcation diagrams calculated using δ = − 10 4 with the ones using δ = 0 reveals a most interesting phenomenon: the complete disappearance of the RPPs in the asymmetric case. The domains of parameter space where RPPs take place are evidenced by the regions marked in green in the bifurcation diagrams in figures 3(a)-(c). These RPP regions emerge interspersed with regions of aperiodic oscillations and are vastly distributed over parameter space, specially seen under the variation of the feedback rate η in figure 3(c). When the optical angular frequencies' difference is tuned to δ = − 10 4 , we observe the complete disappearance of the RPPs over the entire observed range of τ, P, and η, as can be clearly seen in figures 3(d)-(f), respectively.
The disappearance of RPPs for a detuning as small as 10 −4 suggests that the dynamical mechanism behind this property is the destabilization of the RPP attractors by the symmetry-breaking related to the difference in optical angular frequencies of the lasers. In the case of external cavity lasers, the RPPs are known to present a phase space trajectory which visits a defined sequence of ECMs [29,32]. As discussed above, we believe the complex influence of the detuning in CLMs' properties to be such as to destabilize the RPP attractors, leading to much more irregular solutions, such as the ones observed in figures 2(e) and (f).

Behavior of RPPs over wide control parameter ranges
The bifurcation diagrams in figure 3 describe how the laser output in a delay-coupled configuration changes under the variation of just a single of the fundamental parameters τ, P, and η, while holding all other parameters constant. We also showed that the RPPs, which are abundant in the symmetric configuration (δ = 0), disappear when the difference in optical angular frequencies is not zero, even for an amount as small as δ = − 10 4 . A natural question now is to ask, what happens under much larger variations of δ and of other parameters? Is the vanishing of RPPs a general property in parameter space, observable also when more than one parameter is varied simultaneously? To clarify these questions, we performed a detailed numerical investigation, calculating highresolution phase diagrams over wide parameter ranges. The construction of such diagrams follows the procedure described in [29,33,34]. We analyze the laser self-pulsations over large control parameter windows, classifying them dichotomically as periodic or non-periodic oscillations and, subsequently, counting the number of local maxima per period of the periodic pulsations. Figure 4 shows isospike stability diagrams [35] obtained when simultaneously varying τ and P, and τ and η over wide ranges. The diagrams in figure 4 depict the number of spikes (local maxima) contained in one period of the periodic lasers pulsations, recorded using 17 basic colors, as indicated by the colorbar in the figure and explained in the appendix. The ranges of parameters representing the bifurcation diagrams in figure 3 are indicated by the labeled lines over the diagrams. The analysis of laser 1 presented the same diagrams of laser 2 for all set of parameters considered.
The τ × P diagram shown in figure 4(a) was calculated considering the symmetric configuration (δ = 0). As expected, this panel coincides with one obtained in [29] for a single laser with delayed feedback (and using the same set of parameters and initial configuration). A sequence of 'hills' of CW, one-peak oscillations and quasiperiodic oscillations is seen, where the height of the hills decreases as τ is increased. In the region between the hills one finds an alternation of colorful stripes corresponding to RPPs and to non-periodic oscillations. This behavior has been described and explained in detail in [29]. Figure 4(b) shows a stability diagram in the τ η × control plane for the symmetric configuration, evidencing two types of behavior. For low values of the delay, stripes of CW, one-peak oscillations, and quasiperiodic oscillations alternate with stripes of RPPs and nonperiodic oscillations. These stripes follow nearly hyperbolic curves η τ ∼ 1 in the diagram. In contrast, for larger values of the delay, one finds more irregular stripes alternating RPPs and non-periodic oscillations. These two types of behavior can be appreciated along the vertical line 'C' in figure 3(c), a line that clearly crosses both regions despite the fact of crossing the diagram vertically.
The stability diagrams in figures 4(a) and (b) contain several multicolored stripes representing oscillations having many spikes per period, as indicated by the colorbar. The large majority of the oscillations having more than two spikes can be shown to be RPPs, as may be seen in the green regions in the bifurcation diagrams A, B, and C of figure 3. Figures 4(a) and (b) display, respectively, 7.38% and 20.83% solutions with more than two spikes, i.e., solutions that are very likely to be RPPs. When the symmetry is broken, e.g., for δ = − 10 4 , the diagrams of figures 4(c) and (d) reveal a complete absence of the regions with more than two peaks per period over the entire windows, suggesting this absence to be a general property of the symmetric-asymmetric transition. In addition, the hills in figure 4(c) and the stripes in figure 4(d), both formed by CW, one-peak oscillations, two-peak oscillations, and quasiperiodic oscillations, are thinner than the ones seen when δ = 0. As mentioned above, this is due to the role of δ in the creation/destabilization of CLMs.
The nonzero value of the detuning used so far is obviously quite small, since it is much smaller than the reference value of the frequency of the relaxation oscillations of both lasers (ω ∼ × − 3.6 10 RO 2 ). So, it is natural to ask what happens for larger detuning. To answer this, figure 5 illustrates detunings for two much larger values, namely, for δ = × − 2 10 2 , a value close to ω RO , and for δ = 0.2, a value considerably higher than ω RO . In both cases there is a distribution of periodic and non-periodic phases very similar to what was observed previously in figure 4(d) for very small detuning. Therefore, figure 5 shows that the main impact of the detuning in the system performance comes from the transition between the symmetric configuration (δ = 0) to the asymmetric configuration (δ ≠ 0), the relative magnitude of the detuning being relatively unimportant. Despite that, an increase of the detuning seems to gradually reduce the size of the stripes with transitions between CW, one-peak oscillations, and quasiperiodic oscillations. On the other hand, for a large detuning, we can observe a big phase of one peak and quasiperiodic oscillations for small values of the coupling η. Also, some very small regions of more complex periodic solutions appear in the central region of this diagram. These regions suggest that RPP solutions may, eventually, reappear in the system for specific values of the detuning. This tendency evidences the need of a more comprehensive and embracing analysis of detuning effects in the system.

Impact of delay on the coupling versus detuning control plane
So far, we have described the different types of solutions that can be obtained for the laser intensity when considering specific values of the detuning between the lasers. We found that even a very small detuning δ is able to make the RPPs, normally quite abundant in the symmetric configuration, practically disappear over the considered parameter ranges. We also showed that an increase of the detuning tends to enlarge the regions of non-periodic solutions in the system. But what happens when we consider the influence of the detuning over a much wider range? Does a strongly detuning between the lasers necessarily means more non-periodic phases over the parameter space? Analyzing the bifurcation diagram shown in figure 6, which considers the detuning δ between the lasers as the bifurcation parameter, it is possible to recognize that this is not always true. In this diagram, for values of the detuning up to δ ∼ 0.308, the laser presents a CW output. Above this value, the system enters a wide region of δ where the solution displays non-periodic oscillations. Nevertheless, some windows of RPPs can be seen over this phase of non-periodic outputs, showing that periodic solutions can be achieved for a specific set of parameters in the asymmetric, δ ≠ 0, configuration of the lasers.
Thus, despite of the destabilization of RPP attractors observed above, triggered by the symmetry-breaking provided by even very small differences in the lasers' optical angular frequencies, the influence of the detuning on the system dynamics can be such as to promote recurrent returns of these solutions, for suitable values of the detuning. A more specific and detailed explanation on how δ influences the appearance and disappearance of the RPP windows in parameter space would require an in-depth analysis of the complicated influence of this parameter on CLM properties. As already mentioned, such analysis is beyond the scope of the present paper.
Taking into account the relatively ease of experimentally varying the coupling and the detuning intensities in this system, Rogister and Blondel [26] performed a detailed numerical analysis of the η δ × control parameter space for a fixed value of the delay time (τ = 20). They presented stability charts dividing the behavior of the laser intensity in parameter space into four types: stationary, periodic P1 (one peak per period), periodic P2 (two peaks per period), and other regimes. When increasing the coupling, their periodic phases presented transitions stationary/P1 oscillations in almost triangular phases centered in regions with small detuning of the charts.
Their findings motivate us to investigate what happens for other values of the delay τ. More specifically, we wish to understand how the periodic phases self-organize for values of τ spanning a macroscopic range of the short-cavity regime. So, in order to show how the separation between the lasers affects the overall picture in control parameter space, we analyze the η δ × diagrams, constructed like the ones of figure 4, for several values of the delay-time τ, as shown in figure 7.
We start with the case without delay (τ = 0). In this case, the diagram is almost symmetric with respect to δ = 0 in a way that, for τ δ ≳ | | 2 , the laser operates in the CW regime, and for τ δ ≲ | | 2 the laser oscillates with one peak per period, apart from small regions near the curve τ δ = | | 2 , where more complex oscillations appear. For τ = 1 the region of one-peak oscillations diminishes, and a thick phase of non-periodic oscillations emerges between CW and one-peak regions. The symmetry of the diagram also starts to vanish. At τ = 5 the CW regions begin to self-organize, forming a sequence of almost triangular phases. These phases present an internal structure where an increase in η produces a transition CW, one-peak oscillations, and quasiperiodic oscillations. Further increase of the delay gradually reduces the size and increases the number of the triangular phases, which are located along the δ = 0 horizontal line in the diagram, and get deformed for bigger values of the feedback rate η. This effect is shown in the diagrams with τ ranging from 10 up to 50. As the delay-time is further increased, these phases suffer a narrowing, turning into small, almost vertical stripes for τ = 50, placed in a growing background of non-periodic oscillations. Note that our diagram for τ = 20 is very similar to the ones constructed considering the same delay value, shown in [26]. The small differences are due to the different choices of the frequency ω and the fact that we show a much broader range of the coupling and the detuning parameters. The initial history used to produce the phase diagrams of [26] was not mentioned. Our results are robust against slight variations of the initial history. The sequence of transitions observed when increasing the coupling parameter along the triangular phases in figure 7 is the same sequence observed along the hills in figures 4(a) and (c) when the delay-time is increased: CW, one-peak oscillations, and quasiperiodic oscillations (for some parameter values, a two-peak oscillation phase also appears between the one-peak and quasiperiodic phases). The same sequence is also observed when the simultaneous variation of these two parameters is considered, as can be seen in figures 4(b) and (d). When starting from the CW operation, the above scenario shows that an increase in either the delay-time or in the coupling coefficient produces very similar effects in the overall performance, leading to the same sequence of transitions in the oscillatory pattern for a considerably large range of parameter values.
The major effect of the delay time is to compress the periodic phases into regions characterized by smaller values of the magnitude of the detuning. This result contributes to the understanding of how the almost triangular structures of periodic solutions seen previously in the diagrams for τ = 20 arise, here and in [26]. We also note that increased delays tend to enlarge the size of the aperiodic (chaotic) phases for large detuning in the system.

Conclusions and outlook
A remarkable and unexpected result of our simulations is that frequency differences much smaller than the relaxation oscillation frequencies of the lasers are enough to significantly alter their dynamical performance. . Isospike stability diagrams [35] illustrating the impact of increasing the delay time τ on the coupling (η) × detuning frequency (δ) charts. Orange: CW, white: quasiperiodic oscillations, black: aperiodic oscillations.
The major difference observed is that the RPP phases, abundant in the symmetric case, disappear in the asymmetric configuration, giving birth to phases of highly aperiodic oscillations. This effect was observed over large windows of the τ × P and τ η × control parameter spaces. Hohl and Gavrielides [14] had previously considered a symmetry-breaking of two mutually delay-coupled diode lasers for different pump currents and coupling strength. They discovered an interesting type of localized synchronization characterized by intensity oscillations with very different amplitudes. Our stability diagrams serve now as a guide to select control parameters for enhanced synchronization of this type over wide parameter ranges in the small-cavity regime.
In addition, we analyzed how an increase of the delay time alters the distribution of regular and chaotic pulsations in the η δ × parameter plane. We explained how the periodic structures seen in the stability diagrams obtained previously for a single value of the delay are born. We find that the system starts with large regions of CW and one-peak oscillations for the case without delay. As the delay is increased, sequences of almost triangular structures, composed by CW, one-peak oscillations, and quasiperiodic oscillations are arranged along the regions around δ = 0, while aperiodic phases dominate regions of large detuning. We detected a great resemblance between the effects induced in oscillatory patterns when increasing the delay time and the coupling coefficient, leading us to believe that both parameters destabilize the CW solutions in a similar way. Obviously, a large number of questions still remains to be answered: for instance, the characterization of multistability and its impact in the lasers, as well as the robustness of the coupled lasers against unavoidable parameter mismatch between them. From a theoretical point of view, the lack of a mathematical analysis of the questions discussed here draws attention to a pressing reality: there is absolutely no theoretical framework yet to address the systematic unfolding of solutions beyond the traditional analysis of fixed points and the many approximations derived for the tame oscillations that emerge immediately after Hopf bifurcations [37].
The isospike stability diagrams [35,36] in figures 4, 5, and 7 were obtained by numerically solving the laser equations equations (1)-(4) using the standard fourth-order Runge-Kutta algorithm with fixed-step, h = 0.005, over a mesh consisting of 500 × 500 equally spaced points. For every mesh point, we started numerical integrations from a fixed and arbitrarily chosen initial configuration: = N (0) 1 . The first 10 7 time-steps were discarded as transient time needed to approach the final attractor. The subsequent 10 7 iterations were then used to compute the number of spikes contained in one period of the oscillations. As indicated by the colorbar in the figures, a palette of 17 colors is used to represent periodic pulses with more than 17 spikes per period of the oscillations of | | E t ( ) 1 2 . We 'recycled colors modulo 17' meaning that the color index used is obtained as the remainder of the integer division of the number of spikes by 17. Solutions with 18 spikes per period are marked with the same color associated with one spike, 19 spikes are marked with same color of two spikes, and so on. Multiples of 17 were given the index 17. In this way all periodic pulses can be accommodated with a palette of 17 colors. In addition, non-zero fixed points, representing constant output, are plotted in orange, quasiperiodic solutions are plotted in white, and a lack of numerically detectable periodicity is plotted in black. The isospike diagrams were obtained, after discarding the transient, by recording the local maxima of the time series, together with the instant when they occur, and recording repetitions of the maxima. The computation of stability diagrams is a quite demanding numerical task that we performed with the help of 1536 processors of a SGI Altix cluster having a theoretical peak performance of 16 Tflops.