Semiconductor laser with optical feedback: from excitable to deterministic low-frequency fluctuations.

Semiconductor lasers with optical feedback present a regime in which power dropouts are observed. Although this regime has been extensively studied, there is no agreement about whether the dropouts are deterministically or stochastically generated. In this paper we will study the statistics of time intervals between dropouts assuming noise-driven simple excitable models. We explain the appearance of characteristic times in the first return maps.


I. INTRODUCTION
Semiconductor lasers with optical feedback display a rich variety of behaviors that are interesting both from laser theory and from a nonlinear dynamics point of view. A widely studied but not yet fully understood phenomenon is the so-called low-frequency fluctuation ͑LFF͒ regime, which occurs for moderate to strong feedback levels. In this regime, the laser output intensity suddenly drops to almost zero and then gradually recovers to its original value. These dropout events occur at irregular time intervals and become more frequent as the pump rate is increased.
Because both size and shape are quite similar from one pulse to another, the time interval between them is a most interesting observable. Statistical analysis of these time spacings has been performed ͓1-4͔. Specifically, histograms and first return maps present a remarkable feature: when the dropouts are very rare the return map shows a cloud of points and the histogram is a single peak with an exponential tail but, for higher currents values, certain characteristic times arise ͓1,2͔.
There are several approaches to understanding this behavior. Most of them developed from the Lang-Kobayashi equations, which explains the LFF as a form of high-dimensional chaos ͓5,6͔. Recently, a different approach to the problem was proposed by Giudici et al. ͓1͔ based on experimental observations. According to this interpretation, the laser behaves as an excitable system and the noise plays a dominant role in the appearance of the LFF.
The confrontation between these two paradigms highlights a central problem in nonlinear dynamics. Once one acknowledges that deterministic rules are able to provide complex dynamical behavior, how can we recognize whether a dynamically complex system is purely deterministic? The construction of refutable quantities is by no means trivial ͓7͔. In this paper, we investigate the behavior of different noise-driven excitable systems, and analyze the interspike time distribution that each model predicts. We compare the behavior of these observables as parameters are changed, and compare them with the experimental data reported in Refs. ͓1,2͔.
The definition of an excitable behavior is not uniform throughout the literature. Here we assume that a system is excitable when it displays the following features: ͑a͒ it has a stable equilibrium state, ͑b͒ if the system is perturbed from its equilibrium state beyond a certain threshold, it relaxes after a long excursion in phase space, and ͑c͒ the size of this wide relaxation becomes independent of the size of the perturbation provided the latter exceeds the threshold value. In the case of the LFF regime, wide relaxations correspond to dropout events. In this paper, we will analyze two simple noise-driven excitable dynamical models. The first one is a two-dimensional system in which the threshold is given by the distance ͑in phase space͒, between a stable point and a saddle point. We chose to begin our study with this model for historical reasons, following the dynamical scenario proposed in Ref. ͓1͔. Later, we will move on to build a threedimensional model in which the distance between the attractor and an unstable periodic solution acts as a threshold.
For both models, dropouts will be obtained as noise perturbs the equilibrium states. We will show that, as long as we are interested in the longer time scales of the system ͑typically, the time between dropouts͒, the two-dimensional model is enough to reproduce the statistics observed in the experiments. However, the interspike histograms can display some structure in the shorter time scales ͑multipeaked distri-bution͒, as the injection current is raised. In order to explain this feature, we need to extend our dynamical model one more dimension. We will be careful to check that the threedimensional model also reproduces the statistical features observed in the two-dimensional one.
This paper is organized as follows. First, in Sec. II, we describe a simple model for an excitable system ͑following ͓1͔, a system close to an Andronov bifurcation͒, and study its interspike time distribution once noise is added. In Secs. III and IV, we address the issue of multipeaked histograms. We progressively build a three-dimensional excitable model that allows for the shorter time scales of the system and reproduces the interspike histograms and first return maps reported in the literature. The relationship between these two models is discussed in the last section, where we include our conclusions.

