Romain Martinenghi

Sergei Rybalko

Maxime Jacquot

Kouomou Yanne

Laurent Chembo

Yanne K Chembo

Laurent Larger

Photonic nonlinear transient computing with multiple-delay wavelength dynamics

Keywords: numbers: 42.54.Sf, 05.45.-a, 89.70.Eg

The brain research and neural network computing communities proposed independently in the early 2000 novel computational principles [1], which are suspected to mimic actual calculation and processing tasks that have been observed and studied in the brain. This computational principles referred to as Echo State Network (ESN [2]), or Liquid State Machine (LSM, [3]), and also with the generic term Reservoir Computing (RC), are definitely different with respect to the standard Turing Machine principles widely implemented in electronic digital processors. Instead of processing the calculation tasks step by step with static states stored in memories, this new principle is based on computational power developed by complex nonlinear transient motion developed in the high dimensional phase space of a nonlinear system excited by an input signal representing the information to be processed. The complex dynamics is usually materialized by a network of neurons (as in the brain), or by any spatially extended network of coupled nonlinear dynamical nodes.

The corresponding generic architecture is depicted in Fig. 1(a), where strong similarities can be seen compared to standard recurrent neural networks (RNN): an input layer is dedicated to the injection of the input information (input connectivity matrix W I ) into a complex interconnected network of dynamical nodes (internal network connectivity matrix W D ); an output layer (readout matrix W R ) is dedicated to the extraction of the result, computed from the nonlinear transient developed by the network dynamics consequently to the injected input signal. Since one of our aim is to transpose these concepts into Physics and into a real world experimental demonstrator, our system will be referred in the remaining part of the paper, as NTC (Nonlinear Transient Computing [4]). NTC is suggested with the intention to reflect more clearly the actual physical origin of the approach, in a way which is expected to be more meaningful for physicists and for the nonlinear dynamics community, although historically the computer or brain science community are referring to ESN, LSM or RC. There are at least two strong differences in NTC compared to RNN: (i) the internal network structure W D is assumed to be fixed, and not to be optimized according to a learning procedure as it is for RNN. This feature is particularly suited to a physical implementation, since reconfiguration flexibility typically brought by a numerical implementation of the RNN, is not required anymore. (ii) A second difference compared to RNN consists in the fact that the learning procedure in NTC is strongly simplified. It is reduced to the determination of Read-Out coefficients (matrix W R ), which can be easily calculated from a simple regression technique.

A significant innovation in NTC was recently proposed and implemented electronically [5] via the introduction of a delay dynamics instead of the classical network of dynamical nodes. This implemented original theoretical solution appeared also as a convenient technical one, since it provides a simpler way for real-world physical implementation, compared to the usually adopted topology of a network [6]. Such an approach can be justified by a known analogy between a delay dynamics and a spatio-temporal one (as a RNN). Both are infinite dimensional, and the space-time representation of a delay dynamics was already proposed 20 years ago [7], introducing on the one hand a discrete time variable corresponding to one delay step forward, and on the other hand a virtual continuous space variable corresponding to the short time-scale fluctuations within a time delay interval. As illustrated in Fig. 1(b), virtual nodes in a delay dynamics can be defined as temporal positions within a time delay interval τ D , which are separated by a "node distance" δτ [5]. This fixes the number of virtual nodes within a time delay intervalle, N = τ D /δτ (150 in our case). The connectivity of the resulting virtual spatio-temporal network is achieved via two dynamical mechanisms. The neighboring nodes are linearly coupled via the characteristic time τ of the oscillator impulse response h(t) (δτ ≃ τ /5 for an optimal adjacent node spacing, as found in [5]). Additionnally to this "short distance" coupling, long distance ones are also present. The proposed multiple delay topology results in a denser internal dynamics connectivity of the equivalent network, with however still a sparse connectivity involving only N I nodes among the N possible ones. Each node is thus subject to a nonlinear transformation (function f [x]) applied to a linear combination, with random weights w D i , of a few N I "previous nodes" randomly defined among the N possible ones. Each corresponding delay is defined as

k i δτ D , k i ∈ [1, . . . , N ] for i = 1 to N I (= N/10).
Compared to the very recent NTC demonstrations with delay dynamics [5,8,9], we propose an original photonic implementation of RC where the dynamical variable is the wavelength of a tunable laser, and also a multiple delayed feedback topology intended to enhance the internal connectivity of the equivalent dynamical network. While evaluating the processing capability of our approach through its performance on a standard benchmark test of spoken digit classification, we found that our multiple delay photonic system with a reduced number of nodes, exhibits a comparable computational efficiency, with word error rate (WER) of the same order as the best results achieved so far for the same test.

