Stability and Hopf Bifurcation Analysis of Coupled Optoelectronic Feedback Loops

Abstract

The dynamics of a coupled optoelectronic feedback loops are investigated. Depending on the coupling parameters and the feedback strength, the system exhibits synchronized asymptotically stable equilibrium and Hopf bifurcation. Employing the center manifold theorem and normal form method introduced by Hassard et al. (1981), we give an algorithm for determining the Hopf bifurcation properties.

1. Introduction

In recent research [1–5], it is found that even if several individual systems behave chaotically, in the case where the systems are identical, by proper coupling, the systems can be made to evolve toward a situation of exact isochronal synchronism. Synchronization phenomena are common in coupled semiconductor systems, and they are important examples of oscillators in general, and many works are concerned with coupled semiconductor systems [6–15].

We consider a feedback loop comprises a semiconductor laser that serves as the optical source, a Mach-Zehnder electrooptic modulator, a photoreceiver, an electronic filter, and an amplifier. The dynamics of the feedback loop can be modeled by the delay differential equations [14, 15]: d𝑥₁(𝑡)𝛾d𝑡=−₁+𝛾₂𝑥₁(𝑡)−𝛾₂𝑦₁(𝑡)−𝛽𝛾₂cos²𝑥₁(𝑡−𝜏)+𝜑₀,d𝑦₁(𝑡)d𝑡=𝛾₁𝑥₁(𝑡).(1) Here, 𝑥₁(𝑡) is the normalized voltage signal applied to the electrooptic modulator, 𝜏 is the feedback time delay, 𝛾₁ and 𝛾₂ are the filter low-pass and high-pass corner frequencies, 𝛽 is the dimensionless feedback strength, they are all positive constants, and 𝜑₀ is the bias point of the modulator.

Depending on the value of the feedback strength 𝛽 and delay 𝜏, the loop, which is modeled by system (1), is capable of producing dynamics ranging from periodic oscillations to high-dimensional chaos  [1, 14, 15].

We couple two nominally identical optoelectronic feedback loops unidirectionally, that is, the transmitter affects the dynamics of the receiver but not vice versa. Thus, the equations of motion describing the coupled system are given by (1) for the transmitter and d𝑥₂(𝑡)𝛾d𝑡=−₁+𝛾₂𝑥₂(𝑡)−𝛾₂𝑦₂(𝑡)−𝛽𝛾₂cos²𝑘𝑥₁(𝑡−𝜏)+(1−𝑘)𝑥₂(𝑡−𝜏)+𝜑₀,d𝑦₂(𝑡)d𝑡=𝛾₁𝑥₂(𝑡),(2) for the receiver. In (2), 𝑘>0 denotes the coupling strength. We will find that with the variety of 𝑘, the dynamical behavior of the coupled system can be different, while the feedback strength 𝛽 keeps the same value.

The paper is organized as follows. In Section 2, using the method presented in [16], we study the stability, and the local Hopf bifurcation of the equilibrium of the coupled system (1) and (2) by analyzing the distribution of the roots of the associated characteristic equation. In Section 3, we use the normal form method and the center manifold theory introduced by Hassard et al. [17] to analyze the direction, stability and the period of the bifurcating periodic solutions at critical values of 𝛽. In Section 4, some numerical simulations are carried out to illustrate the results obtained from the analysis. In Section 5, we come to some conclusion about the effect caused by the variety of parameters.

2. Stability Analysis

In this section, we consider the linear stability of the nonlinear coupled system d𝑥₁(𝑡)𝛾d𝑡=−₁+𝛾₂𝑥₁(𝑡)−𝛾₂𝑦₁(𝑡)−𝛽𝛾₂cos²𝑥₁(𝑡−𝜏)+𝜑₀,d𝑦₁(𝑡)d𝑡=𝛾₁𝑥₁(𝑡),d𝑥₂(𝑡)𝛾d𝑡=−₁+𝛾₂𝑥₂(𝑡)−𝛾₂𝑦₂(𝑡)−𝛽𝛾₂cos²𝑘𝑥₁(𝑡−𝜏)+(1−𝑘)𝑥₂(𝑡−𝜏)+𝜑₀,d𝑦₂(𝑡)d𝑡=𝛾₁𝑥₂(𝑡).(3)

