Double Hopf bifurcation induces coexistence of periodic oscillations in a diffusive Ginzburg–Landau model

Physics Letters A 383 (2019) 630–639

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Double Hopf bifurcation induces coexistence of periodic oscillations in a diffusive Ginzburg–Landau model Yanfei Du a,b , Ben Niu b,∗ , Yuxiao Guo b , Jianquan Li a a b

College of Arts and Sciences, Shaanxi University of Science and Technology, Xi’an 710021, China Department of Mathematics, Harbin Institute of Technology, Weihai 264209, China

a r t i c l e

i n f o

Article history: Received 15 August 2018 Received in revised form 2 January 2019 Accepted 8 January 2019 Available online 16 January 2019 Communicated by C.R. Doering Keywords: Complex Ginzburg–Landau system Diffusion Time delay feedback Double Hopf bifurcation Coexistence

a b s t r a c t We perform bifurcation analysis in a complex Ginzburg–Landau system with delayed feedback under the homogeneous Neumann boundary condition. We calculate the amplitude death region, and it turns out that the boundary of the amplitude death region consists of two Hopf bifurcation curves with wave number zero. The existence conditions for double Hopf bifurcations are established. Taking the feedback strength and time delay as bifurcation parameters, normal forms truncated to the third order at double Hopf singularity are derived, and the unfolding near the critical points is given. The bifurcation diagram near the double Hopf bifurcation is drawn in the two-parameter plane. The phenomena of amplitude death, the existence of stable bifurcating periodic solutions, and the coexistence of two stable periodic solutions with fast oscillation and slow oscillation respectively are simulated. © 2019 Elsevier B.V. All rights reserved.

1. Introduction The complex Ginzburg–Landau equation (CGLE), which is one of the most popular reaction–diffusion systems, has been widely used in the field of nonlinear science [1], since the CGLE describes dynamics around the onset of a Hopf bifurcation. The model of CGLE has been extensively used for the study of spatiotemporal chaos, where uniform oscillations, spiral waves, coherent structures and turbulence are found. Considerable efforts have been made to understand this type of chaotic behavior, and time-delay feedback scheme has been one of the efficient way to control turbulence. One of the most important time-delay feedback schemes was introduced by Pyragas, who proposed to generate a control signal proportional to the difference between the time-delayed state of the system and its current state [2]. Different kinds of feedback schemes have been applied to CGLE, such as local feedback [3,4], global feedback [5–7], and time delay feedback with both local and global terms [8]. The CGLE has been extensively studied from the point of view of mathematical analysis for studies related to the global attractors, inertial manifolds, and soft and hard turbulences etc., and we refer to [9–14] and the references therein. Amplitude death (AD), referring to a situation where oscillations cease, has been investigated by many researchers (Yam-

*

Corresponding author. E-mail address: [email protected] (B. Niu).

https://doi.org/10.1016/j.physleta.2019.01.016 0375-9601/© 2019 Elsevier B.V. All rights reserved.

aguchi and Shimizu [15], Aronson et al. [16]). The investigation of AD is useful in considering the complete suppression of oscillations [17–19]. The transition from oscillations to AD is obtained via some bifurcations. In this paper, we find out that the boundary of amplitude death region consists of two Hopf bifurcation curves. During the process of the stability switching, double Hopf bifurcation may arise, which is usually accompanied by coexisting of several oscillations, quasi-periodic oscillations, and even chaos. Double Hopf bifurcation analysis has been discussed widely. Guckenheimer and Holmes [20] gave the dynamics in the neighborhood of a double Hopf bifurcation point in an ordinary differential equation. Campbell and Bélair consider resonant double Hopf bifurcation in the damped harmonic oscillator with delayed forcing. The center manifold method was adopted to reduce delay differential equations into finite-dimensional systems, and then normal forms can be calculated on the center manifold [21–25]. Yu [26] proposed a perturbation method based on multiple scales for computing the normal forms of nonlinear dynamical systems. Through carrying out bifurcation analysis (mainly Hopf and double Hopf bifurcation analysis) of the equilibrium of system, we can figure out the structure of AD region, the effect of parameters on AD, and the critical values of parameters corresponding to the spatiotemporal behavior of dynamics of system. Thus we can control the state of system through restricting parameters. CGLE can be associated with Dirichlet [27,28], Neumann [17, 27,29,30], or periodic [31] boundary conditions depending on the physical situation. The Cauchy problem for the 1-D CGL with ei-

Y. Du et al. / Physics Letters A 383 (2019) 630–639

ther zero Dirichlet or zero Neumann boundary data is shown to be well-posed via classical techniques of nonlinear parabolic equations [32]. Ermentrout and Troy [17] considered phaselocking in a particular reaction–diffusion system with a continuous linear gradient of frequencies. Sirovich [29] considered CGLE in the case of stressless boundary conditions in Rayleigh–Bédnard convection, where the velocity parallel to a boundary has a zero normal derivative, but the velocity itself need not be zero so that the fluid can slip. The influence of boundary conditions was highlighted recently for a two-dimensional CGLE [27,33]. The complex Ginzburg–Landau equation (CGLE) [34]

∂ W (x, t ) = {(1 + i ω0 ) − (1 + i α )| W (x, t )|2 } W (x, t ) ∂t ∂2 + (1 + i β) 2 W (x, t ), ∂x

2. Amplitude death region and double Hopf bifurcation In this section, we investigate the stability of the equilibrium and the existence of Hopf bifurcation and double Hopf bifurcation for delayed Ginzburg–Landau system (2) by analyzing characteristic equations at the trivial steady state. It is clear that W = 0 is always a trivial steady state of system (2). The linearization of system (2) at the trivial steady state is

∂2 ∂ W (x, t ) = (1 + i β) 2 W (x, t ) + (1 − ε + i ω0 ) W (x, t ) ∂t ∂x + ε W (x, t − τ ).

(1)

− ϕ = μϕ , x ∈ (0, lπ ), ϕ  (0) = ϕ  (lπ ) = 0, has eigenvalues

