Canard resonance: on noise-induced ordering of trajectories in heterogeneous networks of slow-fast systems

The experiment, inspired by [34, 35], is based on a semiconductor laser array subject to optoelectronic feedback, as shown on figure 1. The laser array consists of 451 Vertical Cavity Surface Emitting Lasers (thus single longitudinal mode), all biased by a single power supply. The average laser threshold of the population is 183.3 mA with a standard deviation of 5.8 mA (more details are available in [63, 64]). Due to the spectral distribution of the lasers (the array's spectrum is 2 nm broad), coherent coupling via field interaction can be neglected. In fact, far field measurements show no sign of interference between the individual lasers. The beam emitted by the laser array is collected by a short focal length lens and split in a detection and a feedback branch. In an image plane along the feedback branch (marked as 'NF' in figure 1), an iris whose aperture and position can be controlled allows to select the part of the array which will be guided further towards the feedback detector. On this path, a clippable mirror can send the beam towards a rotating screen for image analysis in absence of feedback. The selected part of the beam is focussed onto a photodetector (photodetector 1 in figure 1). The voltage produced by this detector is proportional to the total intensity emitted by the region of the array selected by the iris. This voltage is logarithmically amplified and its continuous component is suppressed. After a final linear amplification stage it is injected as a controller for the power supply driving the whole array. This leads to a global coupling through the optoelectronic reinjection. The detection branch is split in two by a prism situated in an image plane (labelled 'NF') on figure 1 and parts of each of these beams can be selected by two additional iris. Each of the selected beams are then focussed on two independent detectors which enable monitoring of parts of the whole network down to the single element. Additional details can be found in [64].

Zoom In Zoom Out Reset image size

We consider the following adimensional model for a population of N coupled lasers, subject to optoelectronic feedback. The model in physical variables together with the scaling is fully described in [63, 64]

$$\begin{align} \dot{x}_i(t) & = x_i(t)(y_i(t)-1)\nonumber\\ \dot{y}_i(t) & = \gamma\left[\delta_i-y_i(t)+k\left(w(t)+A\ln(1+\alpha X(t))\right)-x_i(t)y_i(t)\right]\nonumber\\ \dot{w}(t) & = -\epsilon\left[w(t)+A\ln(1+\alpha X(t))\right]\,, \end{align} \tag{ 1 }$$

where x_i (t) is the normalised intensity and y_i (t) is the normalised carrier density of each laser, and time is normalised to the photon lifetime. The variable w(t) represents the scaled high-pass filtered feedback current, and is common for the whole population. This feedback depends nonlinearly on the mean intensity $X(t) = \frac{1}{N} \sum_i x_i (t)$: the lasers are coupled in an all to all configuration through the feedback. The lasers differ by their lasing thresholds, and thus by their relative pump currents δ_i . We order the lasers in descending order of their pump currents; i.e. ${\delta _1} \ge {\delta _2} \ge {\delta _3} \ge \ldots$ The parameter k represents the coupling strength and A and α are the parameters of the non-linear saturable absorber in the high-pass filtered feedback current. The timescales of the system depend on the ratio of photon and carrier lifetimes $\gamma = \mathcal{O}(10^{-3})$, and $\epsilon = \mathcal{O} (10^{-4}-10^{-5})$, which relates to the feedback filter.

We integrate noise in equations (1) in a way similar to spontaneous emission noise, as an additive noise term added to the evolution of the complex laser electric field E_i (t), $|E_i(t)|^2 = x_i(t)$,

$$\begin{align} \dot{E}_i(t) = \frac{1}{2}(1+i\alpha_H) E_i(t) (y_i(t)-1) + \sigma \xi_i(t)\,, \end{align} \tag{ 2 }$$

where ξ_i (t) are independent white gaussian complex noise terms with zero mean and unit variance, and σ is the noise strength. We use 2 whose deterministic part is equivalent to that of the equation for x in 1 because in 2 the simple zero-mean complex noise term adequately captures the spontaneous emission process, which is less convenient to express in terms of the intensity x. The linewidth enhancement factor α_H does not play a role in the dynamics of this system, but is added for completeness. We chose α_H  = 4 in all simulations with noise. Of course, a proper description of all noise sources in semiconductor lasers leading to accurate modeling of spectral lineshapes or relative intensity noise spectra would require additional stochastic terms (see e.g. [65–67]). However, here we only aim at introducing a minimal yet physically meaningful stochastic term and at observing its impact on the canard trajectories.

We use an Improved Euler algorithm incorporating additive noise, with a time step of dt = 0.1. In our simulations, we have always used the timescale parameters γ = 4 · 10⁻³ and  = 10⁻⁴, and feedback parameters $k = 0.7, A = 1/k$ and α = 2, unless otherwise mentioned, as these parameters reproduce the experimental results qualitatively. Sample checks performed using the stochastic Runge–Kutta algorithm described in [68] implemented by [69] did not show any inconsistencies.

The experimental setup is built to mimic neural dynamics; the neuromorphic laser can show excitable, oscillating and chaotic dynamics, depending on the parameters. We briefly review [63, 64] the dynamics of a single node (N = 1, X(t) = x₁(t)≡x(t)) in order to facilitate comparison with the dynamics of the population. The system is a slow-fast system with being the slowest timescale in the system. The critical manifold, defined by  = 0, has two branches. The OFF-branch is defined by

$$\begin{equation} x = 0 \,, y = \delta + kw \,, \end{equation} \tag{ 3 }$$

the motion on the OFF-branch is towards the non-lasing fixed point at $w = 0, \, x = 0, \, y = \delta$. The OFF-branch is attracting for $w\lt(1-\delta)/k$, the non-lasing fixed point is stable if $\delta\lt1$. The ON-branch is implicitly defined by

$$\begin{equation} x(w) = \delta - 1 + k\left[w+ A\ln\left(1+\alpha x(w)\right)\right] \,, y = 1 \,, \end{equation} \tag{ 4 }$$

the motion is towards the lasing fixed point at $x = \delta-1, y = 1, w = -A\ln\left(1+\alpha(\delta-1)\right)$. The ON-branch emerges at $x = 0, \, w = (1-\delta)/k$ and is attracting for $x\gt kA-\alpha^{-1}$ or, equivalently, $\partial_x w(x)\gt0$. The critical manifold is shown in figure 2(c).

Zoom In Zoom Out Reset image size

Thus, the dynamics is excitable if $\delta\lt1$: The non-lasing fixed point is stable but the system responds to perturbations either by relaxing linearly back to the original state or, if the perturbation overcomes a certain threshold, by undergoing a large excursion (mostly independent of the details of the perturbation) in phase space due to the very different time scales. In particular, noise may drive the laser away from the OFF-branch; if the system reaches the ON-branch, a spike is emitted. The excitable dynamics is shown in figure 2(a), where noise triggers the random emission of excitable pulses.

The lasing fixed point is stable for $\delta\gt1.5136...$, and the dynamics are excitable as well. The fixed point undergoes a Hopf bifurcation as δ decreases, followed by a period doubling cascade and the emergence of a strange attractor with small amplitude. As this attractor grows in size, the critical manifold becomes part of it, and the system emits a large spike on top of the small amplitude chaotic dynamics. This chaotic region only exists for a small parameter range around $\delta \approx 1+kA-\alpha^{-1}$.

For $1\lt\delta\lt1.5136...$ none of the fixed points is stable. The neuromorphic laser exhibits, depending on the pump current δ, chaotic or regular spiking oscillations, as the system alternates between the two branches. Importantly, while the periodic alternation between the two branches is common to many planar slow-fast system, the chaotic dynamics and oscillatory jumps between these branches are enabled by the presence of a third intermediate time scale. The regular spiking regime, which is the focus of our numerical study, is illustrated in figures 2(b) and (c).

From the analysis above, one expects two parameter regions leading to excitable behaviour. We concentrate experimentally on the first one, where the laser is in a stable non-lasing state.

To check this behaviour, we close the population control iris such that only the part of the beam emitted by a single laser is actually transmitted towards the feedback detector. In this configuration, the other lasers are slaved to the dynamics of the selected one and the system can be analysed as a single node. While additional coupling between lasers (mostly through gain medium or thermal dynamics) could exist, we have no evidence of these effects taking place here and therefore analyse our results in the framework set by the model where the only coupling is through the optoelectronic reinjection. We focus on the response of the system to external perturbations when it is biased below its coherent emission threshold and close to the first instability. We direct towards the feedback detector a beam coming from another semiconductor laser emitting at 970 nm. The current of that laser is pulsed so that optical pulses of 40 µs duration are emitted at a repetition rate of 10 Hz. We apply stronger and stronger light pulses from that laser onto the feedback detector and monitor the response of the system. The results are shown in figure 3.

Zoom In Zoom Out Reset image size

For small perturbations, below 0.2 mW approximately, the system's response is almost always linear with only very rarely a large spike response. On panel (a) of figure 3, we show a superposition of 10³ responses to perturbations of 0.22 mW. The perturbations are applied at a repetition rate of 10 Hz, low enough that all realizations can be considered as independent. At this power value, about half of the perturbations elicit a large spike response. The latter arise after a small and variable delay as show by the dispersion on the rising edge. To confirm the only two possible responses, we plot on panel (b) of figure 3 the mean intensity emitted in a short time window of duration 0.1 ms at a time Δ_t  = 0.1 ms after the perturbation signal has been emitted, for increasing perturbation power. Only two regions are often visited by the system. The bottom one, close to 0.02 arb.units corresponds to the linear response and the top one, close to 0.25 arb.units corresponds to the excitable orbit. These two regions are clearly separated by the large empty area in the diagram, which shows that either a linear response or a large response can be emitted. The former is most probable for low perturbation power and the latter is strongly dominant for large perturbation power. We quantify the efficiency in the generation of these excitable responses on panel (c) of figure 3. The efficiency is the ratio between the number of large responses over the number of perturbations. The excitability threshold is easy to identify and it coincides with a strong decrease in the dispersion of the responses arrival time.

While for a single laser the critical manifold has an OFF and an ON branch, for a population of N lasers the critical manifold has 2^N branches: each laser can be 'on', or 'off':

$$\begin{eqnarray} x_i = 0\,, & y_i = \delta_i + k\left[w+A\ln(1+\alpha X(w))\right] & \mbox{(OFF)} \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} y_i = 1\,, & x_i = \delta_i-1 + k\left[w+A\ln(1+\alpha X(w))\right] & \mbox{(ON)}\,. \end{eqnarray} \tag{ 6 }$$

