Semiconductor lasers with optical feedback are well-known nonlinear dynamical systems. Under appropriate feedback conditions, these lasers emit optical pulses that resemble neural spikes. Influenced by feedback delay and various noise sources, including quantum spontaneous emission noise, the dynamics are highly stochastic. A good understanding of the spike timing statistics is needed to develop photonic systems capable of using the fast-spiking laser output for novel applications, such as information processing or random number generation. Here we analyze experimental sequences of inter-spike intervals (ISIs) recorded when a sinusoidal signal was applied to the laser current. Different combinations of the DC value and frequency of the signal applied to the laser lead to ISI sequences with distinct statistical properties. This variability prompts an investigation into the relationship between experimental parameters and ISI sequence statistics, aiming to uncover potential encoding methods for optical spikes, since this can open a new way of encoding and decoding information in sequences of optical spikes. By using ordinal analysis and machine learning, we show that the ISI sequences have statistical ordinal properties that are similar to Flicker noise signals, characterized by a parameter $\alpha $ that varies with the signal that was applied to the laser current when the ISIs were recorded. We also show that for this dataset, the ( $\alpha $, permutation entropy) plane is more informative than the (complexity, permutation entropy) plane because it allows better differentiation of ISI sequences recorded under different experimental conditions, as well as better differentiation of original and surrogate ISI sequences.

Semiconductor lasers, known for being fast, energy efficient, and low cost, exhibit, under the influence of optical feedback, a rich variety of complex dynamic behaviors. This diversity arises from the interplay of three factors: nonlinear light-matter interactions in the laser cavity, high dimensionality of the phase space, introduced by feedback delay time and stochasticity, due to the presence of quantum spontaneous emission noise and other noise sources. Exploiting these characteristics, semiconductor lasers with optical feedback serve as ideal systems for conducting controlled experiments aimed at unraveling nonlinear, stochastic, and time delay-induced phenomena. Moreover, they present an appealing prospect for innovative photonic applications, including random number generation and information processing. One extensively studied feedback-induced dynamical regime is the low-frequency fluctuations regime, characterized by spiking laser output. In this study, we investigate the statistical properties of the timing of optical spikes, analyzing experimental sequences of inter-spike intervals (ISIs). Employing a recently proposed technique that integrates nonlinear time series analysis and machine learning, we demonstrate that the ISI sequences exhibit statistical ordinal properties akin to Flicker noise.

## I. INTRODUCTION

Artificial intelligence (AI) systems are revolutionizing all fields of science and technology; however, their remarkable advancements come with notable drawbacks. The intensive computational requirements of AI, particularly during the training and operation of complex models, contribute to a significant energy footprint. The substantial energy consumption, often associated with very intensive computational implementation, prompts the exploration of alternative hardware architectures that allow more efficient ways to implement AI strategies. The need for efficient alternative solutions is motivating the development of novel photonic information processing systems such as specialized photonic neural networks.^{1–4}

With the aim of exploring fundamental mechanisms for implementing photonic neurons, here we analyze the output of a semiconductor laser. Semiconductor lasers are fast, energy efficient, low cost, and, most importantly, under appropriate conditions, they emit pulses of light (optical spikes) that resemble neural spikes.^{5–9} We analyze experimental data recorded from a semiconductor laser with optical feedback when a sinusoidal signal is applied to the laser current.^{10} We study feedback-induced spiking dynamics known as low-frequency fluctuations, LFFs. Although the LFFs have been extensively studied,^{11} the spike timing statistics remain poorly understood. Using the Fano factor technique, Tiana-Alsina and Masoller^{12} recently found that in the LFF regime, current modulation can generate long-range regularity in the timing of the spikes, which lock to the modulation and, despite the presence of noise, remain in phase over thousands of modulation cycles. Contrasting the behavior under pulsed and sinusoidal modulation,^{13} harmonic and subharmonic locking regions were found under pulsed modulation, but only subharmonic locking was found under sinusoidal modulation, which was interpreted as a signature of time-crystal behavior.

