The effects of periodic forcing on the temporally correlated spikes of a semiconductor laser with feedback

Optical excitable devices that mimic neuronal behavior can be building-blocks of novel, brain-inspired information processing systems. A relevant issue is to understand how such systems represent, via correlated spikes, the information of a weak external input. Semiconductor lasers with optical feedback operating in the low frequency fluctuations regime have been shown to display optical spikes with intrinsic temporal correlations similar to those of biological neurons. Here we investigate how the spiking laser output represents a weak periodic input that is implemented via direct modulation of the laser pump current. We focus on understanding the influence of the modulation frequency. Experimental sequences of inter-spike-intervals (ISIs) are recorded and analyzed by using the ordinal symbolic methodology that identifies and characterizes serial correlations in datasets. The change in the statistics of the various symbols with the modulation frequency is empirically shown to be related to specific changes in the ISI distribution, which arise due to different phase-locking regimes. A good qualitative agreement is also found between simulations of the Lang and Kobayashi model and observations. This methodology is an efficient way to detect subtle changes in noisy correlated ISI sequences and may be applied to investigate other optical excitable devices. © 2014 Optical Society of America OCIS codes: 190.3100 Instabilities and chaos, 140.1540 Chaos, 140.2020 Diode lasers, 140.5960 Semiconductor lasers, 250.5960 Semiconductor lasers. References and links 1. B. Lindner, J. Garcı́a-Ojalvo, A. Neiman and L. Schimansky-Geier, “Effects of noise in excitable systems,” Phys. Rep. 392, 321–424 (2004). 2. E. Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting (The MIT press, 2007). 3. W. Coomans, L. Gelens, L. Mashal, S. Beri, J. Danckaert, and G. Van der Sande, “Solitary and coupled semiconductor ring lasers as optical spiking neurons,” Phys. Rev. E, 84, 036209 (2011). 4. A. Hurtado, K. Schires, I. D. Henning, and M. J. Adams, “Investigation of vertical cavity surface emitting laser dynamics for neuromorphic photonic systems,” Appl. Phys. Lett. 100, 103703 (2012). 5. M. Turconi, B. Garbin, M. Feyereisen, M. Giudici, and S. Barland, “Control of excitable pulses in an injectionlocked semiconductor laser,” Phys. Rev. E 88, 022923 (2013). 6. A. Aragoneses, S. Perrone, T. Sorrentino, M. C. Torrent, and C. Masoller, “Unveiling the complex organization of recurrent patterns in spiking dynamical systems,” Sci. Rep. 4, 4696 (2014). 7. A. Aragoneses, T. Sorrentino, S. Perrone, D. J. Gauthier, M. C. Torrent, and C. Masoller, “Experimental and numerical study of the symbolic dynamics of a modulated external-cavity semiconductor laser,” Opt. Express 22, 4705–4713 (2014). 8. S. Barbay, R. Kuszelewicz, and A. M. Yacomotti, “Excitability in a semiconductor laser with saturable absorber,” Opt. Lett. 36, 4476–4478 (2011). 9. M. A. Nahmias, B. J. Shastri, A. N. Tait, and P. R. Prucnal, “A leaky integrate-and-fire laser neuron for ultrafast cognitive computing,” IEEE J. Sel. Top. Quantum Electron. 19, 1800212 (2013). 10. F. Selmi, R. Braive, G. Breaudoin, I. Sagnes, R. Kuszelewicz, and S. Barbay, “Relative refractory period in an excitable semiconductor laser,” Phys. Rev. Lett. 112, 183902 (2014). 11. B. J. Shastri, M. A. Nahmias, A. N. Tait, B. Wu, and P. R. Prucnal, “SIMPEL: Circuit model for photonic spike processing laser neurons,” arXiv:1409.7030v1 [physics.optics] (2014). 12. J. C. Alexander, E. J. Doedel, and H. G. Othmer, “On the resonance structure in a forced excitable system,” SIAM J. Appl. Math. 50, 1373–1418 (1990). 13. J. M. Buldú, J. Garcı́a-Ojalvo, C. R. Mirasso, and M. C. Torrent, “Stochastic entrainment of optical power dropouts,” Phys. Rev. E 66, 021106 (2002). 14. J. M. Buldú, D. R. Chialvo, C. R. Mirasso, M. C. Torrent, and J. Garcı́a-Ojalvo, “Ghost resonance in a semiconductor laser with optical feedback,” Europhys. Lett. 64, 178 (2003). 15. M. Giudici, C. Green, G. Giacomelli, U. Nespolo, and J. R. Tredicce, “Andronov bifurcation and excitability in semiconductor lasers with optical feedback,” Phys. Rev. E 55, 6414 (1997). 16. D. W. Sukow, J. R. Gardner, and D. J. Gauthier, “Statistics of power-dropout events in semiconductor lasers with time-delayed optical feedback,” Phys. Rev. A. 56, R3370 (1997). 17. J. Mulet and C. R. Mirasso, “Numerical statistics of power dropouts based on the LangKobayashi model,” Phys. Rev. E 59, 5400 (1999). 18. M. Sciamanna, C. Masoller, N. B. Abraham, F. Rogister, P. Mégret, and M. Blondel, “Different regimes of lowfrequency fluctuations in vertical-cavity surface-emitting lasers” J. Opt. Soc. Am. B 20, 37 (2003). 19. J. F. Martı́nez Avila, H. L. D. de S. Cavalcante, and J. R. Rios Leite, “Experimental deterministic coherent resonance,” Phys. Rev. Lett. 93, 144101 (2004). 20. Y. Hong and K. A. Shore, “Statistical measures of the power dropout ratio in semiconductor lasers subject to optical feedback,” Opt. Lett. 30, 3332 (2005). 21. A. Torcini, S. Barland, G. Giacomelli, and F. Marin, “Low-frequency fluctuations in vertical cavity lasers: experiments versus Lang–Kobayashi dynamics,” Phys. Rev. A 74, 063801 (2006). 22. J. Zamora-Munt, C. Masoller, and J. Garcı́a-Ojalvo, “Transient low-frequency fluctuations in semiconductor lasers with optical feedback,” Phys. Rev. A 81, 033820 (2010). 23. K. Hicke, X. Porte, and I. Fischer, “Characterizing the deterministic nature of individual power dropouts in semiconductor lasers subject to delayed feedback,” Phys. Rev. E 88, 052904 (2013). 24. Y. Liu, N. Kikuchi, and J. Ohtsubo, “Controlling dynamical behavior of a semiconductor laser with external optical feedback,” Phys. Rev. E 51, R2697-R2700 (1995). 25. W-S Lam, N. Parvez, and R. Roy, “Effect of spontaneous emission noise and modulation on semiconductor lasers near threshold with optical feedback,” Int. J. of Modern Phys. B 17, 4123-4138 (2003). 26. J. P. Toomey, D. M. Kane, M. W. Lee, and K. A. Shore, “Nonlinear dynamics of semiconductor lasers with feedback and modulation,” Opt. Express 18, 16955–16972 (2010). 27. Y. Takiguchi, Y. Liu, and J. Ohtsubo, “Low-frequency fluctuation induced by injection-current modulation in semiconductor lasers with optical feedback,” Opt. Lett. 23, 1369–1371 (1998). 28. D. W. Sukow and D. J. Gauthier, “Entraining power-dropout events in an external-cavity semiconductor laser using weak modulation of the injection current,” IEEE J. Quantum Electron. 36, 175 (2000). 29. J. M. Mendez, R. Laje, M. Giudici, J. Aliaga, and G. B. Mindlin, “Dynamics of periodically forced semiconductor laser with optical feedback,” Phys. Rev. E 63, 066218 (2001). 30. F. Marino, M. Giudici, S. Barland, and S. Balle, “Experimental evidence of stochastic resonance in an excitable optical system,” Phys. Rev. Lett. 88, 040601 (2002). 31. T. Schwalger, J. Tiana-Alsina, M. C. Torrent, J. Garcı́a-Ojalvo, and B. Lindner, “Interspike-interval correlations induced by two-state switching in an excitable system,” Europhys. Lett. 99, 10004 (2012). 32. S. Thorpe, A. Delorme, and R. Van Rullen, “Spike-based strategies for rapid processing,” Neural Netw. 14, 715–725 (2001). 33. D. Nikolić, P. Fries, and W. Singer, “Gamma oscillations: precise temporal coordination without a metronome,” Trends in Cognitive Sciences 17, 54 (2013). 34. C. Bandt and B. Pompe, “Permutation entropy: a natural complexity measure for time series,” Phys. Rev. Lett. 88, 174102 (2002). 35. J. Tiana-Alsina, M. C. Torrent, O. A. Rosso, C. Masoller, and J. Garcia-Ojalvo, “Quantifying the statistical complexity of low-frequency fluctuations in semiconductor lasers with optical feedback,” Phys. Rev. A 82, 013189 (2010). 36. M. C. Soriano, L. Zunino, O. A. Rosso, I. Fischer and C. R. Mirasso, “Time scales of a chaotic semiconductor laser with optical feedback under the lens of a permutation information analysis,” IEEE J. Quantum Electron. 47, 252-261 (2011). 37. J. P. Toomey and D. M. Kane, “Mapping the dynamic complexity of a semiconductor laser with optical feedback using permutation entropy,” Opt. Express 22, 1713-1725 (2014). 38. J. P. Toomey, D. M. Kane, and T. Ackemann, “Complexity in pulsed nonlinear laser systems interrogated by permutation entropy,” Opt. Express 22, 17840-17853 (2014). 39. A. B. Neiman and D. F. Russell, “Models of stochastic biperiodic oscillations and extended serial correlations in electroreceptors of paddlefish,” Phys. Rev. E 71, 061915 (2005). 40. R. Lang and K. Kobayashi, “External optical feedback effects on semiconductor injection laser properties,” IEEE J. Quantum Electron. 16, 347 (1980).