From these 2^N branches, two branches are symmetrical: the OFF-branch, with all lasers 'off' (x_i  = 0 , ∀i), and the ON-branch, with all lasers having a nonzero intensity. The other 2^N −2 branches are asymmetrical, with $0\lt N^+ \lt N$ with nonzero intensity and N − N⁺ lasers having zero intensity.

The OFF-branch is simply characterized as $x_i(w) = X(w) = 0, y_i(w) = \delta_i+kw$. To characterize the nonzero branches with $N^+\gt0$ lasers with nonzero intensity, we average over the intensities x_i of the lasing elements (6), and find an implicit equation for the average intensity X(w),

$$\begin{equation} X(w) = \frac{N^+}{N}\left[\delta_{\textrm{ON}}-1+k\left(w+A\ln\left(1+\alpha X(w)\right)\right)\right]\,, \end{equation} \tag{ 7 }$$

where the suffix 'ON' denotes an average over the lasing elements, $X_{\textrm{ON}} = \frac{1}{N^+}\sum_{i \in \textrm{ON}}x_i = \frac{N}{N^+} X$ and $\delta_{\textrm{ON}} = \frac{1}{N^+}\sum_{i \in \textrm{ON}}\delta_i.$ The individual intensities of the lasing elements are given by $x_i(w) = X_{\textrm{ON}}(w) +\delta_i - \delta_{\textrm{ON}}$.