∂ W (x, t ) = {(1 + i ω0 ) − (1 + i α )| W (x, t )|2 } W (x, t ) ∂t ∂2 + (1 + i β) 2 W (x, t ) ∂x + ε { W (x, t − τ ) − W (x, t )},

When tem.

(2)

with the homogeneous Neumann boundary condition and the initial condition:

μn =

n2 l2

(n = 0, 1, 2, · · · ), with corresponding

eigenfunctions cos Substituting W (x, t ) = n∞=0 ( pe λt cos nl x) into linear system (5) leads to the characteristic function n x. l

f n (λ) = 1 − ε + i ω0 − n2 /l2 (1 + i β) − λ + ε e −λτ = 0, n = 0, 1 , 2 , · · · .

(6)

ω0 = 0, β = 0 and α = 0, system (2) becomes RCGL sys-

Lemma 2.1. There does not exist amplitude death at the trivial steady state in the RCGL system, i.e. system (2) with ω0 = 0, β = 0 and α = 0. Proof. See Appendix A.

2

Now we consider the case of ω0 = 0. The trivial steady state is stable if and only if all roots of (6) has negative real part.

(3)

To analysis the structure of amplitude death region and the oscillation states clear, we first consider the case of α = 0, and β = 0, i.e.

∂2 ∂ W (x, t ) = {(1 + i ω0 ) − | W (x, t )|2 } W (x, t ) + 2 W (x, t ) ∂t ∂x + ε { W (x, t − τ ) − W (x, t )}.

(5)

It is well known that the eigenvalue problem

have been investigated widely, where W (x, t ) ∈ C is a complex state at time t ≥ 0 and location x ∈ [0, lπ ], ω0 > 0 is frequency, ε ≥ 0 is feedback strength, τ is time delay, α is the nonlinear frequency parameter, β is the linear dispersion coefficient. Motivated by [8,35], we incorporate the time delay feedback due to Pyragas [2] into the complex Ginzburg–Landau system,

∂W ∂W (0, t ) = (lπ , t ) = 0, t ≥ 0, ∂x ∂x W (x, t ) = ϕ (x, t ), − τ ≤ t ≤ 0.

631

(4)

After that, we will consider the effect of α and β . We focus on the joint effect of the feedback strength and time delay on the Ginzburg–Landau model from the view of bifurcation analysis. Taking the feedback strength and time delay as bifurcation parameters, the condition of the existence of double Hopf bifurcation is obtained. After some calculations on the Fourier expansion of the partial functional differential equations corresponding to the original system, the normal form truncated to third order is derived. The bifurcation diagram near the double Hopf bifurcation is drawn in the two-parameter plane. Simulations show that there are two stable homogeneous periodic solutions coexisting around the amplitude death region. The paper is organized as follows. In section 2, by analyzing characteristic equations at the equilibrium, the amplitude death region is found, and the existence conditions of Hopf bifurcation and double Hopf bifurcation are established. In section 3, explicit formulas of normal form truncated to the third order at double Hopf singularity are derived. In section 4, some numerical simulations are presented to illustrate the theoretical results. The paper ends with some conclusions in section 5.

Lemma 2.2. There does not exist amplitude death at the trivial steady state in the system (2) without time delay feedback. Proof. See Appendix B.

2

In order to stabilize the trivial steady state through time delay feedback, there must be roots of (6) on the axis and then crossing it into the left plane. It is obvious is not a root of (6). For τ > 0, we are searching for an root λ = i ξ (ξ = 0) of (6), then we have

increasing imaginary that λ = 0 imaginary

1 − ε − n2 /l2 = −ε cos ξ τ ,

ω0 − β n2 /l2 − ξ = ε sin ξ τ .

(7)

Squaring and adding both sides of the equations, we get

ξ 2 − 2(ω0 −β n2 /l2 )ξ +(ω0 −β n2 /l2 )2 +(1 − ε − n2 /l2 )2 − ε 2 = 0. (8) The roots of Eq. (8) are

⎧  ⎪ ξn± = ω0 − β n2 /l2 ± ε 2 − (1 − n2 /l2 − ε )2 , ⎪ ⎪ ⎨ 1−n2 /l2 if n < l, ε > , 2 ⎪ ξ = ω − β n2 /l2 , if n = l, n 0 ⎪  ⎪ ⎩ ± ξn = ω0 − β n2 /l2 ± i (1 − n2 /l2 − ε )2 − ε 2 , if n > l.

(9)

We can conclude that when n > l, the roots of (6) never intersect the imaginary axis for any ε and τ , thus the roots of (6) always have negative real part when n > l. When n = l, we have

632

Y. Du et al. / Physics Letters A 383 (2019) 630–639

Fig. 1. a) When π2 < ω0 ≤ 32π , β = 0, the amplitude death region A D = S 00 . b) When colors in the figure(s), the reader is referred to the web version of this article.) 2

Re ddτλ |λ=i ω0 = 0, Re ddτ λ2 |λ=i ω0 < 0, which means that the real part of roots of (6) reaches the maximum value 0 at τ = 0. Therefore, the roots of (6) have negative real part when n = l for any τ > 0. Thus the stability of trivial equilibrium depends on the real part of roots for (6) in the case of n < l. By Eq. (7), we have sin ξ + τ < 0, and sin ξ − τ > 0, and thus we obtain the following two sets of critical curves,

2 j π + arccos(n2 /l2 + ε − 1)/ε

τnj+ =

2 j π + 2π − arccos(n /l + ε − 1)/ε 2

,

|ξn+ |

=

{(ε , τ ) | ∅,

j−

j

1 2

(see in Fig. 1a)), moreover, when

ε

ω0 = 0 and

j+

2 j π +π 2

,

τ

j



,

j = 0, 1, · · · , J − 1 and J is the smallest integer such that

τ0J + (ε) holds for any ε .

τ0J − (ε) ≥

· · · < τn0− < τn0−−1 < · · · < τ10− < τ00− < τ00+ < τ10+ < · · · < τn0−+1 < τn0+ < · · ·