II. TWO-DIMENSIONAL EXCITABLE MODEL
Several groups have reported extensive experimental studies of interspike time distributions for semiconductor la-sers with feedback in the LFF regime ͑see ͓1,2͔ and references therein͒. In ͓2͔, a comparison between the predictions of a model ͑derived from the Lang-Kobayashi equations ͓8͔͒ and experimental results is performed. The experimental interspike distributions display, in a parameter range where LFF is observed, a bimodal structure. This is the first issue that we will address.
We are interested in whether this statistical behavior can be reproduced within the framework of the hypothesis of Ref. ͓1͔, i.e., a noise-driven excitable system close to an Andronov bifurcation. This bifurcation is locally a saddlenode one in which the unstable manifold of the saddle is the stable manifold of the node ͑see region II of Fig. 1͒. After the bifurcation, a stable limit cycle remains as a successor of the previous heteroclinic connection ͑region III of Fig. 1͒. A third fixed point is needed ͑unstable focus͒ to feed the saddle point and the stable limit cycle before and after the Andronov bifurcation, respectively.
The equations we are studying are ͓9͔ xЈϭy, ͑1͒ with (x,y)R 2 , and ⑀ 1 , ⑀ 2 R ϩ . A partial bifurcation diagram for this system is displayed in Fig. 1. Within the excitable regime ͑region II of Fig. 1͒, an initial condition close to the node, subject to the action of noise, may cross the stable manifold of the saddle to relax after a long excursion in phase space. We associate this pulse with a dropout event of the laser intensity. In order to obtain excitable pulses in our model, we added a Gaussian white noise term in Eq. ͑1͒ with zero mean and variance 2Dϭ10 Ϫ3 .
The rate of escape from the node is ruled only by local properties, and can be described in terms of Kramers' formula ͓10͔. The theory also predicts a wide peak in the probability distribution, with a fast rise for short times and an exponential tail for the longer ones. The most probable escape time can be viewed as a first characteristic time of the system.
Usually, for excitable systems it is assumed that spikes are only triggered from a quiescent state. However, in our model it is possible to have an early-triggered spike if a trajectory crosses the stable manifold of the saddle at some point before reaching the node. In such a case the interspike time equals the time of flight of a trajectory following the unstable manifold of the saddle. This second characteristic time superimposes another peak in the histogram. Figure 2 shows a typical two-peaked histogram obtained from numerical simulations of the system ͑1͒ and ͑2͒ with parameter values ⑀ 1 ϭ0.26 and ⑀ 2 ϭ0.45. The leftmost peak is related to the early-triggered spikes and its importance decreases as we move far from the global bifurcation that limits regions I and II in Fig. 1 ͑for sufficiently high values of ⑀ 2 we recover the single-peaked distribution predicted by the usual escape rate theory͒. The rightmost peak corresponds to escapes from the node and has an exponential tail. The analysis of how interspike time distributions are changed as we move the parameter values will be reported elsewhere ͓11͔.
Then the excitability scenario is, in one of its simplest incarnations, highly compatible with the existence of two characteristic times in the interspike time distribution. From a geometrical point of view, the appearance of the twopeaked distribution is related to the existence of an invariant manifold that splits the flow. Specifically, in this model the stable manifold of the saddle acts as a separatrix between the early-triggered spikes and the spikes triggered from a quiescent state.

