The statistics of chaos in a bursting laser

TL;DR

This paper experimentally studies how double optical feedback in semiconductor lasers induces a transition from Gaussian to power-law distributed chaotic bursting, with implications for controlling complex laser dynamics.

Contribution

It reveals the statistical transition mechanisms in laser chaos due to feedback ratio changes and links them to Farey fractions and pump current variations.

Findings

Chaotic spiking statistics shift from Gaussian to power-law distributions.

Transitions to bursting are near low-order Farey fractions.

Numerical models confirm the experimental chaos control possibilities.

Abstract

We demonstrate experimentally how semiconductor lasers subjected to double optical feedback change the statistics of their chaotic spiking dynamics from Gaussian to long-tail Power Law distributions associated to the emergency of bursting. These chaotic regimes, which are features of excitable complex systems, are quantified by the tail exponent $\alpha$ and appear by changing the ratio between the feedback times. Transitions to bursting occur in the neighbourhood of low order Farey fractions. The physics behind these transitions is related to the variation of threshold pump current in the compound system as obtained from a deterministic set of rate equations. Numerical integration also verifies the observed chaos transitions indicating the possibility of controlling the bursting chaotic statistics.

Figures

b presents a tail index $\alpha=1.8$, which is close to the Gaussian value. Numerical-theoretical results from Fig. LABEL:sub@Fig2d for the histogram present a tail index $\alpha=2.0$. The similarity between the experiments and the theory is therefore verified in the statistics of spikes. Scanning the $L_{2}$ cavity length with a single feedback over the range from $0.90$ m to $8.02$ m gave us a large collection of time series from which we extracted the average $\langle T\rangle$ and standard deviation $\sigma_{T}$ of the inter-spike time as functions of the feedback time. Smooth, nearly cons

b presents a tail index $\alpha=1.8$, which is close to the Gaussian value. Numerical-theoretical results from Fig. LABEL:sub@Fig2d for the histogram present a tail index $\alpha=2.0$. The similarity between the experiments and the theory is therefore verified in the statistics of spikes. Scanning the $L_{2}$ cavity length with a single feedback over the range from $0.90$ m to $8.02$ m gave us a large collection of time series from which we extracted the average $\langle T\rangle$ and standard deviation $\sigma_{T}$ of the inter-spike time as functions of the feedback time. Smooth, nearly cons

b presents a tail index $\alpha=1.8$, which is close to the Gaussian value. Numerical-theoretical results from Fig. LABEL:sub@Fig2d for the histogram present a tail index $\alpha=2.0$. The similarity between the experiments and the theory is therefore verified in the statistics of spikes. Scanning the $L_{2}$ cavity length with a single feedback over the range from $0.90$ m to $8.02$ m gave us a large collection of time series from which we extracted the average $\langle T\rangle$ and standard deviation $\sigma_{T}$ of the inter-spike time as functions of the feedback time. Smooth, nearly cons

b presents a tail index $\alpha=1.8$, which is close to the Gaussian value. Numerical-theoretical results from Fig. LABEL:sub@Fig2d for the histogram present a tail index $\alpha=2.0$. The similarity between the experiments and the theory is therefore verified in the statistics of spikes. Scanning the $L_{2}$ cavity length with a single feedback over the range from $0.90$ m to $8.02$ m gave us a large collection of time series from which we extracted the average $\langle T\rangle$ and standard deviation $\sigma_{T}$ of the inter-spike time as functions of the feedback time. Smooth, nearly cons

b presents a tail index $\alpha=1.8$, which is close to the Gaussian value. Numerical-theoretical results from Fig. LABEL:sub@Fig2d for the histogram present a tail index $\alpha=2.0$. The similarity between the experiments and the theory is therefore verified in the statistics of spikes. Scanning the $L_{2}$ cavity length with a single feedback over the range from $0.90$ m to $8.02$ m gave us a large collection of time series from which we extracted the average $\langle T\rangle$ and standard deviation $\sigma_{T}$ of the inter-spike time as functions of the feedback time. Smooth, nearly cons

b presents a tail index $\alpha=1.8$, which is close to the Gaussian value. Numerical-theoretical results from Fig. LABEL:sub@Fig2d for the histogram present a tail index $\alpha=2.0$. The similarity between the experiments and the theory is therefore verified in the statistics of spikes. Scanning the $L_{2}$ cavity length with a single feedback over the range from $0.90$ m to $8.02$ m gave us a large collection of time series from which we extracted the average $\langle T\rangle$ and standard deviation $\sigma_{T}$ of the inter-spike time as functions of the feedback time. Smooth, nearly cons

c and LABEL:sub@Fig3e. Varying $R$ further, $\langle T\rangle$ becomes larger as the quiescent time between bursts increases. It may even reach the condition where the intensity time series presents no LFF drops. Apart from LFF dynamics and statistics, the laser goes through many different forms of oscillation including local chaos around the MGM, periodic oscillation and continuous wave operation. None of these regimes is inspect in our study. To obtain $\langle T\rangle$ we have used a LFF pulse counting procedure and therefore, any information about the real dynamical state of the intervals

c and LABEL:sub@Fig3e. Varying $R$ further, $\langle T\rangle$ becomes larger as the quiescent time between bursts increases. It may even reach the condition where the intensity time series presents no LFF drops. Apart from LFF dynamics and statistics, the laser goes through many different forms of oscillation including local chaos around the MGM, periodic oscillation and continuous wave operation. None of these regimes is inspect in our study. To obtain $\langle T\rangle$ we have used a LFF pulse counting procedure and therefore, any information about the real dynamical state of the intervals

c and LABEL:sub@Fig3e. Varying $R$ further, $\langle T\rangle$ becomes larger as the quiescent time between bursts increases. It may even reach the condition where the intensity time series presents no LFF drops. Apart from LFF dynamics and statistics, the laser goes through many different forms of oscillation including local chaos around the MGM, periodic oscillation and continuous wave operation. None of these regimes is inspect in our study. To obtain $\langle T\rangle$ we have used a LFF pulse counting procedure and therefore, any information about the real dynamical state of the intervals