The operating principles of our NTC can be summarized as follow [see Fig. 1(b)]. An input information to be processed is exciting each virtual spatial node of the NTC dynamics; this node addressing procedure is practically performed by temporal division multiplexing, consisting in a random spreading of each input sample over all the nodes within a time delay interval. The spreading is ruled by a sparse and random input connectivity matrix W I applied to each of the original information samples. Each sequence of a spread input information undergoes: (i) a nonlinear transformation, operating on a randomly weighted linear combination of N I delayed feedback; (ii) a linear filtering via an impulse response h(t), which is limiting the rate of change in the feedback loop. The computed output of the NTC is obtained by a read-out procedure, corresponding to a linear separation of the transient motion in the N D -nodes dynamics phase The experimental setup is depicted in Fig. 2. It is based on an optoelectronic architecture originally designed to display an Ikeda-like nonlinear delay dynamics [10]. The nonlinear function (f [x] = sin 2 (x + Φ 0 )) is provided by a tunable interference phenomena: an imbalanced birefringent interferometer is seeded by the light emitted at ∼ 1.5 µm by a tunable two-electrode Distributed Bragg Reflector (DBR) laser diode. The wavelength is varied through an injection current I DBR thus changing the interference condition, while another electrode (current I) serves as the usual input for setting the output power. The interferometer output is detected by a photodiode, and an electronic feedback circuit performs the bandpass filtering ruling the dynamics of the oscillator (impulse response h(t)). This filtered signal is the one used for the NTC read-out. The input information is added to the filtered signal, the resulting sum being nonlinearly transformed through f [x], and multipledelayed by an FPGA (Field Programmable Gate Array) board. The FPGA is programmed to perform a flexible and reconfigurable multiple delay line, in which several elementary FIFO (First In First Out) memories are implemented together with weights w D ki for each delay. The weighted sum of the multiple delayed signal is amplified, combined with an offset, and serves finally as the input drive I DBR of the laser wavelength tuning electrode. The normalized dynamics can be written as follows:

1 θ t t0 x(ξ)dξ + x(t) + τ dx dt (t) = x 0 + ( 1
)
β sin 2 NI i=1 w D ki [x(t -k i δτ )] + u(t) + Φ 0
where x is the output of the bandpass filter (also the signal used for the read-out), u is the input data, Φ 0 is an offset phase ruling the operating point along the nonlinear sin 2 -function, β is the loop gain set by the feedback amplification, and x 0 is a formal integration constant which garanties a zero mean value for the bandpass filter output. The delayed feedback coefficients w D ki are programmed in the FPGA, they are randomly defined from a uniform distribution; as already stated, they are introduced for an enhanced internal connectivity compared to the single delay case [5,8,9]. The response times θ = 1.59 ms and τ = 7.95 µs are associated to the low and high cut-off frequencies (100 Hz and 20 kHz) of the bandpass filter respectively, and δτ ≃ 1.59 µs is the elementary time delay determining the spacing between two adjacent virtual nodes. There are N = 150 different nodes, and N I = 15 randomly distributed delayed feedbacks, the largest delay being N δτ ≃ 238 µs.

In order to ensure a stable fixed point solution in an input-free operation, the feedback gain β has to be set below the oscillation threshold defined as 1, calibrated experimentally as the minimum gain leading to a rising oscillation, while the nonlinear function is operating around its maximum linear slope. In the experiment β is set to ca. 0.5. The amplitude of the input information is set so that it induces a large scan of the nonlinear function (around 1-2 extrema, typically ∆u ≃ π).