III. THREE-DIMENSIONAL EXCITABLE MODEL
In the same paper in which Giudici and co-workers proposed the excitability scenario, histograms and return maps of the time intervals between successive pulses were displayed. When dropouts are very rare the return map shows a cloud of points and the histogram is a single peak with an exponential tail but, for higher current values, certain char-FIG. 1. Partial bifurcation diagram and phase portraits for the system described by Eqs. ͑1͒ and ͑2͒. In regions I and II there are three fixed points: a node, a saddle, and a repulsor. Crossing the separatrix to region III, the saddle and the attractor collapse ͑Andronov bifurcation͒. The two lower regions display qualitatively different behavior. In region I the unstable manifold of the saddle approaches a limit cycle. In region II the unstable manifold of the saddle is the stable manifold of the node and the system behaves as an excitable one.
FIG. 2. Probability distribution of time between dropout events in the model described by Eqs. ͑1͒ and ͑2͒ for parameter values ⑀ 1 ϭ0.26 and ⑀ 2 ϭ0.45 and with Gaussian noise added in the x variable with zero mean and variance 2Dϭ10 Ϫ3 . acteristic times arise, leading to a multipeaked histogram and a symmetric grid in the return map. From the experimental observations reported in Ref. ͓1͔ it is clear also that the time between the spots in the grid can be associated to a characteristic time of the experimental system: the external cavity round trip. Recently, ͓12͔ a high-resolution inspection of the dynamical evolution of the laser intensity was reported. The recovery after a dropout presents fast intensity pseudoperiodic pulses at the round-trip time of the external cavity. Thus, in order to reproduce the detailed statistics, we must make allowance for this fast-pulsing phenomenon.
In this section, we want to extend our dynamical picture to reproduce the appearance of a grid structure in the return map of the interspike times within the framework of excitability plus noise. Since we have to account for a characteristic time, which appears to be associated with fast pulsing oscillations, we have to move beyond our two-dimensional model. We present the derivation of an excitable system in three dimensions. This model also simulates the typical shape of averaged pulses in the LFF regime.
We first analyze the main features of the dynamical behavior of the laser intensity. Typical time series for this variable in the LFF regime was presented in ͓1,3,4͔. In order to clarify this behavior, we split the temporal variations of the laser intensity into three parts: ͑a͒ a slow recovery of the signal up to an almost constant value ͑buildup͒, ͑b͒ a variable period of time during which the intensity is nearly stationary in that value, and ͑c͒ the dropout process in which the signal drops in a short period of time. In a dynamical reconstructed system the buildup can be associated with an orbit approaching a fixed point. This orbit stays in a neighborhood of the fixed point due to a certain trapping process ͓part ͑b͔͒, and then rapidly escapes from this region ͓part ͑c͔͒ and is reinjected again through the buildup process. One can also note an oscillatory behavior near the fixed point, which suggests its rotational nature. These will be the building blocks of our dynamical sketch.
One feature that is clearly displayed in the LFF series is the alternation between fast and slow dynamics. This is very common to a wide variety of systems, from chemical reactions to electronic devices. The simplest example of a system with such an alternation of different time scales is a relaxation oscillator. Typically, there is a fast variable that relaxes to a z-shaped slow manifold, where the upper and lower branches are stable. In these branches the motion is ruled by the slow variables and the orbits can eventually reach the turning point of the z shape and jump to the other stable branch. If we choose a folded two-dimensional slow manifold as shown in Fig. 4, we can combine two planar flows ͑one at each stable branch of the z shape͒ with the fast switching mechanism. This idea was fully developed by Deng ͓13͔ to construct homoclinic orbits and chaotic attractors in three-dimensional flows. Here we use the simplest z-switching mechanism: a singular perturbation of a transverse crossing of two planes at zϭϮ1 and a diagonal plane at zϭxϩd. The set of equations reads where ⑀ is the singular parameter and d is a fixed parameter controlling the reinjection of the orbits. When ⑀ is small enough, the nullcline zЈϭ0 consists of the roots of a cubic polynomial in z which has a z shaped hysteresis cycle in the middle, as displayed in Fig. 3. The first two equations, Eqs. ͑3͒ and ͑4͒, govern the motion in the stable branches lying near the perturbed planes zϭϮ1. For a sufficiently weak ⑀ value, the motion in the upper ͑lower͒ stable branch is approximated by the vector field f(x,y) ͓g(x,y)͔.
In order to map our problem to this dynamical sketch, we relate the dropout and buildup processes to the two stable branches in the z-shaped manifold and take as our x variable a properly scaled variable proportional to the laser intensity. Then, as illustrated in Fig. 4, we have the following temporal ͑6͒ and ͑7͒ and on the lower stable branch the planar system ͑8͒ and ͑9͔͒. For a suitable parameter choice the whole system behaves as an excitable one: if we perturb an initial condition in P by an amount greater than the distance to the saddle S, the orbit relaxes after a long excursion in phase space following the heteroclinic connection R.
behavior ͑a͒ In the upper plane there is a slow recovery of the x variable until the orbit reaches the turning point on the slow manifold (xϭ1Ϫd). ͑b͒ The orbit falls in a vicinity of the fixed point and spends a variable amount of time before crossing towards the x negative direction. ͑c͒ The x variable rapidly drops until the orbit reaches the turning point at x ϭϪ1Ϫd and then turns back to the upper plane. Thus we have recovered our initial analysis for the temporal evolution of the laser intensity in terms of the system described by Eqs. ͑3͒-͑5͒. Now, we need an explicit form for f(x,y) and g(x,y) in order to simulate our caricature of the LFF process.
We choose the simpler form for the buildup of a variable: a linear flow converging to a stable node at (x,y)ϭ(b,0). f 1͑x,y ͒ϭϪa͑ xϪb ͒, ͑6͒ f 2͑x,y ͒ϭϪcy. ͑7͒ If we set bϾ1Ϫd, the fixed point is located beyond the turning point of the z switch. Thus, for a sufficiently small value of ⑀, every trajectory near the half plane ͓xϽ(1 Ϫd),zϭ1͔ moves towards the positive x direction until it reaches the edge of the z switch.
Deciding which is the better form for the planar flow in the lower branch ͓g(x,y)͔ is not a simple task. Here one must decide if the dropout process is a deterministic or a noise-driven one. There is no clear evidence of which bifurcations occur with the appearance of the LFF and with the qualitative change in the return map.
In the excitability scenario we only need one saddle point and one stable focus ͑located at a higher value of x) in the lower stable branch of our dynamical sketch, as shown in Fig. 4. Actually, with these two elements the requirements for the excitable system specified in the Introduction are satisfied ͑a͒ A stable equilibrium point exists ͑the focus P). ͑b͒ If we perturb ͑in the x variable͒ an orbit near the equilibrium point by an amount greater than the distance between the fixed points, this orbit may cross the stable manifold of the saddle (S), escape towards the turning point of the lower branch of the z switch, and return to the equilibrium state after a long excursion in phase space. ͑c͒ Further, as the great relaxation is originated by the z-switch mechanism it becomes independent of the size of the perturbation in the slow variable.
The experimental observations reported show that the rate of dropout events increases as the control parameter is raised until the oscillations becomes almost periodic. In our scenario this can be achieved going from an excitable regime to relaxation oscillations in the hysteresis cycle. This transition can be made through two different local bifurcations of the planar flow g͑x,y͒: a saddle-node bifurcation ͑the stable focus becomes a node before his collapse with the saddle, as discussed in the previous section͒ and a subcritical Hopf bifurcation ͑the focus loses its stability͒. In both cases the orbits lying on the lower stable branch are ejected towards the turning point after the bifurcation.
Actually, both bifurcations are global, due to the heteroclinic connection ͑R in Fig. 4͒. In the saddle-node case, after the merging of the fixed points, a limit cycle remains and we recover the Andronov bifurcation.
If we choose the subcritical Hopf bifurcation instead, the main features of the statistics of time intervals studied in the preceding section are retained but, as we will show in the next section, our dynamical sketch predicts also the grid in the return map. We first present a simple planar flow that undergoes a subcritical Hopf bifurcation and then justify this particular choice.