c and LABEL:sub@Fig3e. Varying $R$ further, $\langle T\rangle$ becomes larger as the quiescent time between bursts increases. It may even reach the condition where the intensity time series presents no LFF drops. Apart from LFF dynamics and statistics, the laser goes through many different forms of oscillation including local chaos around the MGM, periodic oscillation and continuous wave operation. None of these regimes is inspect in our study. To obtain $\langle T\rangle$ we have used a LFF pulse counting procedure and therefore, any information about the real dynamical state of the intervals

c and LABEL:sub@Fig3e. Varying $R$ further, $\langle T\rangle$ becomes larger as the quiescent time between bursts increases. It may even reach the condition where the intensity time series presents no LFF drops. Apart from LFF dynamics and statistics, the laser goes through many different forms of oscillation including local chaos around the MGM, periodic oscillation and continuous wave operation. None of these regimes is inspect in our study. To obtain $\langle T\rangle$ we have used a LFF pulse counting procedure and therefore, any information about the real dynamical state of the intervals

c and LABEL:sub@Fig3e. Varying $R$ further, $\langle T\rangle$ becomes larger as the quiescent time between bursts increases. It may even reach the condition where the intensity time series presents no LFF drops. Apart from LFF dynamics and statistics, the laser goes through many different forms of oscillation including local chaos around the MGM, periodic oscillation and continuous wave operation. None of these regimes is inspect in our study. To obtain $\langle T\rangle$ we have used a LFF pulse counting procedure and therefore, any information about the real dynamical state of the intervals

c for $R=0.95$ and $\kappa=0.5\times 10^{9}$ s*-1*. In the calculations we used $\kappa_{1}=\kappa_{2}=\kappa$, which means that both arms of optical feedback are sending light back to the gain medium with the same strength. Experimentally, to provide that, equal laser threshold reduction due to each optical feedback alone was monitored and adjusted through mirrors alignment. It is important to notice that this pattern is not a result of optical frequency interferometry, once $\Delta\tau$ is very large if compared to the optical period as in our experimental and theoretical scan of $R$. Figure

c for $R=0.95$ and $\kappa=0.5\times 10^{9}$ s*-1*. In the calculations we used $\kappa_{1}=\kappa_{2}=\kappa$, which means that both arms of optical feedback are sending light back to the gain medium with the same strength. Experimentally, to provide that, equal laser threshold reduction due to each optical feedback alone was monitored and adjusted through mirrors alignment. It is important to notice that this pattern is not a result of optical frequency interferometry, once $\Delta\tau$ is very large if compared to the optical period as in our experimental and theoretical scan of $R$. Figure

c for $R=0.95$ and $\kappa=0.5\times 10^{9}$ s*-1*. In the calculations we used $\kappa_{1}=\kappa_{2}=\kappa$, which means that both arms of optical feedback are sending light back to the gain medium with the same strength. Experimentally, to provide that, equal laser threshold reduction due to each optical feedback alone was monitored and adjusted through mirrors alignment. It is important to notice that this pattern is not a result of optical frequency interferometry, once $\Delta\tau$ is very large if compared to the optical period as in our experimental and theoretical scan of $R$. Figure

c for $R=0.95$ and $\kappa=0.5\times 10^{9}$ s*-1*. In the calculations we used $\kappa_{1}=\kappa_{2}=\kappa$, which means that both arms of optical feedback are sending light back to the gain medium with the same strength. Experimentally, to provide that, equal laser threshold reduction due to each optical feedback alone was monitored and adjusted through mirrors alignment. It is important to notice that this pattern is not a result of optical frequency interferometry, once $\Delta\tau$ is very large if compared to the optical period as in our experimental and theoretical scan of $R$. Figure

b, gives a structure to the MGM threshold current as function of $R$. Its reduction has a base line equal to the single feedback case apart from the regions where $R$ is near a low order Farey number. The changes in $J_{th,MGM}$ is more significant with peaks according to the Farey order $p/q$ in $R$. Throughout the range $0<R<1$ equal Farey denominators $q$ give similar peaks for the threshold reduction. For the purpose of comparison with Fig. LABEL:sub@Fig5b, the sixth order Farey sequence (all the fractions p/q with $q\leq 6$) is shown in Fig. LABEL:sub@Fig5c. The vertical axis of Fig. LABE

b, gives a structure to the MGM threshold current as function of $R$. Its reduction has a base line equal to the single feedback case apart from the regions where $R$ is near a low order Farey number. The changes in $J_{th,MGM}$ is more significant with peaks according to the Farey order $p/q$ in $R$. Throughout the range $0<R<1$ equal Farey denominators $q$ give similar peaks for the threshold reduction. For the purpose of comparison with Fig. LABEL:sub@Fig5b, the sixth order Farey sequence (all the fractions p/q with $q\leq 6$) is shown in Fig. LABEL:sub@Fig5c. The vertical axis of Fig. LABE

b, gives a structure to the MGM threshold current as function of $R$. Its reduction has a base line equal to the single feedback case apart from the regions where $R$ is near a low order Farey number. The changes in $J_{th,MGM}$ is more significant with peaks according to the Farey order $p/q$ in $R$. Throughout the range $0<R<1$ equal Farey denominators $q$ give similar peaks for the threshold reduction. For the purpose of comparison with Fig. LABEL:sub@Fig5b, the sixth order Farey sequence (all the fractions p/q with $q\leq 6$) is shown in Fig. LABEL:sub@Fig5c. The vertical axis of Fig. LABE

a. Those behavior do not correspond to optical interferometric oscillations. They have to do with the difference between the two feedback time delays as discussed in the previous sections. Such oscillations have an impact in the chaotic statistics. Experimental and numerically calculated mean inter-spike times $\langle T\rangle$ through the same range of Fig. LABEL:sub@Fig6a are shown in Figs. LABEL:sub@Fig6b and LABEL:sub@Fig6c, respectively.

a. Those behavior do not correspond to optical interferometric oscillations. They have to do with the difference between the two feedback time delays as discussed in the previous sections. Such oscillations have an impact in the chaotic statistics. Experimental and numerically calculated mean inter-spike times $\langle T\rangle$ through the same range of Fig. LABEL:sub@Fig6a are shown in Figs. LABEL:sub@Fig6b and LABEL:sub@Fig6c, respectively.

a. Those behavior do not correspond to optical interferometric oscillations. They have to do with the difference between the two feedback time delays as discussed in the previous sections. Such oscillations have an impact in the chaotic statistics. Experimental and numerically calculated mean inter-spike times $\langle T\rangle$ through the same range of Fig. LABEL:sub@Fig6a are shown in Figs. LABEL:sub@Fig6b and LABEL:sub@Fig6c, respectively.

