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

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.

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 Section 2, 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−𝑢−𝑣)−𝛽𝑢−𝛾𝑢𝑡−𝜏₁𝑣𝑡−𝜏₂,𝑣̇(𝑡)=𝛿(1−𝜅)(1−𝑢−𝑣)³−𝛾𝑢𝑣, (1) where 𝛾=𝛾_𝑐+𝛾 is the bifurcation parameter. It is obvious that there is a unique equilibrium point 𝐸(0,1) in system (1).

Let 𝑥₁=𝑢−0,𝑥₂=𝑣−1,and𝑋=[𝑥₁,𝑥₂]^⊤, then by substituting the corresponding variables in (1), we have 𝑥̇₁(𝑡)=−𝛼𝜅𝑥₁+𝑥₂−𝛽𝑥₁−𝛾𝑥₁𝑡−𝜏₁𝑥₂𝑡−𝜏₂,𝑥̇+1₂(𝑡)=−𝛿(1−𝜅)𝑥₁+𝑥₂³−𝛾𝑥₁𝑥₂+1. (2)

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) in 𝒞∶=𝒞([−𝜏,0],ℜ^𝑛), and its adjoint form 𝑥̇₁𝑡̂=𝛼𝜅+𝛽𝑥₁𝑡̂+𝛼𝜅𝑥₂𝑡̂+𝛾𝑥₁𝑡̂+𝜏₁,𝑥̇₂𝑡̂=𝛾𝑥₁𝑡̂, (4) in 𝒞∶=𝒞([−𝜏,0],ℜ^𝑛) with respect to the associated bilinear relation 𝜓_𝑗(𝜗),𝜙_𝑘(𝜃)=𝜓_𝑗(0),𝜙_𝑘(0)−𝛾⁰_−𝜏1𝜓_𝑗𝜁+𝜏₁𝜙_𝑘(𝜁)𝑑𝜁𝑗,𝑘=1,2,…, (5) has eigenvalues 𝜆 satisfying the transcendental characteristic equation Δ(𝜆,𝑁)∶=𝜆²+𝜆𝛼𝜅+𝛽+𝛾𝑒^−𝜆𝜏1−𝛾𝛼𝜅=0. (6)

By the Hopf bifurcation, we let 𝜆_1,2=±𝑖𝜔(𝛾), with, 𝜔(𝛾)≠0 and ℜ𝑒{𝑑Δ(𝜆,𝛾)/𝑑𝛾}≠0, be the eigenvalue of (6). All the remaining eigenvalues of Δ(𝜆,𝛾)=0 have negative real parts. Then, by substituting 𝜆₁=𝑖𝜔(𝛾) into Δ(𝜆,𝛾)=0 and setting the real and imaginary parts of the resulting algebraic equations to zero, we have −𝜔²−𝛾𝛼𝜅+𝜔𝛾sin𝜔𝜏₁+𝜔𝛼𝜅+𝛽+𝜔𝛾cos𝜔𝜏₁𝑖=0. (7)

Separating the real and imaginary parts of (7) gives the following equations: −𝜔²−𝛾𝛼𝜅+𝜔𝛾sin𝜔𝜏₁𝜔=0,𝛼𝜅+𝛽+𝜔𝛾cos𝜔𝜏₁=0, (8) or equivalently sin𝜔𝜏₁=𝜔²+𝛾𝛼𝜅𝛾𝜔>0,cos𝜔𝜏₁=𝛼𝜅+𝛽𝛾𝜔>0. (9) Possible values of 𝜔𝜏₁ are of the form 𝜔𝜏₁∈(2𝑘𝜋,2𝑘𝜋+𝜋/2), 𝑘=0,1,2,…. Thus, we obtain the bifurcation parameter 𝛾 as 𝛾=−2𝛼𝜅𝜔²−4𝛼²𝜅²𝜔⁴−4𝛼²𝜅²−𝜔²𝛽²+2𝛼𝛽𝜅+𝛼²𝜅²+𝜔⁴×2𝛼²𝜅²−𝜔²⁻¹(𝛼𝜅<𝜔). (10)

Now, we compute Hopf’s transversality condition by differentiating implicitly (6) with respect to 𝛾, obtaining 𝑑𝜆𝑑𝛾_{𝜆=𝑖𝜔}=−𝑑Δ𝜆,𝛾𝑑𝛾𝑑Δ𝜆,𝛾𝑑𝜆⁻¹_{𝜆=𝑖𝜔}=−𝜆𝑒^−𝜆𝜏1−𝛼𝜅2𝜆+(𝛼𝜅+𝛽)+𝛾𝑒^−𝜆𝜏1−𝛾𝜏₁𝑒^−𝜆𝜏1_{𝜆=𝑖𝜔}𝛾𝛾𝐺(𝑅)=𝐺+𝑀𝑖,=−𝑎𝑐+𝑏𝑑𝑎²+𝑏²,𝑀𝛾=−𝑎𝑑−𝑏𝑐𝑎²+𝑏²,𝑎=𝜔sin𝜔𝜏₁−𝛼𝜅,𝑏=𝜔cos𝜔𝜏₁,𝑐=𝛼𝜅+𝛾cos𝜔𝜏₁−𝛾𝜏₁𝜔sin𝜔𝜏₁,𝑑=2𝜔−𝛾sin𝜔𝜏₁−𝛾𝜏₁𝜔cos𝜔𝜏₁. (11)

Lemma 1. The transversality conditions ℜ𝑒(𝑑𝜆/𝑑𝛾)_{𝜆=𝜔𝑖}=−𝑎𝑐−𝑏𝑑≠0 hold. If the choice of values for the model parameters is such that the inequality 𝜔sin𝜔𝜏₁−𝛼𝜅𝛼𝜅+𝛾cos𝜔𝜏₁−𝛾𝜏₁𝜔sin𝜔𝜏₁+𝜔cos𝜔𝜏₁2𝜔−𝛾sin𝜔𝜏₁−𝛾𝜏₁𝜔cos𝜔𝜏₁<0 (12) holds.

This means that the pair of eigenvalues 𝜆_1,2=𝑣(𝛾)+𝑖𝜔(𝛾) of Δ(𝜆,𝛾)=0 with 𝑣(𝛾_𝑐)=0, 𝜔(𝛾_𝑐)≠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 𝛾_𝑐. Also, it ensures that the steady-state stability of the system is lost at 𝛾=𝛾_𝑐, 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 𝜔² and letting 𝜔²=𝜔, we have 𝜔²+𝛼₁𝜔+𝛼₂=0,𝛼₁=2𝛼𝛾𝜅−𝛾²,𝛼₂=𝛼²𝛾²𝜅²+𝛼𝜅+𝛽², (13) whose real and positive solutions are of the form 𝜔₁=12−𝛼₁+𝛼²₁−4𝛼₂,𝜔₂=12−𝛼₁−𝛼²₁−4𝛼₂, (14) if the following holds 𝛼₁𝛼<0∣2𝛼𝜅<𝛾,₂₁−4𝛼₂>0∣𝛾⁴−2𝛼𝜅𝛾³−4𝛼𝜅+𝛽²𝛼>0,₁𝛼<0∣2𝛼𝜅<𝛾,₂₁−4𝛼₂=0∣𝛾⁴=2𝛼𝜅𝛾³+4𝛼𝜅+𝛽². (15) In addition, one has to show that the substitution of 𝜔₁ and 𝜔₂ into (10) will ensure Hopf’s transversality conditions, that is, ℜ𝑒(𝑑𝜆/𝑑𝛾)_{𝜆=𝑖𝜔}. Therefore, with 𝜔₁=𝜔² and ℜ𝑒(𝑑𝜆/𝑑𝛾)_{𝜆=𝑖𝜔}, we obtain 2𝜔²+𝛼₁=0∣𝜔=𝛾²−2𝛼𝛾𝜅2,𝛼²₁−4𝛼₂=0∣𝛾⁴=2𝛼𝜅𝛾³+4𝛼𝜅+𝛽². (16) By setting 𝜔_min=𝜔 and 𝛾_min=𝛾, we have 𝜔_min=𝛾²_min−2𝛼𝛾_min𝜅2,𝛾⁴_min=2𝛼𝜅𝛾³_min+4𝛼𝜅+𝛽². (17)