xЈϭy, ͑8͒
yЈϭxϩyϪ␦x 2 Ϫ␥xy. ͑9͒ Here, ␥ and ␦ are fixed real positive values ͓␦ controls the distance between the two fixed points, (x,y)ϭ(0,0) and (x,y)ϭ(1/␦,0)]. Our control parameter will be , which varies over all real values. Let us analyze the bifurcation points in our linear parameter space (). As we require, the fixed point at the origin is always a saddle while the other fixed point loses its stability at ϭ␥/␦. At this parameter value we find, after a suitable coordinate change, a Hopf bifurcation with a positive leading term. This means that there is some interval in parameter space ͑whose upper limit is ␥/␦) where an unstable limit cycle encircles the stable focus. This cycle collapses at the bifurcation with the fixed point and when Ͼ␥/␦, an unstable focus remains.
As one can check numerically, the two fixed points and the unstable limit cycle are the only limit sets of the planar system ͑8͒ and ͑9͒. Then we have at least three topologically inequivalent flows. These are shown in Fig. 5.
For low values of we have no limit cycle and one branch of the unstable manifold of the saddle feeds the stable focus ͑region I͒. In the three-dimensional model ͑3͒-͑5͒ associated with this flow, every orbit reinjected is captured by the focus. Only the noise can drive the orbit towards negative x values, so we are in an excitable regime. Note that in this case the stable manifold of the saddle represents the excitability threshold. Region II instead has an unstable limit cycle, which is connected with the stable manifold of the saddle. Within this region, the behavior of the threedimensional system is also excitable, but now the limit cycle acts as the threshold. As we will see later, the branch of the unstable manifold of the saddle lying in the xϾ0 half-plane embraces the unstable cycle, thus hindering the perturbed orbits and driving them towards the turning point of the z switch. This is a very interesting point because we are near a subcritical Hopf but there is no bistability ͑we do not need a stable limit cycle͒. Finally, after the Hopf we enter region III where the unstable focus ejects all reinjected orbits and we obtain sustained oscillations in the three-dimensional system.
It remains a last bifurcation point separating regions I and II. Since we have studied all local bifurcations, this must be a global one. Regarding the topologically inequivalent phase portraits associated with each region, it appears that the stable and unstable manifolds of the saddle point must cross at some parameter value. When this happens, we obtain a saddle loop; a homoclinic connection encloses the stable focus. Furthermore, as we cross the global bifurcation point from region II to region I, the unstable limit cycle grows in size and period until it reaches the saddle connection ͑which has infinite period͒. This also explains why the unstable manifold of the saddle embraces the limit cycle in region II.
One can explicitly calculate the approximate parameter value at which the global bifurcation takes place by means of the Melnikov function. In fact, it is possible to rewrite Eqs. ͑8͒ and ͑9͒ as a time-independent perturbation of a Hamiltonian problem. When the Melnikov function ͑written as a function of the parameters͒ vanishes, we would have a homoclinic connection ͓14͔. Hence, we obtain the value ϭ(6/7)(␥/␦). This is in good agreement with the real value at which the global bifurcation occurs. Now we return to our three-dimensional sketch. It can be checked that for a small enough ⑀ the equilibrium points and the main features of the planar system described above persist in a neighborhood of the lower branch of the z-shaped slow manifold. Indeed, inserting Eqs. ͑6͒-͑9͒ into Eqs. ͑3͒-͑5͒ we can derive analytical expressions for the nullclines and graphically demonstrate the persistence of the fixed points.