Each single input information consists of a sequence of samples, each of which has to be spanned over all the nodes within a time delay. The spanning is performed by an input connectivity matrix W I , randomly but uniquely defined for each input information. The temporal waveform to be injected in the setup, is derived from a computed matrix product between W I and a M c , the latter being a 2D representation of the input information (the digit cochleagram in the case of the spoken digit recognition task, see below). The resulting matrix, after being horizontally unfolded, results in a 1D waveform u(t), which is injected into the multiple delay dynamics via the programming of an arbitrary waveform generator (AWG, Lecroy ArbStudio 1102). For each input information signal, the full transient response is recorded by a digital scope. An off-line post-processing is then performed for both the training and testing stages. The training consists in the resolution of a regression problem, and results in the definition of N read-out coefficients w R k that are leading to an optimally correct output, for each input information belonging to a training subset. The efficiency of such an optimal read-out is finally evaluated on a complementary subset of input/output pairs. In principle, the FPGA could also be programmed to implement a direct on-line read-out, as soon as the coefficients are known. This configuration would result in a real-time processing, for which the processing time would be only limited by the analogue bandwidth of the delay dynamical system. For the sake of simplicity, this testing phase was fully processed off-line, right after the training phase.

In order to evaluate the processing efficiency, we performed a standard task typically used in other RNN or NTC reports [5,8,9], and consisting in spoken digit recognition. The main goal of this standard classification test is to recognize a pronounced digit among the ten possible ones from 0 to 9. The spoken digit data base corresponds to 500 digits extracted from the TI46 speech corpus [11]. The digits are pronounced by 10 different female speakers uttering 5 times the 10 digits. Following a standard pre-processing task typically performed in many similar acoustic speech recognition task, the 1D acoustic waveform sampled at 12.5 kHz is transformed into a 2D frequency-time representation (matrix M c ), a so-called cochleagram (the Lyon Ear Model), which provides a monitoring of the average acoustic spectrum (N f = 86 frequency channels) evolving during the spoken digit pronounciation time (N s is typ. 80 samples, varying between 32 and 130). The cochleagram represents the input information signal to be spread over the nodes of the NTC dynamics. This spreading is achieved via a sparsely and randomly determined connectivity matrix W I with dimensions N × N f (sparsity 0.1, and non-zero elements being ±1 randomly distributed). The 2D input data to be injected in the NTC thus consists of an N × N s matrix M u (see Fig. 3), with a number of columns representing the digit duration (variable number of samples N s ). Each value of the N elements (the index of a virtual node within τ D ) of a column, is thus built according to W I . The action of W I on the cochleagram can be interpreted as a random contribution of some selected frequency components at some selected time. The resulting input matrix parameter λ = 10 -3 ),