a. This is shown if Fig. LABEL:sub@Fig7b. The numerical range of $J$ variation was 5% of the solitary threshold current $J_{th,sol}$ with steps of 0.25$\%$ of this value. A corresponding experimental situation of a laser with solitary threshold current of 20 mA would fall in the fairly controllable 50 $\mu$A precision in current tuning. The long-tail index $\alpha$ drops from 1.7 (near Gaussian distribution) to 0.6, which signals very long tail statistics, and then recovers back to near the Gaussian statistics value. Graphics (i), (ii) and (iii) in Fig. LABEL:sub@Fig7c illustrate the nature of

a. This is shown if Fig. LABEL:sub@Fig7b. The numerical range of $J$ variation was 5% of the solitary threshold current $J_{th,sol}$ with steps of 0.25$\%$ of this value. A corresponding experimental situation of a laser with solitary threshold current of 20 mA would fall in the fairly controllable 50 $\mu$A precision in current tuning. The long-tail index $\alpha$ drops from 1.7 (near Gaussian distribution) to 0.6, which signals very long tail statistics, and then recovers back to near the Gaussian statistics value. Graphics (i), (ii) and (iii) in Fig. LABEL:sub@Fig7c illustrate the nature of

a. This is shown if Fig. LABEL:sub@Fig7b. The numerical range of $J$ variation was 5% of the solitary threshold current $J_{th,sol}$ with steps of 0.25$\%$ of this value. A corresponding experimental situation of a laser with solitary threshold current of 20 mA would fall in the fairly controllable 50 $\mu$A precision in current tuning. The long-tail index $\alpha$ drops from 1.7 (near Gaussian distribution) to 0.6, which signals very long tail statistics, and then recovers back to near the Gaussian statistics value. Graphics (i), (ii) and (iii) in Fig. LABEL:sub@Fig7c illustrate the nature of

Equations

$$\frac{dE(t)}{dt}$$

$$\frac{d\Phi (t)}{dt}$$

$$\frac{dN(t)}{dt}$$

$$\varphi (t)=exp(it\delta -|\gamma t{|}^{\alpha }(1-i\beta sgn(t)\Phi ))$$

$$f(x)=$$

$${\omega }_s-{\omega }_0=-{\kappa }_1[\overline{\alpha }\cos ({\omega }_s{\tau }_1)+\sin ({\omega }_s{\tau }_1)]$$

$$-{\kappa }_2[\overline{\alpha }\cos ({\omega }_s{\tau }_2)+\sin ({\omega }_s{\tau }_2)]$$

$$N_s$$

$$E_s$$

Full text

I Introduction

Within optics, the non-linear dynamics of lasers has been studied for a long time, both for fundamental Ott (2002) and applied physics Sciamanna and Shore (2015) interests. Measuring a laser output power is a well developed procedure to determine the amplitude squared of the dynamical variable which is the emitted electric field. Furthermore, the ease of changing parameters in the laboratory permits the observation of large dynamical diversity. Another advantage of dynamical studies in lasers arrangements is the existence of model equations to compare and predict experimental results. So, beyond their interest for themselves they offer examples for analogies in the understanding of others complex systems. Among relevant examples are the dynamical systems with time delayed feedback in biology, as described by Mackey-Glass equations Mackey and Glass (1977).

Excitability, which is an important concept in the study of complex transitory states like in neurons Rinzel and Ermentrout (1989); Izhikevich (2000); Laing et al. (2003); Krahe and Gabbiani (2004), was also addressed in laser dynamics Giudici et al. (1997); Mos et al. (2000); Romariz and Wagner (2007); Hurtado et al. (2010); Coomans et al. (2011); Selmi et al. (2014); Méndez et al. (2005); Martinez Avila et al. (2008). Delayed optical feedback has proved to induce excitable Giudici et al. (1997); Eguia and Mindlin (1999) and chaotic behavior Kim et al. (2015); Masoller and Abraham (1998) on semiconductor lasers. It is well known that these lasers with optical feedback generate chaos which, after an averaging process due to detector bandwidth, manifests with irregular sharp drops in its intensity, usually called low frequency fluctuation (LFF) Risch and Voumard (1977). The statistical properties of these LFF pulses have been extensively investigated Hohl et al. (1995); Sukow et al. (1997); Mulet and Mirasso (1999); Aragoneses et al. (2013).

Two simultaneous different feedbacks are simple schemes to get bursting of power drops with statistics on-demand as we will show. Indeed, the addition of a second time delay feedback on a laser was already proposed and studied, using two external optical cavities Rogister et al. (1999, 2000); Ruiz-Oliveras and Pisarchik (2006); Wu et al. (2009). High dimensional chaotic dynamic was attained experimentally using such configuration Fischer et al. (1994). With this scheme and varying the lengths of the cavities and the intensity of the feedbacks, chaos suppression was also achieved Liu and Ohtsubo (1997). Regions with stable emission and periodical behavior were observed Liu and Ohtsubo (1997), demonstrating the system high sensitivity to double feedback. More than two feedbacks has also been addressed in the literature Többens and Parlitz (2008) aiming at laser stabilization.

Our experiments reported here were done in a semiconductor laser, chaotic by the effect of double optical feedback. Measured laser intensity time series are the essential data in our work. The statistical characterization of the chaotic dynamics is done extracting the inter-spike times represented in histograms and obtaining the so-called tail index ${\bf\alpha}$, for alpha-stable distributions. It is shown how parameter tuning in the system makes the inter spike time go from chaotic with Gaussian statistics to chaos with heavy tail Lévy distributions among pulse bursting. To corroborate the experiments we also solved numerically the rate equations for a mono-mode semiconductor laser including more than one feedback time. We not only verify semi-quantitative experimental aspects of our results but we get evidence that deterministic non-linear equations with multiple feedback times can account for the production of bursting spike dynamics followed by switching between chaotic, periodic, and stable states.

II Experimental Setup