clude that there is a

and n < l.

1 0 εmax ,0 ∈ ( 2 ,

τ

ω02 +1 2

1 2

ε τ

τ

ε

τ

τ

τ

τ

τ

τ

tiating the two sides of Eq. (6) with respective to

 −1 dλ Re dτ

τ =τnj±

=−

1

ξn± (ω0 − ξn± )

that when

τ , we obtain

1

εξ ± sin ξ ± τ 1

 . ±ξn± ε 2 − (1 − n2 /l2 − ε )2

=

ω0 > 0, thus Re ddτλ |τ =τ j+ > 0. Re ddτλ |τ =τ j− < 0 n

0 < ε < εmax ,0 . We can con-

), such that τ00− = τ00+ , i.e.

ω0 >

3π 2

ε<

ω02 +(1−n2 /l2 )2 2(1−n2 /l2 )

n

. We can conclude

,

1− 1− n < n −1 < · · · 1+ 1+ n −1 < n < · · ·

··· < τ

holds

τ

< τ11− < τ01− < τ01+ < τ11+ · · ·

τ

for

1 2

1 < ε < εmax ,0 ,

1 1 2π +arccos(εmax ,0 −1)/εmax,0



1 ω0 − 2εmax ,0 −1

=

where

1 εmax ,0

is

1 1 2π +2π −arccos(εmax ,0 −1)/εmax,0



1 ω0 + 2εmax ,0 −1

defined

by

. Then there is

1− 1+ 1 < ε < εmax ,0 , τ ∈ (τ0 , τ0 )}, 3π 5π which is illustrated in Fig. 1b). So when 2 < ω0 < 2 , A D =  1 0

another stable region S 01 = {(ε , τ ) | S0

τ00− < τ00+ for

τ

τ

<τ

Proof. Since the roots of (6) always has negative real part when τ > 0 if n ≥ l, to find the amplitude death region, we only need to consider the case when n < l. In the following, we prove that when ω0 > π2

<ε<ε

τ

that the trivial equilibrium of (2) is asymptotically stable in region S 00 when π2 < ω0 < 32π . The stable region S 00 is illustrated in Fig. 1a). Moreover, if τ01− < τ01+ holds for some ε , we can also prove

j

j ω0 + 2εmax ,0 − 1

Firstly we investigate

0− 0+ 0 < ε < εmax (see ,0 , τ0 < τ0

τ

holds if ξn− > 0, which means

−1

0 max,0



0 ω0 + 2εmax ,0 −1

ε τ

τ

τ

Obviously, ξn+ >

2 j π + 2π − arccos(εmax,0 − 1)/εmax,0

holds for

ε

=− j

ω0 − 2ε

ε τ



2 j π + arccos(εmax,0 − 1)/εmax,0

1 2

ε

τ

j

=

1 and 2 0 max,0 ,

otherwise,

j max,0

0 0 2π −arccos(εmax ,0 −1)/εmax,0

ε

< ε < εmax,0 , τ ∈ (τ0 , τ0 )}, if ω0 >

j

=

0− and 00+ make a circuit 0 0 0 between max,0 , inside which there is a region S 0 = {( , ) | 0− 0+ 1 < < ∈ ( 0 , 0 )}, which is illustrated in Fig. 1a). 2 1−n2 /l2 0 Similarly, we can also prove that there is an max , ,n ∈ ( 2 ω02 +(1−n2 /l2 )2 1−n2 /l2 0− 0+ 0 ), such that n = n . Denote S n = {( , ) | 2 < 2(1−n2 /l2 ) 0 < max , ∈ ( n0− , n0+ )}. ,n ω2 +1 Secondly, we can prove that when 12 < < 02 , 00+ < 10+ < · · · < n0−+1 < n0+ < · · · (see Appendix D). ω2 +1 − Thirdly, we prove that when 12 < < 02 , · · · < n0− < n0− 1 < 0− 0− · · · < 1 < 0 (see Appendix E). From the previous discussion, we have S 00 ⊂ S 10 ⊂ · · · . Differen-

where εmax,0 is defined by



1 2

Appendix C). It means that the curve

ε

,

Theorem 2.3. The amplitude death region of Eq. (2) with J −1 β = 0 is A D = S 00 ∪ S 01 ∪ · · · ∪ S 0 with





0 ω0 − 2εmax ,0 −1

τ

j = 1, 2, · · · .

2

In the following, to character the amplitude death region of Eq. (2) in detail, in which the trivial steady state is asymptotically stable, we consider the case with ω0 = 0 and β = 0.

j S0

0 0 arccos(εmax ,0 −1)/εmax,0

0 εmax ,0 is defined by

ε

τnj− =

|ξn− |

 < ω0 ≤ 52π , the amplitude death region A D = S 00 S 01 . (For interpretation of the

3π 2

1 2

S0. Since ξn− < ξn+ , there must exist the smallest integer J such

that

τ0J − > τ0J + for any ε . Thus, the amplitude death region A D = J −1

S 00 ∪ S 01 ∪ · · · ∪ S 0 . The proof is complete.

2

Y. Du et al. / Physics Letters A 383 (2019) 630–639

Fig. 2. a) Let l = 3,

633

ω = 2. When we choose a small β = 0.5, the amplitude death region A D = S 00 is bounded by τ00+ and τ00− . Thus, the small β have no inference on AD τ00+ and τ20− .

region. b) When β increases to 1, the amplitude death region A D = S 00 becomes smaller, which is bounded by

From the above discussion, we know that effect on the amplitude death region.

ω0 have important

a

ω0 = 0 and β = 0 in system (2). When ω0 ≤ π2 , the AD region of system (2) does not exist. When π2 < ω0 ≤ 32π , the AD region of system (2) in the ε − τ plane is given by A D = S 00 . 2 j π +π When < ω0 ≤ 2( j +12)π +π , the AD region of system (2) in the 2   ε − τ plane is given by A D = S 00 · · · S 0J −1 . It means that, the

