Tuning the period of square-wave oscillations for delay-coupled optoelectronic systems

We analyze the response of two delay-coupled optoelectronic oscillators. Each oscillator operates under its own delayed feedback. We show that the system can display square-wave periodic solutions that can be synchronized in phase or out of phase depending on the ratio between selfand cross-delay times. Furthermore, we show that multiple periodic synchronized solutions can coexist for the same values of the fixed parameters. As a consequence, it is possible to generate square-wave oscillations with different periods by just changing the initial conditions.


I. INTRODUCTION
Time delays in physical, biological, or chemical systems are known for their oscillatory instabilities [1].For optical and optoelectronic systems with feedback, the time scale associated with the feedback is generally much larger than the intrinsic time scales of the dynamical system.This large delay of the feedback can be useful for applications [2,3].
An interesting dynamical regime that results from a large delayed feedback [4] is square-wave switching [5].The generation of tunable pulsating dynamics has been studied over the past few years [6][7][8][9][10][11][12][13][14], motivated by new applications such as optical clocks [7] and other binary logical applications, generation of stable microwave signals, or optical sensing [8].In particular, stable square waves oscillating in antiphase have been observed in the intensities emitted in each of the polarization directions in edge-emitting diode lasers (EELs) subject to crossed-polarization reinjection (XPR).In this setup, the natural polarization TE mode is rotated 90 • and delayed coupled to the normally unsupported TM mode.Asymmetric square waves with a period close to but longer than twice the coupling delay time have also been reported in mutually coupled EELs, in a scheme where the TE mode of each laser is injected into the TM mode of the other laser [10,11].Moreover, square-wave oscillations have been investigated for vertical-cavity surface-emitting lasers [6,12,15] subject to XPR, semiconductor ring lasers subject to a delayed optical feedback [13], optoelectronic oscillators (OEOs) [16], and mode-locked fiber lasers [14].Other studies on polarized optical feedback have been proposed in [17,18].
One fundamental question raised by all these studies is the possibility to generate square-wave oscillations with a desired period.This question had already been discussed as researchers analyzed the Ikeda delay differential equation [19].It was found that different periodic solutions exhibiting frequencies that are multiples of a basic frequency appear through successive Hopf bifurcations.Ikeda and Matsumoto then realized that such multiple periodic regimes could be used to encode information in high-capacity memory devices [20].Experiments showed that an electro-optical hybrid bistable system with a fiber delay loop could indeed sustain a large number of oscillatory states [21].However, the experiments also produced evidence that such systems were very sensitive to spurious resonances.Particularly disturbing was the fact that a large number of harmonics predicted by numerical simulations were not observed in physical implementations, although they could be recovered by subjecting the system to a periodic modulation with the same frequency as the missing harmonics.This raised the question of the stability of all the bifurcating time-periodic modes and their robustness with respect to external perturbations.In this paper we examine the emergence of stable square waves in a system of two delay-coupled OEOs.Each oscillator operates on its own delayed feedback and the presence of two distinct delays allows a large number of stable square-wave time-periodic regimes.By controlling the ratio of the two delays, we may generate a square wave with a desired period as well as multiple periodic states with different periods coexisting for the same values of the parameters.
The effect of the ratio between different delay times on the synchronization properties of two mutually coupled chaotic lasers was recently investigated in [22].Here we concentrate on the multiplicity of stable periodic regimes generated by the mutually coupled OEOs rather than their synchronization efficiency as they are chaotic.By using asymptotic methods based on the relatively large values of the two delays, we propose a systematic analytical study of their bifurcation mechanisms.The validity of all our results is tested by solving numerically the original evolution equations for the OEOs.Because of the two distinct delays, the bifurcation possibilities are rich, but their derivations are relatively simple because the analysis essentially relies on the solutions of coupled equations for maps.In this sense, we expect that our analysis can be applied to other two-delayed coupled systems and lead to similar results.
Our choice of two mutually coupled OEOs is motivated by the large variety of dynamical regimes [16,23] that are generated by single optoelectronic systems.They have been used as chaos generators for secure chaos-based communications [24][25][26].These devices have also been proposed and studied to produce efficient ultrapure microwaves in the periodic regime [27][28][29][30].In particular, their robustness to noise has been studied both theoretically and experimentally [31].Finally, OEOs operating in the steady-state regime have recently been implemented in an experimental demonstration of a photonic liquid-state machine performing as a kind of neuromorphic computer [32].As compared to optical feedback systems and optical injection systems, in which the dynamics depends on the frequency, phase, and amplitude of the field and for which a frequency detuning of a few hundred MHz between the two systems can lead to a large degradation of the synchronization, optoelectronic systems are more flexible due to their insensitivity to optical phase variations.Another advantage of OEOs is that they can be electrically driven [33].
The paper is organized as follows.In Sec.II we describe the system we are considering and its dynamical model.In Sec.III we develop a theoretical approximation to determine parameter conditions for which synchronized square-wave periodic solutions appear.In Sec.IV we obtain a theoretical prediction for the amplitude of the periodic oscillations.In Sec.V we study the secondary instabilities of the periodic square-wave oscillations.In Sec.VI we compare the theoretical predictions with numerical simulations of the full dynamical model described in Sec.II.In Sec.VII we analyze the effect of a small mismatch in the delay times.Finally, in Sec.VIII we summarize our results.

