## Abstract

Square-wave oscillations exhibiting different plateau lengths have been observed experimentally by investigating an electro-optic oscillator. In a previous study, we analysed the model delay differential equations and determined an asymptotic approximation of the two plateaus. In this paper, we concentrate on the fast transition layers between plateaus and show how they contribute to the total period. We also investigate the bifurcation diagram of all possible stable solutions. We show that the square waves emerge from the first Hopf bifurcation of the basic steady state and that they may coexist with stable low-frequency periodic oscillations for the same value of the control parameter.

## 1. Introduction

Relaxation oscillations with alternate fast and slow phases appear in several areas of science from electronics to neural modelling. They are described mathematically as the solution of two or more nonlinear ordinary differential equations that exhibit different time scales. Over the years, reliable asymptotic techniques such as the method of matched asymptotic expansions [1–3] have been developed and successfully used to determine analytical expressions of physical interest (amplitude and period). The van der Pol equation in the large damping case is the reference problem for the analysis of relaxation oscillations [1–4] but other problems have also emerged in the field of chemical and biological oscillations [5–7].

Analytical studies of relaxation oscillations that are solutions of delay differential equations (DDEs) are however very rare. A notable exception is the analysis of a model for haematological stem cell regulation by Fowler & Mackey [8] and Fowler [9]. The method of matched asymptotic expansions is difficult to implement for DDEs because we often need to anticipate the response of both the state and delayed variables. Much of the mathematical work that has been done [10–14] is concerned with scalar nonlinear DDEs of the form
1.1where *x*^{′} denotes the derivative of *x* with respect to the dimensionless time *s* (*s*≡*t*/*t*_{D} where *t* is the real time and *t*_{D} is the delay of the feedback). *f*(*x*,*λ*) is a nonlinear function of *x* and *λ* is a control parameter. *ε*≡*t*_{0}/*t*_{D}>0, where *t*_{0} is the linear decay time of *x* in the absence of feedback. Equation (1.1) arises in a variety of applications, for example, physiological control systems [15], the transmission of light through a ring cavity [16–18] and population biology [19]. Under particular conditions on *f*(*x*,*λ*), equation (1.1) may exhibit nearly 2-periodic square-wave oscillations provided *ε* is sufficiently small (figure 1). More precisely, these oscillations consist of sharp transition layers of size proportional to *ε* connecting plateaus that are close to the period 2 fixed points of the map
1.2where *x*_{n}≡*x*(*t*) and *x*_{n−1}≡*x*(*t*−1). The period 2 fixed points of the map provide excellent approximations of the extrema of the oscillations. The description of the fast transition layers and the determination of the correction to the period are however much more delicate. Significant contributions to the asymptotic relations between the solutions of the map (1.2) and the solutions of the DDE (1.1) have been made by Chow & Mallet-Paret [10], Mallet-Paret & Nussbaum [11], Chow *et al.* [12] and Hale & Huang [13,14]. In particular, the Hopf bifurcation to the 2-periodic square-wave solutions has carefully been analysed. As the bifurcation parameter deviates from its Hopf bifurcation value, the oscillations quickly change their shape from sinusoidal to square waves [20].

Does equation (1.1) exhibit other types of square-wave oscillations? An analysis of the possible Hopf bifurcation points of equation (1.1) indicates that nearly 1-periodic square-wave solutions are possible but are unstable because they emerge from an unstable steady state [20]. Moreover, transient asymmetric square waves exhibiting different plateau lengths can be initiated by choosing particular initial conditions but they disappear at finite time [21].

In the study of Weicker *et al.* [22], we addressed the question whether stable periodic square-wave oscillations exhibiting different plateau lengths (called duty cycles) are possible for problems modelled by second-order DDEs. The question has been raised by experiments performed on electro-optic oscillators (EOOs), which are modelled mathematically in terms of second-order DDEs [23,24]. An EOO typically incorporates a nonlinear (intensity) modulator, an optical-fibre delay line, and an optical detector 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, as used in radar applications [25–27], or a broadband chaotic carrier typically intended for physical data encryption in high bit rate optical communications [28]. For a specific range of values of the parameters, periodic square-wave oscillations (figure 2) were found exhibiting a period *P* close to one delay *t*_{D} as well as different plateau lengths. It motivated an asymptotic analysis of the EOO equations in the limit of large delays. We obtained a good approximation of the plateaus and were able to explain how their respective lengths depend on the control parameters [22].

In this paper, we concentrate on three different issues that were omitted in [22]. First, We analyse the fast transition layers and show how they contribute to the total period. Second, we numerically investigate the bifurcation diagram of the square-wave oscillations and show how they emerge from a particular Hopf bifurcation. Third, we numerically found another stable time-periodic solution exhibiting a low frequency that may coexist with the square-wave solution. The plan of the paper is as follows. In §2, we introduce the EOO equations and propose a complete asymptotic description of the square-wave oscillations. In §2*a*, the approximation of the slowly varying plateaus is described in more detail than in [22]. The two fast transition layers are examined in §2*b*. We show that they are described by the same equation which we analyse. In §3, we numerically investigate the bifurcation diagram of the stable solutions using two different methods. The main results are summarized in §4.

## 2. Asymptotic analysis

In dimensional form, the evolution equation for an EOO are (eqns (35) and (36) in [23] or eqns (3) and (4) in [22])
2.1and
2.2where prime means differentiation with respect to *s* and *s* is time measured in units of the delay. The parameters *ε*, *δ* and *Φ* are fixed and given by [22]
2.3The feedback amplitude *β* is our bifurcation parameter.

Equations (2.1) and (2.2) admit nearly 1-periodic square-wave oscillations exhibiting different plateau lengths (figure 3). The slow–fast time behaviour of the solution is due to the small value of *ε*. As we shall later demonstrate, the relatively small change of *y* compared with *x* (figure 3*b*) is the result of the small value of *δ*. Furthermore, the asymmetry of the square-wave oscillations ( in figure 3*a*) is related to the deviation *Φ*+*π*/4. Experimentally, we may explore different ranges of values of these parameters. In this paper, we shall keep *δ* and *Φ* fixed as given in (2.3) and consider different small values of *ε* whenever it becomes appropriate for our numerical illustrations or analysis.

We next propose to construct the square-wave solution in the limit . Specifically, we seek a *P*-periodic solution satisfying the condition
2.4where the period *P* is given by
2.5As shown in figure 3*a*, the solution consists of two slowly varying plateaus connected by fast transition layers. We anticipate the analysis of the transition layers (see §2*b*) by assuming that the contribution from these layers to the period *P* is the same (*εr*). We analyse the slow and fast parts of the solution, separately.

### (a) Slowly varying plateaus

The leading approximation is obtained by setting *ε*=0 in equations (2.1) and (2.2). The reduced equations with (2.4) and (2.5) are
2.6
and
2.7
2.8
From equation (2.7), we determine *y*=*y*(*x*) as
2.9
The function (2.9) is represented in figure 3*b* and exhibits three branches provided *β*>1. The evolution of *x* and *y* along the left and right branches corresponds to the evolution along the plateaus of the square-wave periodic solution. They can be determined by inserting (2.9) into the left-hand side of equation (2.6) and by solving the resulting first-order equation for *x*. However, this solution is complicated and we may find simple analytical expressions by taking advantage of the small value of *δ*. Specifically, we seek a perturbation solution of equations (2.6) and (2.7) of the form
2.10
and
2.11
where *j*=1 or 2 refer to the time domains 0<*s*<*s*_{0} and *s*_{0}<*s*<1, respectively (figure 4).

Inserting (2.10) and (2.11) into equations (2.6) and (2.7) and equating to zero the coefficients of each power of *δ* leads to a sequence of problems for the unknowns functions *y*_{0}, *y*_{1j}, *x*_{0j} and *x*_{1j}. The leading-order problem is *O*(1) and is given by
2.12
and
2.13
Equation (2.12) implies that *y*_{0} is a constant. We already know that for a finite range of values of *y*_{0}, equation (2.13) admits more than one root (figure 3*b*). The solutions corresponding to the left and right branches are denoted by *x*_{01}<0 and *x*_{02}>0, respectively. We do not know the values of *y*_{0} and analyse the *O*(*δ*) problem for *y*_{1j}(*s*) and *x*_{1j}(*s*). It is given by
2.14
2.15
2.16
and
2.17
Figure 4 suggests the following initial conditions for *y*_{11} and *y*_{12}:
2.18
where *y*_{1M} and *y*_{1m} correspond to the maximum of *y*_{11} and the minimum of *y*_{12}, respectively. The solution of equations (2.14)–(2.18) is then
2.19
2.20
2.21
and
2.22
Continuity of *y*_{11} and *y*_{12} at times *s*=*s*_{0} and 1 leads to the conditions
2.23
and
2.24
which are two equations for *y*_{1M}−*y*_{1m}. A solution of equations (2.23) and (2.24) is possible only if
2.25
As for *y*_{11} and *y*_{12}, we next assume that the corrections *x*_{11} and *x*_{12} are equal at *s*=*s*_{0} and *s*=1 (figure 4*a*). From (2.21) and (2.22), we then obtain the condition
2.26
or equivalently,
2.27
Equation (2.27) admits multiple solutions. We specifically look for a solution of equation (2.27) which satisfies the perfect square-wave condition *x*_{01}=−*x*_{02} if *Φ*=−*π*/4. This solution is given by
2.28
Using (2.13), (2.28) allows one to determine *y*_{0}, *x*_{01} and *x*_{02}. Substracting equation (2.13) with *x*_{01} and equation (2.13) with *x*_{02} gives
2.29
Using (2.28) then allows one to eliminate *x*_{02} in equation (2.29). We find
2.30
Equation (2.30) provides the solution for *x*_{01}=*x*_{01}(*β*) in the implicit form
2.31
We obtain *x*_{02} and *s*_{0} by using (2.28) and (2.25):
2.32
and
2.33
In figure 5, we compare our approximations with the numerical solution obtained for *β*=1.2. The expression for *y*_{0} as well as for *x*_{03}, defined as the third root of equation (2.13), are documented in the appendix. In §3, we numerically analyse the bifurcation diagram of the possible stable solutions and show that the square-wave oscillations emerge from a Hopf bifurcation.

### (b) The fast transition layers

The plateaus of the square wave are connected by fast transition layers on time intervals proportional to *ε* (figure 3*a*).

#### (i) Jump down at *s*=0

We first consider the fast transition layer at *s*=0 and introduce the inner variable *ζ*_{1}≡*sε*^{−1}. The leading-order transition layer equations for *y*=*Y* _{1}(*ζ*_{1}) and *x*=*X*_{1}(*ζ*_{1}) are then given by
2.34
and
2.35
where we have used the periodicity condition
2.36
Equation (2.34) implies that *Y* _{1} is a constant. It needs to match the constant determined in our analysis of the slowly varying plateaus, i.e. *Y* _{1}=*y*_{0}*δ*^{−1}. Using the expression of *y*_{0} given by (A.2), equation (2.35) can be rewritten as
2.37
This equation can be reformulated in a simpler form by introducing the deviation *z*_{1}≡*X*_{1}−*x*_{03}=*X*_{1}+*Φ*+*π*/4. From equation (2.37), we obtain
2.38
The boundary conditions for the jump down transition are and . In terms of *z*_{1}, they take the simpler form
2.39
where
2.40

#### (ii) Jump up at *s*=*s*_{0}+*εr*

We next consider the transition layer near *s*=*s*_{0}+*εr* and introduce the inner variable *ζ*_{2}≡(*s*−*s*_{0}−*εr*)*ε*^{−1}. The leading-order transition layer equations for *y*=*Y* _{2}(*ζ*_{2}) and *x*=*X*_{2}(*ζ*_{2}) are given by
2.41
and
2.42
where we have used the periodicity condition
2.43
The constant solution for *Y* _{2} is again matching the value obtained from the analysis of slowly varying plateaus, i.e. *Y* _{2}=*y*_{0}*δ*^{−1}. Using the expression of *y*_{0} given by (A.2), equation (2.42) simplifies as
2.44
Introducing the deviation *z*_{2}≡*X*_{2}−*x*_{03}=*X*_{2}+*Φ*+*π*/4, equation (2.44) becomes
2.45
We next note the following relations between the two inner variables:
2.46
Inserting (2.46) into equation (2.45), we formulate an equation for *z*_{2}(*ζ*_{1}) of the form
2.47
The boundary conditions for the second transition layer are now
2.48
where *a* is defined by (2.40). We realize that equations (2.47) and (2.48) are the same as equations (2.38) and (2.39) except that the boundary conditions have been interchanged. This implies that the solution of equations (2.47) and (2.48) is related to the solution of equations (2.38) and (2.39) by
2.49

In conclusion, we found the same DDE for the two fast transition layers. It is given by
2.50
and
2.51
where we have omitted the subscript 1 for *z*_{1} and *ζ*_{1}. We next proceed as in [10]. We note that by rescaling time *ζ* as *ξ*≡−*ζ*/2*r*, equation (2.50) can be rewritten as a DDE with delay 1 and parameter *r*
2.52
and
2.53
*z*=±*a* are both critical points of equation (2.52). This means that we are looking for a heteroclinic orbit for some value of *r*, that is, a trajectory joining these critical points as . The delay parameter *r* is unknown *a priori*, and must be determined as part of the solution. We cannot solve the problem analytically for arbitrary *β* (it is a nonlinear DDE).

#### (iii) Correction to the period

In this section, we solve equations (2.52) and (2.53) for *β* close to 1. Our objective is to demonstrate that there is indeed a unique value of *r* such that equations (2.52) and (2.53) admit a solution. To this end, we introduce a small parameter *μ* defined by
2.54
where *b*=±1 if . We then expand the solution *z* and parameter *r* in power series of *μ*
2.55
and
2.56
where *ν*≡*μξ*. The motivation for introducing (2.54) comes from the fact that in first approximation as , which implies that the amplitude of the solution scales like . After introducing (2.54)–(2.56) into (2.50), we equate to zero the coefficients of each power of *μ*. The leading-order problem is *O*(*μ*) and is given by
2.57
In order to have a non-constant solution for *z*_{1}, we require that . The next problem is *O*(*μ*^{2}) and is given by
2.58
with the boundary conditions
2.59
We choose *b*=1 and note that the damped Hamiltonian equation (2.58) has a unique solution if *r*_{1}=0. We conclude that we have found an analytical expression for the transition layer solution provided
2.60
We have determined numerically the period of the square-wave oscillations with a high precision. The values of *δ* and *Φ* are documented in (2.3), *ε*=5×10^{−2} and *β*=1.2. We find *P*≃1.047. From (2.5), we then compute 2*εr*=4.7×10^{−2}, which implies *r*=0.47. The numerical value of *r* is close to the analytical value *r*=0.5 given in (2.60).

## 3. Numerical bifurcation diagrams

We consider *β* as our bifurcation parameter. All other parameters are documented in (2.3). A linear stability analysis of the steady state (*x*,*y*)=(0,0) allows us to determine the primary Hopf bifurcation points and Hopf frequencies. They satisfy the following equations:
3.1
and
3.2
The first Hopf bifurcation is located at *β*=*β*_{1}≃1.020 and exhibits a frequency close to 2*π* (*σ*_{1}=6.28). Using a continuation method, we find a 1-periodic branch of periodic solutions that connects the asymmetric square waves (figure 6*a*). More precisely, the Hopf bifurcation branch is first subcritical and unstable and then folds back to a branch of stable square-wave oscillations. There are many more Hopf bifurcation points as we further increase *β* from *β*_{1}. Using different initial conditions, we have integrated numerically equations (2.1) and (2.2), and found another branch of stable periodic solutions. By contrast to the square-wave oscillations, the new oscillations exhibit a low frequency. We have found that it emerges from the primary Hopf bifurcation point *β*=*β*_{2}≃1.025 as an unstable branch (as expected since this bifurcation is from an unstable steady state) and then stabilizes as *β*≥1.029. The frequency at the Hopf bifurcation point is *σ*_{2}=0.09 meaning a period *P*=69.81 (numerically, we found *P*=69.12; figure 7*b*).

The bifurcation diagrams shown in figure 6 illustrate the results of our simulations. The first Hopf bifurcation leads to the asymmetric square-wave oscillations that we investigated analytically in §2. Specifically, the extrema of *x* as a function of *β* are given by (2.31) and (2.33), where *x*_{01}≥−*π*/4−*Φ* (full red line in figure 6*a* and full black line in figure 6*b*). In figure 6*b*, the square dots are the solutions obtained numerically from simulating the full equations (2.1) and (2.2). For each point, the initial conditions were *x*=−1 , *x*=1 and *y*(0)=0. The long-time solution was then analysed when *s*>10 000. For *β*<1.009, the system jumps to the zero solution. The stability of the zero solution was also tested by using the initial conditions *x*=0 , *x*=10^{−3} and *y*(0)=0. The long-time solution was again analysed when *s*>10 000. For *β*=1.020, *x*=0 is stable. For *β*=1.021, *x*=0 is unstable and the system jumps to the 1-periodic asymmetric square wave. In addition to the 1-periodic square-wave solution, a stable low-frequency periodic solution was determined as soon as *β*≥1.029. The initial conditions were (−1≤*s*<0) and *y*(0)=0. At *β*=1.029, the frequency of the oscillations is *σ*=0.095 which is close to the Hopf frequency *σ*_{2}. The parabolic lines given by are curve fitting the numerical data and strongly suggest that the unstable branch of periodic solutions emerging at *β*=*β*_{2}=1.025 stabilizes as soon as *β*≥1.029. Similar responses (square-wave or low-frequency oscillations) have been found previously [23] but not for the same values of the bifurcation parameter. Here, the two distinct regimes may coexist (figure 7).

## 4. Discussion

In this paper, we investigated several issues that were missing in the study of Weicker *et al.* [22]. First, we concentrate on the fast transition layers between the plateaus of the square waves and showed how they contribute to the correction of the total period. Second, we show numerically that the square-wave oscillations are the result of a first Hopf bifurcation from the basic steady state. The bifurcation is subcritical and allows for the coexistence of stable square waves with a stable steady state. Experiments done on an EOO using quite different values of the parameters [24] suggest that the same mechanism could be responsible for the onset of asymmetric square waves. There are many other primary Hopf bifurcation points but we found only one leading to stable oscillations. The new periodic solution exhibits a large period and smooth oscillations. An asymptotic description of this solution is also possible [23]. Both the square-wave and the large period oscillations are the result of the large delay. They are dominant attractors in our EOO problem and motivate the investigation of other second-order nonlinear DDEs experiencing a large delay.

## Funding statement

T.E. and L.W. acknowledge the support of the FNRS (Belgium) and the Belgian FRIA for a PhD scholarship, respectively. This work benefited from the support of the Belgian Science Policy Office under grant no. IAP-7/35 ‘photonics@be’. It was also supported by the European project PHOCUS (FP7 grant no. 240763). L.L. thanks the support of the Institut Universitaire de France.

## Acknowledgements

We thank D. P. Rosin, D. J. Gauthier and E. Schöll for fruitful discussions at the conference in Palma. We also thank the referees for valuable questions.

## Appendix A

The plateaus of the square wave are *x*=*x*_{01}<0 and *x*=*x*_{02}>0, in first approximation. They are defined as two roots of equation (2.13) for a fixed *y*_{0}. Figure 3*b* suggests that there is a third root. In this appendix, we determine this third root and formulate an expression for *y*_{0}.

Equations for *x*_{01} and *x*_{02} are given by equations (2.28) and (2.30). From equation (2.30), we determine as
A1
From (2.13) with *j*=1, we formulate an expression for *y*_{0} given by
Using (A.1)
A2
In order to find the third root of equation (2.13), we introduce (A.2) into equation (2.13) and obtain
A3
This equation admits the solution
A4
Using equation (2.28), we then obtain the relation
A5
or equivalently,
A6
The two extreme roots are at equal distance from the central root *x*_{03}. This symmetry property has important consequences. In particular, the two fast transition layers admit the same equation and they contribute in the same way to the correction of the period.

## Footnotes

One contribution of 15 to a Theme Issue ‘Dynamics, control and information in delay-coupled systems’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.