# Deterministic and Stochastic Bifurcations of the Catalytic CO Oxidation on Ir(111) Surfaces with Multiple Delays

Authors: Zaitang Huang; Weihua Lei
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

The main purpose is to investigate both deterministic and stochastic bifurcations of the catalytic CO oxidation. Firstly, super- and subcritical bifurcations are determined by the signs of the Poincaré-Lyapunov coefficients of the center manifold scalar bifurcation equations. Secondly, we explore the stochastic bifurcation of the catalytic CO oxidation on Ir(111) surfaces with multiple delays according to the qualitative changes in the invariant measure, the Lyapunov exponent, and the stationary probability density of system response. Some new criteria ensuring stability and stochastic bifurcation are obtained.


---

## Body

## 1. Introduction

For most of the practical cases, the dynamical systems will be disturbed by some stochastic perturbation. There are real benefits  to be gained in using stochastic rather than deterministic models. It could be intuitively thought that the external noise would smear out the finer details of the nonlinear system, such as bifurcations between well-defined states. However, it could also happen that noise forces the system far away from its deterministic attractors to explore regions of phase space that are otherwise not visited or interact with the nonlinearities giving rise to new phenomena [1–10].Now, a great care should be taken when reducing the internal/external sources of noise and uncertainties such as thermal noise, observation errors, and model misspecification. The CO oxidation on platinum group crystals offers an experimental setup, where the nonlinear mechanisms have been captured with high accuracy, and sources of internal noise are controlled. The CO oxidation on crystals shows, in the absence of noise, a wide variety of phenomena, for example, oscillations, bistability, bifurcation, excitability, spatiotemporal patterns, and turbulence. These reactions have been receiving an increasing attention as a laboratory analogy of catalysis used at the industrial level to manufacture chemical products and process harmful species. Understanding the effects of external noise on this oxidation reaction can be considered as an effort to bridge the quality gap. Besides material and pressure gap, which have been discussed in the literature before, the quality gap separates small experiments in perfectly controlled settings and the industrial processes that are always under the influence of uncertainties.In this paper, we investigate the deterministic and stochastic bifurcations of the catalytic CO oxidation on Ir(111) surfaces with multiple delays. Mechanisms and parameters are known with a high degree of accuracy, and the relevant parameters can be controlled. These ideal properties make it possible to go back and forth between theory and experiments. Our paper focuses on homogeneous states and builds on the work of Hayase et al. [11–18] who studied experimentally and numerically the effect of noise on CO oxidation on Ir(111). They verified that the mathematical model was correct and observed probability distributions that are consistent with the presence of multiplicative noise.The structure of the paper is as follows. In Section2, we analyzed the single and double Hopf bifurcation of deterministic system. In Section 3, we analyzed the stochastic stability and bifurcation of stochastic system.


## 2. Deterministic System

In this section, we will investigate the delayed CO oxidation on Ir(111) surfaces as follows:(1)u˙(t)=ακ(1-u-v)-βu-γu(t-τ1)v(t-τ2),v˙(t)=δ(1-κ)(1-u-v)3-γuv,
where γ=γc+γ~ is the bifurcation parameter. It is obvious that there is a unique equilibrium point E(0,1) in system (1).Letx1=u-0,x2=v-1,andX=[x1,x2]⊤, then by substituting the corresponding variables in (1), we have
(2)x˙1(t)=-ακ(x1+x2)-βx1-γx1(t-τ1)[x2(t-τ2)+1],x˙2(t)=-δ(1-κ)(x1+x2)3-γx1(x2+1).Linearized System. As the first step, we analyze the stability of the trivial solution of the linearized system of (1), which is given by
(3)x˙1(t)=-(ακ+β)x1-ακx2-γx1(t-τ1),x˙2(t)=-γx1,
in 𝒞:=𝒞([-τ,0],ℜn), and its adjoint form
(4)x˙1(t^)=(ακ+β)x1(t^)+ακx2(t^)+γx1(t^+τ1),x˙2(t^)=γx1(t^),
in 𝒞^:=𝒞^([-τ,0],ℜn) with respect to the associated bilinear relation
(5)(ψj(ϑ),ϕk(θ))=(ψj(0),ϕk(0))-γ∫-τ10ψj(ζ+τ1)ϕk(ζ)dζj,k=1,2,…,
has eigenvalues λ satisfying the transcendental characteristic equation
(6)Δ(λ,N):=λ2+λ[(ακ+β)+γe-λτ1]-γακ=0.By the Hopf bifurcation, we letλ1,2=±iω(γ), with, ω(γ)≠0 and ℜe{dΔ(λ,γ)/dγ}≠0, be the eigenvalue of (6). All the remaining eigenvalues of Δ(λ,γ)=0 have negative real parts. Then, by substituting λ1=iω(γ) into Δ(λ,γ)=0 and setting the real and imaginary parts of the resulting algebraic equations to zero, we have
(7)-ω2-γακ+ωγsinωτ1+ω[(ακ+β)+ωγcosωτ1]i=0.Separating the real and imaginary parts of (7) gives the following equations:
(8)-ω2-γακ+ωγsinωτ1=0,ω[(ακ+β)+ωγcosωτ1]=0,
or equivalently
(9)sinωτ1=ω2+γακγω>0,cosωτ1=ακ+βγω>0.
Possible values of ωτ1 are of the form ωτ1∈(2kπ,2kπ+π/2), k=0,1,2,…. Thus, we obtain the bifurcation parameter γ as
(10)γ=(4α2κ2ω4-4(α2κ2-ω2)(β2+2αβκ+α2κ2+ω4)-2ακω2c-4α2κ2ω4-4(α2κ2-ω2)(β2+2αβκ+α2κ2+ω4))×(2(α2κ2-ω2))-1(ακ<ω).Now, we compute Hopf’s transversality condition by differentiating implicitly (6) with respect to γ, obtaining
(11){dλdγ}λ=iω=-{(dΔ(λ,γ)dγ)(dΔ(λ,γ)dλ)-1}λ=iω=-{λe-λτ1-ακ2λ+(ακ+β)+γe-λτ1-γτ1e-λτ1}λ=iω=G(γ)+M(γ)i,G(R)=-ac+bda2+b2,M(γ)=-ad-bca2+b2,a=ωsinωτ1-ακ,b=ωcosωτ1,c=ακ+γcosωτ1-γτ1ωsinωτ1,d=2ω-γsinωτ1-γτ1ωcosωτ1.Lemma 1.
The transversality conditionsℜe(dλ/dγ)λ=ωi=-ac-bd≠0 hold. If the choice of values for the model parameters is such that the inequality
(12)(ωsinωτ1-ακ)(ακ+γcosωτ1-γτ1ωsinωτ1)+ωcosωτ1(2ω-γsinωτ1-γτ1ωcosωτ1)<0
holds.This means that the pair of eigenvaluesλ1,2=v(γ)+iω(γ) of Δ(λ,γ)=0 with v(γc)=0, ω(γc)≠0 will cross the imaginary axis from left to right in the complex plane with a nonzero speed as the bifurcation parameter γ varies near its critical value γc. Also, it ensures that the steady-state stability of the system is lost at γ=γc, and the instability in the sense of supercritical and subcritical bifurcations will persist for larger values of γ. The exact nature of such bifurcations satisfying Hopf’s conditions ((6)–(12)) will depend on the nonlinearity in the delay system.Rewriting (10) as a quadratic in ω2 and letting ω2=ω~, we have
(13)ω~2+α1ω~+α2=0,α1=2αγκ-γ2,α2=α2γ2κ2+(ακ+β)2,
whose real and positive solutions are of the form
(14)ω~1=12(-α1+α12-4α2),ω~2=12(-α1-α12-4α2),
if the following holds
(15)α1<0∣2ακ<γ,α12-4α2>0∣γ4-2ακγ3-4(ακ+β)2>0,α1<0∣2ακ<γ,α12-4α2=0∣γ4=2ακγ3+4(ακ+β)2.
In addition, one has to show that the substitution of ω~1 and ω~2 into (10) will ensure Hopf’s transversality conditions, that is, ℜe(dλ/dγ)λ=iω. Therefore, with ω~1=ω2 and ℜe(dλ/dγ)λ=iω, we obtain
(16)2ω2+α1=0∣ω=γ2-2αγκ2,α12-4α2=0∣γ4=2ακγ3+4(ακ+β)2.
By setting ωmin=ω and γmin=γ, we have
(17)ωmin=γmin2-2αγminκ2,γmin4=2ακγmin3+4(ακ+β)2.These are the minimum values separating the regimes of stable and unstable solutions of the linearized delay equation at the Hopf bifurcation pointγ=γc, λ=±ωi(γmin)=±γmin2-2αγminκ/2. Below and equal to these minimum values, all the infinite numbers of eigenvalues of the transcendental characteristic equation (6) will lie in the left-hand side of the complex plane.Next, deriving an expression forτ1 begins by further rewiring (8) as follows:
(18)tanωτ1=ω2+γακακ+β,
thus, we obtain
(19)τ1,k=1ω(arctan(ω2+γακακ+β)+2kπ)k=0,1,2,…,n,…,
and by using the minimum value for ωmin=γmin2-2αγminκ/2, we get
(20)τ1,k,min=1ωmin(arctan(ωmin2+γακακ+β)+2kπ)ck=0,1,2,…,n,….
With these minimum values for the gain ωmin, time delay τ1,k,min, and the response frequency ωmin=γmin2-2αγminκ/2, one can describe all the critical settings of the model parameters in the delay equation (1) at which the unstable-free oscillation persists. Should unstable oscillations appear at some desired frequency such that ωmin=γmin2-2αγminκ/2, fine-tuning the model parameters in the vicinity of their critical settings can control such oscillations. To find out whether the bifurcations are supercritical and/or subcritical, we carry out a local center manifold analysis of the nonlinear delay equation (1) as presented in the preceding section.


### 2.1. The Single Hopf Bifurcation