These are the minimum values separating the regimes of stable and unstable solutions of the linearized delay equation at the Hopf bifurcation point 𝛾=𝛾_𝑐, 𝜆=±𝜔𝑖(𝛾_min)=±𝛾²_min−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 𝜏₁ begins by further rewiring (8) as follows: tan𝜔𝜏₁=𝜔²+𝛾𝛼𝜅,𝛼𝜅+𝛽 (18) thus, we obtain 𝜏_1,𝑘=1𝜔𝜔arctan²+𝛾𝛼𝜅𝛼𝜅+𝛽+2𝑘𝜋𝑘=0,1,2,…,𝑛,…, (19) and by using the minimum value for 𝜔_min=𝛾²_min−2𝛼𝛾_min𝜅/2, we get 𝜏_1,𝑘,min=1𝜔_min𝜔arctan^2min+𝛾𝛼𝜅𝛼𝜅+𝛽+2𝑘𝜋𝑘=0,1,2,…,𝑛,…. (20) With these minimum values for the gain 𝜔_min, time delay 𝜏_1,𝑘,min, and the response frequency 𝜔_min=𝛾²_min−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=𝛾²_min−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 𝜙^𝑃(𝜃),𝜙^𝑄(𝜃) onto the center and stable subspaces 𝑃,𝑄∈𝒞. For the eigenvalues 𝜆=±𝑖𝜔, we have 𝜙^𝑃(𝜃)=Φ(𝜃)𝑏 with 𝜙^𝑃_𝑘(𝜃)=𝑒^𝜆𝑘𝜃, −𝜏≤𝜃≤0 of the adjoint-generalized eigenspace 𝑃∈𝒞 and 𝜓𝑃(𝑠)=Ψ𝑏̂, 𝜓𝑃_𝑗(𝑠)𝑒^{−𝜆𝑗𝑠}, 0≤𝑠≤𝜏,𝑗,𝑘=1,2 of the adjoint-generalized eigenspace 𝑃∈𝒞. The values of the basis Φ(𝜃)∈𝒞 for 𝑃 are thus given by Φ(𝜃)=𝑏𝑏cos𝜔𝜃−sin𝜔𝜃𝜔sin𝜔𝜃𝜔cos𝜔𝜃, (21) and, similarly, the basis Ψ(𝑠)∈𝒞 for 𝑃 is of the form Ψ(𝜃)=𝜓₁,𝜓₂=cos𝜔𝜃𝜔sin𝜔𝜃−sin𝜔𝜃𝜔cos𝜔𝜃,0≤𝑠≤𝑅. (22)

These bases form the inner product matrix (Ψ(𝑠)(𝜃=𝜓,Φ))¹(𝑠)𝜓₂(𝑠)𝜙₁𝜙₂=𝜓₁(𝑠)𝜙₁(𝜃)𝜓₁(𝑠)𝜙₂(𝜃𝜓₂)(𝑠)𝜙₁(𝜃)𝜓₂(𝑠)𝜙₂(𝜃,) (23) in which 𝜓₁(𝑠)𝜙₁(𝜃)=cos𝜔𝑠cos𝜔𝜃+𝜔²,𝜓sin𝜔𝑠sin𝜔𝜃₁(𝑠)𝜙₂(𝜃)=−cos𝜔𝑠sin𝜔𝜃−𝜔²,𝜓sin𝜔𝑠cos𝜔𝜃₂(𝑠)𝜙₁(𝜃)=−sin𝜔𝑠cos𝜔𝜃−𝜔²,𝜓cos𝜔𝑠sin𝜔𝜃₂(𝑠)𝜙₂(𝜃)=sin𝜔𝑠sin𝜔𝜃+𝜔²cos𝜔𝑠cos𝜔𝜃. (24)

The substitution of the elements of (Ψ,Φ)_𝑛𝑠𝑔≠0(Ψ,Φ)_𝑛𝑠𝑔=𝜓₁,𝜙₁𝜓₁,𝜙₂𝜓₂,𝜙₁𝜓₂,𝜙₂, (25) where 𝜓₁,𝜙₁1=1−2𝛾1𝜔−𝜔sin𝜔𝜏₁+𝜏₁1+𝜔²cos𝜔𝜏₁,𝜓₁,𝜙₂=𝛾1+𝜔²sin𝜔𝜏₁,𝜓₂(𝑠)𝜙₁(𝜃)=−𝛾1+𝜔²sin𝜔𝜏₁,𝜓₂(𝑠)𝜙₂(𝜃)=𝜔²1+𝛾𝜔−𝜔sin𝜔𝜏₁−𝜏₁1+𝜔²cos𝜔𝜏₁. (26)

Then, the basis Ψ(𝑠) for 𝑃∈𝒞 of the adjoint equation (2) is normalized to Ψ=[𝜓₁,𝜓₂]^⊤∈𝒞, computing Ψ=(Ψ,Φ)⁻¹_𝑛𝑠𝑔Ψ(𝑠) to yield Ψ(𝑠)=𝜓₁₁(𝑠)𝜓₁₂(𝑠)𝜓₂₁(𝑠)𝜓₂₂(𝑠)Ψ,0≤𝑠≤𝑅,₁₁(𝑠)=1det(Ψ,Φ)⁻¹_𝑛𝑠𝑔𝜓₂,𝜙₂cos𝜔𝑠+𝜓₁,𝜙₂,Ψsin𝜔𝑠₁₂(𝑠)=𝜔det(Ψ,Φ)⁻¹_𝑛𝑠𝑔𝜓₂,𝜙₂sin𝜔𝑠−𝜓₁,𝜙₂,Ψcos𝜔𝑠₂₁(𝑠)=𝜔det(Ψ,Φ)⁻¹_𝑛𝑠𝑔𝜓₂,𝜙₁cos𝜔𝑠+𝜓₁,𝜙₁,Ψcos𝜔𝑠₂₂(𝑠)=𝜔det(Ψ,Φ)⁻¹_𝑛𝑠𝑔𝜓₂,𝜙₁sin𝜔𝑠−𝜓₁,𝜙₁,()cos𝜔𝑠detΨ,Φ⁻¹_𝑛𝑠𝑔=𝜓₁,𝜙₁𝜓₂,𝜙₂−𝜓₁,𝜙₂𝜓₂,𝜙₁, (27) where the substitution of the new elements (𝜓_𝑗(𝑠),𝜙_𝑘(𝜃)), 𝑗,𝑘=1,2 into (3) will lead to the identity matrix (Ψ,Φ)_id=1det(Ψ,Φ)⁻¹_𝑛𝑠𝑔det(Ψ,Φ)⁻¹_𝑛𝑠𝑔00det(Ψ,Φ)⁻¹_𝑛𝑠𝑔=.1001 (28)

Consequently, the constant matrices 𝐵∈𝒞 and 𝐵∈𝒞 are equivalent, and the elements of 𝐵 at the Hopf bifurcation are 𝐵=[[0,𝜔]^⊤,[−𝜔,0]^⊤].

With the algebraic simplifications 𝑒^𝐵𝜃=cos𝜔𝜃−sin𝜔𝜃sin𝜔𝜃cos𝜔𝜃, (29) one can easily see that Φ(𝜃)=Φ(0)𝑒^𝐵𝜃=cos𝜔𝜃−sin𝜔𝜃𝜔sin𝜔𝜃𝜔cos𝜔𝜃,−𝑅≤𝜃≤0, (30) and it follows that 𝐽(𝑡,𝜇)Φ(𝜃)=Φ(0)𝑒^{𝐵(𝜃+𝑡)} is indeed the solution operator of the linearized delay equations in 𝒞. Similarly, we have Ψ(𝑠)=Ψ(0)𝑒⁻𝐵_𝑠, 0≤𝑠≤𝑅, and 𝐽(̂𝑡,̂𝑅)Ψ(𝑠)=Φ(0)𝑒^{𝐵(𝜃+𝑡)} is the solution operator for the adjoint delay equations in 𝒞.

By the change of variable 𝑥^𝑃_𝑡(𝜃)=Φ(𝜃)𝑧(𝑡)+𝑥^𝑄_𝑡(𝜃) with 𝑧(𝑡)∈ℜ²,𝑧(𝑡)=(Ψ(𝑠),𝜙^𝑃(𝜃)), we obtain the following as a first-order approximation in 𝜀, for 𝜃=−𝜏₁,−𝜏₂: 𝑥₁𝑡−𝜏₁𝑥₂𝑡−𝜏₁=𝑧₁(𝑡)cos𝜔𝜏+𝑧₂(𝑡)sin𝜔𝜏₁−𝑧₁(𝑡)𝜔sin𝜔𝜏₁+𝑧₂(𝑡)𝜔cos𝜔𝜏₁,𝑥₁𝑡−𝜏₂𝑥₂𝑡−𝜏₂=𝑧₁(𝑡)cos𝜔𝜏+𝑧₂(𝑡)sin𝜔𝜏₂−𝑧₁(𝑡)𝜔sin𝜔𝜏₂+𝑧₂(𝑡)𝜔cos𝜔𝜏₂, (31) and for 𝜃=0, 𝑥₁(𝑡)𝑥₂(𝑡)=𝑧₁(𝑡)𝜔𝑧₂(𝑡). (32) Therefore, we obtain (1) of the center manifold as follows: 𝑧̇₁(𝑡)=−𝜔𝑧₂+Δ𝑓⁽¹⁾𝑧₁,𝑧₂,𝑧̇,𝛾₂(𝑡)=𝜔𝑧₁+Δ𝑓⁽²⁾𝑧₁,𝑧₂,𝛾, (33) where the perturbation functions Δ𝑓⁽¹⁾(𝑧₁,𝑧₂,𝛾) and Δ𝑓⁽²⁾(𝑧₁,𝑧₂,𝛾) denote the following: Δ𝑓⁽¹⁾𝑧₁,𝑧₂,𝛾=𝛾𝑐₁₀𝑧₁+𝑐₁₀₀𝑧₂+𝑐₁₁𝑧³₁+𝑐₁₂𝑧²₁𝑧₂+𝑐₁₃𝑧₁𝑧²₂+𝑐₁₄𝑧³₂+𝑐₁₅𝑧²₁+𝑐₁₆𝑧₁𝑧₂+𝑐₁₇𝑧²₂,Δ𝑓⁽²⁾𝑧₁,𝑧₂,𝛾=𝛾𝑐₂₀𝑧₁+𝑐₂₀₀𝑧₂+𝑐₂₁𝑧³₁+𝑐₂₂𝑧²₁𝑧₂+𝑐₂₃𝑧₁𝑧²₂+𝑐₂₄𝑧³₂+𝑐₂₅𝑧²₁+𝑐₂₆𝑧₁𝑧₂+𝑐₂₇𝑧²₂. (34)