Our experimental data was obtained in the form of long time series for the fluctuating laser intensity obtained using the setup shown in Fig. 1. The experiments were done with commercial single mode diode lasers stabilized to 0.01 K, emitting in the infra-red wavelength and with solitary threshold current $J_{th}$ around 20 mA giving output power of few mW. Two specific units were SDL 5401 GaAlAs at 800nm - $J_{th}=26$ mA and Thorlabs L850P010 at 850nm - $J_{th}=10$ mA. As indicated in Fig. 1 the beam splitter BS1 allows the laser output intensity to be measured by a 2.0 GHz bandwidth photodiode PD and the data series collected with a Tektronix DPO 7104, 1 GHz oscilloscope. The two external cavities are created by the $50/50$ beam splitter BS2 and two mirrors $M_{1}$ and $M_{2}$. The fixed mirror $M_{1}$ create the first external cavity with a constant laser light pathway $L_{1}=8.02$ m corresponding to a fixed feedback time delay, $\tau_{1}=$53.5 ns, due to light round trip time. The second external cavity, created by the variable position mirror $M_{2}$, has a variable light pathway $L_{2}$, giving the second time delay, $\tau_{2}$, which is scanned during the experiments.

Experimental cavity alignment conditions were maintained using an aspheric collimating lens $l_{1}$ for beam divergence reduction. Operation current was kept constant close to the solitary laser threshold, but up to 3% current threshold reduction due to feedback could be achieved in both arms. Optimum alignment condition in the two external arms was imposed by inserting AR coated lenses $l_{2}$ and $l_{3}$ in front of each feedback mirror, in order to get mode matching. Consequently, maximum output intensity for lengths below 0.90 m and with one single feedback, remains in stationary conditions (without spikes). In this configuration, semiconductor laser with optical feedback exhibits physical properties of excitable systems when noise or any perturbation is applied Giudici et al. (1997). The best alignment was maintained during the variation of $L_{2}$ for lengths over 0.90 m.

III Theoretical Model

The dynamics of a semiconductor laser under feedback operation is well described by the so called Lang-Kobayashi equations Lang and Kobayashi (1980). This model has been previously used in order to include additional delayed term Fischer et al. (1994); Liu and Ohtsubo (1997); Rogister et al. (1999); Ruiz-Oliveras and Pisarchik (2006). Here, we use the deterministic form of these equations for a system composed by two external cavities which provides a double feedback:

[TABLE]

To be noticed is the absence of stochastic terms in the equations. Thus all random type dynamics results from the deterministic nonlinear nature of the model which reproduce the main dynamics of these kind of systems even to obtain chaotic behavior. Implementation of a fourth order Runge-Kutta routine with an integration step of $\delta t=10^{-13}$ s were done in C++. This value for $\delta t$ was shown to be small enough for the equations integrations during tests varying its value. In the equations $i=1,2$ is the number of external cavities and $G_{N}=G_{0}\left[N(t)-N_{0}\right]/\left[1+\epsilon|E(t)|^{2}\right]$ is the medium gain. The optical frequency when the laser is solitary (with no feedback) is represented by $\omega_{0}$ and $\Theta_{i}=\omega_{0}\tau_{i}+\phi(t)-\phi(t-\tau_{i})$ is the field phase excursion after one external cavity round trip. The dynamical variables are the slow (with respect to the optical period) envelope for the electric field, having $E(t)$ as adimensional amplitude, $\Phi$(t) as the associated phase and $N(t)$ the carrier density. The parameter in the equations are: $\tau_{p}=1/\gamma_{p}$ and $\tau_{s}=1/\gamma_{s}$ corresponding to photon and carrier lifetime, respectively; $G_{0}$ the small signal gain coefficient, $N_{0}$ representing the carrier density at the transparency of the medium, $J$ is the injected pump current ($J_{th}$ represents threshold current), $\epsilon$ the saturation coefficient, $\kappa_{i}$ the feedback coefficients, $\tau_{i}$ the feedback delay times and $\overline{\alpha}$ the line width enhancement factor (not to be confused with the statistical index $\bf{\alpha}$). The parameters values used were $\tau_{s}=2.0$ ns, $\tau_{p}=3.55$ ps, $\epsilon=5\times 10^{-7}$, $G_{0}=3.2\times 10^{3}$ s*-1*, $N_{0}=1.5\times 10^{8}$, $\kappa_{1}=\kappa_{2}=11\times 10^{9}$ s*-1* and $\overline{\alpha}=3.0$.

IV Statistics of Long-Tail Alpha-Stable Lévy Distribution

Observation of very different statistical distributions in our system naturally lead us to examine the well known alpha-stable type of distribution which present a good formalism to characterize shapes of distributions with respect to their location, scale, distortion and tail Nolan (2018); McCulloch (1986). The experimental and numerical-theoretical data in our work were assumed to have stable distribution.

A general stable probability distribution, also known as Lévy alpha-stable, stable Paretian or $\alpha$-stable distribution, is defined by the invariance property of a sum between two independent random variables with the index of stable distribution equal to $\alpha$ having as the result an $\alpha$-distribution with the same index. Often, the probability density function (PDF) $f(x)$ for a stable distribution does not have a general formula and some times it is not even analytical. The description of those distributions is better done by the characteristic function $\phi(t)$. The PDF, when it exists, is given by the Fourier transform of $\phi(t)$ which for stable distributions can be written as:

[TABLE]

where i is the imaginary number, $sgn(t)$ is the signal of $t$ and $\Phi=-\frac{2}{\pi}\log{|t|}$ if $\alpha=1$, whereas $\Phi=\tan{\left(\frac{\pi\alpha}{2}\right)}$ if $\alpha\neq 1$. The stable parameter $\alpha$ is, for our data analysis, the most important among the four that characterize these distributions. It is also known as the tail index, index of stability or the characteristic exponent. Stable distributions have $0<\alpha\leq 2$ and it reveals the concentration of the distribution and its tail character. The upper limit, $\alpha=2$, correspond to the Gaussian case while long tail distributions have $\alpha<2$. Cauchy or Lorentz distributions have $\alpha=1$. Variance is undefined for $\alpha<2$ and mean does not exist for $\alpha\leq 1$. The others three parameters that caracterize an $\alpha$-stable distribution are related to the location $\delta$ $\in[-\infty,+\infty]$, to the scale $\gamma$ $\in[0,\infty]$ and to the skewness $\beta$ $\in[-1,+1]$. For the Gaussian case, i.e., $\alpha=2$ and $\beta=0$ due to the symmetry of this distribution, $f(x)$ can be obtained as it follows:

[TABLE]

A comparison with the well known PDF formula for Gaussian distribution allows us to relate its mean $\mu$ and variance $\sigma$ to the location and scale parameters $\delta=\mu$ and $\gamma=\sigma/\sqrt{2}$, respectively. Once we had histograms of probability distributions for time interval between intensity spikes, we estimated $\alpha$ directly from this data using the software developed by Nolan Nolan . It calculates the four parameters for stable-like distribution with an algorithm is based on the quantile fit process of McCulloch McCulloch (1986). In our data, the occurrence of spike-bursting causes the deviation of $\alpha$ from the Gaussian condition to the long-tail distributions with $\alpha$ much close to 1 than 2.

V Experimental and Theoretical Results

V.1 Chaos with Single Feedback

Chaotic spiking power drops in semiconductor laser with single feedback has been extensively studied Sciamanna and Shore (2015). We present its statistical properties here for the sake of comparison with the case of double feedback to be discussed later on. A typical segment of an experimental series with a single feedback is shown in Fig. LABEL:sub@Fig2a. Power drops with a gaussian distribution is revealed in the corresponding histogram of Fig. LABEL:sub@Fig2b. Corresponding numerical-theoretical calculations are given in the time series segment of Fig. LABEL:sub@Fig2c. The related histogram, in Fig. LABEL:sub@Fig2d, also looks like a Gaussian probability distribution as in the experimental result.

From the histograms with more than $10^{4}$ drop events, we calculate the tail index $\alpha$ that represents the tail decay of the distributions McCulloch (1986). To remind, pure Gaussian distributions have $\alpha=2$ while the long tail Lorentz distribution has $\alpha=1$. The experimental probability distribution of Fig. LABEL:sub@Fig2b presents a tail index $\alpha=1.8$, which is close to the Gaussian value. Numerical-theoretical results from Fig. LABEL:sub@Fig2d for the histogram present a tail index $\alpha=2.0$. The similarity between the experiments and the theory is therefore verified in the statistics of spikes. Scanning the $L_{2}$ cavity length with a single feedback over the range from $0.90$ m to $8.02$ m gave us a large collection of time series from which we extracted the average $\langle T\rangle$ and standard deviation $\sigma_{T}$ of the inter-spike time as functions of the feedback time. Smooth, nearly constant, dependence for $\langle T\rangle$ was obtained collecting $120$ equally separated runs moving $L_{2}$ from approximately $0.90$ m to $8.02$ m while $L_{1}$ was kept constant and equal to $8.02$ m during the entire scan as shown in Fig. LABEL:sub@Fig2e. The numerical solution of the rate equations resulted with the same behavior and is given in Fig. LABEL:sub@Fig2f. The points in these graphics represent the mean time between consecutive LFF drops for time series of intensity with more than $10^{4}$ events each. The vertical bars represent the standard deviation. In both experimental and numerical data the $\alpha$ index was always near $1.9$, confirming Gaussian shape for the histograms, independent of feedback time in the studied range. Small feedback times ($L_{2}<0.50$ m) are known to make the system depart from LFF dynamics Heil et al. (2001). The parameters of the equations were taken to produce the observed quantitative coincidence with the experimental values and are shown in the methods section. These results for a single feedback dynamics were taken for comparison with the striking new results from the double feedback case, to be presented next.

V.2 Chaos with Double Feedback

The inclusion of the second feedback in the laser changes drastically its dynamics as compared to what was reported experimental and theoretically for one delay. For the double feedback experiments the fixed length arm $L_{1}$, with feedback time $\tau_{1}$, is unblocked and give the same (as good as possible) amount of field return established for the variable arm $L_{2}$. Repeating the $120$ experimental steps for the values of $\tau_{2}$, like was done in the one feedback case presented in Fig. LABEL:sub@Fig2e, completely different scenarios are found. This is shown in Fig. LABEL:sub@Fig3a which gives the measured dependence of the average time $\langle T\rangle$ as function of the feedback lengths ratio $R$. The length increments $\Delta L_{2}$ correspond to $\approx$ 6 cm. The theoretical counterpart, calculated with steps of 1.5 cm, is shown in Fig. LABEL:sub@Fig3b.

One clearly see the sensitivity of $\langle T\rangle$ on this parameter. When $R$ is close to rational numbers with small denominators, the low frequency fluctuation drops of intensity become less frequent due to the appearance of long intervals without drops. As a consequence, the mean time between drops $\langle T\rangle$ increases and the data variance diverges. Notice that in both, experimental and numerical, the effect observed occurs with external cavity length steps $\approx$ $\times 10^{4}$ optical wavelengths. Farey sequences, defined as the set of rational numbers written as the ratio $p/q$ of two natural numbers, are known to manifest in frequency locking or entrainment between oscillators. These physical phenomena deal mostly with periodic, quasiperiodic motions and the onset of chaos Xu et al. (1990). In our results the Farey hierarchy plays a critical role in changing chaos statistics. Bunching of spikes appear around the small denominator ratio, also called low order Farey numbers, in a very sensitive way with respect to the ratio $R$. Those transitions from statistic behavior with Gaussian shape to spiking burst long-tail distributions are characterized by peaks of $\langle T\rangle$ near Farey fractions during the scan of $R$. The most significant peak in value of $\langle T\rangle$ happen around $R=1/2$. Outside those Farey ratios we have LFF chaos with $\langle T\rangle$ and $\sigma_{T}$ quite insensitive to the changes in $R$. This indicates Gaussian profiles for the inter-spikes time histogram. As $R$ increase towards $1/2$ (similar to what happens near 1/1 and higher order Farey ratios) $\langle T\rangle$ and the respective variance increase drastically indicating transitions to bursting regime of spikes. Crossing the value $R=1/2$ these statistics parameters decrease again returning to the Gaussian chaos. As will be shown latter on, details of those transitions, including stationary regime and chaos without LFF drops, also occur for values of $R$ very close to these special ratios. Figure LABEL:sub@Fig3c presents an interval of the experimental time series in the case of spiking burst behavior for a length ratio $R=L_{2}/L_{1}\approx 0.498$, i. e., near $1/2$. The corresponding histogram of the $10^{4}$ inter-spike times is shown in Fig. LABEL:sub@Fig3d. Unlike the Gaussian distribution presented in the single feedback, and also in the double feedback case with $R$ far from low order Farey numbers (1/2, 1/3, 2/3, 1/1, etc), in this case a heavy tail distribution is observed. For this experimental data the calculated long tail index is $\alpha=1.1$, confirming the presence of an expressive long tail in the probability distribution $P(T)$ due to the long time interval between bunches of spikes. Figs. LABEL:sub@Fig3e and LABEL:sub@Fig3f present the numerical calculation for $R=0.4946$ showing the same non-Gaussian distribution and spiking burst behavior as in the experimental results. The associated histogram, also composed by approximately $10^{4}$ intervals, also has a long tail index $\alpha=1.1$.