Abstract: Optical excitable devices that mimic neuronal behavior can be building-blocks of novel, brain-inspired information processing systems.A relevant issue is to understand how such systems represent, via correlated spikes, the information of a weak external input.Semiconductor lasers with optical feedback operating in the low frequency fluctuations regime have been shown to display optical spikes with intrinsic temporal correlations similar to those of biological neurons.Here we investigate how the spiking laser output represents a weak periodic input that is implemented via direct modulation of the laser pump current.We focus on understanding the influence of the modulation frequency.Experimental sequences of inter-spike-intervals (ISIs) are recorded and analyzed by using the ordinal symbolic methodology that identifies and characterizes serial correlations in datasets.The change in the statistics of the various symbols with the modulation frequency is empirically shown to be related to specific changes in the ISI distribution, which arise due to different phase-locking regimes.A good qualitative agreement is also found between simulations of the Lang and Kobayashi model and observations.This methodology is an efficient way to detect subtle changes in noisy correlated ISI sequences and may be applied to investigate other optical excitable devices.

Introduction
Excitable behavior is observed in many nonlinear, natural systems [1].Such excitable systems respond to external perturbations in an all-or-none manner.For perturbations below the excitable threshold the response is of small amplitude and linear, unable to drive the system far away from its rest state.For perturbations above the threshold the response is nonlinear and corresponds to a large excursion of the system's variables in phase space, manifested by a spike.Because information processing in the brain is associated with the excitable action potentials that propagate through neuronal axons [2], many efforts are being done to develop excitable devices that could mimic neuronal behavior in bio-inspired information processing networks.Particularly, semiconductor lasers with optical injection [3][4][5], optical feedback [6,7] or saturated absorber [8][9][10][11] have been investigated with this purpose.
A relevant issue is to understand how such neuron-like optical systems represent, in the sequence of spikes, an external input signal.When an excitable system is periodically forced with a stimulus above the threshold it can respond periodically, with a period that is related to the forcing period, a behavior that is called entrainment, or phase locking, in analogy with the behavior observed in periodically forced oscillators [12].Phase locking behavior is referred to as n : m, with n spikes for m forcing cycles.When the stimulus is below the threshold, intrinsic noise leads to stochastic entrainment [13,14] and complex spike patterns which contain information about the subthreshold input.
To study entrainment, external forcing can be implemented via a direct modulation of the laser pump current [24][25][26][27].
The inter-spike-intervals (ISIs) are strongly affected by the current modulation, as the spikes tend to occur in the same phase in the modulation cycle and the ISIs become integer multiples of the modulation period, T mod .For increasing modulation amplitude and fixed modulation frequency f mod , if f mod is high enough, the ISIs become progressively smaller multiples of T mod and if the modulation is strong enough, there is a spike for each drive cycle (1 : 1 locking) [28,29].While the role of the forcing amplitude in the statistical distribution of ISIs has been studied in detail by several groups [14,[27][28][29][30][31], less is known about the influence of the modulation frequency and the role of weak forcing in the ISI sequence (i.e., in the timing of the spikes).Taking into account that neuronal systems encode information in correlated spike sequences [32,33], it is crucial to understand how optical spike encoding works and how simple signals, such as a weak periodic modulation, can be encoded in a sequence of optical spikes.
Recently we studied experimental and simulated ISI sequences [6,7] employing a symbolic method of time-series analysis referred to as ordinal analysis [34].This methodology is capable of identifying the presence of serial correlations in datasets, and by analyzing the statistics of the different symbols (ordinal patterns, OPs), it can capture subtle variations in temporally correlated data [35][36][37][38].We unveiled serial correlations between consecutive spikes both, in the unforced laser (without modulation) and in the forced laser (with current modulation).We identified a minimal model that reproduces the correlations in both cases, and this model is known to describe temporal correlations in spike sequences of sensory neurons [39].Thus, our findings suggest that semiconductor lasers with optical feedback and current modulation, operating with suitable parameters, can be building-blocks to develop optical neurons that emulate the spike activity of biological neurons.
In [6,7] we analyzed the role of external forcing on the sequence of optical spikes by gradually increasing the modulation amplitude at a fixed modulation frequency.Here we address the role of the modulation frequency.We aim to understand how weak periodic signals with different frequencies are encoded in spike sequences.In spite of the fact that modulated LFFs have been extensively investigated, to the best of our knowledge, no study has addressed the influence of the weak modulation over a wide frequency range.Modulation frequencies close to the average LFF frequency were studied in [29,30], and encompassing the free spectral range of the external cavity in [27], in a regime of operation such that LFF spikes did not occur without modulation.
Here we study the laser spiking output when varying the modulation frequency and demonstrate significant changes in the statistics of the symbolic ordinal patterns.These changes are accompanied by clear modifications in the shape and peak of the ISI distribution, which reveals noisy phase-locking.Thus, we show that we can track subtle changes in the correlated spiking laser output, due to stochastic entrainment, by means of the OPs' statistics.This method for detecting noisy phase-locking can be useful to investigate temporal correlations and entrainment in other neuron-like optical systems.In addition, we show that for appropriate parameters simulations of the Lang and Kobayashi model [40] are in good qualitative agreement with the experimental results.

Experimental setup
The experimental setup used is as in Fig. 1a in [7].A semiconductor laser (Sony SLD1137VS), with a solitary threshold current I th =28.4 mA, temperature-and current-stabilized with an accuracy of 0.01 C and 0.01 mA, respectively, using a combi controller (Thorlabs ITC501), emitting at 650 nm, has part of its output power fed back to the laser cavity by a mirror 70 cm apart (the external cavity round-trip time is 4.7 ns).A 50/50 beamsplitter in the external cavity sends light to the detection branch consisting of a photo-detector (Thorlabs DET210) connected to an amplifier (FEMTO HSA-Y-2-40) and a 1 GHz digital storage oscilloscope (Agilent Technologies Infiniium DSO9104A).A neutral density filter in the external cavity allows to control the feedback power.The DC pump current is I DC =29.10 mA, the laser is operated at 19.00 C and a threshold reduction due to feedback of 7% is observed.
Through a bias-tee in the laser mount, the pump current is modulated with a sinusoidal signal provided by a waveform generator (Agilent 33250A), with frequency varying from 1 to 50 MHz in steps of 1 MHz, and a peak-to-peak amplitude varying from 0.8% to 2% of I DC , in steps of 0.4%.For these modulation amplitudes the laser current is always above the solitary threshold.The experiment is controlled by a LabVIEW program that acquires the time series, detects the spikes, and calculates the number of inter-spike-intervals (ISIs) until 40,000 ISIs are detected.Then, the program changes the modulation frequency and/or amplitude, waits a few seconds to let transients die away, and the process is repeated.

Lang and Kobayashi model
The Lang and Kobayashi rate equations [40], in adimensional form, for the slowly varying complex electric field E and the carrier density N are where α is the linewidth enhacement factor, τ p and τ N are the photon and carrier lifetimes respectively, G = N/(1 + ε|E| 2 ) is the optical gain (with ε a saturation coefficient), µ is the pump current parameter, η is the feedback coupling coefficient, τ is the feedback delay time, ω 0 τ is the feedback phase, β sp is the noise strength, representing spontaneous emission, and ξ is a Gaussian distribution with zero mean and unit variance.The current modulation is simulated as µ = µ 0 + asin(2π f mod t), where a is the modulation amplitude, f mod is the modulation frequency and µ 0 is the DC current.
Extensive simulations of the Lang and Kobayashi model were performed.While a detailed analysis is out of the scope of this work, we report in Section 6 a reasonable good qualitative agreement for the following model parameters: µ 0 = 1.01, a = 0.004 (corre- sponding to a peak-to-peak amplitude of 0.8% of µ 0 ), ε = 0.01, τ p = 0.00167 ns, τ N = 1 ns, β sp = 5 × 10 −5 ns −1 , η = 10 ns −1 , and α = 4.For each modulation frequency we simulated 2 ms and averaged the intensity time series over a sliding window of 1 ns to reproduce the bandwidth of the detection system.The averaged series contained between 22800 and 35000 dropouts for low and high frequency, respectively.

Ordinal symbolic analysis
We analyze the experimental and numerical ISI sequences by using ordinal analysis [34], as in [6,7].Each ISI sequence, {∆T i }, is transformed into a sequence of ordinal patterns (OPs), which are defined by considering the relative length of D consecutive ISIs and assigning them a symbol that indicates their relative length, in the same order as they appear in the sequence.The shortest interval is assigned 0 and the longest interval is assigned D − 1.For D = 2 the only two possibilities are: ∆T i > ∆T i+1 that gives the '10' OP, and ∆T i < ∆T i+1 that gives the '01' OP.For D = 3 there are six possibilities: ∆T i < ∆T i+1 < ∆T i+2 gives '012', ∆T i+1 < ∆T i < ∆T i+2 gives '102', and so on.The OP probabilities are then calculated by counting their frequency of occurrence in the sequence.
The symbolic transformation employed here has the drawback that it disregards the information about the precise duration of the ISIs, but it has the advantage that it keeps the information about temporal correlations among them, i.e., about correlations in the timing of the optical spikes.Specifically, in the next section we analyze correlations among 3 spikes (by using D = 2 OPs), 4 spikes (by using D = 3 OPs) and 5 spikes (by using D = 2 TPs).As the number of possible OPs increases as D!, the TP method has the advantage of allowing to infer the presence of correlations among 5 consecutive spikes by computing only 2 TPs (TP 01→01 and TP 10→10 ) instead of computing the probabilities of 4! = 24 OPs.

Experimental Results
Figure 1 displays the measured intensity time series and the ISIs distribution for six values of the modulation frequency, f mod , when the modulation amplitude is 1.2% of I DC .For each frequency 30 modulation cycles are shown.The ISI distribution is computed with bins centered at integer multiples of T mod (with the exception of the first bin, centered in 0).The modulation frequencies displayed are chosen to highlight different behaviors: either n : 1 locking predominates (revealed by a high peak in the ISI distribution at nT mod ), or there is a transition from n : 1 to n + 1 : 1 locking, revealed by the peaks at nT mod and (n + 1)T mod having nearly the same heights.
For f mod = 7 MHz (first row) the ISI distribution peaks at T mod .The time series reveals that the ISIs are in fact heterogeneous, as in this case the bin centered in T mod is about 143 ns wide.As the modulation frequency increases, the peak in the ISI distribution shifts to higher multiples of T mod and the ISIs become more homogeneous.At 26 MHz (third row) phase locking 2:1 occurs with 3:1 intermittency.In the time series one can notice that, after a dropout occurs in a modulation cycle, the next cycle takes place during the intensity recovery time, and the consecutive spike is separated in time by 2T mod .Similar observation holds for higher frequencies, now the modulation cycle being clearly visible in the intensity oscillations between consecutive spikes.For 39 MHz (fifth row) the 3:1 pattern is dominant and for 49 MHz (sixth row) we can see intermittent switching between 3:1 and 4:1.As the frequency increases and the modulation becomes faster, the ISIs become larger multiples of T mod as the dropouts are spaced by an increasing number of cycles.Figure 2 displays the variation of the mean ISI, < ∆T >, with the forcing frequency.In panel 2a we plot < ∆T > as a function of f mod for four modulation amplitudes (indicated as percentages of I DC ).A decreasing trend can be observed, interrupted by "plateaus" where < ∆T > oscillates (at about 30 MHz) and remains nearly constant (at about 45 MHz) for the two lower amplitudes, or continues to decrease for the higher amplitudes.
The origin of this behavior can be identified in Fig. 2(b), where the ratio < ∆T > /T mod is plotted vs. the modulation frequency.After an almost linear increase for low modulation frequencies (where < ∆T > varies very little compared to T mod ), two plateaus occur at < ∆T > /T mod ∼ 2 and < ∆T > /T mod ∼ 3.For strong modulation amplitude the shape of the plateaus become more clear, while for weak modulation, only signatures of the plateaus are evident.To investigate if the changes in the ISI distribution induced by the variation of the modulation frequency occur smoothly or rather abruptly, the probabilities of the first 5 bins [i.e., p n with n = 0 . . .4, p n being the probability of an ISI interval being in the bin (nT mod − T mod /2, nT mod + T mod /2)] are plotted vs. f mod .Figure 3 displays the results for the same modulation amplitude as in Fig. 1 (1.2%).A smooth variation of the probabilities is observed over the entire frequency range.For clarity we indicate with vertical lines the six frequencies corresponding to the panels in Fig. 1.It can be observed that p 1 displays a maximum (close to 1) at 7 MHz, p 2 , at 26 MHz and p 3 , at 39 MHz, while p 1 ∼ p 2 at 14 MHz, p 2 ∼ p 3 at 31 MHz and p 3 ∼ p 4 at 49 MHz. Figure 4 displays the results of ordinal analysis applied to the ISI sequences, for the same modulation amplitude.We present results for D = 2 and D = 3.For D = 2 [Fig.4(a)] we plot simultaneously three probabilities: the probability of one OPs ['01', as the probability of '10' is 1-P('01')] and the probabilities of two transitions (from one OP to the same OP, as the other two TPs can also be readily calculated from the normalization conditions: TP 01→01 +TP 01→10 = 1 and TP 10→01 +TP 10→10 = 1).For D = 3 [Fig.4(b)] we plot simultaneously the probabilities of the 6 OPs.
We can see smooth changes in these probabilities as the modulation frequency varies.The same type of smooth variation that was observed in the p n probabilities (Fig. 3), is seen here, more clear in the TP 01→01 and in the OP '210' probability (red curves).These probabilities are the ones that depart more from equiprobability and are anti-correlated.
To demonstrate that these changes are indeed significant, Figs.4(c,d) display the same probabilities, but calculated after the ISIs series have been shuffled (surrogate data).We can see that in these panels all the OPs and TPs are practically equiprobable, as expected, as no correlations exist in the surrogate data.
Comparing Fig. 3 with Figs.4(a),(b) we see that at 7 MHz (when p 1 is maximum) and at 14 MHz (when p 1 ∼ p 2 ) nothing remarkable occurs in the symbols' statistics; however, for higher modulation frequencies, changes in the ISI distributions manifest also in changes in the statistics of the OP probabilities and transitions between OPs: the maximum of p 2 and p 3 (occurring at 26 MHz and at 39 MHz respectively) are located just after the local maxima (minima) of the TP 01→01 (OP '210' probability), and the "equilibrium" situations (p 2 ∼ p 3 at 31 MHz and p 3 ∼ p 4 and 49 MHz) occur just after the local minima (maxima) of TP 01→01 (OP '210' probability).To demonstrate that the above observations are robust, in Fig. 4 we plot the probabilities p 0 . . .p 4 , as well as the OPs and TPs for weaker and for stronger modulation amplitudes.For the different amplitudes, above 20 MHz, the maxima and minima of the TP 01→01 curve precede the maxima and the "equilibria" of the p n probabilities.For strong amplitude we can see that, as p 4 (pink curve) vanishes (and therefore the "equilibrium point" between p 3 and p 4 also disappears), the local maximum and minimum in the transition probability curve also disappear.One remarkable feature is the structure that appears in the TPs in the low frequencies (< 5 MHz) with the increasing amplitude.This structure does not appear to be linked with changes in the ISI distribution, at least for the wide bins (centered at nT mod ) used in this analysis.To check this hypothesis we used narrower bins and plotted all probability density functions for lower frequencies for all amplitudes but up to now the origin of this structure is unclear.We speculate that it could be related to 1 : m phase-locking (m spikes per modulation cycle).A detailed investigation of the behavior is ongoing and the results will be reported elsewhere.
In the experiment we also used different external cavity lengths corresponding to feedback delays of 2.5, 7.5 and 10 ns.Qualitative similar behavior was observed, with similar structures also appearing in the lower frequencies for delays of 2.5 ns and 7.5 ns.1.6% (e,g); 2.0% (f,h).Legends as in Fig. 3.

Comparison with Lang and Kobayashi simulations
Figure 6 shows both experimental and numerical OP and transition probabilities for lower modulation amplitude.The modulation amplitude is 0.8% of the DC current for the experimental data, panels 6(a) and 6(c).In the simulations, panels 6(b) and 6(d), the modulation amplitude parameter is a = 0.004, corresponding to a peak-to-peak amplitude of 0.8% of µ 0 , the DC current parameter.For the numerical data the probabilities change slower with the varying frequency than in the experimental data, but present the same general trends.As the model is a quite simple one (it takes into account only one reflection in the external cavity and neglects multi-mode emission, spatial and thermal effects) only a qualitative agreement could be expected.We find it remarkable that no re-scaling or major changes in the parameters of the model are needed to reproduce the general behavior of the probabilities.The agreement is not good for stronger modulation amplitude.More detailed numerical investigations are ongoing and should be reported elsewhere.
It is important to note that for D = 3 OPs, two clusters of OPs, formed by '021-102' and '120-201', where the OPs occur with the same probability, appear in the entire frequency range both for experimental and numerical data, being more clearly visible for higher frequencies.These clusters were reported in [6,7] and were explained in terms of a simple model that was previously used to model temporal correlations in sensory neurons.Thus, we demonstrate here that the two clusters are robust and persist over a wide range of modulation frequencies.

Conclusions
We have experimentally investigated the spiking output of a semiconductor laser with optical feedback in the LFF regime, under weak current modulation.In this regime the laser behaves as a weakly forced excitable system.We focused on the effects of varying the modulation frequency in the sequence of inter-spike-intervals (ISIs).With increasing modulation frequency the ISIs become larger multiples of the modulation period.We found that the mean ISI does not decrease monotonically as the modulation frequency increases, but displays smooth oscillations and plateau-like behavior due to noisy phase-locking.By using a symbolic method of analysis capable of detecting temporal correlations in datasets, we identified subtle changes in the correlations present in the ISI sequence (revealed by variations in the probabilities of the ordinal patterns and transitions), that complement the information extracted from the ISI distribution.The smooth variations in the symbolic probabilities were shown to be related to changes seen in the ISI distribution.For increasing modulation amplitude we observed that the phase-locking regions migrate to higher frequencies, became wider and the locking became more clear.We have also shown that, for appropriated parameters, simulations of the Lang and Kobayashi model are in good qualitative agreement with the experimental results.The symbolic methodology used is an efficient way to detect subtle changes in noisy correlated datasets and may be applied to investigate other optical excitable devices.A detailed comparison with model predictions, as well as an investigation of 1 : m phase locking behavior at very low frequencies is in progress and will be reported elsewhere.

Fig. 1 .
Fig. 1.Experimental time series of LFF intensity spikes (left column) and the corresponding inter-spike-interval (ISI) distribution (right column).The modulation amplitude is 1.2% of I DC and the modulation frequency is 7 MHz (a-b); 14 MHz (c-d); 26 MHz (e-f); 31 MHz (g-h); 39 MHz (i-j); 49 MHz (k-l).In the left panels only 30 modulation cycles are shown, but the ISI distributions in the right panels are computed from 40000 ISIs.

Fig. 2 .
Fig. 2. (a) Mean ISI as a function of the modulation frequency when the modulation amplitude is 1.2% of I DC .(b) Ratio between the mean ISI and the modulation period, < ∆T > /T mod , versus modulation frequency for several modulation amplitudes.

Fig. 3 .
Fig. 3. Probability, p n , that an ISI value is with the interval nT mod − T mod /2, nT mod + T mod /2.n values in the legend.The first five probabilities are shown, for the same modulation amplitude as Fig. 1 (1.2%).