Here, the lengthy expressions of 𝑐_𝑖𝑗 are omitted here, and which in polar coordinate, 𝑧₁=𝑎sin𝜃, 𝑧₂=−𝑎cos𝜃, 𝜃=𝜔𝑡+𝜑 is equivalent to 𝑎̇(𝑡)=𝜈𝛾𝑎+Π⁽¹⁾𝜏₁,𝜏₂,𝛾𝑎³+○𝑎⁵,𝜏₁,𝜏₂,𝜑̇(𝑡),𝛾=𝜔+Π⁽²⁾𝜏₁,𝜏₂,𝛾𝑎²+○𝑎⁴,𝜏₁,𝜏₂,𝛾. (35) The so-called Poincaré-Lyapunov coefficient Π⁽¹⁾(𝜏₁,𝜏₂,𝛾) is independently determined by the following formula: Π⁽¹⁾𝑎,𝜏₁,𝜏₂=1,𝛾16Δ𝑓_𝑧⁽¹⁾_1𝑧1𝑧1𝑧₁,𝑧₂,𝛾+Δ𝑓_𝑧⁽¹⁾_1𝑧2𝑧2𝑧₁,𝑧₂,𝛾+Δ𝑓_𝑧⁽²⁾_1𝑧1𝑧2𝑧₁,𝑧₂,𝛾+Δ𝑓_𝑧⁽²⁾_2𝑧2𝑧2𝑧₁,𝑧₂+1,𝛾16𝜔Δ𝑓_𝑧⁽¹⁾_1𝑧2𝑧₁,𝑧₂×,𝛾Δ𝑓_𝑧⁽¹⁾_1𝑧1𝑧1𝑧₁,𝑧₂,𝛾+Δ𝑓_𝑧⁽¹⁾_2𝑧2𝑧₁,𝑧₂,𝛾−Δ𝑓_𝑧⁽²⁾_1𝑧2𝑧₁,𝑧₂×,𝛾Δ𝑓_𝑧⁽²⁾_1𝑧1𝑧₁,𝑧₂,𝛾+Δ𝑓_𝑧⁽²⁾_1𝑧2𝑧₁,𝑧₂,𝛾−Δ𝑓_𝑧⁽¹⁾_1𝑧1𝑧₁,𝑧₂,𝛾Δ𝑓_𝑧⁽²⁾_1𝑧1𝑧₁,𝑧₂,𝛾+Δ𝑓_𝑧⁽¹⁾_2𝑧2𝑧₁,𝑧₂,𝛾Δ𝑓_𝑧⁽²⁾_2𝑧2𝑧₁,𝑧₂,𝛾. (36) 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 Δ𝑓⁽¹⁾(𝑧₁,𝑧₂,𝛾), Δ𝑓⁽²⁾(𝑧₁,𝑧₂,𝛾), for example, Δ𝑓_𝑧⁽²⁾_1𝑧2={𝜕(𝜕Δ𝑓_𝑧⁽²⁾_1𝑧2/𝜕𝑧₁)/𝜕𝑧₂}_{(0,0,𝛾𝑐)}, are evaluated at the bifurcation point (𝑧₁,𝑧₂)→(0,0). Thus, the partial derivatives are given as follows: 𝜕Δ𝑓⁽¹⁾𝑧₁,𝑧₂,𝛾𝜕𝑧₁=𝑐₁₀𝛾+3𝑐₁₁𝑧²₁+2𝑐₁₅+𝑐₁₂𝑧₂𝑧₁+𝑐₁₆+𝑐₁₃𝑧₂𝑧₂,𝜕Δ𝑓⁽¹⁾𝑧₁,𝑧₂,𝛾𝜕𝑧₂=𝑐₁₀₀𝛾+𝑐₁₃𝑧²₁+𝑐₁₆+2𝑐₁₃𝑧₂𝑧₁+2𝑐₁₇+3𝑐₁₄𝑧₂𝑧₂,𝜕Δ𝑓⁽²⁾𝑧₁,𝑧₂,𝛾𝜕𝑧₁=𝑐₂₀𝛾+3𝑐₂₁𝑧²₁+2𝑐₂₅+𝑐₂₂𝑧₂𝑧₁+𝑐₂₆+𝑐₂₃𝑧₂𝑧₂,𝜕Δ𝑓⁽²⁾𝑧₁,𝑧₂,𝛾𝜕𝑧₂=𝑐₂₀₀𝛾+𝑐₂₃𝑧²₁+𝑐₂₆+2𝑐₂₃𝑧₂𝑧₁+2𝑐₂₇+3𝑐₂₄𝑧₂𝑧₂, (37)(38) from which we readily obtain the required values for the formula in (36) as 𝜕Δ³𝑓⁽¹⁾𝑧₁,𝑧₂,𝛾𝜕𝑧₁𝜕𝑧₁𝜕𝑧₁=6𝑐₁₁,𝜕Δ³𝑓⁽¹⁾𝑧₁,𝑧₂,𝛾𝜕𝑧₁𝜕𝑧₂𝜕𝑧₂=2𝑐₁₃,𝜕Δ²𝑓⁽¹⁾𝑧₁,𝑧₂,𝛾𝜕𝑧₁𝜕𝑧₂=2𝑐₁₂𝑧₁+𝑐₁₃𝑧₂+𝑐₁₆,𝜕Δ𝑓⁽¹⁾𝑧₁,𝑧₁,𝛾𝜕𝑧₁𝜕𝑧₁=23𝑐₁₁𝑧₁+𝑐₁₂𝑧₂+𝑐₁₅,𝜕Δ𝑓⁽¹⁾𝑧₂,𝑧₂,𝛾𝜕𝑧₂𝜕𝑧₂=2𝑐₁₃𝑧₁+3𝑐₁₄𝑧₂+𝑐₁₇,𝜕Δ³𝑓⁽²⁾𝑧₁,𝑧₂,𝛾𝜕𝑧₁𝜕𝑧₁𝜕𝑧₁=6𝑐₂₄,𝜕Δ³𝑓⁽²⁾𝑧₁,𝑧₂,𝛾𝜕𝑧₁𝜕𝑧₂𝜕𝑧₂=2𝑐₂₃,𝜕Δ²𝑓⁽²⁾𝑧₁,𝑧₂,𝛾𝜕𝑧₁𝜕𝑧₂=2𝑐₂₂𝑧₁+𝑐₂₃𝑧₂+𝑐₂₆,𝜕Δ𝑓⁽²⁾𝑧₁,𝑧₁,𝛾𝜕𝑧₁𝜕𝑧₁=23𝑐₂₁𝑧₁+𝑐₂₂𝑧₂+𝑐₂₅,𝜕Δ𝑓⁽²⁾𝑧₁,𝑧₂,𝛾𝜕𝑧₂𝜕𝑧₂=2𝑐₂₃𝑧₁+3𝑐₂₄𝑧₂+𝑐₂₇. (39) Substituting these quantities into the formulate in (36) will yield Π⁽¹⁾𝑎,𝜏₁,𝜏₂=1,𝛾83𝑐₁₁+𝑐₂₄+𝑐₁₃+𝑐₂₃+1𝑐8𝜔¹⁶𝑐₁₅+𝑐₁₇−𝑐₂₆𝑐₂₅+𝑐₂₇−𝑐₁₅𝑐₂₅+𝑐₁₇𝑐₂₇. (40) This coefficient together with the amplitude bifurcation equation (35) produces the nonlinear bifurcations for multiple delays as shown in the following: (a) Π⁽¹⁾(𝑎,𝜏₁,𝜏₂,𝛾)<0→ supercritical bifurcation; (b) Π⁽¹⁾(𝑎,𝜏₁,𝜏₂,𝛾)>0→ subcritical bifurcation.

Theorem 2. Let 𝑎𝑐+𝑏𝑑>0.(a) If Π⁽¹⁾(𝑎,𝜏₁,𝜏₂,𝛾)<0, then system (1) occurs supercritical bifurcation.(b) If Π⁽¹⁾(𝑎,𝜏₁,𝜏₂,𝛾)>0, then system (1) occurs subcritical bifurcation.

2.2. The Double Hopf Bifurcation