It is easy to see that 𝐸(0,−𝛽cos²𝜑₀,0,−𝛽cos²𝜑₀) is the only equilibrium of system (3). Linearizing system (3) around 𝐸 and denote 𝛿=sin2𝜑₀, we get the linearization system d𝑥₁(𝑡)𝛾d𝑡=−₁+𝛾₂𝑥₁(𝑡)−𝛾₂𝑦₁(𝑡)+𝛽𝛿𝛾₂𝑥₁(𝑡−𝜏),d𝑦₁(𝑡)d𝑡=𝛾₁𝑥₁(𝑡),d𝑥₂(𝑡)𝛾d𝑡=−₁+𝛾₂𝑥₂(𝑡)−𝛾₂𝑦₂(𝑡)+𝑘𝛽𝛿𝛾₂𝑥₁(𝑡−𝜏)+(1−𝑘)𝛽𝛿𝛾₂𝑥₂(𝑡−𝜏),d𝑦₂(𝑡)d𝑡=𝛾₁𝑥₂(𝑡),(4) and the characteristic equation of system (4) 𝜆²+𝛾₁+𝛾₂𝜆+𝛾₁𝛾₂−(1−𝑘)𝛽𝛿𝛾₂𝜆𝑒^−𝜆𝜏×𝜆²+𝛾₁+𝛾₂𝜆+𝛾₁𝛾₂−𝛽𝛿𝛾₂𝜆𝑒^−𝜆𝜏=0,(5) which is equivalent to 𝜆²+𝛾₁+𝛾₂𝜆+𝛾₁𝛾₂−𝛽𝛿𝛾₂𝜆𝑒^−𝜆𝜏𝜆=0,(6)₂+𝛾₁+𝛾₂𝜆+𝛾₁𝛾₂−(1−𝑘)𝛽𝛿𝛾₂𝜆𝑒^−𝜆𝜏=0.(7)

Notice that when 𝛽=0, (5) becomes 𝜆²+𝛾₁+𝛾₂𝜆+𝛾₁𝛾₂²=0,(8) whose roots are 𝜆_1,2=−𝛾₁,𝜆_3,4=−𝛾₂.(9) So, we have the following lemma.

Lemma 1. The equilibrium 𝐸₀(0,−𝛽cos²𝜑₀,0,−𝛽cos²𝜑₀) is asymptotically stable when 𝛽=0.

Next, we regard 𝛽 as the bifurcation parameter to investigate the distribution of roots of (6) and (7).

Let 𝜆=i𝜔  (𝜔>0) be a root of (6) and substituting 𝜆=i𝜔 into (6), separating the real and imaginary parts yields −𝜔²+𝛾₁𝛾₂=𝛽𝛿𝛾₂𝜔𝛾𝜔sin𝜔𝜏,₁+𝛾₂=𝛽𝛿𝛾₂𝜔cos𝜔𝜏.(10) Then, we can get tan𝜔𝜏=−𝜔²+𝛾₁𝛾₂𝜔𝛾₁+𝛾₂.(11) Hence, (11) has a sequence of roots {𝜔_𝑗}_𝑗≥0 (see Figure 1), and 𝜔_𝑗∈2𝑗𝜋𝜏,2𝑗𝜋+𝜋/2𝜏,𝜔²_𝑗<𝛾₁𝛾₂,2𝑗𝜋+3𝜋/2𝜏,2𝑗+1𝜋𝜏,𝜔²_𝑗>𝛾₁𝛾₂,𝑗=0,1,2,….(12)

Figure 1 

The points of intersection of 𝑓1=tan𝜔𝜏 and 𝑓2=(−𝜔2+𝛾1𝛾2)/𝜔(𝛾1+𝛾2), when 𝛾1=0.1, 𝛾2=2.5,  𝜏=1.5.

Define 𝛽_𝑗=𝛾₁+𝛾₂𝛿𝛾₂cos𝜔_𝑗𝜏.(13) Then, (𝜔_𝑗,𝛽_𝑗) is the solution of (10).

From (10), we know that 𝛽²=1𝛿²𝛾²₂𝛾₁𝛾₂𝜔−𝜔²+𝛾₁+𝛾₂²,(14) which gives that d𝛽=1d𝜔𝛽𝛿²𝛾²₂𝜔³𝜔²+𝛾₁𝛾₂𝜔²−𝛾₁𝛾₂.(15) From Figure 1, we know that 𝜔_𝑗→∞ when 𝑗→∞, which means that 𝜔²_𝑗>𝛾₁𝛾₂; furthermore, 𝛽 is increasing with respect to 𝜔, when 𝑗 is sufficiently big.

Reorder the set {𝛽_𝑗} such that 𝛽₀=min{𝛽_𝑗} and 𝜔_𝑗 is correspondent of 𝛽_𝑗(𝑗=0,1,2,…). Then, we have the following lemma.

Lemma 2. There exists a sequence values of 𝛽 denoted by 0<𝛽₀<𝛽₁<⋯,(16) such that (6) has a pair of imaginary roots ±𝑖𝜔_𝑗 when 𝛽=𝛽_𝑗(𝑗=0,1,2,…), where 𝛽_𝑗 is defined by (13), and 𝜔_𝑗 is the root of (11).

Let 𝜆(𝜏)=𝛼𝛽+𝑖𝜔𝛽(17) be the root of (6) satisfying 𝛼(𝛽_𝑗)=0 and 𝜔(𝛽_𝑗)=𝜔_𝑗. We have the following conclusion.

Lemma 3. 𝛼^′(𝛽_𝑗)>0.