Here we analyze the ISI sequences focusing on the role of two physical parameters: the DC value of the laser current and the frequency of the modulation. We use a technique that combines ordinal analysis and machine learning.^{14,15} Ordinal analysis^{16} is a very popular^{17} symbolic time-series analysis method that captures order relations between values in a sequence of data points. The technique proposed in Refs. 14 and 15 consists of two main steps. First, an artificial neural network (ANN) is trained to infer the parameter $\alpha $ that characterizes the strength of correlations of Flicker noise (FN) time series, using as input features the probabilities of ordinal patterns computed from the FN time series. Then, the ordinal probabilities computed from an ISI sequence are used as input features to the trained neural network, which returns an $\alpha $ value that encapsulates information about the correlations present in the ISI sequence.

The results presented in Refs. 14 and 15 showed that the value of $\alpha $ returned by the ANN was useful to characterize a variety of signals. Deterministic and stochastic synthetic time series were analyzed (chaotic signals generated with the Bernoulli map, the Logistic map, the Schuster map, the Skew Tent map, and the Lorentz map, and also, fully stochastic signals generated by flicker noise, fractional Brownian motion, and fractional Gaussian noise). In addition, six empirical datasets were analyzed (fluctuations generated by a chaotic circuit, the output of a chaotic laser, the light intensity of a variable dwarf star, time series of sunspot numbers, time series of RR intervals of healthy subjects, and human gait data separated into three groups, according to the subjects’ ages). Figures 3 and 6 in Refs. 14 and 15 summarize the results and show that the $\alpha $ value returned by the machine learning algorithm carries information because the difference between the permutation entropy (PE) of a time series and the PE of Flicker noise with the same $\alpha $ is small when the time series is stochastic, but it is large when the time series is not fully stochastic. In a follow-up work,^{15} chaotic and periodic signals contaminated with noise were analyzed, and it was found that, as noise contamination increases, the difference in the two PE values (the PE of the signal and the PE of Flicker noise with the same $\alpha $) does not identify underlying determinism, but still uncover the presence of temporal correlations.

Motivated by these results show that the value of $\alpha $ returned by the ANN provides useful information to characterize different types of signals, here we analyze whether it can be used to characterize the experimental ISI sequences. In particular, we study if it is possible to recover information about physical parameters from the analysis of the ISI sequences, since this can open a new way of encoding and decode information in sequences of optical spikes. We show that the ISI sequences have $\alpha $ values that vary with the DC value and frequency of the sinusoidal signal applied to the laser current. We also show that the ( $\alpha $, permutation entropy) plane is more informative than the (complexity, permutation entropy)^{18–21} plane because it allows a better differentiation of ISI sequences recorded under different conditions, as well as a better differentiation of original and surrogate ISI sequences.

## II. DATA

We analyze experimental data consisting of ISI sequences recorded in Ref. 10 under sinusoidal modulation, which are freely available in Ref. 22. Sinusoidal modulation serves as a controlled and versatile method for systematically investigating the response of semiconductor lasers with optical feedback. With only two distinguishable features—amplitude and frequency—sinusoidal modulation offers a simple input signal to the lasers. This modulation provides a platform for studying the system’s dynamics, exploring nonlinear behaviors, and validating theoretical models. We considered the largest modulation amplitude and the range of physical parameters (DC values of the laser current, modulation frequency) for which the laser emits well-defined optical spikes. The DC value of the current, $I$, varies in the range of 25.5–27.5 mA, in a regular grid with 5 values, and the frequency of the applied signal, $\nu $, in the range of 1–70 MHz, in a regular grid with 70 values. In total, we analyzed 350 ISI sequences that contain about 5500–75 000 ISIs each (for the lowest and the highest pump current, respectively).

Figure 1 displays examples of ISI sequences recorded for $I=26.5$ mA and three modulation frequencies, the corresponding return maps, $\Delta T i$ vs $\Delta T i + 1$, and the corresponding distributions of ISI values. The spike times are evaluated using a threshold.^{10} In all panels, the ISIs are normalized to the modulation period. For $\nu =14$ MHz, the ISIs vary between one up to four modulation periods (a–c), for $\nu =36$ MHz the ISIs show small fluctuations around two modulation periods (i.e., the spikes are sub-harmonically locked to the current modulation and two modulation cycles occur in between consecutive spikes) (d–f), while for $\nu =64$ MHz the ISIs are close to three modulation periods but show some large fluctuations (g–i). As illustrated in Fig. 1, different combinations of the DC value of the laser current and modulation frequency lead to ISI sequences with distinct statistical properties. This variability motivates us to study the relation between these experimental parameters and the statistical characteristics of the ISI sequences. Specifically, we will study if it is possible to recover information on physical parameters from the analysis of the ISI sequence since this can open a new way to encode and decode information in sequences of optical spikes.