M u = W I × M c is converted in a 1D signal, u (
W R opt = arg min W R (W R ) t M x -B 2 + λ W R 2 , (2)
for which a training subset of 475 digits is used among the 500 available in the data base. The complementary subset of 25 digits is used to test the performance of the learnt W R opt . This training and testing procedure is repeated for the 20 different possible partitions (crossvalidation). This results in the calculation of a WER, which is statistically limited by the size of the 500 tests actually performed during the cross-validation. Excellent performances (within the limit of statistical significance of 0.2%) have been obtained in terms of WER for such a complex classification test. Figure 4 displays the theoretical and experimental transfert functions of our multiple delay feedback loops with random weights. Their complex frequency dependence confirms that the number of feedback loops, as well as the correct determination of the read-out weighting coefficients is a critical process for the classification task. As shown in Fig. 4, the average WER performance of our setup is 0.6±0.2% with only 15 feedback delay lines, thereby implying that multiple feedbacks enable to achieve performances comparable to state-of-the-art with a relatively limited number of read-out parameters to compute. This result validates the computational efficiency of our photonic NTC setup performed by a multiple delay wavelength dynamics, and it opens the way to the future optimization of delay dynamics NTC processors, via a more accurate conceptual analogy between multiple delay dynamics and classical network dynamics.

We have demonstrated experimentally the efficiency of a photonic neuromorphic processor. The proposed setup is based on a multiple delayed feedback wavelength dynamics providing an enhanced dynamical connectivity. Additionnally to recent success in applications of complex dynamics [12][13][14], photonic nonlinear delay dynamics are confirmed as an efficient and flexible solution for the practical implementation of NTC. Future work will be devoted to ultra-fast photonic versions making use of standard optical telecommunication devices and principles [13]. Processing speed for the recognition of a single spoken digit is around 20 ms, but this can be improved down to 100 ns with telecom grade devices. Photonic NTC based on multiple delay dynamics should offer an efficient and ultra-fast hardware solution for future neuromorphic computers. Many fundamental issues still need to be addressed (quantitative connections between delay dynamics features and computational power, or optimal topology and system architecture for advanced functionalities such as plasticity and integrated learning capability). We also anticipate that information theory and nonlinear dynamics should be explored together in order to provide a theoretical framework capable of describing and optimizing the computational efficiency provided by delay dynamics-based NTC [15].

FIG. 1 :

 1 FIG. 1: (Color online) Principles of an NTC. Left: with a spatially extended dynamical network. Right: with a multiple delay feedback dynamics.

FIG. 2 :

 2 FIG. 2: (Color online) Optoelectronic wavelength setup performing an NTC as a multiple delay feedback dynamics.

FIG. 3 :

 3 FIG.3:(Color online) Illustration of the NTC processing on the spoken digits; example of the spatio-temporal representation for a digit "seven".

FrequencyFIG. 4 :

 4 FIG. 4: (Color online) Experimental results. Left: transfert functions for 15 (red) and 150 (black) delayed feedbacks with random weights (theoretical curves in full lines, experimental data in symbols). Right: WER histograms with 15 delayed feedbacks and a total of 120 cross-validations.

This work was supported by the PHOCUS project (FP7-240763), and the Région Franche-Comté. LL thanks I. Fischer and C. Mirasso for helpful discussions.


Lionel Weicker

Thomas Erneux

Otti D'huys

Jan Danckaert

Maxime Jacquot

Yanne Chembo

Laurent Larger

Strongly asymmetric square-waves in a time delayed system

Keywords: numbers: 02.30.Ks, 02.30.Mv, 02.30.Oz, 05.10.-a, 42.65.Sf

Time-delayed systems are known to exhibit symmetric square-waves oscillating with a period close to twice the delay. Here, we show that strongly asymmetric square-waves of period close to one delay are possible. The plateau lengths can be tuned by changing a control parameter. The problem is investigated experimentally and numerically using a simple bandpass optoelectronic delay oscillator modeled by nonlinear delay integro-differential equations. An asymptotic approximation of the square-wave periodic solution valid in the large delay limit allows an analytical description of its main properties (extrema and square pulse durations). A detailed numerical study of the bifurcation diagram indicates that the asymmetric square-waves emerge from a Hopf bifurcation.

In the last three years, intensified efforts have been developed to explore the physical effects of time delays [START_REF] Erneux | Applied Delay Differential Equations[END_REF]- [START_REF] Smith | An Introduction to Delay Differential Equations with Applications to the Life Sciences[END_REF]. In biology, time delays are called maturation periods, control reflexes, cellular transport times, and they appear in various areas including neural networks, where propagation delays of the electrical signals connecting different neurons are taken into account [START_REF]Delayed complex systems[END_REF]. In the industry, machine tool chatter is a persistent problem in metal cutting applications, where a delay term is typically involved in the applied force [START_REF] Stépan | Dynamics and Chaos in Manufacturing Processes[END_REF]. Optical feedback from a distant reflector is another example where a time delay is causing erratic intensity oscillations in lasers used in optical communication systems [START_REF]Unlocking Dynamical Diversity: Optical Feedback Effects on Semiconductor Lasers[END_REF]. Time-delayed feedbacks are also used for practical applications. They might be deliberately implemented to control neural disturbances, e.g. to suppress undesired synchrony of firing neurons in Parkinson's disease [6]. They are also designed for producing new effects with optoelectronic oscillators (OEO). An OEO typically incorporates a nonlinear (intensity) modulator, an optical fiber delay line, and an optical detector followed by a bandpass filter, in a closed-loop resonating configuration. This hybrid microwave source is capable of generating, within the same optoelectronic cavity, either an ultra-low-jitter (low phase-noise) single tone microwave oscillation (narrow band filtering), as used in radar applications [START_REF] Lasri | [END_REF], or a wide band chaotic carrier (broadband filter) typically intended for physical encryption in high bit rate optical communications [10].

Because of the large variety of applications, there is an increasing interest for problems modeled by delay differential equations (DDEs) and their specific dynamical phenomena. A fundamental property of nonlinear dynamical systems controlled by a delayed feedback is their tendency to exhibit square-wave oscillations if the delay τ D is sufficiently large. These oscillations typically consist of 2τ D -periodic transitions between two or more flat plateaus. The following scalar DDE

εx ′ = -x + f (x(s -1), λ), (1)
has intensively been studied for its square-wave solutions (prime stands for the derivative with respect to the time s). Equation ( 1) arises in a variety of applications, for example, physiological control systems [11], the transmission of light through a ring cavity [12], and population biology [13]. Here the dimensionless time is s ≡ t/τ D and ε ≡ τ /τ D > 0 is a small parameter defined as the ratio of the linear decay time τ to the delay τ D . f (x, λ) denotes a nonlinear function of x, and λ is a control parameter. Setting ε to 0 reduces Eq. (1) to the following map relating x n = x(s) and x n-1 = x(s -1)

x n = f (x n-1 , λ). (2)
The two fixed points of the Period-2 solution of Eq. ( 2) are expected to approach each of the two plateaus of the square-wave solutions of Eq. (1) if ε → 0. Significant contributions to the asymptotic relations between the solutions of the map (2) and the solutions of the DDE (1) have been made by Chow and Mallet-Paret [START_REF] Chow | Coupled Nonlinear Oscillators[END_REF], Mallet-Paret and Nussbaum [START_REF] Mallet-Paret | [END_REF], Chow et al. [16] and Hale and Huang [17]. Close to a Hopf bifurcation point λ = λ 0 (ε) of the DDE (1), the oscillations quickly change their shape from sinusoidal to square-waves as the amplitude of the solution increases [START_REF] Erneux | [END_REF]. The square-wave oscillations then consist of sharp transition layers of size proportional to ε, connecting the two plateaus, of lengths close to unity.

In this letter, we address the question whether periodic square-wave oscillations of different plateau lengths are possible. These types of patterns may appear in systems modeled by Eq. (1) for particular functions f (x, λ) and judicious choices of the initial conditions. But they are unstable for long time [START_REF] Sharkovsky | Difference Equations and their Applications[END_REF] and have been called metastable [START_REF] Erneux | Applied Delay Differential Equations[END_REF]. Some of our recent experiments have revealed the stable occurrence of such asymmetric squarewaves using an optical system that is modeled by a second order DDE instead of Eq. (1). A particular feature of these square-waves that none of earlier investigations did anticipate, is the fact that they are nearly τ D -periodic. Figure 1 exhibits experimental stable τ D -periodic waves that were obtained from a broadband bandpass opto electronic oscillator (OEO). The setup is depicted in the same figure and is a straightforward implementation of the broadband chaos generator described in [10]. It consists of an integrated optics Mach-Zehnder modulator (MZM) used to provide the nonlinear transformation as an electronically tunable two-wave interference. A DC bias can be adjusted externally to fix the zero drive interference condition (parameter Φ). The MZM output is delayed in time by a 4km fiber spool (τ D ≃ 20 µs), and then detected by a photodiode. A simple self-made electronic circuitry performs amplification as well as the dynamical limitation for the signal to be fed back on the MZM. The dynamics of the oscillation loop is thus controlled via its electronic Fourier filtering feature, which is forced to a broadband bandpass filter. The characteristic response time of the corresponding cascaded high pass and low pass first order filters are ca. 0.8 s (θ) and 0.2 µs (τ ), respectively. The normalized gain β of the oscillation loop is adjusted through the injection current of the CW 1.5µm DFB telecom semiconductor laser seeding optically the MZM. The plateau lengths of the waveform in Fig. 1 can be changed gradually by changing a control parameter (Φ or β), but the overall pattern remains τ D -periodic. Numerical simulations of Eqs. ( 3) and (4) for long intervals of time reproduce these asymmetric waves but offer no clue on their nature. In order to have an analytical understanding on their properties, an asymptotic analysis of the OEO equations is proposed.

In dimensionless form, the evolution equations are given by (see [20] for details)

y ′ = x, (3)
εx ′ = -δy -x + β cos 2 (x(s -1) + Φ) -cos 2 Φ , (4)
where s = t/τ D is time measured in units of the delay. ε = τ /τ D and δ = τ D /θ are two small fixed parameters. The two control parameters are the feedback gain β and the feedback offset phase Φ. Under specific initial conditions and parameter values, Eqs. ( 3) and ( 4) admit square-wave solutions for x with different plateau lengths (see Fig. 2). The periodic solution is nearly 1-periodic. We denote the length of the two plateaus by s 0 and 1-s 0 , respectively. x 01 and x 02 mark the extrema of x. The value of Φ is -π/4 + 0.1. If Φ = -π/4, the square-wave is perfectly symmetric (s 0 = 1/2). Taking advantage of the small value of ε, we propose to construct a periodic solution satisfying the condition where P = 1 + 2εr(ε) (r = O(1)). The leading approximation is obtained by setting ε = 0 in Eqs. ( 3)-( 4). The reduced equations with (5) and P = 1 are

x(s -P ) = x(s), (5)
y ′ = x, (6)
y = δ -1 -x + β cos 2 (x + Φ) -cos 2 Φ . (7)
The function ( 7) is represented in Fig. 2(b) and exhibits three branches provided that β > 1. The evolution of x and y along the left and right branches can be determined analytically by inserting ( 7) into ( 6) and solving the resulting first order equation in x. However, simple expressions of the solution can be obtained by taking advantage of the small value of δ. Specifically, we seek a solution of the form

y = δ -1 y 0 + y 1j (s) + ..., (8)
x = x 0j + δx 1j (s) + ..., (9)
where j = 1 or 2 refers to the time domains 0 < s < s 0 and s 0 < s < 1, respectively (see Fig. 2(a)). The leading order problem as δ → 0 is

y ′ 0 = 0, (10)
-x 0j -y 0 + β cos 2 (x 0j + Φ) -cos 2 Φ = 0. ( 11
)
We already know that for a finite interval of y 0 , Eq. ( 11) admits more than one solution (see Fig. 2(b)). The solutions corresponding to the left and right branches are denoted by x 02 < 0 and x 01 > 0, respectively. Their values are unknown because y 0 is still unknown. This motivates the analysis of the next problem for y 1j and x 1j given by

y ′ 1j = x 0j , (12)
-x 1j -y 1j -β sin(2x 0j + 2Φ)x 1j = 0, ( 13
)
where j = 1 or 2. We solve Eq. ( 12) with the initial conditions y 11 (0) = 0 and y 12 (s 0 ) = y 11 (s 0 ). We then determine x 11 and x 12 using Eq. ( 13). We obtain

y 11 = x 01 s, (14)
y 12 = x 01 s 0 + x 02 (s -s 0 ), ( 15
)
x 11 = - x 01 s 1 + β sin(2x 01 + 2Φ) , ( 16
)
x 12 = - (x 01 s 0 + x 02 (s -s 0 )) 1 + β sin(2x 02 + 2Φ) . ( 17
)
The 1-periodicity condition requires that y 12 (1) = 0 and from [START_REF] Mallet-Paret | [END_REF], we obtain an equation relating the unknowns x 01 , x 02 , and s 0

x 01 s 0 + x 02 (1 -s 0 ) = 0. ( 18
)
Moreover, continuity of x 11 and x 12 at s = s 0 requires cos(x 01 + x 02 + 2Φ) sin(x 01 -

x 02 ) = 0. ( 19
)
This is a first condition for the unknowns x 01 and x 02 . Because x 01 = x 02 , the solution of Eq. ( 19) is x 01 + x 02 + 2Φ = ±π/2 + 2nπ. The numerical solution in Fig. (2) verifies the condition

x 01 + x 02 + 2Φ = -π/2. ( 20
)
We need a second equation relating x 01 and x 02 . Substracting Eqs. ( 11) with j = 1 and 2 gives

-(x 01 -x 02 ) -β sin(x 01 + x 02 + 2Φ) sin(x 01 -x 02 ) = 0. ( 21
) Using ( 20) then allows to eliminate x 02 in (21):

-(2x 01 + 2Φ + π/2) + β sin(2x 01 + 2Φ + π/2) = 0. (22)
This provides x 01 = x 01 (β) which can be determined in the implicit form β = β(x 01 ). The values of x 02 and s 0 are then found by using ( 20) and ( 18) (see Fig. 3). Furthermore, using (22) in (11) with j = 1, one gets Equation ( 11) admits a third root x 03 (x 02 < x 03 < x 01 ), which can be verified from ( 11) and ( 23) as being x 03 = -(Φ + π/4). Combining this value with (20) leads to the roots symmetry property

y 0 = Φ + π 4 - β 2 cos(2Φ) (23)
z 0 = x 01 -x 03 = x 03 -x 02 (24)
which will play a crucial role for the transition layers (see below).

In summary, we concentrated on the branch of strongly asymmetric square-waves. As β progressively decreases from a value β larger than one, we note from Fig. 3 (a) and (b) that both the extrema of the square-wave and the plateau length s 0 decrease to zero as β → 1.

In order to complete our asymptotic analysis, we need to investigate the fast transition layers between x 01 and x 02 . To this end, we introduce the transition layer variable ζ 1 ≡ s/ε (jump down) or ζ 2 ≡ (ss 0 -εr)/ε (jump up) into Eqs. ( 3) and (4). The correction term εr represents the contribution to the period of each transition layer (P = 1+2εr) and r needs to be determined from the transition layer equations. From the leading order equations as ε → 0, we find that δy = y 0 where y 0 is defined by (23) and a first order DDE for x. Because of (24), we note that z(ζ 2 ) = -z(ζ 1 ) where z ≡ xx 03 . We then find that the transition layer problems are described by the same equation for z and only differ by their boundary conditions. For the jump down transition, we find

z ′ = -z + β 2 sin(2z(ζ 1 + 2r)), (25)
where prime means differentiation with respect to ζ 1 . The boundary conditions are z(±∞) = ∓z 0 . Equation ( 25) is covered by the general theory of Hale and coworkers [START_REF] Chow | Coupled Nonlinear Oscillators[END_REF]17,21]. Specifically, we seek a heteroclinic orbit for some value of r, that is a trajectory joining z 0 and -z 0 at times ζ = ∓∞. In the limit, β → 1 we find that r → 1/2. For an arbitrary value of β, the transition layer equation needs to be solved numerically. Close to β = 1, the asymptotic analysis assuming s 0 = O(1) and fast transition layers of size ε is no more valid. The smallest value of β where we numerically found the square-wave is close to β = β c ≃ 1.01 (see Fig. 4). Note that the maxima of x are close to the steady state value x = 0. As we further decrease β, we find x = 0 as the only stable solution. A systematic study of the bifurcation diagram by integrating numerically Eqs. ( 3) and ( 4) and by using a continuation method (DDE-BIFTOOL [23]) indicates that the asymmetric square-wave solutions belong to a Hopf bifurcation branch that emerges from x = 0 at β = 1.020 with a period close to 1. It is the first Hopf bifurcation from the basic steady state x = 0. The branch of periodic solutions is subcritical and unstable but quickly folds back to stable asymmetric square-waves as seen in Fig. 4.