For the eigenvalue 𝜆_1,2=±𝑖𝜔₁, 𝜆_1,2=±𝑖𝜔₂ of the equation Δ(𝜆,𝛾)=0 in (6), we have the base function Φ(𝜃)=[𝜙₁(𝜃),𝜙₂(𝜃),𝜙₃(𝜃),𝜙₄(𝜃)], −𝜏≤𝜃≤0, Ψ(𝑠)=[𝜓₁(𝑠),𝜓₂(𝑠),𝜓₃(𝑠),𝜓₄(𝑠)]^⊤, 0≤𝑠≤𝜏 in 𝒞([−𝜏,0],ℜ^𝑛) and 𝒞([−𝜏,0],ℜ^𝑛), accordingly, Φ(𝜃)=cos𝜔₁𝜃−sin𝜔𝜃₁cos𝜔₂𝜃−sin𝜔₂𝜃sin𝜔₁𝜃cos𝜔₁𝜃sin𝜔₂𝜃cos𝜔₂𝜃,Ψ(𝑠)=cos𝜔₁𝑠−sin𝜔𝑠cos𝜔₂𝑠−sin𝜔₂𝑠sin𝜔₁𝑠cos𝜔₁𝑠sin𝜔₂𝑠cos𝜔₂𝑠^⊤, (41) and they form the inner product matrix (Ψ(𝑠),Φ(𝜃))=(𝜓_𝑗(𝑠)𝜙_𝑘(𝜃)),𝑗,𝑘=1,2,3,4, namely, (Ψ(𝑠)(𝜃=𝜓,Φ))₁(𝑠)𝜙₁(𝜃)𝜓₁(𝑠)𝜙₂(𝜃)𝜓₁(𝑠)𝜙₃(𝜃)𝜓₁(𝑠)𝜙₄(𝜃)𝜓₂(𝑠)𝜙₁(𝜃)𝜓₂(𝑠)𝜙₂(𝜃)𝜓₂(𝑠)𝜙₃(𝜃)𝜓₂(𝑠)𝜙₄(𝜃)𝜓₃(𝑠)𝜙₁(𝜃)𝜓₃(𝑠)𝜙₂(𝜃)𝜓₃(𝑠)𝜙₃(𝜃)𝜓₃(𝑠)𝜙₄(𝜃)𝜓₄(𝑠)𝜙₁(𝜃)𝜓₄(𝑠)𝜙₂(𝜃)𝜓₄(𝑠)𝜙₃(𝜃)𝜓₄(𝑠)𝜙₄(𝜃). (42) The substitution of the elements (Ψ(𝑠),Φ(𝜃))=(𝜓_𝑗(𝑠)𝜙_𝑘(𝜃)) into the bilinear pairing (5) and integrating, we obtain the nonsingular matrix (Ψ,Φ)_𝑛𝑠𝑔=𝜓₁₁−𝜓₁₂𝜓00₂₁𝜓₂₂0000𝜓₂₃−𝜓₃₄00𝜓₄₃𝜓₄₄, (43) where the constant values of this matrix are given by 𝜓₁₁=𝜓₂₂=1−𝛾𝜏₁cos𝜔₁𝜏₁,𝜓₁₂=−𝛾𝜏₁sin𝜔₁𝜏₁,𝜓₁₂=−𝜓₂₁,𝜓₃₃=𝜓₄₄=1−𝛾𝜏₁cos𝜔₂𝜏₁,𝜓₃₄=−𝛾𝜏₁sin𝜔₂𝜏₁,𝜓₄₃=−𝜓₃₄. (44) Then, Ψ∈𝒞 is normalized to the new basis Ψ=[𝜓₁(𝑠),𝜓₃(𝑠),𝜓₃(𝑠),𝜓₄(𝑠)], 0≤𝑠≤𝜏, where the elements 𝜓_𝑗(𝑠), 𝑗=1,2,3,4 of Ψ∈𝒞 are given by 𝜓₁=1𝜓²₁₁+𝜓²₁₂𝜓²₂₂cos𝜔₁𝑠+𝜓²₂₁sin𝜔₁𝑠,𝜓²₂₂sin𝜔₁𝑠−𝜓²₂₁cos𝜔₁𝑠,0,0^⊤,𝜓₂=1𝜓²₁₁+𝜓²₁₂𝜓²₁₂cos𝜔₁𝑠−𝜓²₁₁sin𝜔₁𝑠,𝜓²₁₂sin𝜔₁𝑠+𝜓²₁₁cos𝜔₁𝑠,0,0^⊤,𝜓₃=1𝜓²₃₃+𝜓²₃₄0,0,𝜓²₄₄cos𝜔₂𝑠+𝜓²₄₃sin𝜔₂𝑠,𝜓²₄₄sin𝜔₂𝑠−𝜓²₄₃cos𝜔₂𝑠^⊤,𝜓₄=1𝜓²₁₁+𝜓²₁₂0,0,𝜓²₃₄cos𝜔₂𝑠−𝜓²₃₃sin𝜔₂𝑠,𝜓²₃₄sin𝜔₂𝑠+𝜓²₃₃cos𝜔₁𝑠^⊤. (45) Again, the substitution of the elements 𝜓_𝑗(𝑠), 𝜙_𝑘(𝜃), 𝑗,𝑘=1,2,3,4 of (Ψ(𝑠),Φ(𝜃)) will yield the identity matrix, namely, (Ψ,Φ)_𝑛𝑠𝑔=1Λ𝜓Λ0000Λ0000Λ0000Λ=𝐼,Λ=²₁₁+𝜓²₁₂𝜓²₃₃+𝜓²₃₄. (46) Then, the constant matrix 𝐵 is given by 𝐵=0−𝜔₁𝜔00₁000000−𝜔₂00𝜔₂0. (47) The change of variable 𝑥^𝑃_𝑡(𝜙(𝜃),𝛾)=Φ(𝜃)𝑧(𝑡)+𝑥^𝑄_𝑡(𝜙(𝜃),𝛾), 𝑧(𝑡)∈ℜ⁴,𝑧(𝑡)=(Ψ(𝜃),𝜙^𝑃(𝜃)) will yield Φ(𝜃)𝑧(𝑡)=cos𝜔₁𝜃−sin𝜔₁𝜃cos𝜔₂𝜃−sin𝜔₂𝜃sin𝜔₁𝜃cos𝜔₁𝜃sin𝜔₂𝜃cos𝜔₂𝜃×𝑧₁(𝑡)𝑧₂(𝑡)𝑧₃(𝑡)𝑧₄(𝑡), (48) thereby, for −𝜏≤𝜃<0 and 𝜏₁≤𝜏≤𝜏₂, we have Φ−𝜏₁𝑧(𝑡)=cos𝜔₁𝜏₁sin𝜔₁𝜏₁cos𝜔₂𝜏₁sin𝜔₂𝜏₁−sin𝜔₁𝜏₁cos𝜔₁𝜏₁sin𝜔₂𝜏₁cos𝜔₂𝜏₁×𝑧₁(𝑡)𝑧₂(𝑡)𝑧₃(𝑡)𝑧₄(𝑡),Φ−𝜏₂𝑧(𝑡)=cos𝜔₁𝜏₂sin𝜔₁𝜏₂cos𝜔₂𝜏₂sin𝜔₂𝜏₂−sin𝜔₁𝜏₂cos𝜔₁𝜏₂−sin𝜔₂𝜏₂cos𝜔₂𝜏₂×𝑧₁(𝑡)𝑧₂(𝑡)𝑧₃(𝑡)𝑧₄(𝑡),Φ(0)𝑧(𝑡)=𝑧10100101¹(𝑡)𝑧₂(𝑡)𝑧₃(𝑡)𝑧₄(𝑡), (49) and so 𝛾𝑥₁𝑡−𝜏₁=𝛾𝑧₁cos𝜔₁𝜏₁+𝑧₂sin𝜔₁𝜏₁+𝑧₃cos𝜔₂𝜏₁+𝑧₄sin𝜔₂𝜏₁,𝛾𝑥₁𝑡−𝜏₁𝑥₂𝑡−𝜏₂=𝛾𝑧₁cos𝜔₁𝜏₁+𝑧₂sin𝜔₁𝜏₁+𝑧₃cos𝜔₂𝜏₁+𝑧₄sin𝜔₂𝜏₁×𝑧₁cos𝜔₁𝜏₂+𝑧₂sin𝜔₁𝜏₂+𝑧₃cos𝜔₂𝜏₂+𝑧₄sin𝜔₂𝜏₂. (50) The ODEs of the center manifold 𝑀_𝛾∈𝒞([−𝜏,0],ℜ⁴) are given by as follows: 𝑧̇₁(𝑡)=−𝜔₁𝑧₂+Δ𝑔¹𝑧₁,𝑧₂,𝑧₃,𝑧₄,𝑧̇,𝛾₂(𝑡)=𝜔₁𝑧₁+Δ𝑔²𝑧₁,𝑧₂,𝑧₃,𝑧₄,𝑧̇,𝛾₃(𝑡)=−𝜔₂𝑧₂+Δ𝑔³𝑧₁,𝑧₂,𝑧₃,𝑧₄,𝑧̇,𝛾₄(𝑡)=𝜔₂𝑧₂+Δ𝑔⁴𝑧₁,𝑧₂,𝑧₃,𝑧₄,𝛾, (51) where the perturbations in these equations denote Δ𝑔¹𝑧₁,𝑧₂,𝑧₃,𝑧₄,𝛾=𝑐⁽¹⁾₁₁₁𝑧³₁+𝑐⁽¹⁾₁₁₂𝑧²₁𝑧₂+𝑐⁽¹⁾₁₁₃𝑧²₁𝑧₃+𝑐⁽¹⁾₁₁₄𝑧²₁𝑧₂+𝑐⁽¹⁾₁₁₅𝑧₁𝑧²₂+𝑐⁽¹⁾₁₁₆𝑧₁𝑧²₃+𝑐⁽¹⁾₁₁₇𝑧₁𝑧²₄+𝑐⁽¹⁾₂₂₂𝑧³₂+𝑐⁽¹⁾₃₃₃𝑧³₃+𝑐⁽¹⁾₄₄₄𝑧³₄+𝑑⁽¹⁾₁₁₁𝑧²₁+𝑑⁽¹⁾₁₁₂𝑧₁𝑧₂+𝑑⁽¹⁾₁₁₃𝑧₁𝑧₃+𝑑⁽¹⁾₁₁₄𝑧₁𝑧₄+𝑑⁽¹⁾₂₂₂𝑧²₂+𝑑⁽¹⁾₃₃₃𝑧²₃+𝑑⁽¹⁾₄₄₄𝑧²₄+⋯,Δ𝑔²𝑧₁,𝑧₂,𝑧₃,𝑧₄,𝛾=𝑐⁽²⁾₁₁₁𝑧³₁+𝑐⁽²⁾₁₁₂𝑧²₁𝑧₂+𝑐⁽²⁾₁₁₃𝑧²₁𝑧₃+𝑐⁽²⁾₁₁₄𝑧²₁𝑧₂+𝑐⁽²⁾₁₁₅𝑧₁𝑧²₂+𝑐⁽²⁾₁₁₆𝑧₁𝑧²₃+𝑐⁽²⁾₁₁₇𝑧₁𝑧²₄+𝑐⁽²⁾₂₂₂𝑧³₂+𝑐⁽¹⁾₃₃₃𝑧³₃+𝑐⁽¹⁾₄₄₄𝑧³₄+𝑑⁽²⁾₁₁₁𝑧²₁+𝑑⁽²⁾₁₁₂𝑧₂𝑧₂+𝑑⁽²⁾₁₁₃𝑧₁𝑧₃+𝑑⁽²⁾₁₁₄𝑧₁𝑧₄+𝑑⁽²⁾₂₂₂𝑧²₂+𝑑⁽²⁾₃₃₃𝑧²₃+𝑑⁽²⁾₄₄₄𝑧²₄+⋯,Δ𝑔³𝑧₁,𝑧₂,𝑧₃,𝑧₄,𝛾=𝑐⁽³⁾₁₁₁𝑧³₁+𝑐⁽³⁾₁₁₂𝑧²₁𝑧₂+𝑐⁽³⁾₁₁₃𝑧²₁𝑧₃+𝑐⁽³⁾₁₁₄𝑧²₁𝑧₂+𝑐⁽¹⁾₁₁₅𝑧₁𝑧²₂+𝑐⁽³⁾₁₁₆𝑧₁𝑧²₃+𝑐⁽³⁾₁₁₇𝑧₁𝑧²₄+𝑐⁽³⁾₂₂₂𝑧³₂+𝑐⁽³⁾₃₃₃𝑧³₃+𝑐⁽¹⁾₄₄₄𝑧³₄+𝑑⁽³⁾₁₁₁𝑧²₁+𝑑⁽³⁾₁₁₂𝑧₁𝑧₂+𝑑⁽³⁾₁₁₃𝑧₁𝑧₃+𝑑⁽³⁾₁₁₄𝑧₁𝑧₄+𝑑⁽³⁾₂₂₂𝑧²₂+𝑑⁽³⁾₃₃₃𝑧²₃+𝑑⁽³⁾₄₄₄𝑧²₄+⋯,Δ𝑔⁴𝑧₁,𝑧₂,𝑧₃,𝑧₄,𝛾=𝑐⁽⁴⁾₁₁₁𝑧³₁+𝑐⁽⁴⁾₁₁₂𝑧²₁𝑧₂+𝑐⁽⁴⁾₁₁₃𝑧²₁𝑧₃+𝑐⁽⁴⁾₁₁₄𝑧²₁𝑧₂+𝑐⁽⁴⁾₁₁₅𝑧₁𝑧²₂+𝑐⁽⁴⁾₁₁₆𝑧₁𝑧²₃+𝑐⁽⁴⁾₁₁₇𝑧₁𝑧²₄+𝑐⁽⁴⁾₂₂₂𝑧³₂+𝑐⁽⁴⁾₃₃₃𝑧³₃+𝑐⁽⁴⁾₄₄₄𝑧³₄+𝑑⁽⁴⁾₁₁₁𝑧²₁+𝑑⁽⁴⁾₁₁₂𝑧₁𝑧₂+𝑑⁽⁴⁾₁₁₃𝑧₁𝑧₃+𝑑⁽⁴⁾₁₁₄𝑧₁𝑧₄+𝑑⁽⁴⁾₂₂₂𝑧²₂+𝑑⁽⁴⁾₃₃₃𝑧²₃+𝑑⁽⁴⁾₄₄₄𝑧²₄+⋯, (52) and the above coefficients in these equations denote the constant values in (52). Here, the lengthy expressions of 𝑐^𝑘_𝑖𝑗𝑙 and 𝑑^𝑘_𝑖𝑗𝑙 are omitted here. These equations are further written into the desired standard form of amplitude and phase relations using the polar coordinates 𝑧_𝑗=𝑎_𝑗sin𝜃_𝑗, 𝑧_𝑗=𝑎_𝑗cos𝜃_𝑗, 𝜃_𝑗=𝜔_𝑗𝑡+𝜑_𝑗, 𝑗=1,2,3,4. Again, the application of the integral averaging method to the resulting amplitude and phase equations will yield 𝑑𝑎₁𝑑𝑡=𝛽⁽¹¹⁾₁₁₁𝑎₁+Π⁽¹¹⁾𝑎,𝛾𝑎³₁+𝛽⁽¹¹⁾₁₁₂𝑎₁𝑎²₂𝑑+⋯,𝜑₁𝑑𝑡=𝛽⁽¹²⁾₁₁₁𝑎₁+Π⁽¹²⁾𝑎,𝛾𝑎³₁+𝛽⁽¹²⁾₁₁₂𝑎₁𝑎²₂𝑑+⋯,𝑎₂𝑑𝑡=𝛽⁽²¹⁾₁₁₁𝑎₂+Π⁽²¹⁾𝑎,𝛾𝑎³₂+𝛽⁽²¹⁾₁₁₂𝑎₁𝑎²₂𝑑+⋯,𝜑₂𝑑𝑡=𝛽⁽²²⁾₁₁₁𝑎₁+Π⁽²²⁾𝑎,𝛾𝑎³₂+𝛽⁽²²⁾₁₁₂𝑎₁𝑎²₂+⋯, (53) 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 Π⁽¹¹⁾(𝑎,𝛾), Π⁽²¹⁾(𝑎,𝛾), 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 Π(𝑎,𝛾) in (36) for the single Hopf, there are also explicit formulas for computing the double Hopf coefficients Π⁽¹¹⁾(𝑎,𝛾), Π⁽²¹⁾(𝑎,𝛾), 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 𝛼_𝐷 satisfying in any neighborhood of 𝛼_𝐷, there exists another constant 𝛼 and the corresponding invariant measure 𝜈_𝛼≠𝜇_𝛼 satisfying 𝜈_𝛼→𝜇_𝛼 as 𝛼→𝛼_𝐷. Then, the constant 𝛼_𝐷 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 𝛼₀ satisfying in any neighborhood of 𝛼_𝐷, there exist other two constants 𝛼₁, 𝛼₂, and their corresponding stationary density functions 𝑝_𝛼1, 𝑝_𝛼2 satisfying 𝑝_𝛼1 and 𝑝_𝛼2 are not equivalent. Then, the constant 𝛼₀ 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 𝜇₀ that satisfies the following conditions:(i)when 𝜇<𝜇₀, the stochastic differential equation has only one invariant measure 𝜈₀; moreover, 𝜈₀ is stable (i.e., the maximal Lyapunov exponent is negative).(ii)when 𝜇=𝜇₀, the invariant measure 𝜈₀ loses its stability and becomes unstable (i.e., the maximal Lyapunov exponent is positive); moreover, the rotation number is zero.(iii)when 𝜇>𝜇₀, the stochastic differential equation has three invariant measures 𝜈₀, 𝜈₁, and 𝜈₂; both 𝜈₁ and 𝜈₂ 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 𝜇=𝜇₀.

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 solution 𝑥(𝑡,𝑡₀,𝑥₀) 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 𝛿=𝛿(𝜀,𝛿,𝑡₀)>0 such that 𝑃𝑥𝑡,𝑡₀,𝑥₀<𝛼∀𝑡≥𝑡₀≥1−𝜀, (54) whenever |𝑥₀|<𝛿. 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 𝛾→𝛾+𝛼𝜉(𝑡). 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: 𝑢̇(𝑡)=𝛼𝜅(1−𝑢−𝑣)−𝛽𝑢−𝛾𝑢𝑡−𝜏₁𝑣𝑡−𝜏₂+𝜎₁+𝜎₂𝑢+𝜎₃𝑣𝜉(𝑡),𝑣̇(𝑡)=𝛿(1−𝜅)(1−𝑢−𝑣)³+𝜎−𝛾𝑢𝑣₄+𝜎₅𝑢+𝜎₆𝑣𝜉(𝑡). (55)