Proof. Substituting 𝜆(𝛽) into (6) and taking the derivative with respect to 𝛽, it follows that 2𝜆d𝜆+𝛾d𝛽¹+𝛾₂d𝜆d𝛽−𝛿𝛾₂𝜆𝑒^−𝜆𝜏−𝛽𝛿𝛾₂𝑒^−𝜆𝜏d𝜆d𝛽+𝜏𝛽𝛿𝛾₂𝜆𝑒^−𝜆𝜏d𝜆d𝛽=0.(18) Therefore, noting that 𝛽𝛿𝛾₂𝜆𝑒^−𝜆𝜏=𝜆²+(𝛾₁+𝛾₂)𝜆+𝛾₁𝛾₂, we have d𝜆=1d𝛽𝛽𝜆³+𝛾₁+𝛾₂𝜆²+𝛾₁𝛾₂𝜆𝜆²−𝛾₁𝛾₂+𝜏𝜆³+𝛾₁+𝛾₂𝜆²+𝛾₁𝛾₂𝜆,(19) and by a straight computation, we get 𝛼^′𝛽_𝑗=𝜔²_𝑗𝛽_𝑗Δ𝜏𝜔²_𝑗𝛾²₁+𝛾²₂+𝜔²_𝑗+𝛾₁𝛾₂𝛾₁+𝛾₂+𝜏𝜔⁴_𝑗+𝜏𝛾²₁𝛾²₂>0,(20) where Δ=𝜔²_𝑗+𝛾₁𝛾₂+𝜏𝜔²_𝑗𝛾₁+𝛾₂²+𝜏𝜔_𝑗𝜔²_𝑗−𝛾₁𝛾₂².(21)

As to (7), it can be easily found that −𝛾₁,−𝛾₂ are two negative roots when 𝑘=1, so, next, we only focus on (7) with 𝑘≠1.

Let 𝜆=𝑖𝜛(𝜛)>0 be a root of (7). Using the same method above, we get −𝜛²+𝛾₁𝛾₂=(1−𝑘)𝛽𝛿𝛾₂𝛾𝜛sin𝜛𝜏,₁+𝛾₂𝜛=(1−𝑘)𝛽𝛿𝛾₂𝜛cos𝜛𝜏,(22)tan𝜛𝜏=−𝜛²+𝛾₁𝛾₂𝜛𝛾₁+𝛾₂.(23)

Then, when 0<𝑘<1, (23) has a sequence of roots {𝜛_𝑗}_𝑗≥0, which are the same as those of (11).

When 𝑘>1, (23) has a sequence of roots {𝜛_𝑗}_𝑗≥0, and 𝜛_𝑗∈2𝑗+1𝜋𝜏,2𝑗+1𝜋+𝜋/2𝜏,𝜛²_𝑗<𝛾₁𝛾₂,2𝑗𝜋+𝜋/2𝜏,2𝑗+1𝜋𝜏,𝜛²_𝑗>𝛾₁𝛾₂,𝑗=0,1,2,….(24)

Define 𝛽_𝑗=𝛾₁+𝛾₂(1−𝑘)𝛿𝛾₂.cos𝜛𝜏(25) Then, (𝜛_𝑗,𝛽_𝑗) is the solution of (22).

Repeat the previous process, we have d𝛽=1d𝜛(1−𝑘)²𝛽𝛿²𝛾²₂𝜛³𝜛²+𝛾₁𝛾₂𝜛²−𝛾₁𝛾₂.(26)

Reorder the set {𝛽_𝑗} such that 𝛽₀=min{𝛽_𝑗} and 𝜛_𝑗 is correspondent of 𝛽_𝑗(𝑗=0,1,2,…).

Lemma 4. There exists a sequence values of 𝛽 denoted by 0<𝛽₀<𝛽₁<⋯,(27) such that (7) has a pair of imaginary roots ±𝑖𝜛_𝑗 when 𝛽=𝛽_𝑗(𝑗=0,1,2,…), where 𝛽_𝑗 is defined by (25), and 𝜛_𝑗 is the root of (23).

Let 𝜆(𝜏)=𝛼𝛽+i𝜛𝛽(28) be the root of (7) satisfying 𝛼(𝛽_𝑗)=0, 𝜛(𝛽_𝑗)=𝜛_𝑗. Then, similar to the proof of Lemma 3, we have the following conclusion.

Lemma 5. 𝛼^′(𝛽_𝑗)>0.

Compare 𝛽_𝑗,  𝛽_𝑗 and reorder the set {𝛽_𝑗} and {𝛽_𝑗} and remove the “−” of 𝛽_𝑗, such that 0<𝛽₀<𝛽₁<⋯,(29) then from previous lemmas and the Hopf bifurcation theorem for functional differential equations [18], we have the following results on stability and bifurcation to system (3).

Theorem 6. For system(3), the equilibrium 𝐸 is asymptotically stable when 𝛽∈[0,𝛽₀) and unstable when 𝛽∈(𝛽₀,+∞); system (3) undergoes a Hopf bifurcation at 𝐸 when 𝛽=𝛽_𝑗, 𝑗=0,1,2,…, where 𝛽_𝑗 are defined by (13) or (25).