To obtain the explicit analytical expressions for the stability condition of the Hopf bifurcation solutions (limit cycles), system (1) should be reduced to its center manifold [19–21]. While studying the critical infinite dimensional problem on a two-dimensional center manifold, we begin by defining the elements of the initial continuous function ϕ(θ)∈𝒞 and its projections ϕP(θ),ϕQ(θ) onto the center and stable subspaces P,Q∈𝒞. For the eigenvalues λ=±iω, we have ϕP(θ)=Φ(θ)b with ϕkP(θ)=eλkθ, -τ≤θ≤0 of the adjoint-generalized eigenspace P∈𝒞 and ψP^(s)=Ψb^, ψjP^(s)e-λjs, 0≤s≤τ,j,k=1,2 of the adjoint-generalized eigenspace P^∈𝒞^. The values of the basis Φ(θ)∈𝒞 for P are thus given by
(21)Φ(θ)=[cosωθ-sinωθωsinωθωcosωθ][bb],
and, similarly, the basis Ψ(s)∈𝒞^ for P^ is of the form
(22)Ψ(θ)=[ψ1,ψ2]=[cosωθωsinωθ-sinωθωcosωθ],0≤s≤R.These bases form the inner product matrix(23)(Ψ(s),Φ(θ))=[ψ1(s)ψ2(s)][ϕ1ϕ2]=[(ψ1(s)ϕ1(θ))(ψ1(s)ϕ2(θ))(ψ2(s)ϕ1(θ))(ψ2(s)ϕ2(θ))],
in which
(24)(ψ1(s)ϕ1(θ))=(cosωscosωθ+ω2sinωssinωθ),(ψ1(s)ϕ2(θ))=-(cosωssinωθ-ω2sinωscosωθ),(ψ2(s)ϕ1(θ))=-(sinωscosωθ-ω2cosωssinωθ),(ψ2(s)ϕ2(θ))=(sinωssinωθ+ω2cosωscosωθ).The substitution of the elements of(Ψ,Φ)nsg≠0(25)(Ψ,Φ)nsg=[(ψ1,ϕ1)(ψ1,ϕ2)(ψ2,ϕ1)(ψ2,ϕ2)],
where
(26)(ψ1,ϕ1)=1-12γ{(1ω-ω)sinωτ1+τ1(1+ω2)cosωτ1},(ψ1,ϕ2)=γ(1+ω2)sinωτ1,(ψ2(s)ϕ1(θ))=-γ(1+ω2)sinωτ1,(ψ2(s)ϕ2(θ))=ω2+γ{(1ω-ω)sinωτ1ω2+γcc-τ1(1+ω2)cosωτ1(1ω-ω)(1ω-ω)}.Then, the basisΨ(s) for P^∈𝒞^ of the adjoint equation (2) is normalized to Ψ~=[ψ~1,ψ~2]⊤∈𝒞^, computing Ψ~=(Ψ,Φ)nsg-1Ψ(s) to yield
(27)Ψ~(s)=[ψ~11(s)ψ~12(s)ψ~21(s)ψ~22(s)],0≤s≤R,Ψ~11(s)=1det((Ψ,Φ)nsg-1){(ψ2,ϕ2)cosωs+(ψ1,ϕ2)sinωs},Ψ~12(s)=ωdet((Ψ,Φ)nsg-1){(ψ2,ϕ2)sinωs-(ψ1,ϕ2)cosωs},Ψ~21(s)=ωdet((Ψ,Φ)nsg-1){(ψ2,ϕ1)cosωs+(ψ1,ϕ1)cosωs},Ψ~22(s)=ωdet((Ψ,Φ)nsg-1){(ψ2,ϕ1)sinωs-(ψ1,ϕ1)cosωs},det((Ψ,Φ)nsg-1)={(ψ1,ϕ1)(ψ2,ϕ2)-(ψ1,ϕ2)(ψ2,ϕ1)},
where the substitution of the new elements (ψj(s),ϕk(θ)), j,k=1,2 into (3) will lead to the identity matrix
(28)(Ψ,Φ)id=1det((Ψ,Φ)nsg-1)[det((Ψ,Φ)nsg-1)00det((Ψ,Φ)nsg-1)]=[1001].Consequently, the constant matricesB∈𝒞 and B^∈𝒞^ are equivalent, and the elements of B at the Hopf bifurcation are B=[[0,ω]⊤,[-ω,0]⊤].With the algebraic simplifications(29)eBθ=[cosωθ-sinωθsinωθcosωθ],
one can easily see that
(30)Φ(θ)=Φ(0)eBθ=[cosωθ-sinωθωsinωθωcosωθ],-R≤θ≤0,
and it follows that J(t,μ)Φ(θ)=Φ(0)eB(θ+t) is indeed the solution operator of the linearized delay equations in 𝒞. Similarly, we have Ψ-(s)=Ψ-(0)e-B^s, 0≤s≤R, and J^(t^,R)Ψ-(s)=Φ(0)eB(θ+t) is the solution operator for the adjoint delay equations in 𝒞^.By the change of variablextP(θ)=Φ(θ)z(t)+xtQ(θ) with z(t)∈ℜ2,z(t)=(Ψ~(s),ϕP(θ)), we obtain the following as a first-order approximation in ε, for θ=-τ1,-τ2:
(31)[x1(t-τ1)x2(t-τ1)]=[z1(t)cosωτ+z2(t)sinωτ1-z1(t)ωsinωτ1+z2(t)ωcosωτ1],[x1(t-τ2)x2(t-τ2)]=[z1(t)cosωτ+z2(t)sinωτ2-z1(t)ωsinωτ2+z2(t)ωcosωτ2],
and for θ=0,
(32)[x1(t)x2(t)]=[z1(t)ωz2(t)].
Therefore, we obtain (1) of the center manifold as follows:
(33)z˙1(t)=-ωz2+Δf(1)(z1,z2,γ),z˙2(t)=ωz1+Δf(2)(z1,z2,γ),
where the perturbation functions Δf(1)(z1,z2,γ) and Δf(2)(z1,z2,γ) denote the following:
(34)Δf(1)(z1,z2,γ)=γ~(c10z1+c100z2)+c11z13+c12z12z2+c13z1z22+c14z23+c15z12+c16z1z2+c17z22,Δf(2)(z1,z2,γ)=γ~(c20z1+c200z2)+c21z13+c22z12z2+c23z1z22+c24z23+c25z12+c26z1z2+c27z22.Here, the lengthy expressions ofcij are omitted here, and which in polar coordinate, z1=asinθ, z2=-acosθ, θ=ωt+φ is equivalent to
(35)a˙(t)=ν(γ)a+Π(1)(τ1,τ2,γ)a3+◯(a5,τ1,τ2,γ),φ˙(t)=ω+Π(2)(τ1,τ2,γ)a2+◯(a4,τ1,τ2,γ).
The so-called Poincaré-Lyapunov coefficient Π(1)(τ1,τ2,γ) is independently determined by the following formula:
(36)Π(1)(a,τ1,τ2,γ)=116{Δfz1z1z1(1)(z1,z2,γ)+Δfz1z2z2(1)(z1,z2,γ)+Δfz1z1z2(2)(z1,z2,γ)+Δfz2z2z2(2)(z1,z2,γ)}+116ω{Δfz1z2(1)(z1,z2,γ)×(Δfz1z1z1(1)(z1,z2,γ)+Δfz2z2(1)(z1,z2,γ))-Δfz1z2(2)(z1,z2,γ)×(Δfz1z1(2)(z1,z2,γ)+Δfz1z2(2)(z1,z2,γ))-Δfz1z1(1)(z1,z2,γ)Δfz1z1(2)(z1,z2,γ)+Δfz2z2(1)(z1,z2,γ)Δfz2z2(2)(z1,z2,γ)}.
As is well known [19–24], the amplitude equation (35) together with (36) can determine the stability behaviour of the center manifold ODEs (33), where all the partial derivatives of Δf(1)(z1,z2,γ), Δf(2)(z1,z2,γ), for example, Δfz1z2(2)={∂(∂Δfz1z2(2)/∂z1)/∂z2}(0,0,γc), are evaluated at the bifurcation point (z1,z2)→(0,0). Thus, the partial derivatives are given as follows:
(37)∂Δf(1)(z1,z2,γ)∂z1=c10γ~+3c11z12+2(c15+c12z2)z1+(c16+c13z2)z2,∂Δf(1)(z1,z2,γ)∂z2=c100γ~+c13z12+(c16+2c13z2)z1+(2c17+3c14z2)z2,(38)∂Δf(2)(z1,z2,γ)∂z1=c20γ~+3c21z12+2(c25+c22z2)z1+(c26+c23z2)z2,∂Δf(2)(z1,z2,γ)∂z2=c200γ~+c23z12+(c26+2c23z2)z1+(2c27+3c24z2)z2,
from which we readily obtain the required values for the formula in (36) as
(39)∂Δ3f(1)(z1,z2,γ)∂z1∂z1∂z1=6c11,∂Δ3f(1)(z1,z2,γ)∂z1∂z2∂z2=2c13,∂Δ2f(1)(z1,z2,γ)∂z1∂z2=2(c12z1+c13z2)+c16,∂Δf(1)(z1,z1,γ)∂z1∂z1=2(3c11z1+c12z2+c15),∂Δf(1)(z2,z2,γ)∂z2∂z2=2(c13z1+3c14z2+c17),∂Δ3f(2)(z1,z2,γ)∂z1∂z1∂z1=6c24,∂Δ3f(2)(z1,z2,γ)∂z1∂z2∂z2=2c23,∂Δ2f(2)(z1,z2,γ)∂z1∂z2=2(c22z1+c23z2)+c26,∂Δf(2)(z1,z1,γ)∂z1∂z1=2(3c21z1+c22z2+c25),∂Δf(2)(z1,z2,γ)∂z2∂z2=2(c23z1+3c24z2+c27).
Substituting these quantities into the formulate in (36) will yield
(40)Π(1)(a,τ1,τ2,γ)=18{3(c11+c24)+c13+c23}+18ω{c16(c15+c17)-c26(c25+c27)&18ωccc-c15c25+c17c27(c15+c17)}.
This coefficient together with the amplitude bifurcation equation (35) produces the nonlinear bifurcations for multiple delays as shown in the following: (a) Π(1)(a,τ1,τ2,γ)<0→ supercritical bifurcation; (b) Π(1)(a,τ1,τ2,γ)>0→ subcritical bifurcation.Theorem 2.
Letac+bd>0.(a)
IfΠ(1)(a,τ1,τ2,γ)<0, then system (1) occurs supercritical bifurcation.(b)
IfΠ(1)(a,τ1,τ2,γ)>0, then system (1) occurs subcritical bifurcation.


### 2.2. The Double Hopf Bifurcation