Each of these branches emerges as $x_i\gt0 \,\forall i\in \textrm{ON}$, or equivalently, $X_{\textrm{ON}}\gt\delta_{\textrm{ON}}-\min(\delta_i)_{i\in \textrm{ON}}$. This means that the branches with a single laser i 'on' (and thus $\delta_{\textrm{ON}} = \delta_i$ ) emerge from the OFF-branch at X = 0 and kw = 1 − δ_i . The first of these branches to emerge (for minimal w) is the branch with the laser 1 'on', at $w = (1-\delta_1)/k$, as this laser has the highest pump current δ₁.

A branch with two lasers i and j 'on', $\delta_i\gt\delta_j$, emerges at $x_j = 0, x_i = \delta_i-\delta_j, X_{\textrm{ON}} = (\delta_i-\delta_j)/2$ and $w = \left[1-\delta_j-kA\ln\left(1+\alpha/N(\delta_i-\delta_j)\right)\right]/k$ (equation (7)). It is easy to show that this point lies on the branch with only laser i 'on': the branch with two lasers i and j 'on' emerges from the branch with only laser i 'on'. In an analogous way, we find that each branch emerges from a 'parent branch' with the same lasers being 'on', except for the laser with the lowest intensity: A new branch emerges at $x_k = 0, X_{\textrm{ON}} = \delta_{\textrm{ON}}-\delta_k, \delta_k = \min_{\textrm{ON}}(\delta_i)$, from the branch with N⁺−1 lasers 'on' and laser k 'off'.