3. The Direction and Stability of the Hopf Bifurcation

In Section 2 we obtained some conditions under which system (3) undergoes the Hopf bifurcation at some critical values of 𝛽. In this section, we study the direction, stability, and the period of the bifurcating periodic solutions. The method we used is based on the normal form method and the center manifold theory introduced by Hassard et al. [17].

Move 𝐸(0,−𝛽cos²𝜑₀,0,−𝛽cos²𝜑₀) to the origin 𝑂(0,0,0,0) and denote 𝛿=sin2𝜑₀,  𝜌=cos2𝜑₀, then system (3) can be written as the form d𝑥₁(𝑡)𝛾d𝑡=−₁+𝛾₂𝑥₁(𝑡)−𝛾₂𝑦₁(𝑡)+𝛽𝛿𝛾₂𝑥₁(𝑡−𝜏)+𝛽𝛾₂𝜌𝜑₀𝑥²₁()−2𝑡−𝜏3𝛽𝛿𝛾₂𝑥³₁(𝑡−𝜏)+𝑂(4),d𝑦₁(𝑡)d𝑡=𝛾₁𝑥₁(𝑡),d𝑥₂(𝑡)𝛾d𝑡=−₁+𝛾₂𝑥₂(𝑡)−𝛾₂𝑦₂(𝑡)+𝛽𝛾₂𝑘𝛿𝑥₁(𝑡−𝜏)+(1−𝑘)𝛿𝑥₂(𝑡−𝜏)+𝑘²𝜌𝑥²₁()()𝑡−𝜏+2𝑘1−𝑘𝜌𝑥₁(𝑡−𝜏)𝑥₂()+()𝑡−𝜏1−𝑘²𝜌𝑥²₂()−2𝑡−𝜏3𝑘³𝛿𝑥³₁(𝑡−𝜏)−2𝑘²(1−𝑘)𝛿𝑥²₁(𝑡−𝜏)𝑥₂()()𝑡−𝜏−2𝑘1−𝑘²𝛿𝑥₁(𝑡−𝜏)𝑥²₂()−2𝑡−𝜏3(1−𝑘)³𝛿𝑥³₂(𝑡−𝜏)+𝑂(4),d𝑦₂(𝑡)d𝑡=𝛾₁𝑥₂(𝑡).(30)

Clearly, the phase space is 𝒞=𝒞([−𝜏,0],ℝ⁴). For convenience, let 𝛽^∗∈𝛽_𝑗∪𝛽_𝑗,(31) and 𝛽=𝛽^∗+𝜇, 𝜇∈ℝ. From the analysis above we know that 𝜇=0 is the Hopf bifurcation value for system(30). Let i𝜔^∗ be the root of the characteristic equation associate with the linearization of system (30) when 𝛽=𝛽^∗. For 𝜙=(𝜙₁,𝜙₂,𝜙₃,𝜙₄)∈𝒞, let 𝐿_𝜇𝜙=𝐵𝜙(0)+𝐶𝜙(−𝜏),(32) where −𝛾𝐵=¹+𝛾₂−𝛾₂𝛾00₁𝛾00000−₁+𝛾₂−𝛾₂00𝛾₁0,𝐶=𝛽𝛿𝛾₂0000000𝑘𝛽𝛿𝛾₂0(1−𝑘)𝛽𝛿𝛾₂0.0000(33) By the Rieze representation theorem, there exists a 4×4 matrix, 𝜂(𝜃,𝜇)(−𝜏≤𝜃≤0), whose elements are of bounded variation functions such that 𝐿_𝜇𝜙=⁰_−𝜏d𝜂𝜃,𝜇𝜙(𝜃),𝜙∈𝒞.(34)

In fact, we can choose 𝜂𝜃,𝜇=()𝐵,𝜃=0,0,𝜃∈−𝜏,0−𝐶,𝜃=−𝜏.(35) Then, (30) is satisfied.

For 𝜙∈𝒞, define the operator 𝐴(𝜇) as 𝐴𝜇𝜙(𝜃)=d𝜙(𝜃)[),∫d𝜃,𝜃∈−𝜏,0_0−𝜏d𝜂𝑡,𝜇𝜙(𝑡),𝜃=0,(36) and 𝑅(𝜇)𝜙 as 𝑅𝜇𝜙(𝜃)=[),𝑓0,𝜃∈−𝜏,0𝜇,𝜙,𝜃=0,(37) where 𝑓𝜇,𝜙=𝛽^∗𝛾₂𝜌𝜙²₁()−2−𝜏3𝛿𝜙³₁()(4)0𝑘−𝜏+𝑂₂𝜌𝜙²₁(−𝜏)+𝑘(1−𝑘)𝜌𝜙₁(−𝜏)𝜙₃()+()−𝜏1−𝑘²𝜌𝜙²₃()−2−𝜏3𝑘³𝛿𝜙³₁(−𝜏)−2𝑘²(1−𝑘)𝛿𝜙²₁(−𝜏)𝜙₃()()𝑡−𝜏−2𝑘1−𝑘²𝛿𝜙₁(−𝜏)𝜙²₃()−2−𝜏3(1−𝑘)³𝛿𝜙³₃()(4)0−𝜏+𝑂.(38)Then, system (30) is equivalent to the following operator equation: 𝑢̇_𝑡=𝐴𝜇𝑢_𝑡+𝑅𝜇𝑢_𝑡,(39) where 𝑢(𝑡)=(𝑥₁(𝑡),𝑦₁(𝑡),𝑥₂(𝑡),𝑦₂(𝑡))^𝑇,  𝑢_𝑡=𝑢(𝑡+𝜃), for 𝜃∈[−𝜏,0].