II. DYNAMICAL MODEL
We consider two electro-optical delay systems [34] that are mutually coupled as shown in Fig. 1.The light emitted by a cw semiconductor laser diode (LD) with intensity P is split into two beams, each beam feeding an electro-optical delay loop.Each loop consists of a Mach-Zehnder interferometer (MZI), an optical delay line, a photodiode (PD), and an amplifier.We use subindex i, i = 1,2, to identify the variables associated with loop i.For loop i the optical output of MZI i is split into two parts.A fraction α ii is delayed using a fiber loop by a time T ii .A fraction α ij with i,j = 1,2 and j = i is injected from loop i into loop j after a delay T ij .Self-feedback and cross-feedback optical signals are combined and the resulting intensity is detected by the PD.The electrical signal goes through a bandpass amplifier and is finally used to drive the Mach-Zehnder ac electrode.For each loop the dynamics results from a combination of the nonlinear effect due to the MZI plus a linear filtering process associated with the electrical part of the loop.The dynamics of each electro-optical system can be described as follows.As in [34], we assume there is no reflection in the optical path, so the optical electric field can be described as a scalar E. Each arm of the MZI can be considered as a Pockels cell exhibiting a linear dependence of the refractive index n on E, where n 0 is the intrinsic refractive index of the material (without an electric field).Light with wavelength λ 0 crossing the Pockels medium changes its phase depending on the refractive index as well as on the length L of the arm: When a voltage V is applied across an arm of transversal size d, the corresponding electric field is given by E = V /d and thus the phase change along that arm can be written as where ϕ 0 = 2πn 0 /λ 0 and V π is the voltage required for a modulation of π in the phase, Considering that the optical field at the input of MZI i is √ P ρ i e ϕ i and that the voltage V i applied to the MZI i has a dc and a rf component V B i and V rf i (t), respectively, and introducing the normalized voltages x i (t) = πV rf i (t)/2V π i and i = πV B i /2V π i , the electric field at the output of the MZI i is given by where φ 0 i = ϕ 0 i + ϕ i .Without loss of generality and for simplicity, we take φ 0 1 = 0 and define 0 = φ 0 2 .The optical field arriving at photodiode PD i is given by The output of photodiode PD i with sensitivity S i is Finally, the voltage V PD i (t) is amplified and filtered by the linear bandpass amplifier with effective gain G i and low and high cutoff characteristic times θ i = 5 μs and τ i = 25 ps, respectively.The normalized voltage at the output of the amplifier x i (t) (which is used to modulate the MZI i ) is given by where the dot stands for time derivatives.One can consider the term 1 + τ i θ −1 i ≈ 1 due to the order of magnitude of the actual time scales.Nevertheless, the integral part cannot be neglected in general since it is responsible for the zero mean value of V , which is achieved after a slow transient characteristic of the bandpass dynamics that removes any dc component; indeed, large transient dynamics of the order of the slowest time θ are observed [35].Combining Eqs. ( 5)-( 8), one gets that the dynamics of the system is ruled by two delay integro-differential equations.The integral terms are eliminated through the introduction of two additional variables y i (t) = t t 0 x i (t )dt leading to a system of four delay differential equations where i,j = 1,2; j = i; and γ ij = G j S j ρ i ρ j α ij are effective coupling strengths.
In this work we address the case in which the two systems have identical parameters , and T 11 = T 22 = T f .We also consider 0 = 0 and that the self-and cross-feedback parameters are the same: γ 11 = γ 22 = γ 12 = γ 21 = γ .Finally we define the coupling time as T c = (T 21 + T 12 )/2.Then the steady-state solution is given by Introducing Y i (t) = [y i (t) − y i st ]/T c and scaling the time with T c , s = t/T c , we get where prime means differentiation with respect to s and We consider that the delay time T c has a value of the order of tens of nanoseconds.Therefore, ε is of order 10 −3 and δ of order 10 −2 .Making use of the trigonometric identity cos Eqs. ( 11) can be simplified as These equations admit time-periodic square-wave solutions.Note that the terms multiplying ε are only important for the fast transition layers between plateaus of the square waves.Fortunately, we may ignore these layers as we determine the leading-order bifurcation equations.In the next two sections, we first concentrate on the Hopf bifurcations of the basic steady state and then derive nonlinear maps for the square waves valid in the limit ε → 0. For the single optoelectronic oscillator, this approach has been successfully applied and tested [36].