The critical manifold for three lasers is shown in figure 4(a). From left to right, as the mean intensity X increases, we have the OFF-branch (X = 0, coinciding with the ordinate), three branches with one single laser 'on' and the two other lasers 'off', three branches with two lasers 'on' and the third laser 'off', and, as the rightmost branch, the ON-branch with all three lasers lasing. In figure 4(b) we zoom in on the low intensity region where the branches emerge from each other. We observe the phenomenon described above: From branch '1', at $X = (\delta_1-\delta_2)/3$, branch '1+2' emerges, and, at $X = (\delta_1-\delta_3)/3$, branch '1+3' emerges. From branch '1+2', at $X = (\delta_1+\delta_2 -2\delta_3)/3$, branch '1+2+3' emerges, etc.

Zoom In Zoom Out Reset image size

We calculate the linear stability of these branches in appendix A. The OFF branch (X = 0) has a large attracting part $w\lt(1-\delta_1)/k$; at this point branch '1' emerges and the OFF-branch becomes repelling in the direction of the first laser. The other symmetrical branch, for which all lasers are 'on' (the ON-branch), is attracting for high enough intensities $X\gt kA-\alpha^{-1}$, or equivalently, $\partial_X w(X)\gt0$.

The asymmetrical branches, with $N^+\lt N$ lasers 'on', are mostly repelling: along the directions of the non-lasing nodes they are attracting only for a very small range of intensities $\delta_{\textrm{ON}}-\min\{\delta_i\}_{i\in \textrm{ON}}\lt X_{\textrm{ON}}\lt\delta_{\textrm{ON}}-\max\{\delta_i\}_{i\in \textrm{OFF}}$. Consequently, only the branches for which $\max\{\delta_i\}_{i\in \textrm{OFF}}\lt\min\{\delta_i\}_{i\in \textrm{ON}}$ do have a small attracting part. The point where such a branch looses stability, coincides with the emergence of a new branch, with one extra laser 'on'. Along the lasing directions, these branches are stable if $\partial_X w(X)\gt0$. Combining these conditions, depending on parameters, typically only a few branches with less than half the lasers 'on', and only those with the highest thresholds δ_i , have a small attracting part for low intensities.

The stability of the critical manifold for three coupled lasers is indicated in figure 4, with attracting parts in full blue lines and repelling parts in dashdotted red lines. Clearly, the two symmetrical branches, where all lasers are 'on' (ON-branch) and where all lasers are 'off' (OFF-branch), have a large attracting part. The six asymmetrical branches are mostly repelling, with only one of them (the branch with laser 1 'on') having a small attracting part for very low intensities.