For 𝜓∈𝒞¹([0,𝜏],ℝ⁴), define 𝐴^∗𝜓(𝑠)=−d𝜓(𝑠)(],∫d𝑠,𝑠∈0,𝜏_0−𝜏𝜓(−𝜉)d𝜂(𝜉,0),𝑠=0.(40)

For 𝜙∈𝒞[−𝜏,0] and 𝜓∈𝒞[0,𝜏], define the bilinear form 𝜓(𝑠),𝜙(𝜃)=𝜓(0)𝜙(0)−⁰_−𝜏^𝜃₀𝜓(𝜉−𝜃)d𝜂(𝜃)𝜙(𝜉)d𝜉,(41) where 𝜂(𝜃)=𝜂(𝜃,0). Then, 𝐴(0) and 𝐴^∗ are adjoint operators.

Let 𝑞(𝜃) and  𝑞^∗(𝑠) be eigenvectors of 𝐴(0) and 𝐴^∗ associated to 𝑖𝜔^∗ and −𝑖𝜔^∗, respectively. It is not difficult with verify that 𝑞(𝜃)=𝛾1,¹𝑖𝜔^∗𝛾,1,¹𝑖𝜔^∗𝑇𝑒^{𝑖𝜔∗𝜃},𝑞^∗(𝑠)=1𝐷𝛾1,²𝑖𝜔^∗𝛾,1,²𝑖𝜔^∗𝑒^{𝑖𝜔∗𝑠},(42) where 𝐷=2+2𝛾₁𝛾₂𝜔^∗2+2𝛽^∗𝛿𝛾₂𝜏𝑒^{−𝑖𝜔∗𝜏}.(43) Then, ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=1,  ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=0.

Let 𝑢_𝑡 be the solution of (39) and define 𝑧(𝑡)=𝑞^∗,𝑢_𝑡,𝑊(𝑡,𝜃)=𝑢_𝑡(𝜃)−2Re𝑧(𝑡)𝑞(𝜃).(44) On the center manifold 𝒞₀, we have 𝑊(𝑡,𝜃)=𝑊(𝑧(𝑡),𝑧(𝑡),𝜃),(45) where 𝑊(𝑧,𝑧,𝜃)=𝑊₂₀𝑧²2+𝑊₁₁𝑧𝑧+𝑊₀₂𝑧²2+⋯,(46)𝑧 and 𝑧 are local coordinates for center manifold 𝒞₀ in the direction of 𝑞^∗ and 𝑞^∗. Note that 𝑊 is real if 𝑢_𝑡 is real. We only consider real solutions.

For solution 𝑢_𝑡 in 𝒞₀, since 𝜇=0, we have 𝑧̇(𝑡)=𝑖𝜔^∗𝑧+𝑞^∗(𝜃),𝑓0,𝑊+2Re𝑧(𝑡)𝑞(𝜃)=𝑖𝜔^∗𝑧+𝑞^∗(0),𝑓0,𝑊(𝑧,𝑧,0)+2Re𝑧(𝑡)𝑞(0)=𝑖𝜔^∗𝑧+𝑞^∗(0)𝑓₀(𝑧,𝑧).(47) We rewrite this equation as 𝑧̇(𝑡)=𝑖𝜔^∗𝑧+𝑔(𝑧,𝑧),(48) where 𝑔(𝑧,𝑧)=𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+𝑔₂₁𝑧²𝑧2⋯.(49) By (39) and (48), we have 𝑊̇=𝑢̇_𝑡−𝑧𝑞̇−𝑧̇𝑞=𝐴𝑊−2Re𝑞^∗(0)𝑓₀𝑞(𝜃),𝜃∈[−𝜏,0),𝐴𝑊−2Re𝑞^∗(0)𝑓₀𝑞(0)+𝑓₀(,𝜃=0,=𝐴𝑊+𝐻𝑧,𝑧,𝜃),(50) where 𝐻(𝑧,𝑧,𝜃)=𝐻₂₀(𝜃)𝑧²2+𝐻₁₁(𝜃)𝑧𝑧+𝐻₀₂(𝜃)𝑧²2+⋯.(51) Expanding the above series and comparing the coefficients, we obtain 𝐴−2i𝜔^∗𝐼𝑊₂₀(𝜃)=−𝐻₂₀(𝜃),𝐴𝑊₁₁(𝜃)=−𝐻₁₁(𝜃),….(52) Notice that