III. HOPF BIFURCATIONS OF THE STEADY STATE
Because the delays are large compared to the time scales of each individual oscillator, periodic square-wave oscillations are emerging as the dominant solutions [1].We address now the dependence of the period of these square-wave oscillations on the ratio between the two delays.To do this we determine the Hopf bifurcations of the zero solution.As discussed after Eq. ( 8), the main role of δ is to ensure that the time average of x(t) in the stationary regime is zero.While this is relevant for any nonzero solution, the zero solution satisfies this condition automatically.Therefore, to a first approximation we can neglect the role of δ in the stability analysis of the zero solution.By setting δ = 0 Eqs.( 13) reduce to two coupled equations for x 1 and x 2 , where, as before, i,j = 1,2 and j = i.We start by determining the Hopf bifurcations of the zero solution.Specifically, we consider x i (t) = x st i + u i (t) and formulate the linearized equations for the small perturbations u i , where will play the role of an effective bifurcation parameter as seen below.The solutions of the linearized equations are of the form The zero solution is stable if λ < 0 (perturbations decay in time) and unstable if λ > 0 (perturbations grow).At the Hopf bifurcation, λ = 0.The linearized problem with ε = 0 then reduces to a homogeneous system of two linear algebraic equations for c 1 and c 2 : The condition for a nontrivial solution is given by the determinant which leads to the characteristic equation Thus, there are two families of Hopf bifurcations.First, considering the plus sign in (19) and substituting this into Eq.( 17), one obtains c 1 = c 2 .Thus this Hopf bifurcation leads to oscillations where x 1 and x 2 are in phase.From the real and imaginary parts of ( 19) we obtain two equations for χ and ω: sin(ωs 0 ) + sin(ω) = 0. ( 21) Equation (21) implies that either (i) ωs 0 = −ω + 2nπ or (ii) ωs 0 = π + ω + 2nπ , where n ∈ Z.Only case (i) leads to physical solutions if we consider Eq. ( 20).The solution for Eqs. ( 20) and ( 21) is The fact that s 0 > 0 and ω n > 0 implies that n > 0.
Second, the minus sign in (19) leads to solutions for which c 1 = −c 2 .This is a Hopf bifurcation leading to oscillations where x 1 and x 2 are in antiphase.From the real and imaginary parts of ( 19) we obtain sin(ωs 0 ) − sin(ω) = 0.
Equation (25) implies that either (i) ωs 0 = ω + 2nπ or (ii) ωs 0 = π − ω + 2nπ,where n ∈ Z.Using (24), we note that only case (ii) leads to physical solutions.The solution of Eqs. ( 24) and ( 25) is The fact that s 0 > 0 and ω > 0 implies that n 0. We consider χ > 0. Note from either Eqs.( 20) and ( 21) or Eqs. ( 24) and ( 25) that there is no Hopf bifurcation solution if s 0 = 0 (the system has a single delay). 1 We next consider s 0 as a control parameter and wish to determine the Hopf bifurcation points in terms of χ n .There is a family of Hopf bifurcation curves χ n (s 0 ).In Fig. 2 we show the curves with n = 1, 2, and 3 for the in-phase and antiphase solutions.All the curves exhibit a minimum at χ n = 1 and close to the minimum have a parabolic shape.From a physical point of view we are particularly interested in determining the possible Hopf bifurcations that appear exactly at χ n = 1 for specific values of s 0 since this corresponds to the first Hopf bifurcation that is encountered when increasing χ n , proportional to the LD power P , in a system with given delay times.This situation corresponds to cos ω n = −1 in (23) or cos ω n = 1 in (27).For the first case (in-phase oscillations), the condition cos ω n = −1 implies ω n = (1 + 2m)π , where m ∈ Z. From (22) we then obtain the following values of s 0 : where l ≡ n − m − 1.The condition s in 0 > 0 restricts the value of m to the range 0 m < (2n − 1)/2, thus l 0. The period of the Hopf bifurcation oscillations is determined using (22) and is given by From ( 28) one has namely, for in-phase periodic solutions the dimensionless time difference (T c − T f )/T c = 1 − s 0 has to be an even or odd rational number, and using (28) the period can also be written as For the second case (out-of-phase oscillations), the condition cos ω n = 1 implies ω n = 2mπ where m ∈ Z. From (26), we then determine the following values of s 0 , where k ≡ n − m.The condition s out 0 > 0 now restricts the value of m to the range 0 < m < (1 + 2n)/2, thus k 0. The period of the Hopf bifurcation oscillations is found using (26) and is From ( 28) one has that is, the dimensionless time difference 1 − s 0 has to be an odd or even rational number.From ( 32) and ( 34) the period can be written as A few observations are worth pointing out.First, in-phase and out-of-phase Hopf bifurcation points never appear for the same s 0 .This can be verified by equating (28) with n 1 > 0 and m 1 0 and (32) with n 2 0 and m 2 > 0. After simplifying, we find 4m 2 n 1 = (2m 1 + 1)(1 + 2n 2 ), which implies that an even number should be equal to an odd number.Second, several in-phase Hopf bifurcations or several out-of phase solutions may appear for the same value of s 0 with different periods.This means that the Hopf bifurcation at χ n = 1 can be multiple.This degeneracy is not eliminated if we consider ε = 0.For example, if s 0 = 1, Hopf bifurcations to in-phase solutions appear if m = (n − 1)/2 (n > 0, m 0).From the conditions for a Hopf bifurcation with ε = 0, we find χ n = √ 1 + ε 2 and all the frequencies satisfying tan ω n = −ε (ω n nπ − ε).These two results are illustrated in Fig. 3, where the period of Hopf bifurcation points emerging from χ n = 1 is plotted as a function of s 0 .In this figure, points with the same m are located in horizontal lines, points with the same l or the same k are located in straight lines of positive slope starting from the origin (s 0 = 0,T = 0), and points with the same m − l or the same m − k are located in straight lines that start at (s 0 = 1,T = 1).
We now discuss in more detail the coexisting solutions for a given value of s 0 .We note that once s 0 is fixed, not all the values of l and m are possible.For instance, for an in-phase solution with s 0 = 5/7 smaller possible values for l and m are l f = 2 and m f = 3, so the fundamental solution has a period T in f = 2/7.The first harmonic is associated with values of l and m that can be obtained by multiplying both the numerator and denominator of s 0 = 5/7 by 3, namely, 15/21, so that l 1 = 7 and m 1 = 10 and the period is T in 1 = T in f /3 = 2/21.Notice that multiplying by an even number will lead to a ratio of two even numbers, which does not fulfill the condition (28).In general, for a given s 0 the fundamental in-phase solution is given by the minima l and m that fulfill (28).Higher harmonics are obtained by multiplying the numerator and denominator of ( 28) by an odd number.Therefore, the values allowed for m and l are given by where j stands for the order of the harmonic.The period of the harmonic is where we have used (31).For out-of-phase solutions a similar argument applies.For a given s 0 the fundamental solution is given by the minima k and m that fulfill (32), while the harmonic of order j can be obtained by multiplying the numerator and denominator of ( 32) by (2j + 1); therefore, the values allowed for m and k are given by and the period is where we have used (35).Finally, we note that, according to Eqs. ( 28) and (32), inphase and antiphase periodic solutions appear when the ratio between the self-and cross-feedback delays is odd/odd and odd/even, respectively.Each OEO can be seen as having two delay times τ 1 = T f and τ 2 = 2T c .Therefore, in-and out-ofphase solutions exist when the ratio τ 2 /τ 1 is even/odd.This is in agreement with the conditions for synchronization found for coupled chaotic lasers [22].