Definitively the model equations fit the transitions of chaos statistics observed. Better quantitative comparison would ask for further exploration of the equations parameters. In contrast to the double feedback results, experiments and calculations with single feedback (with our selected parameters) never show spiking bursts throughout the whole range of variation of $L_{2}$ as indirectly indicated in Figs. LABEL:sub@Fig2e and LABEL:sub@Fig2f. It is important to reinforce that, different from other previously published works which include stochastic terms in the models, our equations are fully deterministic and so the theoretical results depend only on the system parameters.

VI Discussion

VI.1 Stationary Solutions and Bifurcations

The stationary solutions ($E_{s}$, $N_{s}$, $\omega_{s}$) for the Eqs. (1)-(3) are obtained setting $\frac{dE(t)}{dt}=0$, $\frac{d\Phi(t)}{dt}=(\omega_{s}-\omega_{0})$, $\frac{dN(t)}{dt}=0$ and using the condition $E(t)=E(t-\tau_{1})=E(t-\tau_{2})=E_{s}$. While $\omega_{0}$, the solitary laser frequency, is always constant, $\omega_{s}$, the frequency of the compound system including the two cavities, vary in time and will be constant for the stationary state. In general, multiple constant values for $\omega_{s}$ result from the steady state condition. In each solution, the dynamical variable $\Phi(t)$ reads $\Phi(t)=(\omega_{s}-\omega_{0})t$. A simplification of the equations is possible setting the gain saturation coefficient equal to zero ($\epsilon=0$) as the general features of the system are still preserved. Equations (1)-(3), for the stationary solution, then became:

[TABLE]

[TABLE]

Equation (6) is transcendental and has no analytical solution. Therefore, numerical procedure is necessary to obtain $\omega_{s}$ for a given set of parameters. All values obtained for $\omega_{s}$ are then introduced in Eqs. (7) and (8) producing the steady state values of $N_{s}$ and $E_{s}$. The graphical inspection of the $\omega_{s}$ solutions can be done defining $F(\omega_{s})$ and $G(\omega_{s})$ as the left-hand and right-hand sides of Eq. (6), respectively.

For the single feedback case, i.e. $L_{1}$ blocked (numerically we set $\kappa_{1}=0$), the stationary solutions are shown in Figs. LABEL:sub@Fig4a and LABEL:sub@Fig4b for moderate feedback strength ($\kappa_{2}=1.0\times 10^{9}$ s*-1*) and for $R=1$ which means $L_{2}$ is equal to the length of the fixed arm $L_{1}$. For this case, the graphical solution of Eq. (6) is presented in Fig. LABEL:sub@Fig4a where $F(\omega_{s})$ (dashed black line) and $G(\omega_{s})$ (red line) are showed in terms of $\omega_{s}$ (multiplied by the roundtrip time of the fixed external cavity $\tau_{1}$, here blocked). The crossing points means that $F(\omega_{s})=G(\omega_{s})$. It gives the frequencies solutions $\omega_{s}$ of Eq. (6). The complete stationary solutions $(E_{s},N_{s},\omega_{s})$ are obtained using the Eqs. (7) and (8). These solutions can be represented in a projection of the phase space ($\Delta\omega\tau_{1}$,$\Delta N$) as it is shown in Fig. LABEL:sub@Fig4b, where $\Delta\omega=\omega_{s}-\omega_{0}$ and $\Delta N=N_{s}-N_{th,sol}$. They appear over an ellipse Sano (1994) which analytical equation can be obtained by manipulating Eqs. (6) and (7). Here, the steady solutions were obtained numerically. The term $N_{th,sol}=N_{0}-\gamma_{p}/G_{0}$ is the threshold value for $N$ when the laser is being operated in the solitary mode (with no feedback).

The appearance of the solutions, by changing parameters like the feedback strength or the time delay, happen through a saddle-node bifurcation of cycles when $F(\omega_{s})$ touch $G(\omega_{s})$. The result of this crossing process is the creation of two new values for $\omega_{s}$ in each bifurcation. These solutions are classified as limit cycles, one stable mode (lower points in the ellipse) and the equivalent unstable antimode (upper points in the ellipse). The comb of stable resonant frequencies $\omega_{s}$ provided by the compound cavity are called external cavities modes (ECMs). The center of the ellipse represents the solitary laser mode when the laser has no external feedback ($\kappa_{1}=\kappa_{2}=0$). The mode with the largest gain, is called Maximum Gain Mode (MGM). In Fig. LABEL:sub@Fig4b it is represented by a red square. It is the lower point of the ellipse in the ($\Delta\omega$, $\Delta N$) phase space projection representation. This mode should operate with the lower value of $N$ and the higher value of $E$ (and consequently higher intensity $|E|^{2}$). It is the mode with the lowest pump current threshold $J_{th,MGM}$. The stability of this mode in numerical calculations is highly sensitive with respect to the lasers parameters Zamora-Munt et al. (2010); Martinez Avila et al. (2008).