𝑞(𝜃)=𝛾1,¹𝑖𝜔^∗𝛾,1,¹𝑖𝜔^∗𝑇𝑒^{𝑖𝜔∗𝜃},𝑢_𝑡(𝜃)=𝑧𝑞(𝜃)+𝑧𝑞(𝜃)+𝑊(𝑧,𝑧,𝜃),(53) where 𝑊^(𝑖)(𝑧,𝑧,𝜃)=𝑊^(𝑖)₂₀(𝜃)𝑧²2+𝑊^(𝑖)₁₁(𝜃)𝑧𝑧+𝑊^(𝑖)₀₂(𝜃)𝑧²2+⋯,𝑖=1,2,3,4.(54) Combing (38) and by straightforward computation, we can obtain the coefficients which will be used in determining the important quantities: 𝑔₂₀=2𝛽^∗𝛾₂𝜌𝐷𝑒^{−2𝑖𝜔∗𝜏}𝑘²,𝑔−𝑘+2₁₁=2𝛽^∗𝛾₂𝜌𝐷𝑘²,𝑔−𝑘+2₀₂=2𝛽^∗𝛾₂𝜌𝐷𝑒^{2𝑖𝜔∗𝜏}𝑘²,𝑔−𝑘+2₂₁=2𝛽^∗𝛾₂𝐷𝜌𝑒^{𝑖𝜔∗𝜏}𝑊⁽¹⁾₂₀(−𝜏)+2𝑒^{−𝑖𝜔∗𝜏}𝑊⁽¹⁾₁₁(−𝜏)+𝑘²𝜌𝑒^{𝑖𝜔∗𝜏}𝑊⁽¹⁾₂₀(−𝜏)+2𝑒^{−𝑖𝜔∗𝜏}𝑊⁽¹⁾₁₁()()𝜌×𝑒−𝜏+𝑘1−𝑘_{−𝑖𝜔∗𝜏}𝑊⁽³⁾₁₁(−𝜏)+𝑒^{𝑖𝜔∗𝜏}𝑊⁽³⁾₂₀(−𝜏)2+𝑒^{𝑖𝜔∗𝜏}𝑊⁽¹⁾₂₀(−𝜏)2+𝑒^{−𝑖𝜔∗𝜏}𝑊⁽¹⁾₁₁()+()−𝜏1−𝑘²𝜌2𝑒^{−𝑖𝜔∗𝜏}𝑊⁽³⁾₁₁(−𝜏)+𝑒^{𝑖𝜔∗𝜏}𝑊⁽³⁾₂₀(−𝜏)−4𝛿𝑒^{−𝑖𝜔∗𝜏}.(55) We still need to compute 𝑊₂₀(𝜃) and 𝑊₁₁(𝜃), for 𝜃∈[−𝜏,0). We have 𝐻(𝑧,𝑧,𝜃)=−𝑞^∗(0)𝑓₀𝑞(𝜃)−𝑞^∗(0)𝑓₀𝑞(𝜃)=−𝑔(𝑧,𝑧)𝑞(𝜃)−𝑔(𝑧,𝑧)𝑞(𝜃).(56) Comparing the coefficients about 𝐻(𝑧,𝑧,𝜃) gives that 𝐻₂₀(𝜃)=−𝑔₂₀𝑞(𝜃)−𝑔₀₂𝑞(𝜃),𝐻₁₁=−𝑔₁₁𝑞(𝜃)−𝑔₁₁𝑞(𝜃).(57) Then, from (52), we get 𝑊̇₂₀(𝜃)=2𝑖𝜔^∗𝑊₂₀(𝜃)+𝑔₂₀𝑞(𝜃)+𝑔₀₂𝑞(𝜃),𝑊̇₁₁(𝜃)=𝑔₁₁𝑞(𝜃)+𝑔₁₁𝑞(𝜃),(58) which implies that 𝑊₂₀(𝜃)=𝑔₂₀𝑞(0)−𝑖𝜔^∗𝑒^{𝑖𝜔∗𝜃}+𝑔₀₂𝑞(0)−3𝑖𝜔^∗𝑒^{−𝑖𝜔∗𝑗𝜃}+𝐸𝑒^{2𝑖𝜔∗𝑗𝜃},𝑊₁₁(𝜃)=𝑔₁₁𝑞(0)𝑖𝜔^∗𝑒^i𝜔∗𝜃+𝑔₁₁𝑞(0)−𝑖𝜔^∗𝑒^{−i𝜔∗𝜃}+𝐹.(59) Here, 𝐸 and 𝐹 are both four-dimensional vectors and can be determined by setting 𝜃=0 in 𝐻(𝑧,𝑧,𝜃). In fact, from (38) and 𝐻(𝑧,𝑧,0)=−2Re𝑞^∗(0)𝑓₀𝑞(0)+𝑓₀,(60) we have 𝐻₂₀(0)=−𝑔₂₀𝑞(0)−𝑔₀₂𝑞(0)+2𝛽^∗𝛾₂𝜌𝑒^{−2𝑖𝜔∗𝜏}1,0,𝑘²−𝑘+1,0^𝑇,𝐻₁₁(0)=−𝑔₁₁𝑞(0)−𝑔₁₁𝑞(0)+2𝛽^∗𝛾₂𝜌1,0,𝑘²−𝑘+1,0^𝑇.(61) It follows from (52) and the definition of 𝐴 that 𝛽^∗𝐵𝑊₂₀(0)+𝛽^∗𝐶𝑊₂₀(−𝜏)=2𝑖𝜔^∗𝑊₂₀(0)−𝐻₂₀(0),𝛽^∗𝐵𝑊₁₁(0)+𝛽^∗𝐶𝑊₁₁(−𝜏)=−𝐻₁₁(0),(62) which implies that 𝐸=𝐵+𝑒^{−2𝑖𝜔∗𝜏}𝐶−2𝑖𝜔^∗I⁻¹×𝐵𝑔₂₀𝑞(0)𝑖𝜔^∗+𝑔₀₂𝑞(0)3𝑖𝜔^∗𝑔+𝐶²⁰𝑞(0)𝑖𝜔^∗𝑒^{−𝑖𝜔∗𝜏}+𝑔₀₂𝑞(0)3𝑖𝜔^∗𝑒^{𝑖𝜔∗𝜏}+1𝛽^∗𝑔₂₀𝑞(0)+𝑔₀₂𝑞(0)×2𝛽^∗𝛾₂𝜌𝑒^{−2𝑖𝜔∗𝜏}1,0,𝑘²−𝑘+1,0^𝑇,𝐹=(𝐵+𝐶)⁻¹𝐵𝑔₁₁𝑞(0)−𝑖𝜔^∗+𝑔₁₁𝑞(0)𝑖𝜔^∗𝑔+𝐶¹¹𝑞(0)−𝑖𝜔^∗𝑒^{−𝑖𝜔∗𝜏}+𝑔₁₁𝑞(0)𝑖𝜔^∗𝑒^{𝑖𝜔∗𝜏}+1𝛽^∗𝑔₁₁𝑞(0)+𝑔₁₁𝑞(0)−2𝛽^∗𝛾₂𝜌1,0,𝑘²−𝑘+1,0^𝑇.(63) Consequently, the above 𝑔₂₁ can be expressed by the parameters and delay in system (30). Thus, we can compute the following quantities: 𝑐₁(0)=𝑖2𝜔^∗𝑔₂₀𝑔₁₁−2𝑔₁₁²−13𝑔₂₀²+𝑔₂₁2,𝜇₂=−Re𝑐₁(0)Re𝜆^′𝛽^∗,𝛽₂=2Re𝑐₁(0),𝑇₂=−Im𝑐₁(0)+𝜇₂Im𝜆^′𝛽^∗𝜔^∗,(64) which determine the properties of bifurcating periodic solutions at the critical value 𝜏₀. The direction and stability of the Hopf bifurcation in the center manifold can be determined by 𝜇₂ and 𝛽₂, respectively. In fact, if 𝜇₂>0(𝜇₂<0), then the bifurcating periodic solutions are forward (backward); the bifurcating periodic solutions on the center manifold are stable (unstable) if 𝛽₂<0(𝛽₂>0); and 𝑇₂ determines the period of the bifurcating periodic solutions: the period increases (decreases) if 𝑇₂>0(𝑇₂<0).