IV. NONLINEAR MAPS FOR PRIMARY PERIODIC SQUARE-WAVE OSCILLATIONS
Following the understanding gained in the previous section on the period of the solutions arising when the zero state becomes unstable, we now develop a map to determine the amplitude of these solutions.To this end, we plan to formulate equations for a map relating x 1 (s) to x 1 (s − T /2), where T is the period taking advantage of the small value of ε.Note that we need to take into account the effect of δ to ensure that x(t) has a zero time average in the stationary regime.
We first consider the case of the in-phase oscillations.With the condition x 2 (s) = x 1 (s) ≡ x(s) and Y (s) ≡ Y 1 (s) = Y 2 (s), Eqs.(13) with ε = 0 reduce to The period T in of the in-phase Hopf bifurcation was defined by (29).Starting with (28) and using (29), we may relate the delays s in 0 to T in , Furthermore, a second useful equation is obtained by relating 1 to T in as where we again used (29).Using ( 42) and (43), Eq. ( 40) considerably simplifies as Equation ( 44) is a first equation for a map relating x(s) and x(s − T in 2 ).To obtain a second equation relating Y (s) and Y (s − T in 2 ), we integrate Eq. (41) from s − T in /2 to s.As illustrated in Fig. 4, we consider square-wave solutions in which x(s) remains practically constant over half of the period.Introducing now the notation , with x k located at the beginning of the semiperiod (see Fig. 4), we have Subtraction of Eq. (45) evaluated at x k and at x k−1 gives Then For the out-of-phase solution, we assume that x 1 (s) and x 2 (s) are T out -periodic solutions where T out is given by (33).
Using the fact that x(s) ≡ x 1 (s) = x 2 (s + T out /2) and Y (s) ≡ Y 1 (s) = Y 2 (s + T out /2) we obtain from Eqs. ( 13) with ε = 0 the following equations for x(s) and Y (s): As for the in-phase solutions, we now use ( 32) and ( 33) and determine useful relations between the two delays s out 0 and 1, and T out .Specifically, we find Using ( 51) and (52), Eq. ( 49) then simplifies as which is the same equation as (44) replacing T in by T out . Introducing 2 ), with x s located at the beginning of the semiperiod one obtains exactly the same map as for in-phase solutions, namely, Eq. (48).