The asymmetric square wave regime is not the only stable periodic solution of Eqs. ( 3) and ( 4). Although the number of primary Hopf bifurcation points increase as β is increased from β = 1.020, we found that only two bifurcations lead to stable branches: the branch of asymmetric squares-wave solutions studied in this letter and a branch of low frequency oscillations that emerges from a Hopf bifurcation point at β = 1.025 with a period close to 70 (≃ δ -1 , related physically to θ as already found in [20]). These two branches dominate the bifurcation diagram as we gradually increase β from 1 or decrease β from a large value. No other stable time-dependent regimes have been found in the interval 0 < β < 2. Details on the analysis of the transition layers and on the numerical bifurcation diagrams will be published elsewhere.

A similar asymmetric square-wave has recently been explored both experimentally and numerically for an OEO operating under quite different parameter values [22]. We are currently looking if a similar Hopf bifurcation mechanism is responsible for these oscillations.

T.E. acknowledges the support of the F.N.R.S. (Belgium). L.W. acknowledges the Belgian F.R.I.A. for a PhD scholarship. This work was supported by the European project PHOCUS (FP7 grant 240763). L.L. thanks the support of the Institut Universitaire de France.

FIG. 1 :FIG. 2 :

 12 FIG.1:(Color online) Upper: τD-periodic square-wave oscillations for three values of the feedback phase φ increasing from left to right. The duty cycle is monotonically modified, but the period is unchanged. Lower: experimental setup performing a bandpass electro-optic Ikeda dynamics.

FIG. 3 :

 3 FIG. 3: (Color online) (a) and (b): Analytical bifurcation diagram for the extrema of x (x01 and x02) and the length of the first plateau (s0). The vertical broken lines mark the value of β = 1.2 used in the simulations in Fig. 2(a). The predicted values of the extrema and the length of the first plateau are in excellent quantitative agreement with the values documented from the numerical simulations.

FIG. 4 :

 4 FIG. 4: (Color online) Same values of the fixed parameters and initial conditions as in Fig. 2 except β = 1.01. (a): Numerical square-wave solution exhibiting a small plateau length s0. (b): Periodic solution in the phase-plane. The dot marks the steady state (x, y) = (0, 0).