From the discussion in Section 2, we have known that Re𝜆^′(𝛽_𝑗)>0; therefore; we have the following result.

Theorem 7. The direction of the Hopf bifurcation for system (3) at the equilibrium 𝐸(0,−𝛽cos²𝜑₀,0,−𝛽cos²𝜑₀) when 𝛽=𝛽^∗ is forward (backward), and the bifurcating periodic solutions on the center manifold are stable (unstable) if Re(𝑐₁(0))<0(>0). Particularly, the stability of the bifurcation periodic solutions of system (3) and the reduced equations on the center manifold are coincident at the first bifurcation value 𝛽=𝛽₀.

4. Numerical Simulations

In this section, we will carry out numerical simulations on system (3) at special values of 𝛽. We choose a set of data as follows: 𝛾₁=0.1,𝛾₂=2.5,𝜑₀=𝜋4,𝜏=1.5,(65) which are the same as those in [1]. Then, 𝛿=1,𝜌=0.

Then, we can obtain 𝜔=₀^⋅0.2225,𝜔₁^⋅𝜛=3.5677,…,_0⋅=1.7294,𝜛₁^⋅𝛽=5.5531,…,_0⋅=1.1008,𝛽=₁^⋅1.7434,…,𝛽=₀^⋅1.3534,𝛽=₁^⋅2.5235,…,𝑘=1.9,𝛽=₀^⋅0.6093,𝛽=₁^⋅1.1356,…,𝑘=3.(66) From the analysis in Section 2, we know that 𝛽(𝛽) is increasing with respect to 𝜔(𝜛) when 𝜔(𝜛)>𝛾₁𝛾₂, which means that 𝛽₀=min𝛽_𝑗,𝛽₀=min𝛽_𝑗,𝑗=0,1,2,…,(67) that is, 𝛽₀(𝛽₀) is the first critical value at which system (3) undergoes a Hopf bifurcation.