For the eigenvalueλ1,2=±iω1, λ1,2=±iω2 of the equation Δ(λ,γ)=0 in (6), we have the base function Φ(θ)=[ϕ1(θ),ϕ2(θ),ϕ3(θ),ϕ4(θ)], -τ≤θ≤0, Ψ(s)=[ψ1(s),ψ2(s),ψ3(s),ψ4(s)]⊤, 0≤s≤τ in 𝒞([-τ,0],ℜn) and 𝒞([-τ,0],ℜn), accordingly,
(41)Φ(θ)=[cosω1θ-sinωθ1cosω2θ-sinω2θsinω1θcosω1θsinω2θcosω2θ],Ψ(s)=[cosω1s-sinωscosω2s-sinω2ssinω1scosω1ssinω2scosω2s]⊤,
and they form the inner product matrix (Ψ(s),Φ(θ))=(ψj(s)ϕk(θ)),j,k=1,2,3,4, namely,
(42)(Ψ(s),Φ(θ))=[(ψ1(s)ϕ1(θ))(ψ1(s)ϕ2(θ))(ψ1(s)ϕ3(θ))(ψ1(s)ϕ4(θ))(ψ2(s)ϕ1(θ))(ψ2(s)ϕ2(θ))(ψ2(s)ϕ3(θ))(ψ2(s)ϕ4(θ))(ψ3(s)ϕ1(θ))(ψ3(s)ϕ2(θ))(ψ3(s)ϕ3(θ))(ψ3(s)ϕ4(θ))(ψ4(s)ϕ1(θ))(ψ4(s)ϕ2(θ))(ψ4(s)ϕ3(θ))(ψ4(s)ϕ4(θ))].
The substitution of the elements (Ψ(s),Φ(θ))=(ψj(s)ϕk(θ)) into the bilinear pairing (5) and integrating, we obtain the nonsingular matrix
(43)(Ψ,Φ)nsg=[ψ11-ψ1200ψ21ψ220000ψ23-ψ3400ψ43ψ44],
where the constant values of this matrix are given by
(44)ψ11=ψ22=1-γτ1cosω1τ1,ψ12=-γτ1sinω1τ1,ψ12=-ψ21,ψ33=ψ44=1-γτ1cosω2τ1,ψ34=-γτ1sinω2τ1,ψ43=-ψ34.
Then, Ψ∈𝒞^ is normalized to the new basis Ψ~=[ψ~1(s),ψ~3(s),ψ~3(s),ψ~4(s)], 0≤s≤τ, where the elements ψ~j(s), j=1,2,3,4 of Ψ∈𝒞^ are given by
(45)ψ~1=1ψ112+ψ122{(ψ222sinω1s-ψ212cosω1s),0,0]⊤[ψ222(ψ222cosω1s+ψ212sinω1s),1ψ112+ψ122ccc(ψ222sinω1s-ψ212cosω1s),0,0]⊤},ψ~2=1ψ112+ψ122{(ψ222sinω1s-ψ212cosω1s),0,0]⊤[ψ222(ψ122cosω1s-ψ112sinω1s),1ψ112+ψ122ccc(ψ122sinω1s+ψ112cosω1s),0,0]⊤},ψ~3=1ψ332+ψ342{(ψ442sinω2s-ψ432cosω2s)]⊤[ψ2220,0,(ψ442cosω2s+ψ432sinω2s),1ψ332+ψ342ccc(ψ442sinω2s-ψ432cosω2s)]⊤},ψ~4=1ψ112+ψ122{(ψ442sinω2s-ψ432cosω2s)]⊤[ψ2220,0,(ψ342cosω2s-ψ332sinω2s),1ψ112+ψ122ccc(ψ342sinω2s+ψ332cosω1s)]⊤}.
Again, the substitution of the elements ψ~j(s), ϕk(θ), j,k=1,2,3,4 of (Ψ~(s),Φ(θ)) will yield the identity matrix, namely,
(46)(Ψ,Φ)nsg=1Λ[Λ0000Λ0000Λ0000Λ]=I,Λ=(ψ112+ψ122)(ψ332+ψ342).
Then, the constant matrix B is given by
(47)B=[0-ω100ω1000000-ω200ω20].
The change of variable xtP(ϕ(θ),γ)=Φ(θ)z(t)+xtQ(ϕ(θ),γ), z(t)∈ℜ4,z(t)=(Ψ~(θ),ϕP(θ)) will yield
(48)Φ(θ)z(t)=[cosω1θ-sinω1θcosω2θ-sinω2θsinω1θcosω1θsinω2θcosω2θ]×[z1(t)z2(t)z3(t)z4(t)],
thereby, for -τ≤θ<0 and τ1≤τ≤τ2, we have
(49)Φ(-τ1)z(t)=[cosω1τ1sinω1τ1cosω2τ1sinω2τ1-sinω1τ1cosω1τ1sinω2τ1cosω2τ1]×[z1(t)z2(t)z3(t)z4(t)],Φ(-τ2)z(t)=[cosω1τ2sinω1τ2cosω2τ2sinω2τ2-sinω1τ2cosω1τ2-sinω2τ2cosω2τ2]×[z1(t)z2(t)z3(t)z4(t)],Φ(0)z(t)=[10100101][z1(t)z2(t)z3(t)z4(t)],
and so
(50)γx1(t-τ1)=γ[z1cosω1τ1+z2sinω1τ1γc+z3cosω2τ1+z4sinω2τ1],γx1(t-τ1)x2(t-τ2)=γ[z1cosω1τ1+z2sinω1τ1yc+z3cosω2τ1+z4sinω2τ1]×[z1cosω1τ2+z2sinω1τ2ccc+z3cosω2τ2+z4sinω2τ2].
The ODEs of the center manifold Mγ∈𝒞([-τ,0],ℜ4) are given by as follows:
(51)z˙1(t)=-ω1z2+Δg1(z1,z2,z3,z4,γ),z˙2(t)=ω1z1+Δg2(z1,z2,z3,z4,γ),z˙3(t)=-ω2z2+Δg3(z1,z2,z3,z4,γ),z˙4(t)=ω2z2+Δg4(z1,z2,z3,z4,γ),
where the perturbations in these equations denote
(52)Δg1(z1,z2,z3,z4,γ)=c111(1)z13+c112(1)z12z2+c113(1)z12z3+c114(1)z12z2+c115(1)z1z22+c116(1)z1z32+c117(1)z1z42+c222(1)z23+c333(1)z33+c444(1)z43+d111(1)z12+d112(1)z1z2+d113(1)z1z3+d114(1)z1z4+d222(1)z22+d333(1)z32+d444(1)z42+⋯,Δg2(z1,z2,z3,z4,γ)=c111(2)z13+c112(2)z12z2+c113(2)z12z3+c114(2)z12z2+c115(2)z1z22+c116(2)z1z32+c117(2)z1z42+c222(2)z23+c333(1)z33+c444(1)z43+d111(2)z12+d112(2)z2z2+d113(2)z1z3+d114(2)z1z4+d222(2)z22+d333(2)z32+d444(2)z42+⋯,Δg3(z1,z2,z3,z4,γ)=c111(3)z13+c112(3)z12z2+c113(3)z12z3+c114(3)z12z2+c115(1)z1z22+c116(3)z1z32+c117(3)z1z42+c222(3)z23+c333(3)z33+c444(1)z43+d111(3)z12+d112(3)z1z2+d113(3)z1z3+d114(3)z1z4+d222(3)z22+d333(3)z32+d444(3)z42+⋯,Δg4(z1,z2,z3,z4,γ)=c111(4)z13+c112(4)z12z2+c113(4)z12z3+c114(4)z12z2+c115(4)z1z22+c116(4)z1z32+c117(4)z1z42+c222(4)z23+c333(4)z33+c444(4)z43+d111(4)z12+d112(4)z1z2+d113(4)z1z3+d114(4)z1z4+d222(4)z22+d333(4)z32+d444(4)z42+⋯,
and the above coefficients in these equations denote the constant values in (52). Here, the lengthy expressions of cijlk and dijlk are omitted here. These equations are further written into the desired standard form of amplitude and phase relations using the polar coordinates zj=ajsinθj, zj=ajcosθj, θj=ωjt+φj, j=1,2,3,4. Again, the application of the integral averaging method to the resulting amplitude and phase equations will yield
(53)da-1dt=β111(11)a1+Π(11)(a,γ)a13+β112(11)a1a22+⋯,dφ-1dt=β111(12)a1+Π(12)(a,γ)a13+β112(12)a1a22+⋯,da-2dt=β111(21)a2+Π(21)(a,γ)a23+β112(21)a1a22+⋯,dφ-2dt=β111(22)a1+Π(22)(a,γ)a23+β112(22)a1a22+⋯,
where the coefficients in these equations are the resulting constant values after averaging. Similarly, the coefficients in these averaged equations, in particular the cubic ones Π(11)(a,γ), Π(21)(a,γ), will determine whether the double Hopf interactions are of the form (i) subsupercritical bifurcation type, (ii) supersubcritical bifurcation type, (iii) subsubcritical bifurcation type, and (iv) supersupercritical bifurcation type. Like the Poincaré-Lyapunov coefficient Π(a,γ) in (36) for the single Hopf, there are also explicit formulas for computing the double Hopf coefficients Π(11)(a,γ), Π(21)(a,γ), which are extremely too long to be displayed here. For the account of these formulas for the double Hopf coefficients, again, the publications [22, 23] are excellent references. Classical studies of dynamical systems resulting to equations of the form (53) have demonstrated unprecedented dynamics ranging from global dynamic bifurcations to ultimately chaotic dynamics (see [22, 24]). There exist reliable and robust mathematical tools to examine the long-term behaviour of such dynamics. Campbell et al. [24], in particular, have shown numerically that equations of the form (53) exhibit four double Hopf’s interactions, namely, two super-super-, one super-sub-, and one sub-sub-critical bifurcations. Other interesting dynamics such as large limit cycles, tori, and period doubling were also characterized by the authors.


## 2.1. The Single Hopf Bifurcation