j max,0

ε

∈

( 12 ,

We will see in section 4 that simulations support this analysis result (see Fig. 4 and Fig. 5a)). Remark 2. Now we discuss the effect of β on AD region of sysj± j− j+ tem (2). Since τ0 do not depend on β , τ0 < τ0 still holds when β > 0. We can also prove that β do not change the monotonicity τnj+ , i.e. it is still monotonically increasing on n just like the case j−

of β = 0. When β is a small positive number, τn is still monotonically decreasing on n, which can be broken by the increasing of β . Thus, small β makes no contribution to the AD region. Similarly as the result given in Theorem 2.3, and diagram shown in Fig. 1,   J −1 we have A D = S 00 · · · S 0 . However, when β get larger, AD region becomes smaller (see Fig. 2). When β is large enough, AD region does not exist. Remark 3. We can easily get that AD region of system (2).

α makes no contribution to the

ω02 +1 2

), such

τ

( j = 0, 1, 2, · · · ), then there exists j−

j

= τ0 , denoted by τ0 . System (4) j j max,0 , 0 ).

undergoes a double Hopf bifurcation at (ε

Remark 1. Suppose

AD region becomes larger and the amount of AD regions increases, for larger ω0 .

2 j π +π 2 j+ that 0

Theorem 2.5. Suppose ω0 >

τ

3. Normal form of double Hopf bifurcation From the previous section, system (4) undergoes a double Hopf j j bifurcation at (εmax,0 , τ0 ). To make the dynamics of system (4) near the double Hopf singularity clear, taking μ = (τ , ε ) as bifurcation parameters, we determine the third-order normal forms of double Hopf bifurcation of system (4). The method we use is due to the normal form method of PFDEs proposed by Faria [24], which is based on the center manifold theory introduced in [37]. Let u 1 (x, t ) = Re( W (x, t )), u 2 (x, t ) = Im( W (x, t )), then system (4) can be written as

⎧ ⎪ ∂ u 1 (x, t ) ∂ 2 u 1 (x, t ) ⎪ ⎪ = + (1 − ε )u 1 (x, t ) − ω0 u 2 (x, t ) ⎪ ⎪ ∂t ∂ x2 ⎪ ⎪ ⎪ ⎪ ⎪ + ε u 1 (x, t − τ ) − (u 21 (x, t ) + u 22 (x, t ))u 1 (x, t ), ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ∂ u 2 (x, t ) ∂ 2 u 2 (x, t ) = + ω0 u 1 (x, t ) + (1 − ε )u 2 (x, t ) (10) ∂t ∂ x2 ⎪ ⎪ 2 2 ⎪ + ε u ( x , t − τ ) − ( u ( x , t ) + u ( x , t )) u ( x , t ), ⎪ 2 2 1 2 ⎪ ⎪ ⎪ ⎪ ∂ ∂ ∂ ∂ u u u u ⎪ 1 1 2 2 ⎪ ⎪ (0, t ) = (lπ , t ) = 0, (0, t ) = (lπ , t ) = 0, ⎪ ⎪ ∂x ∂x ∂x ∂x ⎪ ⎪ ⎩ t ≥ 0. It is obvious that E 0 (0, 0) is the trivial equilibrium of system (10). Corresponding to the Neumann boundary condition, define the real-valued Hilbert space

 Obviously, the AD region of Eq. (4) can be obtained by Theorem 2.3. Due to the general Hopf bifurcation theorem [24,36], we have the following theorem.

Theorem 2.4. Suppose that

1−n2 /l2 2

<ε<

ω02 +(1−n2 /l2 )2 2(1−n2 /l2 )

holds, then

system (4) undergoes a Hopf bifurcation at the trivial equilibrium when τ = τnj− or τnj+ , n = 0, 1, · · · , l − 1, j = 0, 1, 2 · · · . From Theorem 2.3, the boundary curves of the stable region j+ j− for the trivial equilibrium are τ0 and τ0 , which means that the double Hopf bifurcation points is the intersection of two Hopf bifurcation curves with the wave number n/l = 0.

X := (u 1 , u 2 ) T ∈ H 2 (0, lπ ) × H 2 (0, lπ ) :

∂ ui ∂ ui (0, t ) = (lπ , t ) = 0, i = 1, 2 . ∂x ∂x Define the corresponding complexification X C := X ⊕ i X = {U 1 + iU 2 : U 1 , U 2 ∈ X } and the complex-value L 2 inner product u , v =

 lπ 0

(u 1 v 1 + u 2 v 2 )dx, for u = (u 1 , u 2 )T , v = ( v 1 , v 2 )T ∈ X C . Let U (t ) = (u 1 (x, t ), u 2 (x, t )) T , system (10) can be written as

dU (t ) dt

= E2

∂2 U (t ) + ( A − ε E 2 )U t (0) + ε G 1 U t (−τ ) + f (U t ), ∂ x2 (11)

where

634

Y. Du et al. / Physics Letters A 383 (2019) 630–639

A=

1

−ω0

ω0

1



 t

, G 1 = E 2 , f (U ) =

−(u 21 + u 22 )u 1 −(u 21

+ u 22 )u 2

 and the adjoint bilinear form

,

0 θ

and E 2 is 2 × 2 identity matrix. For the sake of notation, we denote j j ∗ ∗ εmax ,0 by ε , and τ0 by τ . Carry out the rescaling t → t /τ , and define (τ , ε ) = (τ ∗ + α1 , ε ∗ + α2 ), then we have an equivalent form of (11)

dU (t ) dt

∂2 U (t ) + L (τ ∗ + α1 , ε ∗ + α2 )U t ∂ x2 + F (τ ∗ + α1 , ε ∗ + α2 , U t ), (12)

= D (τ ∗ + α1 , ε ∗ + α2 )

where D (τ ∗ + α1 , ε ∗ + α2 ) = (τ ∗ + α1 ) E 2 , L (τ ∗ + α1 , ε ∗ + α2 ) = (τ ∗ + α1 )( A − (ε ∗ + α2 ) E 2 )U t (0) + (τ ∗ + α1 )(ε ∗ + α2 )G 1 U t (−1), and F (τ ∗ + α1 , ε ∗ + α2 , U t ) = (τ ∗ + α1 ) f (U t ). Then Eq. (12) can be rewritten as

(ψ, φ)k = ψ(0)φ(0) −

ψ(ξ − θ)dηk (θ)φ(ξ )dξ.

−1 0

The calculations of normal form are very long, so we leave them in Supplementary Material. Based on the derivation in Supplementary Material, the normal form truncated to the third order for double Hopf bifurcation reads as

z˙ 1 = i ξ0− τ ∗ z1 + B α1 z1 α1 z1 + B α2 z1 α2 z1 + C 2100 z12 z2

+ C 1011 z1 z3 z4 , z˙ 2 = −i ξ0− τ ∗ z2 + B α1 z1 α1 z2 + B α2 z1 α2 z2 + C 2100 z1 z22

+ C 1011 z2 z3 z4 ,

∂2 U (t ) + L 0 U t +  G (α , U t ), (13) dt ∂ x2 where D 0 = τ ∗ E 2 , L 0 U t = (τ ∗ A − τ ∗ ε ∗ E 2 )U t (0) + τ ∗ ε ∗ G 1 U t (−1),

z˙ 3 = i ξ0+ τ ∗ z3 + B α1 z3 α1 z3 + B α2 z3 α2 z3 + C 0021 z32 z4

 G (α , U ) = ( D (τ ∗

z˙ 4 = −i ξ0+ τ ∗ z4 + B α1 z3 α1 z4 + B α2 z3 α2 z4 + C 0021 z3 z42

dU (t )

= D0

∂2 + α1 + α2 ) − D 0 ) 2 U (t ) + ( L (τ ∗ + α1 , ε ∗ + ∂x α2 ) − L 0 )U t + F (τ ∗ + α1 , ε∗ + α2 , U t ). , ε∗

t

Consider the linearized system of (13)

dU (t ) dt

(14)

A : C01 ∩ BC → BC , A ϕ = ϕ˙ + X 0 [ D 0 ϕ (0) + L 0 (ϕ ) − ϕ˙ (0)], with dom( A ) = {ϕ ∈ C : ϕ˙ ∈ C , ϕ (0) ∈ dom( )}, C = C ([−1, 0], X C ), BC := {ψ : [−1, 0] → X C : ψ is continuous on [−1, 0),  0, −1 ≤θ < 0, ∃ limθ→0− ψ(θ) ∈ X C }, and X 0 is given by X 0 (θ) = I , θ = 0. In the space BC , (13) becomes an abstract ODE

G (α , U ). U = AU + X 0 t

dt

t

t

μn =

It is well known that

(15) 2 − nl2

with n ∈ N0 are eigenvalues

∂2 , with corresponding normalized eigenfunctions ∂ x2

of

cos nl x

(1)

 cos nl x  L 2

γn (x) =

(2)

. Let βn (x) = γn (x)(1, 0) T , βn (x) = γn (x)(0, 1) T . Then

( j)

{βn }n≥0 are eigenfunctions of 2

∂2 with corresponding eigenvalues ∂ x2

− nl2 ( j = 1, 2), which form an orthonormal basis of X . Define Bn   ( j) ( j) the subspace of C , by Bn := span v (·), βn βn | v ∈ C , j = 1, 2 .   (1) (2) T Denote v (·), βn = v (·), βn , v (·), βn for simplification of notations. Then on Bn , the linear equation d dt

dt

=−

n2 l2

− ∗ B α1 z1 = 2D 1 (1 − ε ∗ + ε ∗ e −i ξ0 τ + i ω0 ), − ∗ B α1 z1 = 2D 1 (−τ ∗ + τ ∗ e −i ξ0 τ ), + ∗ B α1 z3 = 2D 2 (1 − ε ∗ + ε ∗ e −i ξ0 τ + i ω0 ), + ∗ B α2 z3 = 2D 2 (−τ ∗ + τ ∗ e −i ξ0 τ ),

C 2100 = − C 1011 = − C 0021 = − C 1110 = −

8 lπ 16 lπ 8 lπ 16 lπ

D1τ ∗,

(18)

D1τ ∗, D2τ ∗, D2τ ∗.

Make the polar coordinate transformation

z1 = ρ1 cos θ1 + i ρ1 sin θ1 , z2 = ρ1 cos θ1 − i ρ1 sin θ1 , z3 = ρ2 cos θ2 + i ρ2 sin θ2 , z4 = ρ2 cos θ2 − i ρ2 sin θ2 ,

ρ1 , ρ2 > 0. Define 1 = Sign(Re C 2100 ) = Sign(Re(− lπ8 D 1 τ ∗ )), 2 = Sign(Re(− lπ8 D 2 τ ∗ )). From√Re D 1 = Re D 2 , we have 1 = 2 . √ Carry out the rescaling  r1 = ρ1 |C 2100 |,  r2 = ρ2 |C 0021 |,  t = t 1 , where

ρ˙1 = ρ1 (c1 + ρ12 + b0 ρ22 ),

(19)

ρ˙2 = ρ2 (c2 + c0 ρ12 + d0 ρ22 ), where

t

D 0 z(t ) + L 0 z .

Define the matrix-valued function of bounded variation B V ([−1, 0], R) such that

n2 − 2k l

+ C 1110 z1 z2 z4 ,

and drop the hats, then system (17) becomes

U (t ) = D 0 U (t ) + L 0 (U t )

can be written as the retarded functional differential equation on C2 :

dz

+ C 1110 z1 z2 z3 ,

where

∂2 = D 0 2 U (t ) + L 0 U t . ∂x

We have that the solution operator of (13) is a C 0 semigroup, and the infinitesimal generator A is given by

d

(17)

0 D 0 ϕ (0) + L 0 (ϕ ) = −1

dηk (θ)ϕ (θ), ϕ ∈ C ,

(16)

c 1 = 1 (Re B α1 z1 α1 + Re B α1 z1 α2 ),

ηk ∈

c 2 = 1 (Re B α1 z3 α1 + Re B α2 z3 α2 ), b0 =

1 2 ReC 1011 Re C 0021

, c0 =

Re C 1110 Re C 2100

, d0 = 1 2 .

From chapter 7.5 in [20], Eq. (19) has twelve kinds of unfoldings (see Table 1).

Y. Du et al. / Physics Letters A 383 (2019) 630–639

635

Table 1 The twelve unfoldings of system (19). Case

Ia

Ib

II

III

IVa

IVb

V

VIa

VIb

VIIa

VIIb

VIII

d0 b0 c0 d0 − b 0 c 0

+1 + + +

+1 + + −

+1 + − +

+1 − + +

+1 − − +

+1 − − −

−1 + + −

−1 + − +

−1 + − −

−1 − + +

−1 − + −

−1 − − −

In fact, we have b0 = 2, c 0 = 2, d0 = 1, and thus d0 − b0 c 0 = −3. It means that the unfolding system is of type Ib in Table 1, thus we draw bifurcation set and phase portraits for the unfolding of case Ib in Fig. 3 [20]. Remark 4. From Theorem 2.2 in Section 6.2 of [36], system (13) has two families of periodic solutions bifurcating from (0, 0) parameterized by small  . When μ = μ( ) and  are near 0, (i.e. when ε is near ε ∗ ), the periodic solutions have the following representations − ∗ U t1 (μ, θ)(x) =  Re φ1 (θ)e i ξ0 τ t + O ( 2 ),

(20)

and + ∗ U t2 (μ, θ)(x) =  Re φ3 (θ)e i ξ0 τ t + O ( 2 ),

(21)

where φ1 (θ) and φ3 (θ) is defined in (1) of Supplementary Material. 4. Simulations

Fig. 3. Phase portraits for the unfoldings of case Ib with

Fig. 4. a) When