## III. METHODS

A time series characterized by a power spectrum $P(f)\u221d1/ f \alpha $ is referred to as Flicker noise (FN) or colored noise, where $\alpha $ serves as a measure of the correlations in the signal.^{23} For $\alpha =0$, the power spectrum is flat, and the noise is white noise. In the case of $\alpha >0$ ( $\alpha <0$), the larger the value of $\alpha $, the faster the power spectrum decreases (increases) with increasing frequency. In practical terms, large $\alpha $ reveals long-range correlations (or long memory).^{23}

We analyzed the ISI sequences with the method proposed in Refs. 14 and 15, which combines ordinal analysis and machine learning. While the method was originally proposed to distinguish noise and deterministic chaos, it can also be used to quantify the strength of the correlations in terms of $\alpha $. The code is freely available in Ref. 24. The method can be described in three main steps:

Calculate the probabilities of occurrence of ordinal patterns of length $D=4$, $p(i)$ with $i=1\u2026D!=24$, in $50000$ FN time series with $N= 2 20$ points generated with different values of $\alpha $, randomly chosen in the interval $[\u22124,4]$, and use them as features to train an ML algorithm to return the (known) value of $\alpha $. The FN series were generated with the open Python library

*colorednoise.py*.^{25}We connect the input layer to a single dense layer with 64 output units connected to a final layer, regressing all the information of the inputs ( $p(i)$) into the $\alpha $ number;Calculate the probabilities of occurrence of ordinal patterns of length $D$ in the ISI time series, $p(i)$ with $i=1\u2026D!$, and use them as features to the trained ML algorithm, which returns the (unknown) value of $\alpha $.

- Calculate the normalized permutation entropy,
^{16}^{,}$H$, and use $H$ and the value of $\alpha $ returned by the ML algorithm, as features to characterize the ISI time series. $H$ is calculated as^{16}where $S(D)$ is Shannon’s entropy, computed from the ordinal probabilities,$H(D)= S ( D ) ln \u2061 D !,$$S(D)=0$ if only one pattern occurs in the ISI sequence ( $p(i)=1$ and $p(j)=0$ $\u2200j\u2260i$) while $S(D)=ln\u2061 D !$ if all the patterns occur in the ISI sequence with the same probability ( $p(i)=1/D!$ $\u2200i$).$S(D)=\u2212 \u2211 i = 1 D !p(i)ln\u2061 p ( i ),$

^{15,26}$C$ considers the distance to both, the most regular and the most random states of the system. It is expressed as

^{26}

To test the ability of the ( $H$, $\alpha $) and ( $H$, $C$) planes to differentiate random from non-random ISI sequences, we also analyze surrogate ISI sequences, generated with the iterative amplitude-adjusted Fourier transform (IAAFT) algorithm that preserves both the distribution of ISI values and the power spectrum.^{27,28}

## IV. RESULTS

In the following, we present results obtained with ordinal patterns of length $D=4$, i.e., the patterns are defined by the relative values of four consecutive ISI intervals that carry information about the relative timing of five consecutive spikes. Similar results were obtained with $D=3$ and $D=5$. Figure 2 compares ISI sequences (black) and FN series (red) with very similar $H$ values. Panels (a) and (b) depict the ISI data and the $p(i)$ distribution for $I=26$ mA and $\nu =38$ GHz, which is associated with $\alpha \u2248\u22122$ (output of the ANN). Panels (e) and (f) depict the ISI data and the $p(i)$ distribution for $I=27$ mA and $\nu =18$ GHz, which produces $\alpha \u22482$. Panels (i) and (j) present the ISI data and the $p(i)$ distribution for $I=26.5$ mA and $\nu =17$ GHz, which produces $\alpha \u22480.3$. In the other panels [(c) and (d), (g) and (h), and (k) and (l)], we generate FN signals with $\alpha =\u22122,2$, and $0.3$, respectively. We observe that despite the absolute values being different the distributions of ordinal patterns are indeed very similar since the ordinal patterns only take into account relative differences of consecutive data values. As expected, FN and ISI sequences that have $\alpha \u22480$ have almost uniform distributions of ordinal patterns. This unveils the existence of ISI sequences that have statistical ordinal properties that are similar to Flicker noise with $\alpha $ that varies with the experimental parameters. It is noted that the way FN series are generated produces practically symmetric ordinal probability distributions [panels (d), (h), and (l)], a characteristic that is not fully reproduced by ISI laser data.