V. NONLINEAR MAPS FOR SQUARE-WAVE OSCILLATIONS GENERATED BY SECONDARY BIFURCATIONS
We now consider a generalization of the map obtained in the previous section in order to describe the instabilities of the primary periodic square-wave solutions seen in the numerical simulations described in the next section.In fact, the map (48) can be extended to a more general class of lagged solutions of the form The lag time 1 − s 0 corresponds physically to the difference T c − T f normalized to T c .In-and out-of-phase solutions are particular cases.For the in-phase solution 1 − s 0 is a multiple of T in [see Eq. ( 31)] and therefore x 2 (s) = x 1 (s).For out-ofphase solutions, from (35) Substituting (54) in the set (13), cos[2x 2 (s − 1) We consider square-wave solutions in which x(s) remains practically constant over a plateau of duration s p ; for consistency we require qs p = s 0 , q being an integer.Then we can introduce and Y k−1 = Y (s − s p ) and choose x k located at the beginning of the plateau.Considering ε = 0, one has from ( 55) and ( 56) Subtracting (57) evaluated at k − 1 from the same equation evaluated at k one obtains We now make a few remarks about this map.First notice that to obtain this map we have not imposed that the solutions are periodic.Therefore, the map is, in principle, useful to determine the secondary instabilities of the primary periodic square waves.Second, for a given set of parameters, namely, for a fixed value of s 0 , one has several coexisting solutions with plateaus of length s p = s 0 /q.The amplitude of all these solutions is given by the map with the corresponding value of q.Third, for the primary in-phase square-wave periodic solutions discussed in the previous section s p = T in /2; as a consequence, from (42) it follows that s 0 = (2l + 1)s p , so q = 2l + 1. Similarly for the primary out-of-phase squarewave periodic solutions s p = T out /2 and from (51) it follows that s 0 = (2k + 1)s p , so q = 2k + 1. Fourth, for solutions whose period is twice the length of the plateau, as is the case of primary square waves, x k−2l−1 = x k−1 ; this is why the maps obtained in the previous section, where we explicitly considered the periodicity of the solution, correspond to q = 1.Fifth, the map for the primary square waves can in fact be further simplified considering that the square wave is centered at zero and defining Then one has Namely, the amplitude of all coexisting primary in-phase periodic solutions is the same and it is given by the fixed points of the sinus map (60) and similarly for all coexisting out-of-phase periodic solutions with different periods.
Figure 5(a) shows the bifurcation diagram of Eq. ( 59) for q = 1 as function of χ , where χ is changed by increasing P γ 2 while the offset phase is kept constant at = π/4.For small values of χ , the only stationary solution of the map is the fixed point x * = 0, which corresponds to the zero solution in (14).Increasing χ , the zero solution becomes unstable and at χ = 1, as expected from the linear stability analysis, one encounters for offset phase = 0.3π and 0.35π , respectively.Again for small χ the system goes to the zero solution, which becomes unstable at χ = 1.Above that value one encounters in-or out-of-phase periodic solutions whose amplitude x * is given by same sinus map (60).Then these solutions become unstable and periodic solutions of higher order appear.The point where the simple periodic solution becomes unstable depends on the offset phase.The higher-order periodic solutions that appear beyond this point are asymmetric, as discussed in the next section.Finally, the map becomes chaotic.

VI. NUMERICAL SIMULATIONS OF PERIODIC SOLUTIONS
We have performed numerical simulations of the dynamical model (11).We encounter that the zero solution becomes unstable at χ = 1 for any value of s 0 leading to an oscillatory solution.In Fig. 6 we show the in-phase oscillatory solution arising when increasing χ for s 0 = 1/3.Already for χ = 1.01 the shape of the oscillation resembles a square wave, although the plateaus are slightly tilted [see Figs.6(a) and 6(b)].As shown in [36] for a single OEO, the tilting of the plateau is an effect of δ.The square-wave form becomes even more clear as χ increases, as shown in Figs.6(c) and 6(d).Similarly, Fig. 7 displays the shape of the antiphase periodic solution obtained when considering s 0 = 3/4.In both cases the period of the oscillations coincides with the predicted one within order ε as expected.Comparing Fig. 6 with Fig. 7, one sees that for a given value of χ the amplitude of both the in-and out-of-phase oscillations is the same.This is in agreement with the fact that the amplitude of both in-and out-of-phase oscillations are all given by (60).We further illustrate this result in Fig. 8, where the solid line shows the amplitude given by the fixed point of the map (60) for P = 1.11 and γ = 0.5 as a function of the offset phase while the symbols show the numerical results obtained for several values of s 0 leading to different in-and out-of-phase solutions.As can be seen, the theoretical prediction for the oscillation amplitude is in excellent agreement with the numerical simulations of the full dynamical model.
As illustrated in the previous section, with increasing P γ 2 the in-and out-of-phase solutions become unstable, leading to higher-order periodic solutions.Figure 9 shows the bifurcation diagram of the map (59) for q = 1, γ = 0.5, and P = 1.8.There is a central region around = π/4 in which the system has periodic solutions such as the in-phase solution shown in Fig. 10 for = 0.18π .Further away from the center, the system has a period-doubling bifurcation in which it shows lagged synchronization [see Figs.10(a) and 10(b)].The values of the plateaus for these period-2 solutions are not symmetrically located around x = 0.This is also noticeable in the bifurcation diagram shown in Fig. 9, which is not symmetrical around the axis x = 0. Instead, it is symmetrical around the point x = 0, = π/4.The map (59) predicts quite accurately the amplitude of these lag synchronized solutions obtained numerically as shown by the red points.
We now address the coexistence of multiple periodic solutions for a fixed set of parameters.As shown in Fig. 3 and discussed in Sec.III, for a given value of s 0 , that is, for a given value for the delay times T f and T c , there are several Hopf bifurcations leading to oscillations with different period T , all of them taking place at χ = 1.Thus the linear stability analysis of the zero solution indicates that multiple periodic solutions with different periods can exist for χ = 1, although it does not give information about the stability of these solutions.It turns out that indeed there are several stable periodic square-wave solutions coexisting for a fixed set of parameters as illustrated in Fig. 11.The different square-wave solutions are obtained by integrating numerically the dynamical equations (11) starting from different initial conditions.More precisely, we take as an initial condition for x 1 (s) within the interval − max(1,s 0 ) < s < 0 a square wave with amplitude given by the fixed point x * of map (45) and with a period T .For in-phase solutions, the initial condition for x 2 is given by x 2 (s) = x 1 (s), while for out-of-phase solutions we take x 2 (s) = −x 1 (s).Regarding the initial condition for Y i , we take Y 1 (0) = x * T /4 and Y 2 (0) = x * T /4.
Figures 11(a)-11(d) show the coexistence of several inphase solutions obtained considering P = 1.5, γ = 0.5, = π/4, T f = 30 ns, and T f = 90 ns so that s 0 = 1/3 and χ = 1.5.For s 0 = 1/3 the fundamental solution corresponds to l f = 0 and m f = 1 as given by Eq. ( 28) and its period is given by T in f = 2/3 while the period of the higher harmonics is given by (37). Figure 11(a) shows x 1 (s), the fundamental solution.The time trace for x 2 (s) (not shown) coincides with the one for x 1 (s).Similarly, Figs.11(b) and 11(c) show x 1 (s) for the first and second harmonics, which have periods 1/3 and 1/5, respectively, as the fundamental period.Higher-order harmonics are also found.As an example, Fig. 11(d) shows the 20th harmonic, which has a period T in 20 = T in f /41 = 2/123 (notice that we have used a different scale on the time axis).
Figures 11(e)-11(h) illustrate the coexistence of several out-of-phase solutions obtained with the same parameters but T c = 60 ns so that s 0 = 1/2.The fundamental solution corresponds to k f = 0 and m f = 1 as given by Eq. ( 32) and has a period T out f = 1.The period of the higher harmonics is given by (39).As in the previous case, we only display the time traces for x 1 (s).The time traces for x 2 (s) are identical, but in opposite phase with respect to x 1 (s).All these solutions are stable against small numerical perturbations.For ε = 0 one could have an infinite number of such square-wave periodic orbits.In practice, in the full model, the transition between the plateaus of the square waves takes a time of order ε and this limits the minimal period for square-wave periodic solutions.Nevertheless, as shown in Fig. 11, it is perfectly feasible in this system to have tens of coexisting square-wave periodic solutions.The fundamental solution is the one with a larger basin of attraction, therefore starting from arbitrary initial conditions one usually ends up in the fundamental solution.However, by setting appropriately the initial condition, one can make the system operate in any of the other higher-order square-wave solutions.

VII. EFFECT OF MISMATCH IN THE DELAY TIMES
In order to analyze the effect of a small mismatch in the time delays we have performed numerical simulations of the dynamical system (9) for filter characteristic times θ = 5 μs and τ = 25 ps, low feedback and coupling rates (P = 1.5, γ = 0.5), offset phases = π/4 and 0 = 0, varying the self-feedback delay T f , and keeping fixed the coupling delay T c .As an initial condition we use an in-phase or antiphase periodic solution as described in the previous section.
The dynamics of periodic solutions in which x 1 and x 2 are in phase is shown in Fig. 12 for the fundamental square-wave solution [Figs.12(a), 12(c), 12(e), and 12(g)] and the first harmonic [Figs.12(b), 12(d), 12(f), and 12(h)].We consider T c = 90 ns, while T f is changed from T f = 30.0 to 29.2 ns.For T f = 30 ns one has a perfect matching condition for in-phase solutions, s 0 = 1/3.The fundamental solution has a period T in f = 2/3, while the period of the first harmonic is T in f /3 = 2/9.In this case one observes periodic square waves as discussed before [see Figs.12(a) and 12(b)].When a small mismatch is introduced the square-wave form is maintained, but the transition layer becomes progressively wider as illustrated in Figs.12(c) and 12(d).For a larger mismatch, of the order of 2% in T in f , the square-wave form starts to degrade and small intermediate plateaus appear.These small plateaus are not symmetrically located around x = 0.In the case illustrated in the figure the jump from the main plateau to the intermediate one at x > 0 is smaller than the equivalent jump at x < 0. We note that, despite the degradation in the shape of the square wave, the amplitude of the solution does not change.Higher-order harmonics are more sensitive to the mismatch since they have shorter periods and therefore will be more affected by widening and progressive degradation of the transition layer.Nevertheless, one can conclude that the square-wave solutions are robust to mismatch in the delay times of the order of a few percent and that several stable in-phase solutions can coexist even in presence of mismatch.Similar results are found for out-of-phase periodic solutions, as shown in Fig. 13.

VIII. CONCLUSION
We have studied the synchronization conditions in the delay times for in-and out-of-phase periodic square-wave solutions in a system with two delay-coupled optoelectronic delay loops.We have demonstrated that in-or out-of-phase synchronization is possible if the ratio between the self-and cross-delay times satisfies a rational relationship.In particular, the synchronization is in phase if the ratio involves two odd numbers, while it is out of phase for ratios involving an odd and an even number.We have derived analytical expressions for the period of the solutions and an approximated map for the amplitude of the square-wave oscillations.Remarkably, the map turns out to be the same for in-phase and out-of-phase solutions.We have also given an extension of this map that allows us to predict the secondary instabilities of the periodic square-wave solutions.The theoretical predictions are confirmed with numerical simulations of the full dynamical model.While the analytical mathematical treatment presented here has been done for identical systems in the ideal case without noise, we expect the results to be relevant for real systems.The effect of noise in the periodic regime of OEOs is small [31] and we have shown that the periodic square-wave solutions are robust to mismatches in the delay times of the order of a few percent.The rich dynamics of the system allows for the coexistence of many in-or out-of-phase stable periodic orbits with different periods for the same values of the fixed parameters.This provides a large degree of flexibility, allowing us to change the period of the square-wave periodic solution without changing any parameter of the system, just varying the initial condition or inducing a suitable perturbation to the system.Such a system turns out to be interesting for applications such as information encoding, which benefits from high-frequency oscillations of controllable period [20].We finally note that the methodology and the results presented here can be generalized to other systems consisting of coupled nonlinear oscillators with multiple different delays.

FIG. 1 .
FIG. 1. (Color online) Setup of two mutually coupled electrooptical delay systems with T f = T ii and T c = T ji,j =i .

FIG. 3 .
FIG. 3. (Color online) Hopf bifurcation points appearing at χ n = 1 as a function of s 0 leading to oscillatory solutions with different periods T .Squares and crosses correspond to in-and out-of-phase Hopf bifurcations, respectively.As T approaches zero, the number of Hopf bifurcations increases dramatically (only the bifurcations verifying T 0.1 are shown).

FIG. 9 .
FIG.9.(Color online) Amplitude of the square-wave oscillations as a function of the offset phase with P = 1.8 and γ = 0.5.The solid line corresponds to the theoretical prediction (59) with q = 1.Symbols correspond to the numerical integration of the full dynamical model with T f = 30 ns and T c = 90 ns.
Figure 11(e) shows the fundamental solution while Figs.11(f) and 11(g) show the first and second harmonics.Finally, Fig. 11(h) displays the twentieth harmonic.

FIG. 12 FIG. 13
FIG. 12. (Color online) Dynamics of in-phase periodic solutions with T c = 90 ns and different values of T f : (a) and (b) T f = 30.0ns, (c) and (d) T f = 29.6 ns, (e) and (f) T f = 29.4ns, and (g) and (h) T f = 29.2ns.