ω0 = 5.5, S 00 = {(ε , τ ) |

Fig. 5. a) When ω0 = 2.5, the partial bifurcation set on the the complete bifurcation sets near HH.

1 2

1 = −1.

In this section, we fix l = 4, β = 0, α = 0, and take ε and τ as bifurcation parameters. From Theorem 2.4, system (4) undergoes

< ε < 11.91, τ ∈ (τ00− , τ00+ )}. b) When ω0 = 5.5, S 01 = {(ε , τ ) |

ε − τ plane. The amplitude death region A D = S 00 = {(ε , τ ) |

1 2

1 2

< ε < 0.78, τ ∈ (τ01− , τ01+ )}.

< ε < 1.8117, τ ∈ (τ00− , τ00+ )}. b) When ω0 = 2.5,

636

Y. Du et al. / Physics Letters A 383 (2019) 630–639 j±

a Hopf bifurcation at the trivial equilibrium when τ = τn (n = 0, 1, 2, 3), and there are no Hopf bifurcations when n ≥ 4. Firstly, we figure out the amplitude death region with different values of ω0 . If we fix ω0 = 5.5 > 32π , we have τ02+ > τ02− for any ε , j

thus S 0 is empty for j ≥ 2. From Theorem 2.3, we know that

1 < < 2 1− 1+ , 0 0 )},