To uncover experimentally indications of the phase space structure described in the previous subsection, we run the experiment in the relaxation oscillations regime of the fully coupled network, where the system periodically switches between two main branches of the slow manifold, the ON branch and the OFF branch. From the theoretical analysis, one expects that this switching actually occurs through a series of steps where each laser leaves its 'off' state, the whole network going through a complex set of branches of the slow manifold. We show an example of this on figure 5.

Zoom In Zoom Out Reset image size

The complex structure of the critical manifold of the fully connected network cannot be observed easily at the network level especially since, as observed in section 4.1, only a small number of the 2^N branches (besides the ON and OFF ones) are stable and only in a very small parameter region. However, one can still expect to observe that the network cycles among many distinct states, each one corresponding to a different number of lasers in their 'off' or 'on' states. Thus, a signature of the phase space structure can be identified by monitoring two very different nodes, as shown on figure 5(a) and (b). In this case, lasers n and q, selected for their different threshold current, leave their unstable 'off' state at different times during the periodic oscillations. However, most of the dynamics clearly takes place between the two main OFF and ON branches, where all lasers are either 'on' or 'off'.

In this oscillatory regime, the logarithmic and linear amplifiers of the feedback loop are set such that the current driving the whole network oscillates between 158.1 and 219.2 mA during the oscillation cycle, so that nearly all lasers are periodically switching 'on' and 'off'. This of course leads to exciting laser relaxation oscillations in each individual node. These correspond to the fast oscillations at the onset of the 'on' phase in figure 2 panel (b). Each laser shows a different radiofrequency power spectrum (see figure 5(c)), since all of them orbit around a different region with respect to their coherent emission threshold. The different local maxima in some peaks also exist even in the absence of optoelectronic reinjection. Their position does not appear to change under modification of any externally accessible parameter. They are therefore an intrinsic feature of the laser array, which does not seem to influence much our experimental observations beyond that. These laser relaxation oscillations (detected here by a Thorlabs PDA8GS detector with 9 GHz bandwidth) can only be observed in the power spectrum due to detection and physical noise. This observation is not unexpected when the slowest variable is orders of magnitude slower than the other two. Moreover, it clearly hints at the relevance of physical noise for a proper description of the system when each node leaves its 'off' state [70, 71].

As underlined in section 4.2, laser relaxation oscillations can be experimentally observed only in the Fourier domain, which hints at noise causing much smaller overshoots than those presented on the numerical traces in figure 7(a) or even figure 2(b). We verify numerically at a basic level that the addition of noise in the simulations leads to a better correspondence between theoretical and experimental data: in figure 6, we show the normalized intensity of two different lasers together with the total intensity in logarithmic scale, in presence and absence of noise, mimicking the experimental figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the numerical simulations without noise (figure 6(a)) the lasers switch on one by one, in order of their pump currents δ_i . The transition is abrupt for all lasers (but smooths out for the mean intensity), and there is a considerable time difference between the switching on of the first and the last laser. When including noise in the simulations (figure 6(b)) we see some experimental features reproduced: the first laser intensity (green) increases first in a gradual way, but the spikes are fairly simultaneous.

The above effect is to some extent not unexpected: from the phase space structure described in section 4.1 the network should leave the OFF branch via a series of asymmetrical canard branches, each laser leaving its 'off' state one after the other. Since canard orbits are known to be extremely sensitive to noise (see e.g. [72, 73]), one can expect noise to strongly influence that part of the dynamics. Beyond that however, as we shall see, noise can also have an unexpected ordering effect across this heterogeneous network. We analyse this first for a small network of three lasers by studying the departure from the OFF manifold for increasing amount of noise.

A typical time trace and phase plot in the deterministic regime is shown in figure 7(a) and (b): Starting on the OFF branch, the system follows this branch as w(t) increases. Although this branch becomes unstable as $w\gt(1-\delta_1)/k$, we find that the system stays close to it for a while longer, until, at some point, the first laser (blue curve) switches on. This constitutes a first canard. After a few laser relaxation oscillations, the system settles on the repelling branch with laser 1 'on', and the two other lasers switched off. The next laser to switch on is laser 2 (red curve), which again follows a piece of canard orbit, the system remaining close to this repelling branch corresponding to lasers 1 and 2 switched on. Then follows laser 3, and the system follows the attracting ON-branch, until it reaches $w_{\min}$ and all lasers switch off.