To obtain the explicit analytical expressions for the stability condition of the Hopf bifurcation solutions (limit cycles), system (1) should be reduced to its center manifold [19–21]. While studying the critical infinite dimensional problem on a two-dimensional center manifold, we begin by defining the elements of the initial continuous function ϕ(θ)∈𝒞 and its projections ϕP(θ),ϕQ(θ) onto the center and stable subspaces P,Q∈𝒞. For the eigenvalues λ=±iω, we have ϕP(θ)=Φ(θ)b with ϕkP(θ)=eλkθ, -τ≤θ≤0 of the adjoint-generalized eigenspace P∈𝒞 and ψP^(s)=Ψb^, ψjP^(s)e-λjs, 0≤s≤τ,j,k=1,2 of the adjoint-generalized eigenspace P^∈𝒞^. The values of the basis Φ(θ)∈𝒞 for P are thus given by
(21)Φ(θ)=[cosωθ-sinωθωsinωθωcosωθ][bb],
and, similarly, the basis Ψ(s)∈𝒞^ for P^ is of the form
(22)Ψ(θ)=[ψ1,ψ2]=[cosωθωsinωθ-sinωθωcosωθ],0≤s≤R.These bases form the inner product matrix(23)(Ψ(s),Φ(θ))=[ψ1(s)ψ2(s)][ϕ1ϕ2]=[(ψ1(s)ϕ1(θ))(ψ1(s)ϕ2(θ))(ψ2(s)ϕ1(θ))(ψ2(s)ϕ2(θ))],
in which
(24)(ψ1(s)ϕ1(θ))=(cosωscosωθ+ω2sinωssinωθ),(ψ1(s)ϕ2(θ))=-(cosωssinωθ-ω2sinωscosωθ),(ψ2(s)ϕ1(θ))=-(sinωscosωθ-ω2cosωssinωθ),(ψ2(s)ϕ2(θ))=(sinωssinωθ+ω2cosωscosωθ).The substitution of the elements of(Ψ,Φ)nsg≠0(25)(Ψ,Φ)nsg=[(ψ1,ϕ1)(ψ1,ϕ2)(ψ2,ϕ1)(ψ2,ϕ2)],
where
(26)(ψ1,ϕ1)=1-12γ{(1ω-ω)sinωτ1+τ1(1+ω2)cosωτ1},(ψ1,ϕ2)=γ(1+ω2)sinωτ1,(ψ2(s)ϕ1(θ))=-γ(1+ω2)sinωτ1,(ψ2(s)ϕ2(θ))=ω2+γ{(1ω-ω)sinωτ1ω2+γcc-τ1(1+ω2)cosωτ1(1ω-ω)(1ω-ω)}.Then, the basisΨ(s) for P^∈𝒞^ of the adjoint equation (2) is normalized to Ψ~=[ψ~1,ψ~2]⊤∈𝒞^, computing Ψ~=(Ψ,Φ)nsg-1Ψ(s) to yield
(27)Ψ~(s)=[ψ~11(s)ψ~12(s)ψ~21(s)ψ~22(s)],0≤s≤R,Ψ~11(s)=1det((Ψ,Φ)nsg-1){(ψ2,ϕ2)cosωs+(ψ1,ϕ2)sinωs},Ψ~12(s)=ωdet((Ψ,Φ)nsg-1){(ψ2,ϕ2)sinωs-(ψ1,ϕ2)cosωs},Ψ~21(s)=ωdet((Ψ,Φ)nsg-1){(ψ2,ϕ1)cosωs+(ψ1,ϕ1)cosωs},Ψ~22(s)=ωdet((Ψ,Φ)nsg-1){(ψ2,ϕ1)sinωs-(ψ1,ϕ1)cosωs},det((Ψ,Φ)nsg-1)={(ψ1,ϕ1)(ψ2,ϕ2)-(ψ1,ϕ2)(ψ2,ϕ1)},
where the substitution of the new elements (ψj(s),ϕk(θ)), j,k=1,2 into (3) will lead to the identity matrix
(28)(Ψ,Φ)id=1det((Ψ,Φ)nsg-1)[det((Ψ,Φ)nsg-1)00det((Ψ,Φ)nsg-1)]=[1001].Consequently, the constant matricesB∈𝒞 and B^∈𝒞^ are equivalent, and the elements of B at the Hopf bifurcation are B=[[0,ω]⊤,[-ω,0]⊤].With the algebraic simplifications(29)eBθ=[cosωθ-sinωθsinωθcosωθ],
one can easily see that
(30)Φ(θ)=Φ(0)eBθ=[cosωθ-sinωθωsinωθωcosωθ],-R≤θ≤0,
and it follows that J(t,μ)Φ(θ)=Φ(0)eB(θ+t) is indeed the solution operator of the linearized delay equations in 𝒞. Similarly, we have Ψ-(s)=Ψ-(0)e-B^s, 0≤s≤R, and J^(t^,R)Ψ-(s)=Φ(0)eB(θ+t) is the solution operator for the adjoint delay equations in 𝒞^.By the change of variablextP(θ)=Φ(θ)z(t)+xtQ(θ) with z(t)∈ℜ2,z(t)=(Ψ~(s),ϕP(θ)), we obtain the following as a first-order approximation in ε, for θ=-τ1,-τ2:
(31)[x1(t-τ1)x2(t-τ1)]=[z1(t)cosωτ+z2(t)sinωτ1-z1(t)ωsinωτ1+z2(t)ωcosωτ1],[x1(t-τ2)x2(t-τ2)]=[z1(t)cosωτ+z2(t)sinωτ2-z1(t)ωsinωτ2+z2(t)ωcosωτ2],
and for θ=0,
(32)[x1(t)x2(t)]=[z1(t)ωz2(t)].
Therefore, we obtain (1) of the center manifold as follows:
(33)z˙1(t)=-ωz2+Δf(1)(z1,z2,γ),z˙2(t)=ωz1+Δf(2)(z1,z2,γ),
where the perturbation functions Δf(1)(z1,z2,γ) and Δf(2)(z1,z2,γ) denote the following:
(34)Δf(1)(z1,z2,γ)=γ~(c10z1+c100z2)+c11z13+c12z12z2+c13z1z22+c14z23+c15z12+c16z1z2+c17z22,Δf(2)(z1,z2,γ)=γ~(c20z1+c200z2)+c21z13+c22z12z2+c23z1z22+c24z23+c25z12+c26z1z2+c27z22.Here, the lengthy expressions ofcij are omitted here, and which in polar coordinate, z1=asinθ, z2=-acosθ, θ=ωt+φ is equivalent to
(35)a˙(t)=ν(γ)a+Π(1)(τ1,τ2,γ)a3+◯(a5,τ1,τ2,γ),φ˙(t)=ω+Π(2)(τ1,τ2,γ)a2+◯(a4,τ1,τ2,γ).
The so-called Poincaré-Lyapunov coefficient Π(1)(τ1,τ2,γ) is independently determined by the following formula:
(36)Π(1)(a,τ1,τ2,γ)=116{Δfz1z1z1(1)(z1,z2,γ)+Δfz1z2z2(1)(z1,z2,γ)+Δfz1z1z2(2)(z1,z2,γ)+Δfz2z2z2(2)(z1,z2,γ)}+116ω{Δfz1z2(1)(z1,z2,γ)×(Δfz1z1z1(1)(z1,z2,γ)+Δfz2z2(1)(z1,z2,γ))-Δfz1z2(2)(z1,z2,γ)×(Δfz1z1(2)(z1,z2,γ)+Δfz1z2(2)(z1,z2,γ))-Δfz1z1(1)(z1,z2,γ)Δfz1z1(2)(z1,z2,γ)+Δfz2z2(1)(z1,z2,γ)Δfz2z2(2)(z1,z2,γ)}.
As is well known [19–24], the amplitude equation (35) together with (36) can determine the stability behaviour of the center manifold ODEs (33), where all the partial derivatives of Δf(1)(z1,z2,γ), Δf(2)(z1,z2,γ), for example, Δfz1z2(2)={∂(∂Δfz1z2(2)/∂z1)/∂z2}(0,0,γc), are evaluated at the bifurcation point (z1,z2)→(0,0). Thus, the partial derivatives are given as follows:
(37)∂Δf(1)(z1,z2,γ)∂z1=c10γ~+3c11z12+2(c15+c12z2)z1+(c16+c13z2)z2,∂Δf(1)(z1,z2,γ)∂z2=c100γ~+c13z12+(c16+2c13z2)z1+(2c17+3c14z2)z2,(38)∂Δf(2)(z1,z2,γ)∂z1=c20γ~+3c21z12+2(c25+c22z2)z1+(c26+c23z2)z2,∂Δf(2)(z1,z2,γ)∂z2=c200γ~+c23z12+(c26+2c23z2)z1+(2c27+3c24z2)z2,
from which we readily obtain the required values for the formula in (36) as
(39)∂Δ3f(1)(z1,z2,γ)∂z1∂z1∂z1=6c11,∂Δ3f(1)(z1,z2,γ)∂z1∂z2∂z2=2c13,∂Δ2f(1)(z1,z2,γ)∂z1∂z2=2(c12z1+c13z2)+c16,∂Δf(1)(z1,z1,γ)∂z1∂z1=2(3c11z1+c12z2+c15),∂Δf(1)(z2,z2,γ)∂z2∂z2=2(c13z1+3c14z2+c17),∂Δ3f(2)(z1,z2,γ)∂z1∂z1∂z1=6c24,∂Δ3f(2)(z1,z2,γ)∂z1∂z2∂z2=2c23,∂Δ2f(2)(z1,z2,γ)∂z1∂z2=2(c22z1+c23z2)+c26,∂Δf(2)(z1,z1,γ)∂z1∂z1=2(3c21z1+c22z2+c25),∂Δf(2)(z1,z2,γ)∂z2∂z2=2(c23z1+3c24z2+c27).
Substituting these quantities into the formulate in (36) will yield
(40)Π(1)(a,τ1,τ2,γ)=18{3(c11+c24)+c13+c23}+18ω{c16(c15+c17)-c26(c25+c27)&18ωccc-c15c25+c17c27(c15+c17)}.
This coefficient together with the amplitude bifurcation equation (35) produces the nonlinear bifurcations for multiple delays as shown in the following: (a) Π(1)(a,τ1,τ2,γ)<0→ supercritical bifurcation; (b) Π(1)(a,τ1,τ2,γ)>0→ subcritical bifurcation.Theorem 2.
Letac+bd>0.(a)
IfΠ(1)(a,τ1,τ2,γ)<0, then system (1) occurs supercritical bifurcation.(b)
IfΠ(1)(a,τ1,τ2,γ)>0, then system (1) occurs subcritical bifurcation.


## 2.2. The Double Hopf Bifurcation

For the eigenvalueλ1,2=±iω1, λ1,2=±iω2 of the equation Δ(λ,γ)=0 in (6), we have the base function Φ(θ)=[ϕ1(θ),ϕ2(θ),ϕ3(θ),ϕ4(θ)], -τ≤θ≤0, Ψ(s)=[ψ1(s),ψ2(s),ψ3(s),ψ4(s)]⊤, 0≤s≤τ in 𝒞([-τ,0],ℜn) and 𝒞([-τ,0],ℜn), accordingly,
(41)Φ(θ)=[cosω1θ-sinωθ1cosω2θ-sinω2θsinω1θcosω1θsinω2θcosω2θ],Ψ(s)=[cosω1s-sinωscosω2s-sinω2ssinω1scosω1ssinω2scosω2s]⊤,
and they form the inner product matrix (Ψ(s),Φ(θ))=(ψj(s)ϕk(θ)),j,k=1,2,3,4, namely,
(42)(Ψ(s),Φ(θ))=[(ψ1(s)ϕ1(θ))(ψ1(s)ϕ2(θ))(ψ1(s)ϕ3(θ))(ψ1(s)ϕ4(θ))(ψ2(s)ϕ1(θ))(ψ2(s)ϕ2(θ))(ψ2(s)ϕ3(θ))(ψ2(s)ϕ4(θ))(ψ3(s)ϕ1(θ))(ψ3(s)ϕ2(θ))(ψ3(s)ϕ3(θ))(ψ3(s)ϕ4(θ))(ψ4(s)ϕ1(θ))(ψ4(s)ϕ2(θ))(ψ4(s)ϕ3(θ))(ψ4(s)ϕ4(θ))].
The substitution of the elements (Ψ(s),Φ(θ))=(ψj(s)ϕk(θ)) into the bilinear pairing (5) and integrating, we obtain the nonsingular matrix
(43)(Ψ,Φ)nsg=[ψ11-ψ1200ψ21ψ220000ψ23-ψ3400ψ43ψ44],
where the constant values of this matrix are given by
(44)ψ11=ψ22=1-γτ1cosω1τ1,ψ12=-γτ1sinω1τ1,ψ12=-ψ21,ψ33=ψ44=1-γτ1cosω2τ1,ψ34=-γτ1sinω2τ1,ψ43=-ψ34.
Then, Ψ∈𝒞^ is normalized to the new basis Ψ~=[ψ~1(s),ψ~3(s),ψ~3(s),ψ~4(s)], 0≤s≤τ, where the elements ψ~j(s), j=1,2,3,4 of Ψ∈𝒞^ are given by
(45)ψ~1=1ψ112+ψ122{(ψ222sinω1s-ψ212cosω1s),0,0]⊤[ψ222(ψ222cosω1s+ψ212sinω1s),1ψ112+ψ122ccc(ψ222sinω1s-ψ212cosω1s),0,0]⊤},ψ~2=1ψ112+ψ122{(ψ222sinω1s-ψ212cosω1s),0,0]⊤[ψ222(ψ122cosω1s-ψ112sinω1s),1ψ112+ψ122ccc(ψ122sinω1s+ψ112cosω1s),0,0]⊤},ψ~3=1ψ332+ψ342{(ψ442sinω2s-ψ432cosω2s)]⊤[ψ2220,0,(ψ442cosω2s+ψ432sinω2s),1ψ332+ψ342ccc(ψ442sinω2s-ψ432cosω2s)]⊤},ψ~4=1ψ112+ψ122{(ψ442sinω2s-ψ432cosω2s)]⊤[ψ2220,0,(ψ342cosω2s-ψ332sinω2s),1ψ112+ψ122ccc(ψ342sinω2s+ψ332cosω1s)]⊤}.
Again, the substitution of the elements ψ~j(s), ϕk(θ), j,k=1,2,3,4 of (Ψ~(s),Φ(θ)) will yield the identity matrix, namely,
(46)(Ψ,Φ)nsg=1Λ[Λ0000Λ0000Λ0000Λ]=I,Λ=(ψ112+ψ122)(ψ332+ψ342).
Then, the constant matrix B is given by
(47)B=[0-ω100ω1000000-ω200ω20].
The change of variable xtP(ϕ(θ),γ)=Φ(θ)z(t)+xtQ(ϕ(θ),γ), z(t)∈ℜ4,z(t)=(Ψ~(θ),ϕP(θ)) will yield
(48)Φ(θ)z(t)=[cosω1θ-sinω1θcosω2θ-sinω2θsinω1θcosω1θsinω2θcosω2θ]×[z1(t)z2(t)z3(t)z4(t)],
thereby, for -τ≤θ<0 and τ1≤τ≤τ2, we have
(49)Φ(-τ1)z(t)=[cosω1τ1sinω1τ1cosω2τ1sinω2τ1-sinω1τ1cosω1τ1sinω2τ1cosω2τ1]×[z1(t)z2(t)z3(t)z4(t)],Φ(-τ2)z(t)=[cosω1τ2sinω1τ2cosω2τ2sinω2τ2-sinω1τ2cosω1τ2-sinω2τ2cosω2τ2]×[z1(t)z2(t)z3(t)z4(t)],Φ(0)z(t)=[10100101][z1(t)z2(t)z3(t)z4(t)],
and so
(50)γx1(t-τ1)=γ[z1cosω1τ1+z2sinω1τ1γc+z3cosω2τ1+z4sinω2τ1],γx1(t-τ1)x2(t-τ2)=γ[z1cosω1τ1+z2sinω1τ1yc+z3cosω2τ1+z4sinω2τ1]×[z1cosω1τ2+z2sinω1τ2ccc+z3cosω2τ2+z4sinω2τ2].
The ODEs of the center manifold Mγ∈𝒞([-τ,0],ℜ4) are given by as follows:
(51)z˙1(t)=-ω1z2+Δg1(z1,z2,z3,z4,γ),z˙2(t)=ω1z1+Δg2(z1,z2,z3,z4,γ),z˙3(t)=-ω2z2+Δg3(z1,z2,z3,z4,γ),z˙4(t)=ω2z2+Δg4(z1,z2,z3,z4,γ),
where the perturbations in these equations denote
(52)Δg1(z1,z2,z3,z4,γ)=c111(1)z13+c112(1)z12z2+c113(1)z12z3+c114(1)z12z2+c115(1)z1z22+c116(1)z1z32+c117(1)z1z42+c222(1)z23+c333(1)z33+c444(1)z43+d111(1)z12+d112(1)z1z2+d113(1)z1z3+d114(1)z1z4+d222(1)z22+d333(1)z32+d444(1)z42+⋯,Δg2(z1,z2,z3,z4,γ)=c111(2)z13+c112(2)z12z2+c113(2)z12z3+c114(2)z12z2+c115(2)z1z22+c116(2)z1z32+c117(2)z1z42+c222(2)z23+c333(1)z33+c444(1)z43+d111(2)z12+d112(2)z2z2+d113(2)z1z3+d114(2)z1z4+d222(2)z22+d333(2)z32+d444(2)z42+⋯,Δg3(z1,z2,z3,z4,γ)=c111(3)z13+c112(3)z12z2+c113(3)z12z3+c114(3)z12z2+c115(1)z1z22+c116(3)z1z32+c117(3)z1z42+c222(3)z23+c333(3)z33+c444(1)z43+d111(3)z12+d112(3)z1z2+d113(3)z1z3+d114(3)z1z4+d222(3)z22+d333(3)z32+d444(3)z42+⋯,Δg4(z1,z2,z3,z4,γ)=c111(4)z13+c112(4)z12z2+c113(4)z12z3+c114(4)z12z2+c115(4)z1z22+c116(4)z1z32+c117(4)z1z42+c222(4)z23+c333(4)z33+c444(4)z43+d111(4)z12+d112(4)z1z2+d113(4)z1z3+d114(4)z1z4+d222(4)z22+d333(4)z32+d444(4)z42+⋯,
and the above coefficients in these equations denote the constant values in (52). Here, the lengthy expressions of cijlk and dijlk are omitted here. These equations are further written into the desired standard form of amplitude and phase relations using the polar coordinates zj=ajsinθj, zj=ajcosθj, θj=ωjt+φj, j=1,2,3,4. Again, the application of the integral averaging method to the resulting amplitude and phase equations will yield
(53)da-1dt=β111(11)a1+Π(11)(a,γ)a13+β112(11)a1a22+⋯,dφ-1dt=β111(12)a1+Π(12)(a,γ)a13+β112(12)a1a22+⋯,da-2dt=β111(21)a2+Π(21)(a,γ)a23+β112(21)a1a22+⋯,dφ-2dt=β111(22)a1+Π(22)(a,γ)a23+β112(22)a1a22+⋯,
where the coefficients in these equations are the resulting constant values after averaging. Similarly, the coefficients in these averaged equations, in particular the cubic ones Π(11)(a,γ), Π(21)(a,γ), will determine whether the double Hopf interactions are of the form (i) subsupercritical bifurcation type, (ii) supersubcritical bifurcation type, (iii) subsubcritical bifurcation type, and (iv) supersupercritical bifurcation type. Like the Poincaré-Lyapunov coefficient Π(a,γ) in (36) for the single Hopf, there are also explicit formulas for computing the double Hopf coefficients Π(11)(a,γ), Π(21)(a,γ), which are extremely too long to be displayed here. For the account of these formulas for the double Hopf coefficients, again, the publications [22, 23] are excellent references. Classical studies of dynamical systems resulting to equations of the form (53) have demonstrated unprecedented dynamics ranging from global dynamic bifurcations to ultimately chaotic dynamics (see [22, 24]). There exist reliable and robust mathematical tools to examine the long-term behaviour of such dynamics. Campbell et al. [24], in particular, have shown numerically that equations of the form (53) exhibit four double Hopf’s interactions, namely, two super-super-, one super-sub-, and one sub-sub-critical bifurcations. Other interesting dynamics such as large limit cycles, tori, and period doubling were also characterized by the authors.