Here, 𝜉(𝑡) is the multiplicative random excitation and the external random excitation directly, and 𝜉(𝑡) is independent, in possession of zero mean value and standard variance Gaussian white noises, that is, ℰ[𝜉(𝑡)]=0,ℰ[𝜉(𝑡)𝜉(𝑡+𝜏)]=𝛿(𝜏). Note that the intensity of the noise 𝜎_𝑖(𝑖=1,2,…,6), 𝛾=𝛾_𝑐+𝛾 is the bifurcation parameter. From (3)–(33), we obtain (55) of the stochastic center manifold: 𝑧̇₁(𝑡)=−𝜔𝑧₂+𝑐₁₀𝑧₁+𝑐₁₀₀𝑧₂+𝑐₁₁𝑧³₁+𝑐₁₂𝑧²₁𝑧₂+𝑐₁₃𝑧₁𝑦²₂+𝑐₁₄𝑧³₂+𝑐₁₅𝑧²₁+𝑐₁₆𝑧₂𝑧₁+𝑐₁₇𝑧²₂+𝑘₁₁𝑧₁+𝑘₁₂𝑧₂+𝑘₁₃𝜉(𝑡),𝑧̇₂(𝑡)=𝜔𝑧₁+𝑐₂₀𝑧₁+𝑐₂₀₀𝑧₂+𝑐₂₁𝑧³₁+𝑐₂₂𝑧²₁𝑧₂+𝑐₂₃𝑧₁𝑧²₂+𝑐₂₄𝑧³₂+𝑐₂₅𝑧²₁+𝑐₂₆𝑧₂𝑧₁+𝑐₂₇𝑧²₂+𝑘₂₁𝑧₁+𝑘₂₂𝑧₂+𝑘₂₃𝜉(𝑡). (56) Here, the lengthy expressions of 𝑐_𝑖𝑗 and 𝑘_𝑖𝑗 are omitted.