A more general framework is presented in Fig. 3, which displays the ( $H$, $\alpha $) and ( $H$, $C$) planes. The colored dots correspond to each frequency $\nu $ for a fixed DC value of the laser current, $I$, with the solid line representing the $H$ and $C$ values obtained from FN signals generated with $\alpha $ uniformly distributed in $[\u22124.4]$ (to obtain an almost-continuous line, $50000$ signals were analyzed). The gray square represents the FN with $\alpha =\u22124$. By comparing the first and the second row, we observe that the ( $H$, $\alpha $) plane can be used for characterizing the ISI sequences recorded under different conditions, while in the ( $H$, $C$) plane, the points collapse in an almost straight line, following the FN results. This can be expected as for highly stochastic, or high-dimensional deterministic time series, $H$ and $C$ are almost linearly related.^{18,21}

In Fig. 4, we analyze the results obtained from the analysis of the IAAFT surrogated ISI sequences, which preserves the distribution of ISI values and the power spectrum. We can see from the first row that $H$ and $\alpha $ plane preserves short temporal correlations (as $H$ is not always close to 1), but long correlations (captured by large $ |\alpha |$ values) are washed out. Hence, $H$ and $\alpha $ are more independent features in comparison with $H$ and $C$.

The values of $\alpha $ and $H$ are displayed in color code in Fig. 5, as a function of the modulation frequency and the DC value of the laser current. Here we see that for low-frequency modulation and high DC values of the pump current, the ISI sequences tend to have large positive values of $\alpha $, while negative $\alpha $ values tend to occur in narrow regions. The examples shown in Fig. 1 for $I=26.5$ mA and $\nu =14$, 36, and 64 MHz correspond to positive, nearly zero, and negative $\alpha $, respectively. Additionally, Fig. 6 displays the characterization of the ISIs in the three-dimensional space defined by $\alpha $, $H$, and the coefficient of variation, $CV$, which is the standard deviation of the ISI divided by the mean value. Here we see that ISI sequences recorded under low or high-frequency current modulation (blue or red symbols, respectively) tend to be located in different regions of this space. This suggests that by combining ordinal analysis with other time series analysis methods and by using machine learning, it may be possible to “decode” information from the ISI sequences. In particular, it may be possible to infer the frequency of the sinusoidal signal applied to the laser current.

## V. DISCUSSION

Instead of estimating $\alpha $ by fitting the power spectrum of the ISI time series to $ P ( f ) \u223c 1 / f \alpha $, we have used an “indirect” way: we estimated $\alpha $ using a machine learning algorithm, trained with features which were the ordinal probabilities of Flicker noise signals generated with different values of $\alpha $. We have used this “indirect” approach because previous studies^{14,15} have shown that, in this way, the machine learning algorithm returned a value of $\alpha $ that carried information useful to characterize different types of signals.

While for Flicker noise signals the value of $\alpha $ returned by the machine learning algorithm is easy to interpret, as it is equal to the value of $\alpha $ obtained by fitting the power spectrum, for the ISI sequences, it is not easy to interpret. We have done preliminary studies and found that the values of $\alpha $ obtained by using the machine learning algorithm and by fitting the power spectrum are in general different. Moreover, the two $\alpha $ values can vary, with the DC value and frequency of the signal applied to the laser current, in a partially correlated or completely uncorrelated way. This can be expected because the $\alpha $ value returned by the machine learning algorithm is estimated by using features such as the set of ordinal probabilities, and ordinal analysis is a nonlinear method that only considers the relative temporal ordering of the ISIs in the sequence—disregarding their actual values; in contrast, the value of $\alpha $ estimated by fitting the power spectrum is based on a linear method, the Fourier transform, which considers both the values of the ISIs and their temporal ordering.