IV. RESULTS
We simulated the system ruled by Eqs. ͑3͒-͑9͒ adding Gaussian white noise in the x variable with zero mean and variance 2D. The parameter values are the following: a ϭ0. 3, bϭ0.3, cϭ1.2, dϭ0.6, ⑀ϭ0.1, ␦ϭ20, ␥ϭ1, and Dϭ10 Ϫ4 , and were chosen to fit the functional form of the experimental dropouts. In Fig. 6 we present the time evolu- FIG. 5. Bifurcation diagram and phase portraits for the planar system described by Eqs. ͑8͒ and ͑9͒. For Ͻ(6/7)(␥/␦) one branch of the unstable manifold of the saddle feeds the stable focus ͑region I͒. As is increased the branches of the stable and unstable manifolds lying on the xϾ0 half-plane approach each other until a saddle loop is formed. In region II an unstable limit cycle is born and it encloses the stable focus. When ϭ␥/␦ a subcritical Hopf bifurcation occurs. In region III the focus is unstable and the stable manifold of the saddle twists around it. tion of the x variable in the excitable region (ϭ0).
Next, we performed a statistical analysis of the time between dropouts, taking up to 10 6 events. As expected, in the excitable regime ͑regions I and II͒ we obtain a peaked histogram and a fuzzy return map ͑see Fig. 7͒. For a wide range of parameter values within region III ͑relaxation oscillations͒ we have obtained a histogram and a return map with welldefined characteristic times. This is displayed in Fig. 8 for ϭ0.1.
The easiest way to understand the appearance of these characteristic times is to look carefully what happens near the fixed point ͑an unstable focus͒. In Fig. 9 we show a close view of the vicinity of the focus for ϭ0.1. It is little after the subcritical Hopf bifurcation; then orbits spiral outwards from the unstable focus in a dense coil. As the noise introduces some uncertainty in the place where the orbit is reinjected, the number of twists that the orbit makes before the dropout is not always the same. Each characteristic time ͑each peak in the histogram͒ corresponds to a defined number of twists around the unstable focus.
It is crucial to note that without added noise we would return to a single-peaked histogram since we would have a fixed number of twists and a constant period. When we add a minute amount of noise we begin to observe different spots in the return map, as the orbits can make a few more or a few less twists before the dropout. The fact that the spot distribution is symmetric and that there are no other structures aside from the characteristic times is a clear indication that noise is still playing a non-negligible role in the time statistics. We are only viewing the old fuzzy spot through a grid of privileged times.
Then we claim that the statistical properties of the time between pulses in the LFF regime can be explained as noise ͑or a form of high-dimensional chaos͒ acting on a very simple deterministic structure, which is excitable in some parameter region. The bifurcation which limits this excitable region, can be a saddle-node one ͑an Andronov bifurcation, globally speaking͒ if we are only interested in the longer times scales, but a subcritical Hopf as we want to reproduce the detailed structure of the interspike histograms.
Finally, we want to point out that the model analyzed in this section does indeed allow us to reproduce also the bimodal structure of the interspike histograms observed in the analysis of the two-dimensional toy model discussed in Sec. II. In Fig. 10 we display a histogram displaying this feature for the parameter values in the caption. Notice that in this case, in the approximation for the planar flow in the lower plane, we are in parameter region I of Fig. 5. Thus, the invariant manifold that splits the flow into the ones associated with the fast and slow pulses, is the stable manifold of the saddle ͑a bidimensional manifold in our three-dimensional system͒. In this parameter region we still have rapid oscillations around the stable focus and, consequently, a collection of characteristic times in the distribution. But, far enough from the Hopf bifurcation, the statistical trace of these oscillations is washed out by the binning in the histogram.