the amplitude death region A D = S 00 ∪ S 01 = {(ε , τ ) |

ε 1 τ ∈ (τ τ ∪ {(ε , τ ) | < ε < εmax , τ ∈ ( τ τ ,0 0 1 where εmax ,0 = 11.91 and εmax,0 = 0.78. We draw the bifurcation 0 max,0 ,

ε

0− 0+ 0 , 0 )}

1 2

diagram in Fig. 4. If we fix ω0 = 2.5 > π2 , we have j S0

τ01+ > τ01− for any ε , thus

is empty for j ≥ 1. From Theorem 2.3, we know that the

amplitude death region A D = S 00 = {(ε , τ ) |

Fig. 6. When

Fig. 7. When region.

1 2

0 < ε < εmax ,0 , τ ∈

ω0 = 2.5, Eq. (9) with β = 0 and n = 0 has two positive roots.

0 (τ00− , τ00+ )}, where εmax ,0 = 1.8117. We draw the bifurcation di-

agram in Fig. 5a). Comparing the AD region when ω0 = 2.5 and the AD region when ω0 = 5.5, it is easy to see that larger ω0 makes the AD region larger and increases the amount of AD regions, which coincides with the analysis result in Remark 1. Now we consider the double Hopf bifurcation. When 1/2 < ε < 3.625, Eq. (9) with β = 0 and n = 0 has two positive roots (see + 0 Fig. 6). When ε = εmax ,0 = 1.8117, two roots are ξ0 = 4.1197,

ξ0− = 0.8803, and we have τ00+ = τ00− = 1.2566, correspondingly. 0 It means that when ε = εmax ,0 = 1.8117, the two Hopf curves

τ00+ and τ00− intersect, and we denote the double Hopf bifur-

cation point HH (see Fig. 5a)). For HH, we can get the coefficients of normal form B α1 z1 = −0.2179 + 0.2162i, B α1 z1 = 0.1077 − 0.4477i, B α1 z3 = 1.0197 + 1.0120i, B α2 z3 = 0.1077 + 0.4477i, C 2100 = −0.0983 − 0.0990i, C 1011 = −0.1965 − 0.1980i, C 0021 = −0.0983 + 0.0990i, C 1110 = −0.1965 + 0.1980i. Furthermore, we can get 1 = −1, 2 = −1, b0 = 2, c 0 = 2, d0 = 1, d0 − b0 c 0 = −3. According to Fig. 3, we have the bifurcation set near HH showing in Fig. 5b), in which the two black lines are two pitchfork bifurcation curves τ = (ε − 1.8117)/(13.5158) + 1.2566, τ = (ε − 1.8117)/(−20.9621) + 1.2566. The bifurcation set is divided into six regions. In D1, the equilibrium E 0 (0, 0) is asymptotically stable, which is shown in Fig. 7. It means that the region D1 is an amplitude death region. When the parameters (τ , ε ) vary and enter the region D2 (or D6), the stable equilibrium bifurcates into a stable periodic solution, which is shown in Fig. 8 (Fig. 9). For parameters in D3 (or D5), periodic solutions appear via secondary