Setting the coordinate transformation 𝑧₁=𝑟cos𝜃, 𝑧₂=𝑟sin𝜃 and by substituting the variable in (34), we obtain 𝑟̇(𝑡)=𝑟𝑐₁₀cos²𝜃+𝑟𝑐₁₀₀+𝑐₂₀sin𝜃cos𝜃+𝑟𝑐₂₀₀sin²𝜃+𝑟³𝑐₁₁cos⁴𝜃+𝑐₁₂+𝑐₂₁cos³+𝑐𝜃sin𝜃₁₃+𝑐₂₂cos²𝜃sin²𝜃+𝑐₁₄+𝑐₂₃cos𝜃sin³𝜃+𝑐₂₄sin⁴𝜃+𝑟²𝑐₁₅cos³𝜃+𝑐₁₆+𝑐₂₅cos²+𝑐𝜃sin𝜃₁₇+𝑐₂₆cos𝜃sin²𝜃+𝑐₂₇sin³𝜃+𝑟𝑘₁₁cos²𝜃+𝑘₁₂+𝑘₂₁cos𝜃sin𝜃+𝑘₂₂sin²𝜃𝜉(𝑡)+𝑘₁₃cos𝜃+𝑘₂₃𝜉(𝑡),𝜃̇(𝑡)sin𝜃=−𝜔+𝑐₂₀cos²𝜃+𝑐₂₀₀−𝑐₁₀sin𝜃cos𝜃−𝑐₁₀₀sin²𝜃+𝑟²𝑐₂₁cos⁴𝜃+𝑐₂₁−𝑐₁₂cos³+𝑐𝜃sin𝜃₂₃+𝑐₁₂cos²𝜃sin²𝜃+𝑐₂₄+𝑐₂₃cos𝜃sin³𝜃+𝑐₁₄sin⁴𝜃+𝑟𝑏₂₅cos𝜃+𝑐₂₆−𝑐₁₅cos²+𝑐𝜃sin𝜃₂₇−𝑐₁₆cos𝜃sin²𝜃−𝑐₁₇sin³𝜃+𝑘₂₁cos²𝜃+𝑘₂₂−𝑘₁₁cos𝜃sin𝜃−𝑘₁₂sin²𝜃𝜉(𝑡)+1𝑟𝑘₂₃cos𝜃−𝑘₁₃sin𝜃𝜉(𝑡). (57)

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 𝜎_𝑖(𝑖=1,…,6) are small enough, the response process {𝑟(𝑡),𝜃(𝑡)} 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: 𝑑𝑟=𝜇₁+𝜇₂8𝜇𝑟+³𝑟+𝜇₇8𝑟³+𝜇𝑑𝑡₃+𝜇₄8𝑟^21/2𝑑𝑊_𝑟+𝑟𝜇₅^1/2𝑑𝑊_𝜃,𝑑𝜃=𝜇₁₀+𝜇₈8𝑟²𝑑𝑡+𝑟𝜇₅^1/2𝑑𝑊_𝑟+𝜇₃𝑟²+𝜇₆8^1/2𝑑𝑊_𝜃, (58) where 𝑊_𝑟(𝑡) and 𝑊_𝜃 are the independent and standard Wiener processes. Set the parameter as follows: 𝜇₁=12𝑐₁₀+𝑐₂₀₀,𝜇₁₀1=−𝜔+2𝑐₂₀+𝑐₁₀₀,𝜇₂=5𝑘²₁₁+5𝑘²₂₂+3𝑘²₂₁+3𝑘²₁₂+6𝑘₁₂𝑘₂₁−2𝑘₁₁𝑘₂₂,𝜇₃=12𝑘²₁₃+𝑘²₂₃,𝜇₄=3𝑘²₁₁+3𝑘²₂₂+𝑘²₁₂+𝑘²₂₁+2𝑘₁₂𝑘₂₁+2𝑘₁₁𝑘₂₂,𝜇₅=14𝑘₁₁+𝑘₂₂𝑘₂₁−𝑘₁₂,𝜇₆=𝑘²₁₁+𝑘²₂₂+3𝑘²₁₂+3𝑘²₂₁−2𝑘₁₂𝑘₂₁−2𝑘₁₁𝑘₂₂,𝜇₇=3𝑐₁₁+3𝑐₂₁+𝑐₁₃+𝑐₂₂,𝜇₈=3𝑐₂₁−3𝑐₁₄+𝑐₂₃−𝑐₁₂. (59)

3.1. Stochastic D-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 𝜇₃=0,𝜇₇=0. Then, system (55) undergoes stochastic D-bifurcation.

Proof. When 𝜇₃=0,𝜇₇=0, Then, system (58) becomes 𝑑𝑟=𝜇₁+𝜇₂8𝑟𝜇𝑑𝑡+⁴8𝑟^21/2𝑑𝑊_𝑟. (60) When 𝜇₄=0, (60) is a deterministic system, and there is no bifurcation phenomenon. Here, we discuss the situation 𝜇₄≠0, let 𝑚(𝑟)=𝜇₁+𝜇₂8−𝜇₄(𝑟)=𝜇16𝑟,𝜎⁴8^1/2𝑟. (61) The continuous random dynamical system generated by (60) is 𝜑(𝑡)𝑥=𝑥+^𝑡₀𝑚𝜑(𝑠)𝑥𝑑𝑠+^𝑡₀𝜎𝜑(𝑠)𝑥∘𝑑𝑊_𝑟, (62) where ∘𝑑𝑊_𝑟 is the differential in the sense of Statonovich, and it is the unique strong solution of (60) with the initial value 𝑥. And 𝑚(0)=0, 𝜎(0)=0, so 0 is a fixed point of 𝜑. Since 𝑚(𝑟) is bounded and for any 𝑟≠0, it satisfies the ellipticity condition 𝜎(𝑟)≠0. This ensures that there is at most one stationary probability density. According to the Itô equation of amplitude 𝑟(𝑡), we obtain its Fokker-Planck-Kolmogorov equation corresponding to (60) as follows: 𝜕𝑝𝜕𝜕𝑡=−𝜇𝜕𝑟¹+𝜇₂8𝑟𝑝+𝜕²𝜕𝑟²𝜇₄8𝑟²𝑝. (63) Let 𝜕𝑝/𝜕𝑡=0, then we obtain the solution of system (63) 𝑝(𝑡)=𝑐𝜎⁻¹(𝑡)exp^𝑡₀2𝑚(𝑢)𝜎²(𝑢)𝑑𝑢. (64)
The above equation (63) has two kinds of equilibrium states: fixed point and nonstationary motion. The invariant measure of the former is 𝛿₀, and its probability density is 𝛿_𝑥. 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): 𝑟(𝑡)=𝑟(0)exp^𝑡₀𝑚^′(0)+𝜎(0)𝜎^′(0)2+𝑑𝑠_𝑡0𝜎^′(0)𝑑𝑊_𝑟. (65) The Lyapunov exponent with regard to 𝜇 of the dynamical system 𝜑 is defined as: 𝜆_𝜑𝜇=lim_𝑡→+∞1𝑡ln‖𝑟(𝑡)‖, (66) substituting (65) into (66), note that 𝜎(0)=0, 𝜎^′(0)=0, we obtain the Lyapunov exponent of the fixed point: 𝜆_𝜑𝛿₀=lim_𝑡→+∞1𝑡ln‖𝑟(0)‖+𝑚^′(0)^𝑡₀𝑑𝑠+𝜎^′(0)^𝑡₀𝑑𝑊_𝑟(𝑠)=𝑚^′(0)+𝜎^′(0)lim_𝑡→+∞𝑊_𝑟(𝑡)𝑡=𝑚^′(0)=𝜇₁+𝜇₂8−𝜇₄.16 (67) For the invariant measure which regard (65) as its density, we obtain the Lyapunov exponent: 𝜆_𝜑(𝜈)=lim_𝑡→+∞1𝑡^𝑡₀𝑚^′(𝑟)+𝜎(𝑟)𝜎^′(𝑟=)𝑑𝑠_𝑅𝑚^′(𝑟)+𝜎(𝑟)𝜎^′(𝑟)2𝑝(𝑟)𝑑𝑟=−2_𝑅𝑚(𝑟)𝜎(𝑟)²𝑝(𝑟)√𝑑𝑟=−322𝜇₄^3/2𝑚(𝑟)²exp16𝜇₄𝑚(𝑟√)=−322𝜇₄^3/2𝜇₁+𝜇₂8−𝜇₄16²×exp16𝜇₄𝜇₁+𝜇₂8−𝜇₄.16 (68)
Let 𝛼=𝜇₁+𝜇₂/8−𝜇₄/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 𝛼=𝛼_𝐷=0 is a point of 𝐷-bifurcation.