V. CONCLUSIONS
In the last few years, the semiconductor laser with feedback has illustrated one of the fundamental problems in nonlinear dynamics: how to refute or build confidence in a model for an experiment displaying complex dynamical behavior. This system presents, for a region of its parameter space, a regime called low frequency fluctuations. Two different scenarios were proposed to explain this regime: noisedriven excitability ͑at least when the regime first arises͒ and high-dimensional chaos ͑eventually contaminated with noise͒. In this paper, we identified as a key observational quantity the interspike time distribution of dropouts, and inspected simple noise-driven excitable models.
In order to explain the distribution of pulses at the time scale of the dropout events ͑bimodal distributions͒, a simple two-dimensional noise-driven excitable model, consistent with hypothesis in Ref. ͓1͔, was able to give satisfactory results. Elsewhere, we will report the changes that this distribution undergoes as parameters are changed, and how it compares with the changes that are observed in the experiments when feedback and injection current are varied ͓11͔.
The explanation of the fine structure of the time distribution of the pulses ͑to account for the structure present at the round trip time scale͒ required us to extend the model to a three-dimensional one. The model could also reproduce the bimodal structure of the interspike time histograms at the scale of the dropout events.
Finally, let us stress that we do not mean to replace any set of equations derived from first principles, but to highlight the topological elements that are needed dynamically to reproduce the reported statistical observables. We found that the bimodal structure of the interspike histograms can be explained in terms of an invariant manifold ͑separatrix͒ that splits the flow. Part of it revisits the attractor ͑assumed in any excitable scenario͒ for long times before being expelled, and part of it rapidly escapes along the unstable manifold of the separatrix. This behavior can be found in both models presented in this paper. To account for the grid observed in the time returns at the scales of the round trip, we needed the occurrence of a subcritical Hopf bifurcation near a saddle point. The folding in the stable manifold of the saddle is a reminder of that bifurcation.