Stability and Hopf bifurcation periodic orbits in delay coupled Lotka-Volterra ring system

1 Introduction

Mathematical models are typically used in ecology to illustrate the basic processes and dynamic mechanisms of ecosystems [1, 2, 3]. By the description and analysis of dynamics models, the essential characteristics of life processes can be understood more impressively. Mathematically, we usually describe the ecological mathematical model depending on the theory of functional differential equation. Among them, Lotka-Volterra system described by the theory of functional differential equation is one of the most famous and important ecological population dynamic models [4, 5].

where g_i represents the linear growth rate of species i, a_ij represents the interaction between species i and j, A = (a_ij) represents the interaction matrix with a_ij. It is easy to see that the model is bidirectional, that is, the growth of the species i depends on self-feedback and feedback from the species i + 1 and i − 1. Among these feedbacks, the interaction between species is not necessarily symmetrical. So in general a_ij ≠ a_ji. Without loss of generality, we assume that all g_i = a_i = 1, which is equivalent to the carrying capacity of each population x_i in the absence of other species and the unit time of the reverse growth rate of each species [5].

Golubitsky and his collaborators [6] proved that some phase relations can be modeled by coupled systems with observing the gait of animals. Under some conditions, coupled systems produce vibration, while uncoupled systems do not produce vibration. Therefore, the rich dynamic characteristics of the coupled oscillator can be understood by discussing the coupled system. Based on the practical significance of coupled Lotka-Volterra ring system, it has been widely used in natural science and ecology. There are many researchers who have conducted deeply studies on the dynamic characteristics of Lotka-Volterra system such as stable, unstable and oscillatory behavior [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. In [20, 21], various continuity theorems play an important role in studying the existence of periodic solutions for Lotka-Volterra systems.

Recently, symmetry has grown up to be an important subject in the study of nonlinear dynamic systems. Generally speaking, symmetry reflects some spatial invariants of dynamical systems. When the system is symmetric, it can exchange with the action of a compact Lie group Γ in Euclidean space [6]. Although symmetry makes the analysis of system more complicated, but it also imposes many special restrictions on the system. Bifurcation can occur to the system with specific symmetry and smaller size, and in some cases, even if the system is symmetric, bifurcation will not occur. Bifurcation phenomenon refers to the qualitative change of some attributes of the object of study. In nature, bifurcation phenomenon is ubiquitous. Therefore, whether in mathematical theory or in practical application, bifurcation theory research has its great significance, especially in symmetric system. So the symmetric Lotka-Volterra system will produce more interesting bifurcation phenomenon, which have been studied initially by some scholars [23, 24]

We mainly study a class of Lotka-Volterra ring systems with coupling [25].

where i represents the number of species from 1 to N, and assume x_N+1 = x₁ for create periodic boundary conditions. The system is a closed coupled Lotka-Volterra system consisting of N identical species. Each species competes with two of the four neighbouring species for limited resources. Splott [25] studied this coupled ring system deeply and found that the system exhibits spatiotemporal chaos in a spatial dimension, and its quasi-periodic paths are chaos, bifurcation, spontaneous symmetry destruction and spatial pattern formation, however, due to its impossible connectivity, it is not a very realistic model, because the author neglected the growth cycle of species and did not consider the effect of time delay on the model (1.1). Because time delay is an important controls parameter, it is imperative to introduce species growth time. Various scholars also incorporate time delays into the symmetric model [26, 27, 28, 29, 30, 31, 32, 33, 34]. Wu [27, 28, 29] used bifurcation theory to study local and global Hopf bifurcations of symmetric functional differential equations. Zheng and Zhang [30, 31, 32] obtained some results of the symmetric neural network model with delay. Hu [34] and Guo [35] discussed Hopf bifurcation periodic orbits and the spatial patterns of periodic orbits respectively in neuron ring systems with delay.

Based on the original model, the dynamic behavior of three identical species connected into a ring system is considered, a new three-dimensional coupled Lotka-volterra ring system with delays is constructed by adding appropriate delay τ. We just introduce the time delay into adjacent species, assuming that any species is coupled with the nearest species and symmetrical. In the model (1.1), a simpler choice of parameters are a₋₂ = a₂ = 0 and a₋₁ = a₁ = b.

where b > 0 and τ > 0. Because the system is symmetry, then the characteristic equation corresponding to the linearization of the Eq. (1.2) has multiple pure imaginary roots at specific parameter values, so the classical Hopf bifurcation theory can not be applied.

In this paper we mainly consider the dynamical properties of Eq. (1.2). The remainder of this paper is organized as follows. In section 2, we proved that a series of Hopf bifurcations will occur when the delay τ increases. In section 3, we obtained the existence and spatial pattern of multiple periodic solutions of Eq. (1.2). In section 4, the detailed calculations of the normal form on center manifold of Eq. (1.2) near to Equivariant-Hopf bifurcation points are determined. We also analysis the direction and the stability conditions of bifurcation nonsynchronous periodic solutions. In section 5, some numerical simulation are given to illustrate the results. Finally, we provides a brief conclusion of our results. The Appendix contains some detailed calculation procedures of coefficients h.

2 Hopf bifurcation

Consider the complex delay Eq. (1.2). In order to study the effect of τ on the stability of equilibrium point, we need to analyze the distribution of roots of eigenvalue equation corresponding to the linear part of the system. It is clear that (11+2b,11+2b,11+2b) is an unique positive equilibrium point of Eq.(1.2). Let c=11+2b, then make an equilibrium transformation of Eq. (1.2), we have

The linearization of Eq. (2.1) at origin as follows:

regarding τ as the bifurcating parameter, the associated characteristic equation of Eq. (2.2) takes the form

then we have Δ₁(λ) = 0 or Δ₂(λ) = 0.

3 Multiple periodic solutions

Next, we consider the symmetric characteristic of the Eq. (1.2). we know that the Eq. (1.2) is D₃-equivariant with

the dihedral group D₃ is generated by the cyclic subgroup subgroup of Z₃ acts with the generator ρ and the flip of κ.

Let T=2πω2, P_T represents the set of all continuous T-periodic function x(t) : R → R³. According to maximum norm, P_T is a Banach space. Apply the action of D₃ × S¹ on P_T with

Let SP_T be a subspace of P_T, which consist of all T-periodic solution of Eq. (1.2) when parameter τ = τj(2) , then for each closed subgroup of D₃ × S¹ is called the isotropy group Σ, and Σ = {(r, θ) ∈ D₃ × S¹; (r, θ)x(t) = x(t)}. Under usual non-resonance and transversality conditions, the Σ-fixed-point subspace of SP_T as follows.

See [6], symmetric delay differential equations have a bifurcation of periodic solution whose spatiotemporal symmetry can be completely characterized by Σ. We consider the following subgroups of D₃ × S¹ to describe the symmetry of periodic solution of system Eq. (1.2). Σ₁ = {(κ, 1)}, Σ₂ = {(κ, −1)}, Σ3={(ρ,ei2π3)}, Σ4={(ρ,e−i2π3)}. Since the Eq. (1.2) is cyclic, we choose

We know ν(θ) is the corresponding eigenvectors of Δ₂ (λ) with ±i ω₂. Generalized eigenspace U_±iω2 is

From [27], we have the following theorem.

4 Normal form for equivariant-Hopf bifurcation

In this part, we only consider the case that the Eq. (2.3) has double characteristic values ±iω₂, where the equivariant-Hopf bifurcation occurs. Center manifold theory and normal form method [36, 37] are used to study Hopf bifurcation. Firstly, rescale the time by t→tτ, Eq. (2.1) can be written as

where

Suppose that the Eq. (4.1) undergoes Equivariant-Hopf bifurcation at τ = τj(2) = τ_*. Choosing the phase space C = C([−1, 0]; R³), where for ϕ = (ϕ₁, ϕ₂, ϕ₃)^T ∈ C. Then Eq. (4.1) can be written as

with i = mod(3).

The linearized equation of Eq. (4.1) at zero as follows

where

with δ=011101110. The characteristic equation of Eq. (4.2) at origin is detΔ2(0,λτ)=0, where detΔ₂(0, λ) = 0 is the characteristic equation of the linearizaion of Eq. (2.1). Since Δ₂(0, i ω₂)ν_j = 0, j = 1, 2, the center space at τ = τ_* and in complex coordinates is P = span(φ₁, φ₂, φ₃, φ₄), where

and

Let Φ = (φ₁, φ₂, φ₃, φ₄) and ν̄₁ = ν₂. Note that νjT ν_i = 3, i ≠ j ∈ 1, 2 and νiT ν_i = 0, i ∈ (1, 2). It is easy to check that a basis of the adjoint space of P^* is

with 〈Ψ,Φ〉 = I_4×4 for the adjoint bilinear form on C^* × C define by

with φ ∈ C, ψ ∈ C^*, and

Introducing new parameter variables μ = τ − τ_*, we can rewrite Eq. (4.1) as

where

Let B = (iω₂τ_*, −iω₂τ_*, iω₂τ_*, −iω₂τ_*) and P is the generalized eigenspace associated with B, P^* is the adjoint space of P. Then C can be decomposed as C = P ⊕ Q where Q = (φ ∈ C :< ψ, φ > = 0, for all ψ ∈ P^*). Using the decomposition z_t = Φ x(t) + y(t), then we have

We can decompose Eq. (4.3) as

with x ∈ C⁴, y ∈ Q¹. We will write the Taylor expansion

and we have

Using the idea of Faria [37], we know that Eq.(4.4) can be written as

where f_i(x, y, μ) is homogeneous polynomials of degree j about (x, y, μ) with coefficients in C⁴. Then the normal form of Eq. (1.2) on the center manifold

where g21,g31 will be calculated in the following part of this section.

Fist of all, we get

These are the second-order terms of (μ, x) in Eq. (4.6). From Faria and Hal[36, 37, 38, 39], we have the second-order terms of (μ, x) in the normal form on center manifold as follows:

Here, define M_j to be the operator in Vj5 (C⁴ × Ker_π) with the range in the same space by

where

In particular,

where j ≥ 2, 1 ≤ k ≤ 4, and {e₁, e₂, e₃, e₄} is the canonical basis for C⁴. Therefore, if ∣q∣ = 1, then

and

To compute g31 (x, 0, μ), we first note that from Eq. (4.7), it follows that

since μ x^qe_j ∈∉ Ker(M31), for ∣q∣ = 2, j = 1, 2, 3, 4. For μ = 0, f21 (x, 0, 0) = g21 (x, 0, 0) = 0, a simplified formula is given.

1. Compute Projker(M31)13!f31(x,0,0)=0.

2. Compute Projker(M31)[14(Dxf21)(x,0,0)U21(x.0)].

   In fact, from the e(±iω2τ∗)=c∓iω2bc, we have

   12f21(x,0,0)=−τ∗a¯−1[2(x2+x3)2+iω2c(x32−x22)]a−1[2(x1+x4)2+iω2c(x12−x42)]a¯−1[2(x1+x4)2+iω2c(x12−x42)]a−1[2(x2+x3)2+iω2c(x32−x22)]

   since

   U21(x,0)=U21(x,μ)|μ=0=(M21)−1ProjIm(M21)f21(x,0,0)=−τ∗(M21)−1a¯−1[2(x2+x3)2+iω2c(x32−x22)]a−1[2(x1+x4)2+iω2c(x12−x42)]a¯−1[2(x1+x4)2+iω2c(x12−x42)]a−1[2(x2+x3)2+iω2c(x32−x22)]=−τ∗iω2a¯−1[−23x22−4x2x3+2x32+iω2c(x32+13x22)]a−1[23x12+4x1x4−2x42+iω2c(13x12+x42)]a¯−1[2x12−4x1x4−23x42+iω2c(x12+13x42)]a−1[−2x22+4x2x3+23x32+iω2c(13x32+x22)]

   we have

   Projker(M31)[14(Dxf21)(x,0,0)U21(x.0)]=12Projker(M31)[12(Dxf21)(x,0,0)U21(x.0)]=b11x1x3x4+b12x12x2b¯11x2x3x4+b¯12x22x1b11x1x2x3+b12x32x4b¯11x1x2x4+b¯12x42x3

   where

   b11=−τ∗2K22ω2+iτ∗2(K1−16p)2ω2b12=τ∗22ω2K5−iτ∗2(K3+K4)2ω2K1=16p2(r12−r22)−2ω2cp2r1r2K2=−ω2cp2(r12−r22)−32p2r1r2K3=p(8c2+2ω22)3c2K4=8p2(r12−r22)+8ω2cp2r1r2K5=−4ω2cp2(r12−r22)−16p2r1r2

3. Compute 14Projker(M31)(Dyf21)(x,0,0)U22(x,0)].

   Define h = h(x)(θ) = U22 (x, 0), and write

   h(θ)=h(1)(θ)h(2)(θ)h(3)(θ)=h2000x12+h0200x22+h0020x32+h0002x42+h1100x1x2+h1010x1x3+h1001x1x4+h0110x2x3+h0101x2x4+h0011x3x4

   where h₂₀₀₀, h₀₂₀₀, h₀₀₂₀, h₀₀₀₂, h₁₁₀₀, h₁₀₁₀, h₁₀₀₁, h₀₁₁₀, h₀₁₀₁, h₀₀₁₁ ∈ Q¹. The coefficients of h are determined by M22h(x)=f22(x,0,0), which is equivalent to

   DxhBx−AQ1(h)=(I−π)X0F2(Φx,0)

   Applying the definition of A_Q¹ and π, we obtain

   h˙−DxhBx=Φ(θ)Ψ(0)F2(Φx,0)h˙(0)−Lh=F2(Φx,0)

   where ḣ denotes the derivative of h(θ) relative to θ. Let

   F2(Φx,0)=A2000x12+A0200x22+A0020x32+A0002x42+A1100x1x2+A1010x1x3+A1001x1x4+A0110x2x3+A0101x2x4+A0011x3x4

   where A_ijmn ∈ C²,   0 ≤ i, j, m, n ≤ 2 and i + j + m + n = 2. Comparing the coefficients of x12,x22,x32,x42, x₁x₂, x₁x₃, x₁x₄, x₂x₃, x₂x₄ and x₃x₄, we have h̄₂₀₀₀ = h₀₂₀₀ = h₀₀₂₀ = h̄₀₀₀₂, h₁₀₁₀ = h̄₀₁₀₁ and h₁₁₀₀ = h₀₀₁₁ = h₁₀₀₁ = h₀₁₁₀ = 0, and h₂₀₀₀, h₁₀₁₀ satisfy the following differential equations respectively,

   h2000˙−2iω2τ∗h2000=Φ(θ)Ψ(0)A2000h2000˙(0)−L(h2000)=A2000 (4.8)

   h1010˙−2iω2τ∗h1010=Φ(θ)Ψ(0)A1010h1010˙(0)−L(h1010)=A1010 (4.9)

   Since

   F2(zt,0)=−τ∗z12(0)z22(0)z32(0)−bτ∗z1(0)(z2(−1)+z3(−1))z2(0)(z1(−1)+z3(−1))z3(0)(z1(−1)+z2(−1))

   and

   12f21(x,y,0)=Ψ(0)F2(Φx+y,0)=−τ∗3a¯−1B1a−1B2a¯−1B2a−1B1

   where

   B1=(c4+y1(0))2+b(c4+y1(0))(−c2−c3+y2(−1)+y3(−1))+e−23iπ[(c1+y2(0))2+b(c1+y2(0))(−c2e−23iπ−c3e23iπ+y1(−1)+y3(−1))]+e23iπ[(c¯1+y3(0))2+b(c¯1+y3(0))(−c2e23iπ−c3e−23iπ+y1(−1)+y2(−1))],B2=(c4+y1(0))2+b(c4+y1(0))(−c2−c3+y2(−1)+y3(−1))+e23iπ[(c1+y2(0))2+b(c1+y2(0))(−c2e−23iπ−c3e23iπ+y1(−1)+y3(−1))]+e−23iπ[(c¯1+y3(0))2+b(c¯1+y3(0))(−c2e23iπ−c3e−23iπ+y1(−1)+y2(−1))],

   Thus

   12(Dyf21)(x,y,0)(h)=−τ∗3a¯−1∂B1∂y1a¯−1∂B1∂y2a¯−1∂B1∂y3a−1∂B2∂y1a−1∂B2∂y2a−1∂B2∂y3a¯−1∂B2∂y1a¯−1∂B2∂y2a¯−1∂B2∂y3a−1∂B1∂y1a−1∂B1∂y2a−1∂B1∂y3h1(θ)h2(θ)h3(θ)

   where

   ∂B1∂y1=2(c4+y1(0))+b(−c2−c3+y2(−1)+y3(−1))+be−23iπ(c1+y2(0))+be23iπ(c¯1+y3(0))∂B1∂y2=e−23iπ[2(c1+y2(0))+b(−c2e−23iπ−c3e23iπ+y1(−1)+y3(−1))]+b(c4+y1(0))+be23iπ(c¯1+y3(0))∂B1∂y3=e23iπ[2(c¯1+y3(0))+b(−c2e23iπ−c3e−23iπ+y1(−1)+y2(−1))]b(c4+y1(0))+be−23iπ(c1+y2(0))∂B2∂y1=2(c4+y1(0))+b(−c2−c3+y2(−1)+y3(−1))+be23iπ(c1+y2(0))+be−23iπ(c¯1+y3(0))∂B2∂y2=e23iπ[2(c1+y2(0))+b(−c2e−23iπ−c3e23iπ+y1(−1)+y3(−1))]+b(c4+y1(0))+be−23iπ(c¯1+y3(0))∂B2∂y3=e−23iπ[2(c¯1+y3(0))+b(−c2e23iπ−c3e−23iπ+y1(−1)+y2(−1))]+b(c4+y1(0))+be23iπ(c1+y2(0))

   Thus

   14ProjKerM31(Dyf21)(x,y,0)U22=d11x1x3x4+d12x12x2d¯11x2x3x4+d¯12x22x1d11x1x2x3+d12x32x4d¯11x1x2x4+d¯12x42x3

   where

   d11=−τ∗2p[r1K6+r2K7−i(r2K6−r1K7)]d12=−τ∗2p[r1K8+r2K9−i(r2K8−r1K9)]K6=(2+2b)Re(h1010(1)(0))−Re(h1010(1)(−1))−ω2cIm(h1010(1)(−1))K7=(2+2b)Im(h1010(1)(0))−Im(h1010(1)(−1))+ω2cRe(h1010(1)(−1))K8=(2−b)Re(h2000(1)(0))−Re(h2000(1)(−1))−ω2cIm(h2000(1)(−1))K9=(2−b)Im(h2000(1)(0))−Im(h2000(1)(−1))+ω2cRe(h2000(1)(−1))

   and h_ijmn will be calculated in Appendix. From Cases I, II, III, we get

   16g31(x,0,μ)=(b11+d11)x1x3x4+(b12+d12)x12x2(b¯11+d¯11)x2x3x4+(b¯12+d¯12)x22x1(b11+d11)x1x2x3+(b12+d12)x32x4(b¯11+d¯11)x1x2x4+(b¯12+d¯12)x42x3

   So, we can express Eq. (4.1) as

   x1˙=iω2τ∗x1+μa¯−1iω2x1+(b11+d11)x1x3x4+(b12+d12)x12x2x2˙=−iω2τ∗x2−μa−1iω2x2+(b¯11+d¯11)x2x3x4+(b¯12+d¯12)x22x1x3˙=iω2τ∗x3+μa¯−1iω2x3+(b11+d11)x1x2x3+(b12+d12)x32x4x4˙=−iω2τ∗x4−μa−1iω2x4+(b¯11+d¯11)x1x2x4+(b¯12+d¯12)x42x3 (4.10)

   Since x₁ = x̄₂, x₃ = x̄₄, through the change of variables x₁ = α₁ − i α₂, x₂ = α₁ + i α₂, x₃ = α₃ − i α₄, x₄ = α₃ + i α₄, we obtain

   α˙1α˙2=ω2τ∗α2−α1+μω2−Im[a¯−1]α1+Re[a¯−1]α2−Re[a¯−1]α1−Im[a¯−1]α2+α1(Re[b11+d11]ρ22+Re[b12+d12]ρ12)+α2(Im[b11+d11]ρ22+Im[b12+d12]ρ12)−α1(Im[b11+d11]ρ22+Im[b12+d12]ρ12)+α2(Re[b11+d11]ρ22+Re[b12+d12]ρ12)α˙3α˙4=ω2τ∗α4−α3+μω2−Im[a¯−1]α3+Re[a¯−1]α4−Re[a¯−1]α3−Im[a¯−1]α4+α3(Re[b11+d11]ρ12+Re[b12+d12]ρ22)+α4(Im[b11+d11]ρ12+Im[b12+d12]ρ22)−α3(Im[b11+d11]ρ12+Im[b12+d12]ρ22)+α4(Re[b11+d11]ρ12+Re[b12+d12]ρ22)

   If we use double polar coordinates α₁ = ρ₁cosχ₁, α₂ = ρ₁sinχ₁, and α₃ = ρ₂cosχ₂, α₄ = ρ₂sinχ₂, then we get

   ρ˙1=ρ1(−μω2Im[a¯−1]+Re[b11+d11]ρ22+Re[b12+d12]ρ12)+o(|μ|2|,(ρ1,ρ2)|)+o(|(ρ1,ρ2)|4)ρ˙2=ρ2(−μω2Im[a¯−1]+Re[b11+d11]ρ12+Re[b12+d12]ρ22)+o(|μ|2|,(ρ1,ρ2)|)+o(|(ρ1,ρ2)|4)χ1˙=−ω2τ∗−μω2Re[a¯−1]−Im[b11+d11]ρ22−Im[b12+d12]ρ12+o(|μ|2|,(ρ1,ρ2)|)+o(|(ρ1,ρ2)|4)χ2˙=−ω2τ∗−μω2Re[a¯−1]−Im[b11+d11]ρ12−Im[b12+d12]ρ22+o(|μ|2|,(ρ1,ρ2)|)+o(|(ρ1,ρ2)|4) (4.11)

   Introducing periodic variable parameters ς, and

   z1(t)=α1(s)+iα2(s),z2(t)=α3(s)+iα4(s),s=t(1+ς)ω2τ∗

   we obtain

   (1+ς)z˙1(t)=α2(s)−iα1(s)+μτ∗[−Im(a¯−1)α1(s)+Re(a¯−1)α2(s)−iRe(a¯−1)α1(s)−iIm(a¯−1)α2(s)]+1ω2τ∗{(Re(b11+d11)|z2|2+Re[b12+d12]|z1|2)α1(s)+i(Re(b11+d11)|z2|2+Re(b12+d12)|z1|2)α2(s)+(Im(b11+d11)|z2|2+Im(b12+d12)|z1|2)α2(s)−i(Im(b11+d11)|z2|2+Im(b12+d12)|z1|2)α1(s)}+o(μ2|z|2)+o(|z|4)=−iz1(t)−iμτ0a−1z1(t)+1ω2τ∗(b11+d11)¯|z2(t)|2z1(t)+1ω2τ∗(b12+d12)¯|z1(t)|2z1(t)+o(μ2|z|2)+o(|z|4)

   Similarly, we get an equation for z₂(t). Thus, ignoring the terms o(μ²∣z∣²)+o(∣z∣⁴), we get the normal form

   (1+ς)z˙1(t)=−iz1(t)−iμτ∗a−1z1(t)+1ω2τ∗(b11+d11)¯|z2(t)|2z1(t)+1ω2τ∗(b12+d12)¯|z1(t)|2z1(t)(1+ς)z˙2(t)=−iz2(t)−iμτ∗a−1z2(t)+1ω2τ∗(b11+d11)¯|z1(t)|2z2(t)+1ω2τ∗(b12+d12)¯|z2(t)|2z2(t) (4.12)

   Let g : C ⊕ C ⊕ R → C ⊕ C be given so that −g(z₁, z₂, μ) is the right-hand side of Eq.(4.12), then Eq.(4.12) can be written as

   (1+ς)z˙+g(z,μ)=0 (4.13)

   Note that

   Dzg(0,0)(z1,z2)=i(z1,z2)z=(z1,z2)∈C⊕C

   Also note that g(., μ) : C ⊕ C → C ⊕ C is D₃ × S¹−equivariant with respect to the following D₃ × S¹−action on C ⊕ C:

   γ(z1,z2)=(ei2π3z1,e−i2π3z2)Z3=〈γ〉≤D3κ(z1,z2)=(z2,z1)Z2=〈κ〉≤D3eiθ(z1,z2)=(eiθz1,eiθz2)eiθ∈S1

   According to [28] and [35], the bifurcations of small-amplitude periodic solutions of Eq.(4.13) are completely determined by the signs of three eigenvalues of

   −i(1+ς)z+g(z,μ)=0 (4.14)

   and their orbital stability is determined by the signs of three eigenvalues of

   Dg(z,0)−i(1+ς)Id (4.15)

   that are not forced to zero by the group action. To be more precise, we note that Eq. (4.13) is equivalent to

   −iςz1+iμτ∗a−1z1−1ω2τ∗(b11+d11)¯|z2(t)|2z1(t)−1ω2τ∗(b12+d12)¯|z1(t)|2z1(t)=0−iςz2+iμτ∗a−1z2−1ω2τ∗(b11+d11)¯|z1(t)|2z2(t)−1ω2τ∗(b12+d12)¯|z2(t)|2z2(t)=0 (4.16)

   It is known that Eq. (4.16) can be written as

   Az1z2+Bz12z¯1z22z¯2

   with

   A=A0+AN(|z1|2+|z2|2)B=B0

   for some complex numbers A₀, A_N,B₀ given by

   A0=iμτ∗a−1−iςAN=−1ω2τ∗(b11+d11)¯B0=1ω2τ∗((b11+d11)¯−(b12+d12)¯)

   By the results of [6, 28] and [35], Re(A_n + b₀) > 0 or Re(A_n+B₀) < 0 determines whether the bifurcation of the phase-locked oscillation occurring in the system is supercritical or subcritical. When Re(A_n + B₀) > 0 and Re(B₀) < 0 these are orbitally asymptotically stable. In addition, Re(2A_N + B₀) > 0 or Re(2A_N + B₀) < 0 determines whether the bifurcation of mirror-reflecting waves and standing waves are supercritical or subcritical. When Re(2A_N + B₀) > 0 and Re(B₀) > 0 these are orbitally asymptotically stable.

   Note that

   H1=Re(B0+AN)=Re(−1ω2τ∗(b12+d12)¯)=12ω22[ω2p(r1K8+r2K9)−τ∗K5]H2=Re(2AN+B0)=Re(−1ω2τ∗((b11+d11)¯+(b12+d12)¯))=12ω22[ω2p(r1(K6+K8)+r2(K7+K9))+τ∗(K2−K5)]H3=Re(B0)=Re(1ω2τ∗((b11+d11)¯−(b12+d12)¯))=12ω22[ω2p(r1(K8−K6)+r2(K9−K7))−τ∗(K2+K5)]