Therefore, the value of $\alpha $ returned by the machine learning algorithm can be considered just a number that formally may serve as a feature. Indeed, we show here that this is a feature that can be used to characterize the ISI sequences. Our main conclusion is that, by combining ordinal analysis and machine learning, we can obtain a new feature, $\alpha $, which carries useful information, because, when combined with two other well-known features (the permutation entropy and the coefficient of variation) allows differentiating in a 3D feature space, the ISI sequences. Additional work is needed to quantify the performance of this feature. We have done preliminary studies that suggest that the value of $\alpha $ returned by the machine learning algorithm is less sensitive to the length of the ISI sequence than the value of $\alpha $ obtained by fitting the spectrum.

The value of $\alpha $ returned by the machine learning algorithm can locate, in Fig. 3, some ISI sequences close to the black line that represents Flicker noise, particularly for low frequencies; however, this does not necessarily mean that the ISIs behave as Flicker noise, as the $\alpha $ value obtained by fitting the spectrum is often different from the one returned by the machine learning algorithm. Therefore, we can only argue that, for low frequencies, the ISIs have statistical *ordinal* properties that are consistent with Flicker noise. Future work is planned to shed light on the relationship between the two $\alpha $ values—the one returned by the machine learning algorithm and the one obtained by fitting the power spectrum—for different types of time series and for different lengths of the time series.

## VI. CONCLUSIONS

We have characterized the statistical ordinal properties of the time intervals between consecutive spikes emitted by a semiconductor laser with optical feedback and current modulation, by analyzing experimental ISI sequences using a technique that combines ordinal analysis and machine learning. We have found that the ISI sequences have statistical ordinal properties that are similar to Flicker noise, as in the plane ( $\alpha $, permutation entropy), the ISIs can be located, depending on the experimental parameters, close to the line that represents Flicker noise. We have also found that the ( $\alpha $, permutation entropy) plane is more informative for characterizing the ISI sequences than the (complexity, permutation entropy) plane. In the three-dimensional feature space defined by $\alpha $, $H$, and $CV$ (the ISI’s coefficient of variation), we have found that ISI sequences recorded under low- or high-frequency modulation tend to be located in different regions. However, additional work is needed to “decode” the information encoded in the spike sequence and recover reliable information about the physical parameters, i.e., the frequency and the DC value of the signal applied to the laser current. Future work will be aimed at testing the performance of other ordinal-based quantifiers, such as the Fischer information, the Fischer–Shannon complexity plane, and quantifiers recently proposed to capture temporal symmetries.^{29,30}

For future work, it will also be interesting to analyze the relationship between the $\alpha $ values obtained by fitting the power spectrum and by using the machine learning algorithm. They can be expected to be different features (that may provide different information) because ordinal analysis only takes into account the relative ordering of the ISIs in the sequence—disregarding their actual values, while Fourier analysis considers both the ISI values and their temporal ordering.

## ACKNOWLEDGMENTS

This work was supported by Ministerio de Ciencia, Innovación y Universidades (No. PID2021-123994NB-C21), Institució Catalana de Recerca i Estudis Avançats (ICREA Academia), and de Gestió d’Ajuts Universitaris i de Recerca (AGAUR, No. 2021 SGR 00606). B.R.R.B. and E.E.N.M. acknowledge the support of the São Paulo Research Foundation (FAPESP), Brazil, Proc. 2018/03211-6, 2021/09839-0, Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), and the support of National Council for Scientific and Technological Development – (CNPq).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Bruno R. R. Boaretto:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). **Elbert E. N. Macau:** Conceptualization (equal); Investigation (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). **Cristina Masoller:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data used in this study are available in Ref. 22. The code that supports the results of this study is openly available in Github at https://github.com/brunorrboaretto/Laser_Spike_Timing_Analysis_Repository, Ref. 24.

## REFERENCES

*Semiconductor Lasers: Stability, Instability and Chaos*