## 3. Stochastic System

In the section, we present some preliminary results to be used in a subsequent section to establish the stochastic stability and stochastic bifurcation. Before proving the main theorem, we give some lemmas and definitions.Definition 3 (see [1] (D-bifurcation)).
Dynamical bifurcation is concerned with a family of random dynamical systems which is differential and has the invariant measureμα. If there exists a constant αD satisfying in any neighborhood of αD, there exists another constant α and the corresponding invariant measure να≠μα satisfying να→μα as α→αD. Then, the constant αD is a point of dynamical bifurcation.Definition 4 (see [1] (P-bifurcation)).
Phenomenological bifurcation is concerned with the change in the shape of density (stationary probability density) of a family of random dynamical systems as the change of the parameter. If there exists a constantα0 satisfying in any neighborhood of αD, there exist other two constants α1, α2, and their corresponding stationary density functions pα1, pα2 satisfying pα1 and pα2 are not equivalent. Then, the constant α0 is a point of phenomenological bifurcation.Definition 5 (see [1] (stochastic pitchfork bifurcation)).
In the viewpoint of phenomenological bifurcation, the stationary solution of the Fokker-Planck-Kolmogorov (FPK) equation which is corresponded with the stochastic differential equation changes from one peak into two peaks.
In the viewpoint of dynamical bifurcation, there exists a constantμ0 that satisfies the following conditions:(i)
whenμ<μ0, the stochastic differential equation has only one invariant measure ν0; moreover, ν0 is stable (i.e., the maximal Lyapunov exponent is negative).(ii)
whenμ=μ0, the invariant measure ν0 loses its stability and becomes unstable (i.e., the maximal Lyapunov exponent is positive); moreover, the rotation number is zero.(iii)
whenμ>μ0, the stochastic differential equation has three invariant measures ν0, ν1, and ν2; both ν1 and ν2 are stable (i.e., their maximal Lyapunov exponents are both negative).(iv)
the global attractors of the stochastic differential equation change from a single-point set into a collection of one-dimensional (the closure of the unstable manifold of the unstable invariant measure).
If the stochastic bifurcation of a stochastic differential equation has the above characters, then the stochastic differential equation admits stochastic pitchfork bifurcation atμ=μ0.Definition 6 (see [1] (the stochastic Hopf bifurcation)).
In the viewpoint of phenomenological bifurcation, the stationary solution of the FPK equation which is corresponded to the stochastic differential equation changes from one peak into a crater.
In the viewpoint of dynamical bifurcation, the following can be noticed:(i)
if one of the invariant measures of the stochastic differential equation loses its stability and becomes unstable (i.e., two Lyapunov exponents are positive), then the rotation number is not zero. Meanwhile, there at least appears one new stable invariant measure,(ii)
the global attractors of the stochastic differential equation change from a single-point set into a random topological disk (the closure of the unstable manifold of the unstable invariant measure).
If the stochastic bifurcation of a stochastic differential equation has the above characters, then the stochastic differential equation admits the stochastic Hopf bifurcation.Definition 7 (see [1] (stochastically stable)).
The trivial solutionx(t,t0,x0) of stochastic differential equation is said to be stochastically stable or stable in probability if for every pair of ε∈(0,1) and α>0, there exists δ=δ(ε,δ,t0)>0 such that
(54)P{|x(t,t0,x0)|<α∀t≥t0}≥1-ε,
whenever |x0|<δ. Otherwise, it is said to be stochastically unstable.In this section, taking some stochastic factors into account, we introduce randomness into the model (1) by replacing the parameters γ by γ→γ+αξ(t). This is only a first step in introducing stochasticity into the model. Ideally, we would also like to introduce stochastic environmental variation into the other parameters such as the transmission coefficient α and δ, but this would make the analysis much too difficult. Hence, we get the following system of stochastic differential equations:
(55)u˙(t)=ακ(1-u-v)-βu-γu(t-τ1)v(t-τ2)+(σ1+σ2u+σ3v)ξ(t),v˙(t)=δ(1-κ)(1-u-v)3-γuv+(σ4+σ5u+σ6v)ξ(t).Here,ξ(t) is the multiplicative random excitation and the external random excitation directly, and ξ(t) is independent, in possession of zero mean value and standard variance Gaussian white noises, that is, ℰ[ξ(t)]=0,ℰ[ξ(t)ξ(t+τ)]=δ(τ). Note that the intensity of the noise σi(i=1,2,…,6), γ=γc+γ~ is the bifurcation parameter. From (3)–(33), we obtain (55) of the stochastic center manifold:
(56)z˙1(t)=-ωz2+c10z1+c100z2+c11z13+c12z12z2+c13z1y22+c14z23+c15z12+c16z2z1+c17z22+(k11z1+k12z2+k13)ξ(t),z˙2(t)=ωz1+c20z1+c200z2+c21z13+c22z12z2+c23z1z22+c24z23+c25z12+c26z2z1+c27z22+(k21z1+k22z2+k23)ξ(t).
Here, the lengthy expressions of cij and kij are omitted.Setting the coordinate transformationz1=rcosθ, z2=rsinθ and by substituting the variable in (34), we obtain
(57)r˙(t)=rc10cos2θ+r(c100+c20)sinθcosθ+rc200sin2θ+r3[c11cos4θ+(c12+c21)cos3θsinθ+(c13+c22)cos2θsin2θ+(c14+c23)cosθsin3θr3c+c24sin4θ]+r2[c15cos3θ+c24sin4θcccccccc+(c16+c25)cos2θsinθ+c24sin4θcccccccccc+(c17+c26)cosθsin2θ+c27sin3θ]+r[k11cos2θ+(k12+k21)cosθsinθcccc+k22sin2θ]ξ(t)+[k13cosθ+k23sinθ]ξ(t),θ˙(t)=-ω+c20cos2θ+(c200-c10)sinθcosθ-c100sin2θ+r2[c21cos4θ+(c21-c12)cos3θsinθ+(c23+c12)cos2θsin2θ+(c24+c23)cosθsin3θ+c14sin4θ]+r[b25cosθ+(c26-c15)cos2θsinθrccc+(c27-c16)cosθsin2θ-c17sin3θ]+[k21cos2θ+(k22-k11)cosθsinθ-k12sin2θ]ξ(t)+1r[k23cosθ-k13sinθ]ξ(t).It is difficult to calculate the exact solution for system (57) today. According to the Khasminskii limit theorem, when the intensities of the white noises σi(i=1,…,6) are small enough, the response process {r(t),θ(t)} weakly converged to the two-dimensional Markov diffusion process [1–5]. Through the stochastic averaging method, we obtained the Ito^ stochastic differential equation (55) as follows:
(58)dr=[(μ1+μ28)r+μ3r+μ78r3]dt+(μ3+μ48r2)1/2dWr+(rμ5)1/2dWθ,dθ=(μ10+μ88r2)dt+(rμ5)1/2dWr+(μ3r2+μ68)1/2dWθ,
where Wr(t) and Wθ are the independent and standard Wiener processes. Set the parameter as follows:
(59)μ1=12(c10+c200),μ10=-ω+12(c20+c100),μ2=5k112+5k222+3k212+3k122+6k12k21-2k11k22,μ3=12(k132+k232),μ4=3k112+3k222+k122+k212+2k12k21+2k11k22,μ5=14(k11+k22)(k21-k12),μ6=k112+k222+3k122+3k212-2k12k21-2k11k22,μ7=3c11+3c21+c13+c22,μ8=3c21-3c14+c23-c12.