Theorem 9. Let 𝜇₃=0, 𝜇₇=0. Then, system (55) does not undergo stochastic P-bifurcation.

Proof. By simplifying (64), we can obtain 𝑝_st(𝑟)=𝑐𝑟^{2(8𝜇1+𝜇2−𝜇4)/𝜇4}, (69) where 𝑐 is a normalization constant; thus, we have 𝑝_st(𝑟)=𝑜𝑟^𝑣𝑟⟶0, (70) where 𝑣=2(8𝜇₁+𝜇₂−𝜇₄)/𝜇₄. Obviously, when 𝑣<−1, that is, 𝜇₁+𝜇₂/8−𝜇₄/16<0, 𝑝_st(𝑟) is a 𝛿 function. When −1<𝑣<0, that is, 𝜇₁+𝜇₂/8−𝜇₄/16>0, 𝑟=0 is a maximum point of 𝑝_st(𝑟) in the state space; thus, the system admits D-bifurcation when 𝑣=−1, that is, 𝜇₁+𝜇₂/8−𝜇₄/16=0, is the critical condition of   D-bifurcation at the equilibrium point. When 𝑣>0, there is no point that makes 𝑝_st(𝑟) 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, 𝜇₃=0,𝜇₇<0. Then, system (55) undergoes stochastic pitchfork bifurcation.

Proof. When 𝜇₃=0,𝜇₇≠0, then (58) can be rewritten as 𝑑𝑟=𝜇₁+𝜇₂8𝜇𝑟+⁷8𝑟³𝜇𝑑𝑡+⁴8𝑟^21/2𝑑𝑊_𝑟. (71) Let 𝛼=𝜇₁+𝜇₂/8, 𝜙=(−𝜇₇/8)𝑟, 𝜇₇<0, then the system (71) becomes 𝑑𝜙=𝜇₁+𝜇₂8𝜙−𝜙³𝜇𝑑𝑡+⁴8^1/2𝜙∘𝑑𝑊_𝑡, (72) which is solved by 𝜙→𝜓_𝛼()𝜙=𝜇𝑡,𝜔𝜙exp₁+𝜇₂8𝜇𝑡+⁴8^1/2𝑊_𝑡(𝜔)×1+2𝜙^2𝑡₀exp2𝜇₁+𝜇₂8𝑡+𝜇₄8^1/2𝑊_𝑠(𝜔)𝑑𝑠^1/2−1. (73)
We now determine the domain 𝐷_𝛼(𝑡,𝜔), where 𝐷_𝛼(𝑡,𝜔)∶={𝜙∈ℜ∶(𝑡,𝜔,𝜙)∈𝐷}(𝐷=ℜ×Ω×𝑋) is (in general possibly empty) the set of initial values 𝜙∈ℜ for which the trajectories still exist at time 𝑡 and the range 𝑅_𝛼(𝑡,𝜔) of 𝜓_𝛼(𝑡,𝜔)∶𝐷_𝛼(𝑡,𝜔)→𝑅_𝛼(𝑡,𝜔).
We have 𝐷_𝛼(𝑡,𝜔)=ℜ,𝑡≥0,−𝑑_𝛼(𝑡,𝜔),𝑑_𝛼(𝑡,𝜔),𝑡<0, (74) where 𝑑_𝛼()2𝑡,𝜔=1×^𝑡₀exp2𝜇₁+𝜇₂8𝑡𝜇+2⁴8^1/2𝑊_𝑠(𝜔)𝑑𝑠^1/2−1𝑅>0,_𝛼(𝑡,𝜔)=𝐷_𝛼((𝑡)𝜔)=−𝑡,𝜗−𝑟_𝛼(𝑡,𝜔),𝑟_𝛼(𝑡,𝜔),𝑡>0,ℜ,𝑡≤0, (75) where 𝑟_𝛼(𝑡,𝜔)=𝑑_𝛼((𝑡)𝜔)=𝜇−𝑡,𝜗exp₁+𝜇₂8𝜇𝑡+⁴8^1/2𝑊_𝑡(𝜔)×2^𝑡₀exp2𝜇₁+𝜇₂8𝑡𝜇+2⁴8^1/2𝑊_𝑠(𝜔)𝑑𝑠^1/2−1>0. (76) We can now determine that 𝐸_𝛼(𝜔)∶=_𝑡∈ℜ𝐷_𝛼(𝑡,𝜔) (77) and obtain 𝐸_𝛼(𝜔)=−𝑑⁻_𝛼(𝑡,𝜔),𝑑⁻_𝛼(𝑡,𝜔),𝛼=𝜇₁+𝜇₂8{0}>0,,𝛼=𝜇₁+𝜇₂8≤0, (78) where 0<𝑑^±_𝛼()2𝑡,𝜔=1×₀^±∞exp2𝜇₁+𝜇₂8𝑡+𝜇₄8^1/2𝑊_𝑠(𝜔)𝑑𝑠^1/2−1<∞. (79) The ergodic invariant measures of system (71) are as follows.(i)For 𝛼≤0, the only invariant measure is 𝜇^𝛼_𝜔=𝛿₀.(ii)For 𝛼>0, we have the three invariant forward Markov measures, 𝜇^𝛼_𝜔=𝛿₀ and 𝜐^𝛼_±,𝜔=𝛿_{±𝑘𝛼(𝜔)}, where 𝑘_𝛼(𝜔)∶=2⁰_−∞exp2𝜇₁+𝜇₂8𝑡𝜇+2⁴8^1/2𝑊_𝑡(𝜔)𝑑𝑠^−1/2. (80)
We have 𝔼𝑘²_𝛼(𝜔)=𝛼. Solving the forward Fokker-Planck-Kolmogorov equation 𝐿^∗𝑝_𝛼=−𝜇₁+𝜇₂8𝜇𝜙+⁴16𝜙−𝜙³𝑃_𝛼𝜙^′+𝜇₄𝜙16²𝑃_𝜈1𝜙^′′=0, (81) yield(i)𝑝_𝛼=𝛿₀ for all 𝛼,(ii)for 𝑝_𝛼>0, 𝑞⁺_𝛼𝜙=𝑁_𝛼𝜙^{(2(𝜇1+𝜇2/8)/𝜇4)−1}exp16𝜙²𝜇₄,𝜙>0,0,𝜙≤0, (82)and 𝑞⁺_𝛼(𝜙)=𝑞⁺_𝛼(−𝜙), where 𝑁⁻_𝛼=Γ(𝜈₁/𝜇₄)(𝜇₄/8)^{(𝜇1+𝜇2/8)/𝜇4}.
Naturally, the invariant measures 𝜐^𝛼_±,𝜔=𝛿_{±𝑘𝛼(𝜔)} are those corresponding to the stationary measures 𝑞⁺_𝛼. Hence, all invariant measures are Markov measures.
We determine all invariant measures (necessarily Dirac measure) of local RDS 𝜒 generated by the SDE 𝑑𝜙=𝜇₁+𝜇₂8𝜙−𝜙³𝜇𝑑𝑡+⁴8^1/2𝜙∘𝑑𝑊, (83) on the state space ℜ, 𝜇₁+𝜇₂/8∈ℜ. We now calculate the Lyapunov exponent for each of these measures.
The linearized RDS 𝜒_𝑡=𝐷Υ(𝑡,𝜔,𝜙)𝜒 satisfies the linearized SDE 𝑑𝜒_𝑡=𝜇₁+𝜇₂8−3Υ𝑡,𝜔,𝜙²𝜒_𝑡𝜇𝑑𝑡+⁴8^1/2𝜒_𝑡∘𝑑𝑊. (84) Hence, 𝐷Υ𝑡,𝜔,𝜙𝜒=𝜒exp𝜇₁+𝜇₂8𝜇𝑡+⁴8^1/2𝑊_𝑡(𝜔)−3^𝑡₀Υ𝑠,𝜔,𝜙²𝑑𝑠. (85) Thus, if 𝜈_𝜔=𝛿_𝜙0(𝜔) is a Υ-invariant measure, its Lyapunov exponent is 𝜆𝜇=lim_𝑡→∞1𝑡𝜒=𝜇log𝐷Υ𝑡,𝜔,𝜙₁+𝜇₂8−3lim_𝑡→∞^𝑡₀Υ𝑠,𝜔,𝜙²=𝜇𝑑𝑠₁+𝜇₂8−3𝔼𝜙²₀, (86) providing the IC 𝜙²₀∈𝐿¹(ℙ) is satisfied.(i) For 𝜇₁+𝜇₂/8∈ℜ, the IC for 𝜈^𝛼_±,𝜔=𝛿₀ is trivially satisfied, and we obtain 𝜆𝜈^𝛼₁=𝜇₁+𝜇₂8. (87) So, 𝜈^𝛼₁ is stable for 𝜇₁+𝜇₂/8<0 and unstable for 𝜇₁+𝜇₂/8>0.(ii) For 𝜇₁+𝜇₂/8>0, 𝜈^𝛼_2,𝜔=𝛿_𝑑^𝛼_𝜔 is ℱ⁰_−∞ measurable; hence, the density 𝑝^𝛼 of 𝜌^𝛼=𝔼𝜈^𝛼₂ satisfies the Fokker-Planck-Kolmogorov equation 𝐿^∗_𝜈1=−𝜇₁+𝜇₂8𝜇𝜙+⁴16𝜙−𝜙³𝑝_𝛼𝜙^′+𝜇₄𝜙16²𝑝_𝛼𝜙^′′=0, (88)which has the unique probability density solution 𝑃^𝛼𝜙=𝑁_𝛼𝜙^{(2(𝜇1+𝜇2/8)/𝜇4)−1}𝜙exp²8𝜇₄,𝜙>0. (89) Since 𝔼_𝜈^𝛼₂𝜙²=𝔼𝑑^𝛼₋²=^∞₀𝜙²𝑝^𝛼𝜙𝑑𝜙<∞, (90) the IC is satisfied. The calculation of the Lyapunov exponent is accomplished by observing that 𝑑^𝛼₋𝜗_𝑡𝜔²=exp2𝜇₁+𝜇₂/8𝑡+2𝜇₄/8^1/2𝑊_𝑡(𝜔)2∫^𝑡_−∞exp2𝜇₁+𝜇₂/8𝑠+2𝜇₄/8^1/2𝑊_𝑠(𝜔)=Ψ𝑑𝑠_′(𝑡),2Ψ (91) where Ψ(𝑡)=^𝑡_−∞exp𝜇₁+𝜇₂8𝜇𝑠+⁴8^1/2𝑊_𝑠(𝜔)𝑑𝑠. (92) Hence, by the ergodic theorem 𝔼𝑑^𝛼₋²=12lim_𝑡→∞1𝑡logΨ(𝑡)=𝜇₁+𝜇₂8, (93) finally 𝜆𝜈^𝛼₂=−2𝜇₁+𝜇₂8<0. (94)(iii) For (𝜇₁+𝜇₂/8)>0, 𝜈^𝛼_2,𝜔=𝛿_𝑑^𝛼_𝜔 is ℱ⁰_−∞ measurable. Since ℒ(𝑑^𝛼₊)=ℒ(𝑑^𝛼₋)𝔼−𝑑^𝛼₋²=𝔼𝑑^𝛼₋²=𝜇₁+𝜇₂8, (95) thus 𝜆𝜈^𝛼₂=−2𝜇₁+𝜇₂8<0. (96)
The two families of densities (𝑞⁺_𝛼)_𝛼>0 clearly undergo a P-bifurcation at the parameter value 𝛼_𝑃=𝜇₄/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 𝛿₀ at 𝛼_𝐷=0 and a P-bifurcation of 𝛼_𝑃=(𝜇₁+𝜇₂/8)²/2. Then, system (55) has stochastic pitchfork bifurcation.