When no mode on the ellipse can stabilizes the laser, the LFF regime in which the emitted intensity of the laser $|E|^{2}$ presents irregular spike drops is observed Zamora-Munt et al. (2010); Martinez Avila et al. (2008); Sacher et al. (1989); Sano (1994); Hohl et al. (1995); Giudici et al. (1997); Sukow et al. (1997); Mulet and Mirasso (1999); Eguia and Mindlin (1999). A well accepted theoretical explanation of this chaotic behaviour was given by Sano Sano (1994). Accordingly, the ECMs are destabilized and became quasi-attractors causing mode hopping between them through wave mixing process. According to Sano, the trajectory of the system in the phase space projection surrounds the quasi-attractors (fixed points represented in the ellipses) in a local chaos transient stepping to the next ECM with lower frequency (to the left in the ellipse). The mode hopping repeat until the system reach the neighborhood of the MGM. The system remains surrounding the MGM until the trajectory crosses the stable manifold of the saddle anti-mode being repulsed from it and making the system experiment a frequency shift backing near to the solitary frequency operation. In this process the variable $N$ increases approaching to its threshold value $N_{th}$ and the intensity $|E|^{2}$ drops in a very short time scale. This entire process repeat in an irregular way characterizing the chaotic LFF regime.

For the double feedback configuration, the right-hand side $G(\omega_{s})$ of Eq. (6) has two contributions giving a beating effect due to the difference $\Delta\tau=\tau_{1}-\tau_{2}$. This is shown in Fig. LABEL:sub@Fig4c for $R=0.95$ and $\kappa=0.5\times 10^{9}$ s*-1*. In the calculations we used $\kappa_{1}=\kappa_{2}=\kappa$, which means that both arms of optical feedback are sending light back to the gain medium with the same strength. Experimentally, to provide that, equal laser threshold reduction due to each optical feedback alone was monitored and adjusted through mirrors alignment. It is important to notice that this pattern is not a result of optical frequency interferometry, once $\Delta\tau$ is very large if compared to the optical period as in our experimental and theoretical scan of $R$. Figure LABEL:sub@Fig4d shows that in this case, due to this beating effect between the external cavities lengths, the attractor skeleton represented by the points in the plane ($\Delta\omega\tau_{1}$,$\Delta N$) now settle over three ellipses. The scan of $R$ can also produce five, seven, nine and more ellipses. These multiple ellipses attractors result from the beating modulation, presented in Fig. LABEL:sub@Fig4c. This effect is directly related to the trajectory in phase space Sano (1994) changing the chaotic dynamics of the systems variables incorporating bursting in the LFF regime.

VI.2 Laser Threshold and Chaos Statistics

A clarifying view of the system physical behavior with double feedback can be reach examining the threshold $J_{th,MGM}$ variation of the MGM as lengths ratio $R$ is varied. The MGM is the most relevant fixed point of the system. Changes in it, through parameters variation, may induce bifurcations in the system. Even in chaotic regime of operation, the MGM remains the fixed point around which the trajectory spends most of the time. Actually, due to the infinity dimension of the system, the size of MGM basin of attraction must be small enough that usual trajectories never penetrate it. This is consistent with our deterministic numerical results and experiments. Futhermore, even if the basin of attraction is visited Zamora-Munt et al. (2010), the injection of external optical signal or noise can excite LFF drops. This justifies the proposition that this system is excitable Giudici et al. (1997). An alternative view point is to argue that the ultrafast deterministic oscillations in the dynamics excite the drops in an effective slow time scale attractor Martinez Avila et al. (2008).

Differently from the single feedback case, where the fixed points ($E_{s}$, $N_{s}$, $\omega_{s}$) when plotted in the phase space projection ($\omega$, $N$) remain over a single ellipse Sano (1994), the double feedback makes these points to be distributed in more than one ellipse changing the attractor in such a way that trajectories reach certain values producing non-gaussian - long tail - time distribution. The MGM threshold current $J_{th,MGM}$ dependence on $R$ reveals the effect of the parameters scan on the system dynamics changes. Figures LABEL:sub@Fig5a and LABEL:sub@Fig5b show how the system change its MGM threshold when it is subjected to a single and a double optical feedback, respectively. It shows analytically how $J_{th,MGM}$ (normalized to the solitary threshold current $J_{th,sol}$) depends on the scan of $R$, which for the single cavity means the variation of the corresponding feedback time.

For the laser to oscillate, with some non-null emission, the pump current $J$ has to be above $J_{th,MGM}$. It can be seen in Fig. LABEL:sub@Fig5a that the threshold current is reduced with respect to the solitary laser ($J_{th,MGM}<J_{th,sol}$ in about 5%). Notice that $J_{th,MGM}$ is constant for most values of $R$ except near $R=0$. Near this value, Figs. LABEL:sub@Fig5a and LABEL:sub@Fig5b are similar. To be reminded, for single feedback (Fig. LABEL:sub@Fig5a) the $R$ near zero means that the only active external cavity becomes comparable to the inner semiconductor laser cavity. It gives feedback time sufficient short to compete with the relaxations dynamics in the laser. In fact this range of values do not contemplate the occurrence of LFF Fischer et al. (1994) and it was not observed in our experiments.

Double feedback, as shown in Fig. LABEL:sub@Fig5b, gives a structure to the MGM threshold current as function of $R$. Its reduction has a base line equal to the single feedback case apart from the regions where $R$ is near a low order Farey number. The changes in $J_{th,MGM}$ is more significant with peaks according to the Farey order $p/q$ in $R$. Throughout the range $0<R<1$ equal Farey denominators $q$ give similar peaks for the threshold reduction. For the purpose of comparison with Fig. LABEL:sub@Fig5b, the sixth order Farey sequence (all the fractions p/q with $q\leq 6$) is shown in Fig. LABEL:sub@Fig5c. The vertical axis of Fig. LABEL:sub@Fig5c represent the Farey order as $1/q$. It is remarkable the comparison of these orders with peaks sizes of $J_{th,MGM}$ from Fig. LABEL:sub@Fig5b. The lower the Farey number order, the stronger the effect on the threshold reduction. Previous work related to our report remotes to the first studies of instabilities in light propagation through ring optical cavities. Ikeda and Mizuno have studied instabilities in double cavity nonlinear optical resonators Ikeda and Mizuno (1984, 1985); Mizuno and Ikeda (1989). They argue with ”frustrated instabilities” in a non-chaotic system with many oscillation modes. In their system, the frustration behavior always happens around the Farey fractions as in ours. However, as we treat an active laser system, our study is different, referring to chaotic dynamics and statistics.