When 𝑘=1.9, by the previous results, it follows that 𝜆^′𝛽=₀^⋅0.2440−0.0491i,𝑐₁(0)^⋅𝜇=−0.5373+0.1082i,_2⋅=2.2020,𝛽=₂^⋅−1.0746,𝑇₂=−0.0018.(68)

Hence, we arrive at the following conclusion: the equilibrium 𝐸 is asymptotically stable when 𝛽∈[0,1.1008) and unstable when 𝛽∈(1.1008,+∞), and, at the first critical value, the bifurcating periodic solutions are asymptotically stable, and the direction of the bifurcation is forward (see Figures 2 and 3).

(a)

(b)

(c)

(d)

Figure 2 

𝛾1=0.1, 𝛾2=2.5, 𝜏=1.5, 𝑘=1.9, which means that condition (𝐻1) holds, and 𝛽=0.7<𝛽0. The initial value is (0.1,−0.5,0.1,−0.5).

(a)

(b)

(c)

(d)

Figure 3 

𝛾1=0.1, 𝛾2=2.5, 𝜏=1.5, 𝑘=1.9, which means that condition (𝐻1) holds, and 𝛽=1.2>𝛽0. The initial value is (0.1,−0.5,0.1,−0.5).

When 𝑘=3, we can get 𝜆^′𝛽=₀^⋅0.9004+0.0924i,𝑐₁(0)^⋅𝜇=−1.9116+0.7930i,_2⋅=2.1231,𝛽=₂^⋅−3.8232,𝑇₂=−0.5720.(69)

Then, we have the following: the equilibrium 𝐸 is asymptotically stable when 𝛽∈[0,0.6093), and unstable when 𝛽∈(0.6093,+∞), and, at the first critical value, the bifurcating periodic solutions are asymptotically stable, and the direction of the bifurcation is forward (see Figures 4, 5, and 6).

(a)

(b)

(c)

(d)

Figure 4 

𝛾1=0.1, 𝛾2=2.5, 𝜏=1.5, 𝑘=3, which means that condition (𝐻2) holds, and 𝛽=0.6<𝛽0. The initial value is (0.1,−0.5,0.1,−0.5).

(a)

(b)

(c)

(d)

Figure 5 

𝛾1=0.1, 𝛾2=2.5, 𝜏=1.5, 𝑘=3, which means that condition (𝐻2) holds, and 𝛽0<𝛽=0.7<𝛽0. The initial value is (0.1,−0.5,0.1,−0.5).

(a)

(b)

(c)

(d)

Figure 6 

𝛾1=0.1, 𝛾2=2.5, 𝜏=1.5, 𝑘=3, which means that condition (𝐻2) holds, and 𝛽=1.2>𝛽0>𝛽0. The initial value is (0.1,−0.5,0.1,−0.5).

5. Conclusion

Ravoori et al. [1] explored an experimental system of two nominally identical optoelectronic feedback loops coupled unidirectionally, which are described by system (3). In the experiment, they found that depending on the value of the feedback strength 𝛽 and delay 𝜏, system (1) is capable of producing dynamics ranging from periodic oscillations to high-dimensional chaos [14, 15].

This paper investigates the stability and the existence of periodic solutions. We find that with the variety of the coupling strength 𝑘, even if all other parameters keep the same, the dynamical behavior can change greatly. In fact, it is clear that the first two equations, 𝑥₁(𝑡) and 𝑦₁(𝑡) are uncoupled with equations 𝑥₂(𝑡) and 𝑦₂(𝑡), so system (1) are independent of (2), which means that coupling strength 𝑘 does not appear in (1). The characteristic equation of (1) has the same form as (6), so the first critical value 𝛽₀ is independent of 𝑘. The analysis of characteristic equation (7) shows that the value of 𝑘 can affect the first critical value 𝛽₀ definitely. And we draw a conclusion that when 𝑘 is in an interval, in which 𝛽₀<𝛽₀ holds, solutions of system (1) and (2) keep synchronous; when 𝑘 belongs to the interval, in which 𝛽₀<𝛽₀ holds, solutions of system (1) and (2) can also keep synchronous with 𝛽<𝛽₀, while they lose their synchronization when 𝛽>𝛽₀, no matter whether 𝛽<𝛽₀ or not.

As a result, the modulation of the coupling strengths 𝑘 together with the feedback strength 𝛽 would be an efficient and an easily implementable method to control the behavior of the coupled chaotic oscillators.

Acknowledgment

This paper is supported by the National Natural Science Foundation of China (no. 11031002) and the Research Fund for the Doctoral Program of Higher Education of China (no. 20122302110044).