This pattern, of the lasers switching on one by one, in order of the pump currents, is found for any population size and parameter set, and becomes more pronounced as the difference between the lasing thresholds δ_i increases. As is typical for canard dynamics, the system remains close to the branches of the slow manifold, after they have become repelling. In particular, for our choice of parameters, the system follows a repelling part of the OFF-branch, and repelling parts of the asymmetric branches. Thus, in the deterministic case, the whole process of switching from the main OFF to the main ON branches occurs via a series of canard orbits.

In the stochastic case instead, noise can drive the system away quicker from the repelling branches, towards attracting branches. Because the attracting part of the critical manifold mostly consists of the symmetric branches (ON-branch with all lasers 'on', and the OFF branch), and the asymmetric branches are mostly repelling, the coherence of the population increases with low noise. This is easily observed on figure 7: noise reduces the time the system spends on the repelling asymmetric branches, as can be seen by comparing the waveforms and phase plots without noise (panels a, b) and with low noise (panels c, d); if noise is too strong the oscillations become less coherent (panels g, h). Between the extreme values, there is an optimal noise level (panels e, f), that balances these two effects, and for which the dynamics of the population is maximally ordered.

This observation indicates that the population across the network is maximally coherent for an optimal noise level, beyond which point noise starts to dominate the dynamics.

We quantify this canard resonance effect in the following two ways: First, we average the variance of the intensities of the population over time. Second, we compute the cross-correlation between the intensities of the lasers in the population with the highest and the lowest threshold, which are expected to have the two waveforms that differ most. The results are shown in figure 8 the variance, shown in panel (a), decreases logarithmically with the noise strength σ, until it reaches a resonance point around σ ≈ 4 × 10⁻³. When increasing the noise level further, the variance rapidly increases again. This resonance effect is more pronounced when the population is less homogenous and the different δ_i are further apart—this is clear when comparing the variances for $\Delta\delta = \max(\delta_i)-\min(\delta_i) = 0.1$ (blue diamonds) and Δδ = 0.01 (red circles). The variance is minimal for slightly lower noise values in a more homogeneous population, but this effect is not large. Moreover, the population size does not seem to have much influence: the variance and correlation are almost identical in populations of 11 and 51 lasers (blue diamonds and black dots, respectively).

Zoom In Zoom Out Reset image size

The effect is dynamical: looking at the amplitude of the oscillations of the slow variable $w_{\max}$, and the mean period T, both decrease linearly with the logarithm of the noise strength lnσ (shown in figures 8(c) and (d)). This is an effect of so-called 'early jumps' [73] away from the repelling part of the canard solution. An intuitive explanation is that the fast variables move away exponentially slow from the repelling branch, $(x(t),y(t))\approx (x_0,y_0)e^{\lambda(w(t))t}$, and thus, a perturbation of $x_0 = \mathcal{O}(\sigma^2))$ results in an approximately logarithmic decrease of both the time and distance travelled along the repelling branch (see also [71]). At resonance, the maximal amplitude point $w_{\max}$ approaches the maximal value of w on the attracting branches (the difference is of the order of the typical variations due to noise). Thus, past the resonance point, the system no longer follows any repelling branch of the critical manifold, and is thus no longer on a canard trajectory.

Past the resonance point, we observe different scaling: both the amplitude $w_{\max}$ and the mean period T show a decrease with the noise strength σ that is closer to linear, as shown in figures 8(b) and (d). The scaling can be understood from the following analysis. Since canard trajectories have essentially disappeared in the noisy system, the switch between two branches involves leaving the still attracting OFF branch towards the attracting ON branch. Thus, the system needs to leave the basin of attraction of the OFF branch, which is determined by the stable manifold of the repelling part of the ON branch. The distance between this manifold and the ON branch is linear with w in first order, and thus the maximal jumping point $w_{\textrm{max}}$ decreases linearly with the noise strength σ.