Concerning the significant regions where the system shows spike-bursts we focus on $R$ values zooming the calculations close to Farey numbers. Fine details in $J_{th,MGM}$ indicate oscillations, as shown in Fig. LABEL:sub@Fig6a. Those behavior do not correspond to optical interferometric oscillations. They have to do with the difference between the two feedback time delays as discussed in the previous sections. Such oscillations have an impact in the chaotic statistics. Experimental and numerically calculated mean inter-spike times $\langle T\rangle$ through the same range of Fig. LABEL:sub@Fig6a are shown in Figs. LABEL:sub@Fig6b and LABEL:sub@Fig6c, respectively.

Approaching the $1/2$ quotient, the system maintains the mean $\langle T\rangle$ and its variance with an almost constant value until $R\approx 0.494$ while the laser operates in the LFF regime with Gaussian statistics of spikes. Proceeding towards $1/2$ makes $\langle T\rangle$ to grow and the variance diverges. The laser then exhibits bursting operation, in which we observe very long intervals without LFF as was shown in Figs. LABEL:sub@Fig3c and LABEL:sub@Fig3e. Varying $R$ further, $\langle T\rangle$ becomes larger as the quiescent time between bursts increases. It may even reach the condition where the intensity time series presents no LFF drops. Apart from LFF dynamics and statistics, the laser goes through many different forms of oscillation including local chaos around the MGM, periodic oscillation and continuous wave operation. None of these regimes is inspect in our study. To obtain $\langle T\rangle$ we have used a LFF pulse counting procedure and therefore, any information about the real dynamical state of the intervals without LFF drops could not be accessed through our procedure. When no LFF pulse was found in an entire time series, we placed $\langle T\rangle$ as the total duration of the series as shown in the top values of $\langle T\rangle$ in Figs. LABEL:sub@Fig6b and LABEL:sub@Fig6c. Continuous increasing of $R$ induces large variations with ups and downs of $\langle T\rangle$ similar to what happens with threshold current $J_{th,MGM}$ for the MGM. This indicates that the succession of transitions between chaotic LFF regime with Gaussian statistics to non-LFF regime are physically produced by the variations of the compound laser threshold in our model. Increasing $R$ beyond $1/2$ shows again the regions of different instabilities until chaotic LFF regime with Gaussian statistics is recovered. Experimentally and numerically, no exact symmetry is observed with respect to $R=1/2$. This is consistent with the analytical result for $J_{th,MGM}$ in Fig. LABEL:sub@Fig6a. It is also observed that on the exact value of the rational Farey fraction $R=1/2$, $J_{th,MGM}$ is equal to the single feedback threshold value (red dashed line). For this value, Figs. LABEL:sub@Fig6b and LABEL:sub@Fig6c show that chaotic LFF regime is recovered with Gaussian statistics as if the system had a single feedback with $\langle T\rangle$ down to the single feedback chaos case. In general similar chaos statistics occur for all Farey numbers $p/q$ where $p$ round trips in one cavity match $q$ roundtrips of the other one. The novelty to be emphasised here is that such quasi-periodicity in feedback times enters the laser dynamics to guarantee Gaussian chaos. When $q$ is big, the two cavities periodicity will only match in such long time, with respect to all relaxations times in the gain medium, so that all field coherence are washed out in between the compound cavities roundtrips. This restates the chaotic behavior again as in a single cavity.

VI.3 Control of the Chaotic Statistics

The long tail spike-bursting behavior is reached within these regions of transitions of ups and downs of $\langle T\rangle$. Considering that the change in $R$ within one of those $\langle T\rangle$ strong variations occurs for over thousands optical wavelengths, it suggests that without interferometric precision, one can adjust the system to have the desired statistics of spikes. Selecting the appropriate values for $R$ creates the possibility to have a system with probability distribution spikes with a tail index $\alpha$ chosen on demand. Figure LABEL:sub@Fig7a shows an even higher resolution numerical result for averaged output intensity of the laser and the number of LFF spikes of one of the chaos statistics transitions near the ratio $R=1/1$.

It is important to notice that the averaged output intensity varies almost linearly with the mean inter-spike time $\langle T\rangle$. So, in Fig. LABEL:sub@Fig7a, the growth of (black squares online) is matched with the decrease of the number of spikes drops (red balls online). The step used for $R$ scan in the calculations was $\delta R=0.00002$ which experimentally would correspond to a variation of 150 $\mu$ m in $L_{2}$. The optical wavelength in our lasers were 0.8 and 0.85 $\mu$m, so optical interferometry is secondary with respect to the observed dynamics as one tunes $R$. The laser suffers a transition from LFF chaos with Gaussian statistics to a regime with no intensity drops near $R=0.9925$. In this region the system presents the long-tail spike-bursting dynamics.

Other parameters are available to control the dynamics in addition to $R$ in the system. The easiest to manipulate is the pump current $J$. Envisaging control of the bursting statistics we calculated a fine tuning of the pump current maintaining $R$ fixed at the center of the transition presented in Fig. LABEL:sub@Fig7a. This is shown if Fig. LABEL:sub@Fig7b. The numerical range of $J$ variation was 5% of the solitary threshold current $J_{th,sol}$ with steps of 0.25$\%$ of this value. A corresponding experimental situation of a laser with solitary threshold current of 20 mA would fall in the fairly controllable 50 $\mu$A precision in current tuning. The long-tail index $\alpha$ drops from 1.7 (near Gaussian distribution) to 0.6, which signals very long tail statistics, and then recovers back to near the Gaussian statistics value. Graphics (i), (ii) and (iii) in Fig. LABEL:sub@Fig7c illustrate the nature of the described statistics with segments of intensity time series for different values of the pump current. This is a clear indication of the potential tuning of the current to control the chaotic bursting statistics.

VII Conclusion and Outlook

In this work, we have demonstrated experimentally, and confirmed with numerical solution of equation model, that a semiconductor laser with double feedback shows chaos statistics that can be controlled by realistic parameters variation. The chaotic behavior was manifested in the time intervals between power drops characteristic of the dynamics in such system. The feedback times ratio and the pump current were shown to play a fundamental role in the transitions from chaos with Gaussian statistics to long-tail distributions characterized by spiking bursts. Our chaotic device with this possibility of control may find application in optical encrypted communications and in simulations of neural network systems.

Acknowledgements.

Work partially supported by Brazilian Agencies Conselho Nacional de Pesquisa e Desenvolvimento (CNPq), Fundação de Ciência de Pernambuco (FACEPE) and CAPES Programa PROEX 534/2018 and by a Brazil-France, CAPES-COFECUB project 456/04.