ω0 = 2.5, ε = 1 and τ = 1.25 in D1, the equilibrium E 0 (0, 0) is asymptotically stable. The amplitude tends to zero, which means that D1 is an amplitude death

Y. Du et al. / Physics Letters A 383 (2019) 630–639

Fig. 8. When ω0 = 2.5, ε = 1.6 and (4) has a periodic solution.

637

τ = 1.29 in D2, the bifurcating periodic solution is stable. The amplitude tends to a positive value, which means the original CGL system

Fig. 9. When

ω0 = 2.5, ε = 1.9 and τ = 1.23 in D6, the bifurcating periodic solution is stable.

638

Y. Du et al. / Physics Letters A 383 (2019) 630–639

Fig. 10. When ω0 = 2.5, ε = 1.85 and τ = 1.257 in D4, two stable periodic solutions with slow oscillation and fast oscillation coexist when we choose different initial conditions. In the left, the initial conditions are u 1 (x, θ) = 5, u 2 (x, θ) = 5, θ ∈ [−τ , 0]; in the right, the initial conditions are u 1 (x, θ) = 0.01, u 2 (x, θ) = 0.01, θ ∈ [−τ , 0].

supercritical Hopf bifurcations, which are type of saddle and only stable on the center manifold. When the parameters enter the region D4, there are two stable periodic solutions coexisting in D4. Moreover, from Remark 4, the two periodic solutions have − ∗

 Re φ1 (θ)e iξ0 τ t + O ( 2 ), and U t2 (μ, θ)(x) =  Re φ3 (θ)e + O ( 2 ), respectively. Since ξ0+ > − ξ0 (see Fig. 6), we know that U t2 (μ, θ)(x) oscillates faster than U t1 (μ, θ)(x). Thus, we can see the phenomena of fast oscillation the representation U t1 (μ, θ)(x) = i ξ0+ τ ∗ t

and slow oscillation coexisting. It shows that for the two coexisting periodic solutions, although the amplitudes of W (x, t ) are the same (see Fig. 10a) and b)), the frequencies are different (see Fig. 10c) and d)).

Acknowledgements The authors wish to express their gratitude to the editors and the reviewers for the helpful comments. Y. Du is supported by Education Department of Shaanxi Province (grant No. 18JK0123) and National Natural Science Foundation of China (No. 61872227). B. Niu and Y. Guo are supported by National Natural Science Foundation of China (grant Nos. 11701120 and No. 11771109) and the Foundation for Innovation at HIT(WH). Appendix A. Proof of Lemma 2.1 In this case, f 0 (λ) = 1 − ε − λ + ε e −λτ = 0, we have f 0 (0) = 1 > 0, and lim f 0 (λ) = −∞, thus, f 0 (λ) = 0 has at least one positive λ→+∞

5. Conclusion

root independent of

ε and τ . Therefore, the origin is unstable.

Appendix B. Proof of Lemma 2.2 In this paper, we investigate the amplitude death region and nonresonant double Hopf bifurcation for a one-dimensional CGLE with time delay feedback. By the characteristic equation analysis, the amplitude death region is found, and the condition of existence of double Hopf bifurcation is given. Applying center manifold theorem and normal form method, the normal form truncated to third order near the bifurcation point is computed, and all coefficients of normal form are expressed by the parameters of the original system. The bifurcation diagram near the double Hopf bifurcation point is drawn in the two-parameter plane (τ , ε ). Simulations also illustrate the phenomena of fast oscillation and slow oscillation coexisting.

2

2

When τ = 0, f n (λ) has a root 1 − nl2 + i (ω0 − nl2 β). When n < l it has a positive real part; when n = l, it has a real part of 0; and when n > l, it has a negative real part, and thus W = 0 is unstable. Appendix C. Proof of τ00− < τ00+ for Theorem 2.3 In fact,

dτ00− dε

√1 [− ε1 ξ0− ξ0−2 2ε −1 dH (ε ) ε−1

=

− 1ε ξ0− + arccos ε . We have

1 2

0 < ε < εmax ,0 in

1 + arccos ε− ε ]. Denote H (ε ) = √ = ε12 (ω0 − 2ε − 1) > 0 when dε

Y. Du et al. / Physics Letters A 383 (2019) 630–639 1 2

ω02 +1

<ε<

arccos ω02 +1 2

ω02 −1 ω02 +1

2

. H ( 12 ) = −2ω0 +

π < 0 if ω0 >

π , H ( ω0 +1 ) = 2 2 2

1 2

> 0. Thus, we know that when ε vary from

to

, H (ε ) vary from a negative value to a positive value. Cor-

respondingly,

τ00− decreases at the beginning and then increases ω2 +1

ε vary from 12 to 02 . Similarly, we can prove that when π ω0 > 2 , τ00+ increases at the beginning and then decreases. Notice that τ00− = τ00+ when ε = 12 , and τ00− >> τ00+ when ε is close when

to

ω02 +1 2

that

ω2 +1

1 0 0 εmax ,0 ∈ ( 2 , 2 ), such 0− 0+ 0 < ε < εmax ,0 , τ0 < τ0 .

. We can conclude that there is a

τ00− = τ00+ , moreover, when

1 2

+ 0+ Appendix D. Proof of τ00+ < τ10+ < ··· < τn0− 1 < τn < ··· for 1 2

<ε<

ω02 +1 2

in Theorem 2.3

To consider the monotonicity of

τn0+ , treat

2π − arccos(n2 /l2 + ε − 1)/ε

τn0+ =

ξn+

as a continuous function of n.

dτn0+ dn