3.2. P-Bifurcation

In the following, we investigate the steady-state probability density 𝑝_st(𝑟) 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 𝜇₃=0,𝜇₇≠0, then (38) can be rewritten as 𝑑𝑟=𝜇₁+𝜇₂8𝜇𝑟+⁷8𝑟³𝜇𝑑𝑡+⁴8𝑟^21/2𝑑𝑊_𝑟. (97)

According to the Itô equation of amplitude 𝑟(𝑡), we obtain its Fokker-Planck-Kolmogorov equation form (97) as follows: 𝜕𝑝𝜕𝜕𝑡=−𝜇𝜕𝑟¹+𝜇₂8𝜇𝑟−⁴2𝜇₇𝑟−𝑟³𝑝(𝑟−𝜇₄)2𝜇₇𝜕²𝜕𝑟²𝑟²𝑝(𝑟), (98) with the initial value condition 𝜇₇=0,𝑝(𝑟,𝑡|𝑟₀,𝑡₀)→𝛿(𝑟−𝑟₀),𝑡→𝑡₀, where 𝑝(𝑟,𝑡|𝑟₀,𝑡₀) is the transition probability density of diffusion process 𝑟(𝑡). The invariant measure of 𝑟(𝑡) is the steady-state probability density 𝑝_st(𝑟) which is the solution of the degenerate system as follows: 𝜕0=−𝜇𝜕𝑟¹+𝜇₂8𝜇𝑟−⁴2𝜇₇𝑟−𝑟³𝑝(𝑟−𝜇₄)2𝜇₇𝜕²𝜕𝑟²𝑟²𝑝(𝑟). (99) Through calculation, we can obtain 𝑝_st(𝑟)=exp𝑟²𝜇₇/𝜇₄𝑟^{−1−(2(𝜇1+𝜇2/8)𝜇7/𝜇4)}Γ−𝜇₁+𝜇₂/8𝜇₇/𝜇₄−𝜇₇/𝜇₄^{(𝜇1+𝜇2/8)𝜇7/𝜇4}. (100)

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 of 𝑝_st(𝑟) approximate to show the behavior of the deterministic system. If the process 𝑟(𝑡) is ergodic, then 𝑝_st(𝑟) can be regarded as the time measurement for staying in the neighborhood of 𝑟(𝑡) according to the Oseledec ergodic theorem.

From the analysis above, if 𝑝_st(𝑟) has a maximum value at 𝑟^∗, the sample trajectory will stay for a longer time in the neighborhood of 𝑟^∗, that is, 𝑟^∗ is stable in the meaning of probability (with a bigger probability). If 𝑝_st(𝑟) has a minimum value (zero), it is just the opposite.

We now calculate the most possible amplitude 𝑟^∗ of system (97), that is, 𝑝_st(𝑟) has a maximum value at 𝑟^∗. So, we have 𝑑𝑝_st(𝑟)𝑑𝑟_𝑟=𝑟^∗𝑑=0,²𝑝_st(𝑟)𝑑𝑟²_𝑟=𝑟∗<0 (101) and the solution 𝑟=𝑟̃=(𝜈₁𝜈₂+4𝜈₃)/8𝜈₂. The probabilities and the positions of the Hopf bifurcation occur with different parameter.

Since 𝑑²𝑝_st(𝑟)𝑑𝑟²₌𝑟̃=2^{4+(3(𝜇1+𝜇2/8)𝜇2/𝜇3)}1exp8𝜇4+⁷𝜇₃×𝜇₇−𝜇₇𝜇₄^{−(𝜇1+𝜇2/8)𝜇7/𝜈4}×8𝜇₁+𝜇₂/8𝜇₇+4𝜇₄𝜇₇^{−2(𝜇1+𝜇2/8)𝜇7/4𝜇4}×Γ−𝜇₁+𝜇₂/8𝜇₇𝜇₄−𝜇₇𝜇₄^{(𝜇1+𝜇2/8)𝜇7/𝜇4−1}<0, (102) thus, what we need is 𝑟^∗=𝑟̃. 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 𝑟=𝑟̃, where 𝑥²₁+𝑥²₂=8𝜇₁+𝜇₂/8𝜇₇+4𝜇₄8𝜇₇. (103)

We now choose some values of the parameters in the equations, draw the graphics of 𝑝_st(𝑟) (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, 𝜎₁=0.1, 𝜎₂=0.1, 𝜎₃=0.2, 𝜎₄=0.3, 𝜎₅=0.1, 𝜎₆=0.5. After further calculations, we obtain 𝜇₁=−0.2,𝜇₂=9.5907,𝜇₃=0,𝜇₄=5.1183,𝜇₇=−20.4732, then 𝑝(𝑟)=0.28209𝑒^−4𝑟2. (104) What is more is that 𝑟̃=0.28209, where 𝑝_st(𝑟) has the maximum value.

(a)

(b)

(c)

Figure 1 

P-Bifurcation of system (55) with 𝜇1=−0.2,𝜇3=0,𝜇4=5.1183,𝜇7=−20.4732; (a) P-bifurcation of 𝑝𝛼(𝑟) at 𝛼𝑝<𝜎2/2,𝜇2=5.6, where 𝛼=𝜇1+𝜇2/8,𝜎2=−𝜇4/𝜇7; (b) P-bifurcation of 𝑝𝛼(𝑟) at 𝛼𝑝=𝜎2/2,𝜇2=9.5907, where 𝛼=(𝜇1+𝜇2/8),𝜎2=−𝜇4/𝜇7; (c) P-bifurcation of 𝑝𝛼(𝑟) at 𝛼𝑝>𝜎2/2,𝜇2=17.6, where 𝛼=𝜈1/8,𝜎2=−𝜈3/𝜈2.

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 𝛾−𝛾_𝑐, it can be readily seen that the addition of the noise gives rise to a shift in the critical bifurcation point 𝛾_𝑐 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.

Acknowledgments

This paper was supported by the National Natural Science Foundation of China (no. 11201089) and the Guangxi Natural Science Foundation (no. 2011GXNSFB018065).