### 3.1. StochasticD-Bifurcation

In the section, we will see how the introduction of randomness changes the stochastic behavior significantly from both the dynamical and phenomenological points of view.Theorem 8 (D-bifurcation).
Letμ3=0,μ7=0. Then, system (55) undergoes stochastic D-bifurcation.Proof.
Whenμ3=0,μ7=0, Then, system (58) becomes
(60)dr=[(μ1+μ28)r]dt+(μ48r2)1/2dWr.
When μ4=0, (60) is a deterministic system, and there is no bifurcation phenomenon. Here, we discuss the situation μ4≠0, let
(61)m(r)=(μ1+μ28-μ416)r,σ(r)=(μ48)1/2r.
The continuous random dynamical system generated by (60) is
(62)φ(t)x=x+∫0t‍m(φ(s)x)ds+∫0t‍σ(φ(s)x)∘dWr,
where ∘dWr is the differential in the sense of Statonovich, and it is the unique strong solution of (60) with the initial value x. And m(0)=0, σ(0)=0, so 0 is a fixed point of φ. Since m(r) is bounded and for any r≠0, it satisfies the ellipticity condition σ(r)≠0. This ensures that there is at most one stationary probability density. According to the Itô equation of amplitude r(t), we obtain its Fokker-Planck-Kolmogorov equation corresponding to (60) as follows:
(63)∂p∂t=-∂∂r{[(μ1+μ28)r]p}+∂2∂r2{[μ48r2]p}.
Let ∂p/∂t=0, then we obtain the solution of system (63)
(64)p(t)=c|σ-1(t)|exp(∫0t‍2m(u)σ2(u)du).
The above equation (63) has two kinds of equilibrium states: fixed point and nonstationary motion. The invariant measure of the former is δ0, and its probability density is δx. The invariant measure of the latter is ν, and its probability density is (64). In the following, we calculate the Lyapunov exponent of the two invariant measures.
Using the solution of linear Itô stochastic differential equation, we obtain the solution of system (60):
(65)r(t)=r(0)exp(∫0t‍[m′(0)+σ(0)σ′(0)2]dsr(0)expcc+∫0t‍σ′(0)dWr).
The Lyapunov exponent with regard to μ of the dynamical system φ is defined as:
(66)λφ(μ)=limt→+∞1tln∥r(t)∥,
substituting (65) into (66), note that σ(0)=0, σ′(0)=0, we obtain the Lyapunov exponent of the fixed point:
(67)λφ(δ0)=limt→+∞1t(ln∥r(0)∥+m′(0)∫0t‍ds+σ′(0)∫0t‍dWr(s))=m′(0)+σ′(0)limt→+∞Wr(t)t=m′(0)=μ1+μ28-μ416.
For the invariant measure which regard (65) as its density, we obtain the Lyapunov exponent:
(68)λφ(ν)=limt→+∞1t∫0t‍[m′(r)+σ(r)σ′(r)]ds=∫R‍[m′(r)+σ(r)σ′(r)2]p(r)dr=-2∫R‍[m(r)σ(r)]2p(r)dr=-322μ43/2m(r)2exp[16μ4m(r)]=-322μ43/2(μ1+μ28-μ416)2×exp[16μ4(μ1+μ28-μ416)].
Letα=μ1+μ2/8-μ4/16. We can obtain that the invariant measure of the fixed point is stable when α<0, but the invariant measure of the nonstationary motion is stable when α>0, so α=αD=0 is a point of D-bifurcation.Theorem 9.
Letμ3=0, μ7=0. Then, system (55) does not undergo stochastic P-bifurcation.Proof.
By simplifying (64), we can obtain
(69)pst(r)=cr2(8μ1+μ2-μ4)/μ4,
where c is a normalization constant; thus, we have
(70)pst(r)=o(rv)r⟶0,
where v=2(8μ1+μ2-μ4)/μ4. Obviously, when v<-1, that is, μ1+μ2/8-μ4/16<0, pst(r) is a δ function. When -1<v<0, that is, μ1+μ2/8-μ4/16>0, r=0 is a maximum point of pst(r) in the state space; thus, the system admits D-bifurcation when v=-1, that is, μ1+μ2/8-μ4/16=0, is the critical condition of   D-bifurcation at the equilibrium point. When v>0, there is no point that makes pst(r) has the maximum value; thus, the system does not admits P-bifurcation.Remark 10.
Unfortunately, these two definitions (D-bifurcation and P-bifurcation) do not necessarily yield the same result.Theorem 11 (stochastic pitchfork bifurcation).
Letγακ=0, μ3=0,μ7<0. Then, system (55) undergoes stochastic pitchfork bifurcation.Proof.
Whenμ3=0,μ7≠0, then (58) can be rewritten as
(71)dr=[(μ1+μ28)r+μ78r3]dt+(μ48r2)1/2dWr.
Let α=μ1+μ2/8, ϕ=(-μ7/8)r, μ7<0, then the system (71) becomes
(72)dϕ=[(μ1+μ28)ϕ-ϕ3]dt+(μ48)1/2ϕ∘dWt,
which is solved by
(73)ϕ→ψα(t,ω)ϕ=(ϕexp((μ1+μ28)t+(μ48)1/2Wt(ω)))×(+(μ48)1/2Ws(ω)))ds∫0t‍)1/2(1+2ϕ2∫0t‍exp((μ48)1/22((μ1+μ28)t+(μ48)1/2Ws(ω)))ds∫0t‍)1/2)-1.
We now determine the domainDα(t,ω), where Dα(t,ω):={ϕ∈ℜ:(t,ω,ϕ)∈D}(D=ℜ×Ω×X) is (in general possibly empty) the set of initial values ϕ∈ℜ for which the trajectories still exist at time t and the range Rα(t,ω) of ψα(t,ω):Dα(t,ω)→Rα(t,ω).
We have(74)Dα(t,ω)={ℜ,t≥0,(-dα(t,ω),dα(t,ω)),t<0,
where
(75)dα(t,ω)=1×(+2(μ48)1/2Ws(ω))ds|)1/2(2|∫0t‍exp((μ48)1/22(μ1+μ28)t+2(μ48)1/2Ws(ω))ds|)1/2)-1>0,Rα(t,ω)=Dα(-t,ϑ(t)ω)Rα(t,ω)={(-rα(t,ω),rα(t,ω)),t>0,ℜ,t≤0,
where
(76)rα(t,ω)=dα(-t,ϑ(t)ω)=(exp((μ1+μ28)t+(μ48)1/2Wt(ω)))×(2(μ48)1/2Ws(ω))ds∫0t‍|)1/2(2|∫0t‍exp((μ48)1/22(μ1+μ28)t+2(μ48)1/2Ws(ω))ds∫0t‍|)1/2)-1>0.
We can now determine that
(77)Eα(ω):=⋂t∈ℜDα(t,ω)
and obtain
(78)Eα(ω)={(-dα-(t,ω),dα-(t,ω)),α=μ1+μ28>0,{0},α=μ1+μ28≤0,
where
(79)0<dα±(t,ω)=1×(+(μ48)1/2Ws(ω))ds∫0±∞‍|)1/2(2|∫0±∞‍exp((μ48)1/22(μ1+μ28)t+(μ48)1/2Ws(ω))ds∫0±∞‍|)1/2)-1<∞.
The ergodic invariant measures of system (71) are as follows.(i)
Forα≤0, the only invariant measure is μωα=δ0.(ii)
Forα>0, we have the three invariant forward Markov measures, μωα=δ0 and υ±,ωα=δ±kα(ω), where
(80)kα(ω):=(2∫-∞0‍exp((μ48)1/22(μ1+μ28)t2∫-∞0‍expccc+2(μ48)1/2Wt(ω))ds)-1/2.
We have𝔼kα2(ω)=α. Solving the forward Fokker-Planck-Kolmogorov equation
(81)L*pα=-(((μ1+μ28)ϕ+μ416ϕ-ϕ3)Pα(ϕ))′+μ416(ϕ2Pν1(ϕ))′′=0,
yield(i)
pα=δ0 for all α,(ii)
forpα>0,
(82)qα+(ϕ)={Nαϕ(2(μ1+μ2/8)/μ4)-1exp(16ϕ2μ4),ϕ>0,0,ϕ≤0,and qα+(ϕ)=qα+(-ϕ), where Nα-=Γ(ν1/μ4)(μ4/8)(μ1+μ2/8)/μ4.
Naturally, the invariant measuresυ±,ωα=δ±kα(ω) are those corresponding to the stationary measures qα+. Hence, all invariant measures are Markov measures.
We determine all invariant measures (necessarily Dirac measure) of local RDSχ generated by the SDE
(83)dϕ=[(μ1+μ28)ϕ-ϕ3]dt+(μ48)1/2ϕ∘dW,
on the state space ℜ, μ1+μ2/8∈ℜ. We now calculate the Lyapunov exponent for each of these measures.
The linearized RDSχt=DΥ(t,ω,ϕ)χ satisfies the linearized SDE
(84)dχt=[(μ1+μ28)-3(Υ(t,ω,ϕ))2χt]dt+(μ48)1/2χt∘dW.
Hence,
(85)DΥ(t,ω,ϕ)χ=χexp((μ1+μ28)t+(μ48)1/2Wt(ω)χexpcc-3∫0t‍(Υ(s,ω,ϕ))2ds(μ48)1/2Wt(ω)).
Thus, if νω=δϕ0(ω) is a Υ-invariant measure, its Lyapunov exponent is
(86)λ(μ)=limt→∞1tlog∥DΥ(t,ω,ϕ)χ∥=(μ1+μ28)-3limt→∞∫0t‍(Υ(s,ω,ϕ))2ds=(μ1+μ28)-3𝔼ϕ02,
providing the IC ϕ02∈L1(ℙ) is satisfied.(i)
Forμ1+μ2/8∈ℜ, the IC for ν±,ωα=δ0 is trivially satisfied, and we obtain
(87)λ(ν1α)=μ1+μ28.
So, ν1α is stable for μ1+μ2/8<0 and unstable for μ1+μ2/8>0.(ii)
Forμ1+μ2/8>0, ν2,ωα=δdωα is ℱ-∞0 measurable; hence, the density pα of ρα=𝔼ν2α satisfies the Fokker-Planck-Kolmogorov equation
(88)Lν1*=-(((μ1+μ28)ϕ+μ416ϕ-ϕ3)pα(ϕ))′+μ416(ϕ2pα(ϕ))′′=0,which has the unique probability density solution
(89)Pα(ϕ)=Nαϕ(2(μ1+μ2/8)/μ4)-1exp(ϕ28μ4),ϕ>0.
Since
(90)𝔼ν2αϕ2=𝔼(d-α)2=∫0∞‍ϕ2pα(ϕ)dϕ<∞,
the IC is satisfied. The calculation of the Lyapunov exponent is accomplished by observing that
(91)d-α(ϑtω)2=exp(2(μ1+μ2/8)t+2(μ4/8)1/2Wt(ω))2∫-∞t‍exp(2(μ1+μ2/8)s+2(μ4/8)1/2Ws(ω))ds=Ψ′(t)2Ψ,
where
(92)Ψ(t)=∫-∞t‍exp((μ1+μ28)s+(μ48)1/2Ws(ω))ds.
Hence, by the ergodic theorem
(93)𝔼(d-α)2=12limt→∞1tlogΨ(t)=μ1+μ28,
finally
(94)λ(ν2α)=-2(μ1+μ28)<0.(iii)
For(μ1+μ2/8)>0, ν2,ωα=δdωα is ℱ-∞0 measurable. Since ℒ(d+α)=ℒ(d-α)(95)𝔼(-d-α)2=𝔼(d-α)2=μ1+μ28,
thus
(96)λ(ν2α)=-2(μ1+μ28)<0.
The two families of densities(qα+)α>0 clearly undergo a P-bifurcation at the parameter value αP=μ4/8, which is the same value as the transcritical case, since the SDE linearized at ϕ=0 is the same in the next section. In both cases, we have a D-bifurcation of the trivial reference measure δ0 at αD=0 and a P-bifurcation of αP=(μ1+μ2/8)2/2. Then, system (55) has stochastic pitchfork bifurcation.


### 3.2.P-Bifurcation

In the following, we investigate the steady-state probability densitypst(r) of the linear Itô stochastic differential equation. Calculating extreme values of the steady-state probability density is one of the most popular efficient methods in studying the bifurcation of a nonlinear dynamical system. The steady-state probability density is an important characteristic value of stochastic bifurcation.Case I.
Whenμ3=0,μ7≠0, then (38) can be rewritten as
(97)dr=[(μ1+μ28)r+μ78r3]dt+(μ48r2)1/2dWr.According to the Itô equation of amplituder(t), we obtain its Fokker-Planck-Kolmogorov equation form (97) as follows:
(98)∂p∂t=-∂∂r{[(μ1+μ28)r-μ42μ7r-r3]p(r)}-μ42μ7∂2∂r2(r2p(r)),
with the initial value condition μ7=0,p(r,t|r0,t0)→δ(r-r0),t→t0, where p(r,t|r0,t0) is the transition probability density of diffusion process r(t). The invariant measure of r(t) is the steady-state probability density pst(r) which is the solution of the degenerate system as follows:
(99)0=-∂∂r{[(μ1+μ28)r-μ42μ7r-r3]p(r)}-μ42μ7∂2∂r2(r2p(r)).
Through calculation, we can obtain
(100)pst(r)=exp(r2μ7/μ4)r-1-(2(μ1+μ2/8)μ7/μ4)Γ[-(μ1+μ2/8)μ7/μ4](-μ7/μ4)(μ1+μ2/8)μ7/μ4.According to Namachivaya’s theory, the extreme value of an steady-state probability density contains the most important essence of the nonlinear stochastic system. In other words, the steady-state probability density can uncover the characteristic information of the steady state. When the intension of the noise tends to zero, the extreme values ofpst(r) approximate to show the behavior of the deterministic system. If the process r(t) is ergodic, then pst(r) can be regarded as the time measurement for staying in the neighborhood of r(t) according to the Oseledec ergodic theorem.From the analysis above, ifpst(r) has a maximum value at r*, the sample trajectory will stay for a longer time in the neighborhood of r*, that is, r* is stable in the meaning of probability (with a bigger probability). If pst(r) has a minimum value (zero), it is just the opposite.We now calculate the most possible amplituder* of system (97), that is, pst(r) has a maximum value at r*. So, we have
(101)dpst(r)dr|r=r*=0,d2pst(r)dr2|r=r*<0
and the solution r=r~=(ν1ν2+4ν3)/8ν2. The probabilities and the positions of the Hopf bifurcation occur with different parameter.Since(102)d2pst(r)dr2|=r~=(((8(μ1+(μ2/8))μ7+4μ4)μ7)-2(μ1+(μ2/8))μ7/4μ424+(3(μ1+μ2/8)μ2/μ3)exp(18(4+μ7μ3))c×μ7(-μ7μ4)-(μ1+μ2/8)μ7/ν4c×((8(μ1+μ2/8)μ7+4μ4)μ7)-2(μ1+μ2/8)μ7/4μ4)c×(Γ[-(μ1+μ2/8)μ7μ4](-μ7μ4)(μ1+μ2/8)μ7/μ4)-1<0,
thus, what we need is r*=r~. The conclusion is to go all the way with what has been obtained by the singular boundary theory. The original nonlinear stochastic system has a stochastic Hopf bifurcation at r=r~, where
(103)x12+x22=8(μ1+μ2/8)μ7+4μ48μ7.We now choose some values of the parameters in the equations, draw the graphics ofpst(r) (Figure 1). It is worth putting forward that calculating the Hopf bifurcation with the parameters in the original system is necessary. If we now have values of the original parameters in system (55), then α=0.778205, β=0.024530, κ=0.08, γ=5.640593,δ=0.200101, σ1=0.1, σ2=0.1, σ3=0.2, σ4=0.3, σ5=0.1, σ6=0.5. After further calculations, we obtain μ1=-0.2,μ2=9.5907,μ3=0,μ4=5.1183,μ7=-20.4732, then
(104)p(r)=0.28209e-4r2.
What is more is that r~=0.28209, where pst(r) has the maximum value.P-Bifurcation of system (55) with μ1=-0.2,μ3=0,μ4=5.1183,μ7=-20.4732; (a) P-bifurcation of pα(r) at αp<σ2/2,μ2=5.6, where α=μ1+μ2/8,σ2=-μ4/μ7; (b) P-bifurcation of pα(r) at αp=σ2/2,μ2=9.5907, where α=(μ1+μ2/8),σ2=-μ4/μ7; (c) P-bifurcation of pα(r) at αp>σ2/2,μ2=17.6, where α=ν1/8,σ2=-ν3/ν2.
(a)
(b)
(c)


## 3.1. StochasticD-Bifurcation

In the section, we will see how the introduction of randomness changes the stochastic behavior significantly from both the dynamical and phenomenological points of view.Theorem 8 (D-bifurcation).
Letμ3=0,μ7=0. Then, system (55) undergoes stochastic D-bifurcation.Proof.
Whenμ3=0,μ7=0, Then, system (58) becomes
(60)dr=[(μ1+μ28)r]dt+(μ48r2)1/2dWr.
When μ4=0, (60) is a deterministic system, and there is no bifurcation phenomenon. Here, we discuss the situation μ4≠0, let
(61)m(r)=(μ1+μ28-μ416)r,σ(r)=(μ48)1/2r.
The continuous random dynamical system generated by (60) is
(62)φ(t)x=x+∫0t‍m(φ(s)x)ds+∫0t‍σ(φ(s)x)∘dWr,
where ∘dWr is the differential in the sense of Statonovich, and it is the unique strong solution of (60) with the initial value x. And m(0)=0, σ(0)=0, so 0 is a fixed point of φ. Since m(r) is bounded and for any r≠0, it satisfies the ellipticity condition σ(r)≠0. This ensures that there is at most one stationary probability density. According to the Itô equation of amplitude r(t), we obtain its Fokker-Planck-Kolmogorov equation corresponding to (60) as follows:
(63)∂p∂t=-∂∂r{[(μ1+μ28)r]p}+∂2∂r2{[μ48r2]p}.
Let ∂p/∂t=0, then we obtain the solution of system (63)
(64)p(t)=c|σ-1(t)|exp(∫0t‍2m(u)σ2(u)du).
The above equation (63) has two kinds of equilibrium states: fixed point and nonstationary motion. The invariant measure of the former is δ0, and its probability density is δx. The invariant measure of the latter is ν, and its probability density is (64). In the following, we calculate the Lyapunov exponent of the two invariant measures.
Using the solution of linear Itô stochastic differential equation, we obtain the solution of system (60):
(65)r(t)=r(0)exp(∫0t‍[m′(0)+σ(0)σ′(0)2]dsr(0)expcc+∫0t‍σ′(0)dWr).
The Lyapunov exponent with regard to μ of the dynamical system φ is defined as:
(66)λφ(μ)=limt→+∞1tln∥r(t)∥,
substituting (65) into (66), note that σ(0)=0, σ′(0)=0, we obtain the Lyapunov exponent of the fixed point:
(67)λφ(δ0)=limt→+∞1t(ln∥r(0)∥+m′(0)∫0t‍ds+σ′(0)∫0t‍dWr(s))=m′(0)+σ′(0)limt→+∞Wr(t)t=m′(0)=μ1+μ28-μ416.
For the invariant measure which regard (65) as its density, we obtain the Lyapunov exponent:
(68)λφ(ν)=limt→+∞1t∫0t‍[m′(r)+σ(r)σ′(r)]ds=∫R‍[m′(r)+σ(r)σ′(r)2]p(r)dr=-2∫R‍[m(r)σ(r)]2p(r)dr=-322μ43/2m(r)2exp[16μ4m(r)]=-322μ43/2(μ1+μ28-μ416)2×exp[16μ4(μ1+μ28-μ416)].
Letα=μ1+μ2/8-μ4/16. We can obtain that the invariant measure of the fixed point is stable when α<0, but the invariant measure of the nonstationary motion is stable when α>0, so α=αD=0 is a point of D-bifurcation.Theorem 9.
Letμ3=0, μ7=0. Then, system (55) does not undergo stochastic P-bifurcation.Proof.
By simplifying (64), we can obtain
(69)pst(r)=cr2(8μ1+μ2-μ4)/μ4,
where c is a normalization constant; thus, we have
(70)pst(r)=o(rv)r⟶0,
where v=2(8μ1+μ2-μ4)/μ4. Obviously, when v<-1, that is, μ1+μ2/8-μ4/16<0, pst(r) is a δ function. When -1<v<0, that is, μ1+μ2/8-μ4/16>0, r=0 is a maximum point of pst(r) in the state space; thus, the system admits D-bifurcation when v=-1, that is, μ1+μ2/8-μ4/16=0, is the critical condition of   D-bifurcation at the equilibrium point. When v>0, there is no point that makes pst(r) has the maximum value; thus, the system does not admits P-bifurcation.Remark 10.
Unfortunately, these two definitions (D-bifurcation and P-bifurcation) do not necessarily yield the same result.Theorem 11 (stochastic pitchfork bifurcation).
Letγακ=0, μ3=0,μ7<0. Then, system (55) undergoes stochastic pitchfork bifurcation.Proof.
Whenμ3=0,μ7≠0, then (58) can be rewritten as
(71)dr=[(μ1+μ28)r+μ78r3]dt+(μ48r2)1/2dWr.
Let α=μ1+μ2/8, ϕ=(-μ7/8)r, μ7<0, then the system (71) becomes
(72)dϕ=[(μ1+μ28)ϕ-ϕ3]dt+(μ48)1/2ϕ∘dWt,
which is solved by
(73)ϕ→ψα(t,ω)ϕ=(ϕexp((μ1+μ28)t+(μ48)1/2Wt(ω)))×(+(μ48)1/2Ws(ω)))ds∫0t‍)1/2(1+2ϕ2∫0t‍exp((μ48)1/22((μ1+μ28)t+(μ48)1/2Ws(ω)))ds∫0t‍)1/2)-1.
We now determine the domainDα(t,ω), where Dα(t,ω):={ϕ∈ℜ:(t,ω,ϕ)∈D}(D=ℜ×Ω×X) is (in general possibly empty) the set of initial values ϕ∈ℜ for which the trajectories still exist at time t and the range Rα(t,ω) of ψα(t,ω):Dα(t,ω)→Rα(t,ω).
We have(74)Dα(t,ω)={ℜ,t≥0,(-dα(t,ω),dα(t,ω)),t<0,
where
(75)dα(t,ω)=1×(+2(μ48)1/2Ws(ω))ds|)1/2(2|∫0t‍exp((μ48)1/22(μ1+μ28)t+2(μ48)1/2Ws(ω))ds|)1/2)-1>0,Rα(t,ω)=Dα(-t,ϑ(t)ω)Rα(t,ω)={(-rα(t,ω),rα(t,ω)),t>0,ℜ,t≤0,
where
(76)rα(t,ω)=dα(-t,ϑ(t)ω)=(exp((μ1+μ28)t+(μ48)1/2Wt(ω)))×(2(μ48)1/2Ws(ω))ds∫0t‍|)1/2(2|∫0t‍exp((μ48)1/22(μ1+μ28)t+2(μ48)1/2Ws(ω))ds∫0t‍|)1/2)-1>0.
We can now determine that
(77)Eα(ω):=⋂t∈ℜDα(t,ω)
and obtain
(78)Eα(ω)={(-dα-(t,ω),dα-(t,ω)),α=μ1+μ28>0,{0},α=μ1+μ28≤0,
where
(79)0<dα±(t,ω)=1×(+(μ48)1/2Ws(ω))ds∫0±∞‍|)1/2(2|∫0±∞‍exp((μ48)1/22(μ1+μ28)t+(μ48)1/2Ws(ω))ds∫0±∞‍|)1/2)-1<∞.
The ergodic invariant measures of system (71) are as follows.(i)
Forα≤0, the only invariant measure is μωα=δ0.(ii)
Forα>0, we have the three invariant forward Markov measures, μωα=δ0 and υ±,ωα=δ±kα(ω), where
(80)kα(ω):=(2∫-∞0‍exp((μ48)1/22(μ1+μ28)t2∫-∞0‍expccc+2(μ48)1/2Wt(ω))ds)-1/2.
We have𝔼kα2(ω)=α. Solving the forward Fokker-Planck-Kolmogorov equation
(81)L*pα=-(((μ1+μ28)ϕ+μ416ϕ-ϕ3)Pα(ϕ))′+μ416(ϕ2Pν1(ϕ))′′=0,
yield(i)
pα=δ0 for all α,(ii)
forpα>0,
(82)qα+(ϕ)={Nαϕ(2(μ1+μ2/8)/μ4)-1exp(16ϕ2μ4),ϕ>0,0,ϕ≤0,and qα+(ϕ)=qα+(-ϕ), where Nα-=Γ(ν1/μ4)(μ4/8)(μ1+μ2/8)/μ4.
Naturally, the invariant measuresυ±,ωα=δ±kα(ω) are those corresponding to the stationary measures qα+. Hence, all invariant measures are Markov measures.
We determine all invariant measures (necessarily Dirac measure) of local RDSχ generated by the SDE
(83)dϕ=[(μ1+μ28)ϕ-ϕ3]dt+(μ48)1/2ϕ∘dW,
on the state space ℜ, μ1+μ2/8∈ℜ. We now calculate the Lyapunov exponent for each of these measures.
The linearized RDSχt=DΥ(t,ω,ϕ)χ satisfies the linearized SDE
(84)dχt=[(μ1+μ28)-3(Υ(t,ω,ϕ))2χt]dt+(μ48)1/2χt∘dW.
Hence,
(85)DΥ(t,ω,ϕ)χ=χexp((μ1+μ28)t+(μ48)1/2Wt(ω)χexpcc-3∫0t‍(Υ(s,ω,ϕ))2ds(μ48)1/2Wt(ω)).
Thus, if νω=δϕ0(ω) is a Υ-invariant measure, its Lyapunov exponent is
(86)λ(μ)=limt→∞1tlog∥DΥ(t,ω,ϕ)χ∥=(μ1+μ28)-3limt→∞∫0t‍(Υ(s,ω,ϕ))2ds=(μ1+μ28)-3𝔼ϕ02,
providing the IC ϕ02∈L1(ℙ) is satisfied.(i)
Forμ1+μ2/8∈ℜ, the IC for ν±,ωα=δ0 is trivially satisfied, and we obtain
(87)λ(ν1α)=μ1+μ28.
So, ν1α is stable for μ1+μ2/8<0 and unstable for μ1+μ2/8>0.(ii)
Forμ1+μ2/8>0, ν2,ωα=δdωα is ℱ-∞0 measurable; hence, the density pα of ρα=𝔼ν2α satisfies the Fokker-Planck-Kolmogorov equation
(88)Lν1*=-(((μ1+μ28)ϕ+μ416ϕ-ϕ3)pα(ϕ))′+μ416(ϕ2pα(ϕ))′′=0,which has the unique probability density solution
(89)Pα(ϕ)=Nαϕ(2(μ1+μ2/8)/μ4)-1exp(ϕ28μ4),ϕ>0.
Since
(90)𝔼ν2αϕ2=𝔼(d-α)2=∫0∞‍ϕ2pα(ϕ)dϕ<∞,
the IC is satisfied. The calculation of the Lyapunov exponent is accomplished by observing that
(91)d-α(ϑtω)2=exp(2(μ1+μ2/8)t+2(μ4/8)1/2Wt(ω))2∫-∞t‍exp(2(μ1+μ2/8)s+2(μ4/8)1/2Ws(ω))ds=Ψ′(t)2Ψ,
where
(92)Ψ(t)=∫-∞t‍exp((μ1+μ28)s+(μ48)1/2Ws(ω))ds.
Hence, by the ergodic theorem
(93)𝔼(d-α)2=12limt→∞1tlogΨ(t)=μ1+μ28,
finally
(94)λ(ν2α)=-2(μ1+μ28)<0.(iii)
For(μ1+μ2/8)>0, ν2,ωα=δdωα is ℱ-∞0 measurable. Since ℒ(d+α)=ℒ(d-α)(95)𝔼(-d-α)2=𝔼(d-α)2=μ1+μ28,
thus
(96)λ(ν2α)=-2(μ1+μ28)<0.
The two families of densities(qα+)α>0 clearly undergo a P-bifurcation at the parameter value αP=μ4/8, which is the same value as the transcritical case, since the SDE linearized at ϕ=0 is the same in the next section. In both cases, we have a D-bifurcation of the trivial reference measure δ0 at αD=0 and a P-bifurcation of αP=(μ1+μ2/8)2/2. Then, system (55) has stochastic pitchfork bifurcation.


## 3.2.P-Bifurcation

In the following, we investigate the steady-state probability densitypst(r) of the linear Itô stochastic differential equation. Calculating extreme values of the steady-state probability density is one of the most popular efficient methods in studying the bifurcation of a nonlinear dynamical system. The steady-state probability density is an important characteristic value of stochastic bifurcation.Case I.
Whenμ3=0,μ7≠0, then (38) can be rewritten as
(97)dr=[(μ1+μ28)r+μ78r3]dt+(μ48r2)1/2dWr.According to the Itô equation of amplituder(t), we obtain its Fokker-Planck-Kolmogorov equation form (97) as follows:
(98)∂p∂t=-∂∂r{[(μ1+μ28)r-μ42μ7r-r3]p(r)}-μ42μ7∂2∂r2(r2p(r)),
with the initial value condition μ7=0,p(r,t|r0,t0)→δ(r-r0),t→t0, where p(r,t|r0,t0) is the transition probability density of diffusion process r(t). The invariant measure of r(t) is the steady-state probability density pst(r) which is the solution of the degenerate system as follows:
(99)0=-∂∂r{[(μ1+μ28)r-μ42μ7r-r3]p(r)}-μ42μ7∂2∂r2(r2p(r)).
Through calculation, we can obtain
(100)pst(r)=exp(r2μ7/μ4)r-1-(2(μ1+μ2/8)μ7/μ4)Γ[-(μ1+μ2/8)μ7/μ4](-μ7/μ4)(μ1+μ2/8)μ7/μ4.According to Namachivaya’s theory, the extreme value of an steady-state probability density contains the most important essence of the nonlinear stochastic system. In other words, the steady-state probability density can uncover the characteristic information of the steady state. When the intension of the noise tends to zero, the extreme values ofpst(r) approximate to show the behavior of the deterministic system. If the process r(t) is ergodic, then pst(r) can be regarded as the time measurement for staying in the neighborhood of r(t) according to the Oseledec ergodic theorem.From the analysis above, ifpst(r) has a maximum value at r*, the sample trajectory will stay for a longer time in the neighborhood of r*, that is, r* is stable in the meaning of probability (with a bigger probability). If pst(r) has a minimum value (zero), it is just the opposite.We now calculate the most possible amplituder* of system (97), that is, pst(r) has a maximum value at r*. So, we have
(101)dpst(r)dr|r=r*=0,d2pst(r)dr2|r=r*<0
and the solution r=r~=(ν1ν2+4ν3)/8ν2. The probabilities and the positions of the Hopf bifurcation occur with different parameter.Since(102)d2pst(r)dr2|=r~=(((8(μ1+(μ2/8))μ7+4μ4)μ7)-2(μ1+(μ2/8))μ7/4μ424+(3(μ1+μ2/8)μ2/μ3)exp(18(4+μ7μ3))c×μ7(-μ7μ4)-(μ1+μ2/8)μ7/ν4c×((8(μ1+μ2/8)μ7+4μ4)μ7)-2(μ1+μ2/8)μ7/4μ4)c×(Γ[-(μ1+μ2/8)μ7μ4](-μ7μ4)(μ1+μ2/8)μ7/μ4)-1<0,
thus, what we need is r*=r~. The conclusion is to go all the way with what has been obtained by the singular boundary theory. The original nonlinear stochastic system has a stochastic Hopf bifurcation at r=r~, where
(103)x12+x22=8(μ1+μ2/8)μ7+4μ48μ7.We now choose some values of the parameters in the equations, draw the graphics ofpst(r) (Figure 1). It is worth putting forward that calculating the Hopf bifurcation with the parameters in the original system is necessary. If we now have values of the original parameters in system (55), then α=0.778205, β=0.024530, κ=0.08, γ=5.640593,δ=0.200101, σ1=0.1, σ2=0.1, σ3=0.2, σ4=0.3, σ5=0.1, σ6=0.5. After further calculations, we obtain μ1=-0.2,μ2=9.5907,μ3=0,μ4=5.1183,μ7=-20.4732, then
(104)p(r)=0.28209e-4r2.
What is more is that r~=0.28209, where pst(r) has the maximum value.P-Bifurcation of system (55) with μ1=-0.2,μ3=0,μ4=5.1183,μ7=-20.4732; (a) P-bifurcation of pα(r) at αp<σ2/2,μ2=5.6, where α=μ1+μ2/8,σ2=-μ4/μ7; (b) P-bifurcation of pα(r) at αp=σ2/2,μ2=9.5907, where α=(μ1+μ2/8),σ2=-μ4/μ7; (c) P-bifurcation of pα(r) at αp>σ2/2,μ2=17.6, where α=ν1/8,σ2=-ν3/ν2.
(a)
(b)
(c)


## 4. Conclusion

In this work, we investigate both the deterministic and stochastic bifurcations of the catalytic CO oxidation on Ir(111) surfaces with multiple delays. By replacingγ~ with γ-γc, it can be readily seen that the addition of the noise gives rise to a shift in the critical bifurcation point γc at the Hopf bifurcation. The results of our computational studies using the Lyapunov exponents and the integral stochastic averaging combined with the Hopf bifurcation produced stability boundaries for structural systems that were much sharper and well defined than those obtained in the earlier stability studies of the same systems. The influence of noise on the reaction-diffusion system model has been studied by the authors [11–18], but explicit and derived conditions for stochastic stability and bifurcation, as the ones displayed here, were not transparent in their publications.



---

*Source: 790946-2013-03-06.xml*