=

1 2n +  [ξn l2 2 ξn+2 ε 2 −(1− n2 −ε )2

−

l

(2π − arccos n arccos

/l2 +ε −1

)(1 −

)(1 −

n2 l2

2

n2 /l2 +ε −1

ε

ε

− ε ),

2n/l2

1

− ε )]. Let G (n) = ξn+ − (2π − dG (n) dn

2 2 ε (1 − n /l − ε ) + (2π

n2 /l2 +ε −1 2 1−( )

ε

n2 l2

1−n2 /l2 −ε 2n − ε2 −(1−n2 /l2 −ε)2 l2 2 2 − arccos n /l ε+ε−1 ) 2n = l2

=



> 0, for n > 0. Denote p (ε ) = G (0) = (2π − arccos n /l ε+ε−1 ) 2n l2 1 ξ0+ − (2π − arccos ε− )( 1 − ε ). By direct calculation, we have ε 2

2

√

1 π = 2εε−1 + 2π − arccos ε− ε > 0. Therefore, when ω0 > 2 , 1 π G (0) = p (ε ) > p ( 2 ) = ω0 − 2 > 0, and G (n) > G (0) > 0, and thus dp (ε ) dε dτn0+ dn

> 0.

− < ··· < τ 0− < τ 0− for Appendix E. Proof of ··· < τn0− < τn0− 1 1 0 1 2

<ε<

ω02 +1 2

in Theorem 2.3

For the monotonicity of dτn0−

τn0− , we can prove that when

1 2

<ε<

1, dn < 0 if ω0 > π2 using the similar method as above. When 2 2 ε > 1, we can easily verify that arccos n /l ε+ε−1 is monotonically decreasing, and ξn− is monotonically increasing, and thus monotonically decreasing.

τn0− is

639

Appendix F. Supplementary material Supplementary material related to this article can be found online at https://doi.org/10.1016/j.physleta.2019.01.016. References [1] I. Aranson, L. Kramer, Rev. Mod. Phys. 74 (2001) 99–143. [2] K. Pyragas, Phys. Lett. A 170 (1999) 421–428. [3] M.E. Bleich, J.E. Socolar, Phys. Rev. E, Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 54 (1996) R17. [4] I. Harrington, J.E. Socolar, Phys. Rev. E, Stat. Nonlinear Soft Matter Phys. 64 (2001) 056206. [5] D. Battogtokh, A. Preusser, A. Mikhailov, Physica D 106 (1997) 327–362. [6] C. Beta, A.S. Mikhailov, Physica D 199 (2004) 173–184. [7] J.B. Gonpe Tafo, L. Nana, T.C. Kofane, Phys. Rev. E, Stat. Nonlinear Soft Matter Phys. 88 (2013) 032911. [8] M. Stich, C. Beta, Physica D 239 (2010) 1681–1691. [9] M. Bartuccelli, P. Constantin, C.R. Doering, J.D. Gibbon, Physica D 44 (1990) 421–444. [10] C.R. Doering, J.D. Gibbon, D.D. Holm, B. Nicolaenko, Nonlinearity 1 (1999) 279–309. [11] H.G. Kaper, P. Takc, Nonlinearity 11 (1998) 291–305. [12] I. Kukavica, Nonlinearity 8 (1995) 113–129. [13] R. Montagne, E. Hernández-García, Phys. Lett. A 273 (2000) 239–244. [14] Q. Tang, S. Wang, Physica D 88 (1995) 139–166. [15] Y. Yamaguchi, H. Shimizu, Physica D 11 (1984) 212–226. [16] D. Aronson, G. Ermentout, N. Kopell, Physica D 41 (1990) 403–449. [17] G.B. Ermentrout, W.C. Troy, SIAM J. Appl. Math. 46 (1986) 359–367. [18] D.V.R. Reddy, A. Sen, G.L. Johnston, Phys. Rev. Lett. 85 (1998) 3381–3384. [19] Y. Guo, B. Niu, Nonlinearity 28 (2015) 1841–1858. [20] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Appl. Math. Sci., vol. 42, Springer, New York, 1983. [21] S.A. Campell, J. Bélair, Can. Appl. Math. Q. 3 (1995) 137–154. [22] C. Elphick, E. Tiraopegui, M.E. Brachet, P. Coullet, G. Iooss, Physica D 29 (1987) 95–127. [23] T. Faria, L.T. Magalhães, J. Differ. Equ. 122 (1995) 181–200. [24] T. Faria, Trans. Am. Math. Soc. 352 (2000) 2217–2238. [25] P.L. Buono, J. Bélair, J. Differ. Equ. 189 (2003) 234–266. [26] P. Yu, Nonlinear Dyn. 27 (2002) 19–53. [27] V.M. Eguiluz, E. Hernandezgarcia, O. Piro, Int. J. Bifurc. Chaos Appl. Sci. Eng. 9 (1998) 2209–2214. [28] L. Nana, E.I. Mutabazi, Proc. R. Soc., Math. Phys. Eng. Sci. 465 (2009) 2251–2265. [29] L. Sirovich, Physica D 37 (1989) 126–145. [30] L.R. Keefe, Stud. Appl. Math. 73 (1985) 91–153. [31] S. Kuksin, A. Shirikyan, J. Phys. A 37 (2003) 3805–3822. [32] J.L. Lions, Quelques méthodes de résolution des problémes aux limites non linéaires, Dunod, Paris, 1969. [33] V.M. Eguiluz, E. Hernáandez-García, O. Piro, Phys. Rev. E 64 (2001) 036205. [34] M. Ipsen, L. Kramer, P.G. Sørensen, Phys. Rep. 337 (2000) 193–235. [35] H. Teki, K. Konishi, N. Hara, Phys. Rev. E 95 (2017) 062220. [36] J. Wu, Theory and Applications of Partial Functional-Differential Equations, Springer, New York, 1996. [37] X. Lin, J.W.H. So, J. Wu, Proc. R. Soc. Edinb. 122 (1992) 237–254.