# Double Hopf Bifurcation in Microbubble Oscillators with Delay Coupling

Authors: Jinbin Wang; Rui Zhang; Lifenq Ma
Journal: Advances in Mathematical Physics
 (2020)

---

## Abstract

Using center manifold reduction methodswe investigate the double Hopf bifurcation in the dynamics of microbubble with delay couplingwith main attention focused on nonresonant double Hopf bifurcation. We obtain the normal form of the system in the vicinity of the double Hopf point and classify the bifurcations in a two-dimensional parameter space near the critical point. Some numerical simulations support the applicability of the theoretical results. In particularwe give the explanation for some physical phenomena of the system using the obtained mathematical results.


---

## Body

## 1. Introduction

At present, owing to the modeling of many technological and biological problems the dynamical systems with delay, delayed differential equations (DDEs), or functional differential equations (FDEs) have been an active area of research [1–6]. Double Hopf refers to same critical value, and there exist two different pure imaginary eigenvalues. We can observe much interesting phenomena through double Hopf bifurcation, such as periodic and quasi-periodic motions, three dimensional invariant torus [7], period doubling [8, 9], homoclinic and heteroclinic connections, and chaos. Recently, the center manifold reduction (CMR) and multiple time scales (MTS) are widely useful techniques which are used by researchers to study bifurcation of DDEs.Bubble dynamics is an important branch of fluid mechanics which can be applied to some research fields. For example, the microbubbles filled with a drug as a carrier for local drug delivery can be improved drug utilization in medical field. The fact is that the dynamics of bubbles in a liquid are influencing each other via acoustic waves, which is a complex process described by time-delay system.Plesset [10], Plesset and Prosperetti [11] set up the famous Rayleigh–Plesset equation (1) which used to describe the microbubble system:(1)a˙−caa¨+32a˙2−a−3r+Δ−a˙3−a−1a2Δ=0here,Δ=ρ−1pa−p0, where ρ is the density of the liquid, and pa=κ4/3πa3−r is the pressure inside the bubble, and r is the adiabatic exponent of the gas. However, Equation (1) merely describes the dynamics of a single micro-bubble system. In the general cases, there are a lot of bubbles (micro-bubbles) with the coupling effect in the liquid. Thus, to better analyze of micro-bubble nature of the system, we need to study an optimum system for describing further. In 2010, Heckman and Rand [12] set up a bubble coupled system with time delay:(2)a˙−caa¨+32a˙2−a−3r+1−a˙3−3r−2a−3ra˙−2a˙=Pb˙t−τ,b˙−cbb¨+32b˙2−b−3r+1−b˙3−3r−2b−3rb˙−2b˙=Pa˙t−τ,wherer is positive-valued parameters, c is sound speed, and τ is the delay and P is a coupling coefficient. a and b are the radius of the microbubbles, respectively. a and b also denote, respectively, the two bubbles. Heckman et al. studied bifurcations of the in-phase manifold given by a=b,a˙=b˙ [13–16]. Recently, Heckman has studied Equation (2) and developed explanations for bifurcation structure of the in-phase manifold, However, these studies were limited in that he investigated only given by a=b,a˙=b˙. This work will extend the previous, by making use of the theory of center manifold [17] and the normal form method [18] to consider double Hopf bifurcation, which occur with more general cases.This article is organized as follows: In Section2, we discuss the occurring conditions for the existence of double Hopf bifurcation on Equation (2). In Section 3, we analyze the nonresonant double Hopf bifurcation for the system using the normal form method and the center manifold theory. We classify the bifurcations in a two-dimensional parameter space near the critical point and some numerical simulations support the applicability of the theoretical results in Section 4. Mathematical results explain the mechanical background of the model in Section 5. The concluding section contains a brief conclusion of this work.


## 2. Stability Analysis

After calculatingae,be=1,1 is an equilibrium point of Equation (2), we set a=xt+1 and b=yt+1, then Equation (2) can be written as(3)Py˙t−τ=x˙t−cxt+1x¨t+32x˙t2−xt+1−3r+1−x˙t3−3r−2xt+1−3rx˙t−2x˙t,Px˙t−τ=y˙t−cyt+1y¨t+32y˙t2−yt+1−3r+1−y˙t3−3r−2yt+1−3ry˙t−2y˙t.Further, we let(4)x˙t=ht,y˙t=gt,then Equation (3) can be written as(5)x˙t=ht,h˙t=Pgt−τ+2ht+ht3+3r−2xt+1−3rhtht−c1+xt+−3/2ht2+xt+1−3r−11+xt,y˙t=gt,g˙t=Pht−τ+2gt+gt3+3r−2yt+1−3rgtgt−c1+yt+−3/2gt2+yt+1−3r−11+yt.Clearly,0,0,0,0 is an equilibrium point of Equation (5). The characteristic equation of its corresponding linear system around the origin 0,0,0,0 is(6)cλ2+3rc+3rλ+Pλe−λτcλ2+3rc+3rλ−Pλe−λτ=0.Whenτ=0, the roots for Equation (6) are(7)λ1,2=12c−3r−P±3r−P2−12c2r,(8)λ3,4=12c−3r+P±3r+P2−12c2r.Clearly(9)Reλ1,2<0,Reλ3,4<0,3r>P,λ1,2=±3ri,Reλ3,4<0,3r=P,Reλ1,2>0,Reλ3,4<0,3r<P,.Whenτ≠0, it is essential to discuss the properties of the roots into two cases.Case 1.
Let±iω1ω1>0 be a root of cλ2+3rc+3rλ+Pλe−λτ=0, then ω1 satisfies(10)Pω1sinω1τ=cω12−3r,Pω1cosω1τ=−3rω1,
and the solutions to Equation (10) are then found to be(11)ω1±=P2−9r2+12c2r±P2−9r22c>0,
under the assumptionP>3r.(12)τk1±=1ω1±θ1+2kπ,k=0,1,2,…,
whereθ1=π−arcsincω1±2−3r/Pω1±. By defining τk1+k=0,1,2,3,… and τk1−k=0,1,2,3,…, we can obtain τk1+<τk+11+ and τk1−<τk+11−. Furthermore, we get τk1+<τk1− and τk+11+−τk1+<τk+11−−τk1− from ω1+<ω1−. Thus, there exists a positive integer j>0, yielding(13)0<τ01+<τ01−<τ11+<τ11−<⋯<τj−11−<τj1+<τj+11+<τj1−.Case 2.
Let±iω2ω2>0 be a root of cλ2+3rc+3rλ−Pλe−λτ=0. Likewise available as a conclusion(14)ω1±=ω2±,(15)τk2±=1ω2±θ2+2kπ,k=0,1,2,…,(16)0<τ02+<τ02−<τ12+<τ12−<⋯<τi−12−<τi2+<τi+12+<τi2−.
here,θ2=arcsin−cω2±2−3r/Pω2±,P>3r.After calculating, the corresponding transverse condition is the following:(17)signdaτdττ=τ0=1,τ0=τk1+,k=0,1,2,3,…,−1,τ0=τk1−,k=0,1,2,3,…,−1,τ0=τk2+,k=0,1,2,3,…,1,τ0=τk2−,k=0,1,2,3,…,From [18, 19] and above analysis, we can obtain the following theorem about the existence of Hopf bifurcation. If P>3r, system Equation (5) undergoes a Hopf bifurcation at its zero equilibrium when τ=τk1± or τ=τk2±, where k=0,1,2,3…In the following, we will discuss the occurring conditions for the existence of double Hopf bifurcation on Equation (5), and choose the P,τ as the parameters plane. Further, we can draw a P,τ diagram as shown in Figure 1. In fact, double Hopf bifurcation occurs at these points.Pc,τcis called a double Hopf bifurcation point, when τk1+=τk2−=τc is satisfied for some k=0,1,2,…Figure 1
The Hopf bifurcation curves for Equation (5) at the trivial equilibrium 0,0,0,0 when r=4/3,c=94, and k=0.From(18)τk1+=1ω1+θ1+2kπ=τk2−=1ω2−θ2+2kπ,we can obtain(19)ω1+:ω2−=P2−9r2+12c2r+P2−9r2:P2−9r2+12c2r−P2−9r2≜b1:b2,the double Hopf bifurcation is called resonance whenb1:b2∈z+ and it is called nonresonance if b1:b2 is  an irrational number [16, 17].


## 3. Normal Form for Double Hopf Bifurcation of in Coupled Microbubble with Nonresonance

In previous section, we have established some results of Hopf and double Hopf bifurcation and drew this conclusion by analyzing the stability of the origin after bifurcation. However, there is a direct way to approach the dynamic of near the non-resonance double Hopf bifurcation critical point by using the center manifold theory and the normal form method.To begin, we will try to transform the infinite dimensional dynamic system into a finite dimensional system with the method of center manifold theory. Rescaling the time byt→t/τ to normalize the delay so that system Equation (5) becomes(20)x˙t=τht,h˙t=τPgt−1+2ht+ht3+3r−2xt+1−3rhtht−c1+xt+−3/2ht2+xt+1−3r−11+xt,y˙t=τgt,g˙t=τPht−1+2gt+gt3+3r−2yt+1−3rgtgt−c1+yt+−3/2gt2+yt+1−3r−11+yt.Suppose that the system Equation (20) undergoes a double Hopf bifurcation at the critical point Pc,τc with two pairs of eigenvalues ±iω1+ and ±iω2−, we can choose(21)ηθ=τcA,θ=0,0,θ∈−1,0,−τcB,θ=−1with(22)A=0100−3r−3rc00000100−3r−3rc,B=0000000−Pc00000−Pc00,.Then, the linearized equation of Equation (20) at the trivial equilibrium is(23)X˙t=L0Xt,whereL0ϕ=∫−10dηθϕθ, ϕ∈C−1,0,R4, and the bilinear form on C∗×C∗standsforadjoint is(24)ψ,φ=ψ0φ0−∫−r0∫ξ=0θψξ−θdηθφξdξ,where(25)φθ=φ1θ,φ2θ,φ3θ,φ4θ∈C,ψs=ψ1s,ψ2s,ψ3s,ψ4sT∈C∗.DesignateΛ=±iτcω1+,±iτcω2−, then as discussed above, we know that the phase space C can be decomposed as C=PΛ⊕QΛ, where PΛ⊂C is the 4-dimensional center subspace spanned by the basis vectors of the linear operator L0 associated with the imaginary characteristic roots, and QΛ is the complement subspace of PΛ. Further, Equation (20) can be rewritten as an ordinary differential equation in the Banach space BC (see 10) of functions bounded; therefore, we can suppose that the bases for PΛ and its adjoint PΛ∗ are given, respectively, by(26)Φθ=q1θ,q1¯θ,q2θ,q2¯θ,(27)ψs=q1∗s,q1∗¯s,q2∗s,q2∗¯sT,where(28)q1θ=1,iω1+τc,−3r+iω1+τc+3r/ciω1+τcP/ce−iω1+τc,−3r+iω1+τc+3r/ciω1+τcP/ce−iω1+τciω1+Teiω1+τcθ≜q110,q120,q130,q140Teiω1+τcθ,(29)q2θ=1,iω2−τc,−3r+iω2−τc+3r/ciω2−τcP/ce−iω2−τc,−3r+iω2−τc+3r/ciω2−τcP/ce−iω2−τciω2−Teiω2−τcθ≜q210,q220,q230,q240Teiω2−τcθ,(30)q1∗s=D11,iω1+τc3r,9r2−3irω1+τciω1+τc−3r/c3iPrω1+τceiω1+τc/c,3rc−iω1+τcciω1+τc−3r3Preiω1+τce−iω1+τcs≜D1q11∗0,q12∗0,q13∗0,q14∗0e−iω1+τcs,(31)q1∗s=D21,iω2−τc3r,9r2−3irω2−τciω2−τc−3r/c3iPrω2−τceiω2−τc/c,3rc−iω2−τcciω2−τc−3r3Preiω2−τce−iω2−τcs≜D2q21∗0,q22∗0,q23∗0,q24∗0e−iω2−τcs,(32)D1=q110q11∗0¯+q120q12∗0¯+q130q13∗0¯+q140q14∗0¯+P/cq120q14∗0¯+q140q12∗0¯τceiω1+τc−1,(33)D2=q210q21∗0¯+q220q22∗0¯+q230q23∗0¯+q240q24∗0¯+Pcq220q24∗0¯+q240q22∗0¯τceiω2−τc−2.We letP=Pc+ε1 and τ=τc+ε2 in Equation (20), where ε1 and ε2 are perturbation parameters, and denote ε=ε1,ε2. Then Equation (20) can be rewritten as(34)X˙t=LεXt+FXt,ε,where(35)LεXt=τc+ε20100−3r−3rc00000100−3r−3rcx0h0y0g0+τc+ε20100000−P+ε1c00000−P+ε1c00x−1h−1y−1g−1and(36)FXt,ε=τc+ε20B1x03+B2x02h0+B3h03+B4ht2g−10B1y03+B2y02g0+B3g03+B4gt2h−1B1=−P2c2,B2=−3r+26c,B3=−1c+3rc3+12,B4=−3Pc.Thus, we can easily get Equation (20) that becomes an abstract ODE in the space BC.(37)ddtμ=Aμ+X0F¯,where,μ∈C, and A is defined by(38)A:C1→BC,Aμ=μ˙+X0L0μ−μ˙0and(39)F¯μ,ε=Lε−L0μ+Fμ,ε.By the continuous projectionπ, we can decompose the enlarged phase space by ∧=±τcω1+,±τcω2− as BC=PΛ⊕kerπ, where kerπ=ϕ+X0c:=πϕ+X0c=0. Let μt=Φzt+vt, Then, Equation (20) is decomposed as(40)z˙=Bz+ψ0F¯Φz+v,ε,v˙=AQ1v+Id−πX0F¯Φz+v,ε,where, we denoteν∈Q1=Q∩C1⊂kerπ, and AQ1 is the restriction of A as an operator from Q1 to the Banach space kerπ.For the convenience of computing the normal form, Equation (40) can be written by neglecting higher order terms with respect to parameters, with(41)z˙1=iτcω1+z1+∑i=14D1q1i∗0fi+h.o.t,z˙2=−iτcω1+z2+∑i=14D¯1q1i∗0¯fi+h.o.t,z˙3=iτcω1−z3+∑i=14D2q2i∗0fi+h.o.t,z˙4=−iτcω1−z4+∑i=14D¯2q2i∗0¯fi+h.o.t,where(42)f1=ε2h0,f2=−3rε2x0−3rε2ch0−τcε1+ε2Pccg−1+B1x30+B2x20h0+B3h30+B4h20g−1,f3=ε2g0,f4=−3rε2y0−3rε2cg0−τcε1+ε2Pcch−1+B1x30+B2y20h0+B3g30+B4g20h−1and(43)x0h0y0g0=Φ0z+v10v20v30v40,x−1h−1y−1g−1=Φ−1z+v1−1v2−1v3−1v4−1.Next, letM21 denote the operator defined in V26C4×kerπ, with(44)M21:V26C4↦V26C4,(45)M21Pz,ε=DzPz,εBz−BPz,ε,whereV26C4 represents the linear space of the second-order homogeneous polynomials in six variables z1,z2,z3,z4,ε1,ε2 with coefficients in C4. It is easy to verify that one may choose the decomposition without considering the strong resonant(46)V26C4=ImM21C4⊕ImM21C4cwith complementary space(47)ImM21C4c=spanz1εi000,0z2εi00,00z3εi0,000z4εii=1,2.Then, the normal form of Equation (37) on the center manifold associated with the equilibrium near ε=0 has the form(48)z˙=Bz+12g21z,0,ε+h.o.t,whereg21 is the function giving the quadratic terms in z,ε for V=0, and is determined by g21z,0,ε=ProjImM21C4c×f21z,0,ε, where f21z,0,ε is the function giving the quadratic terms in z,ε for V=0 defined by the first equation of Equation (41). Thus, the normal form truncated at the quadratic order terms is given by(49)z˙1=iτcω1+z1+λ1z1+h.o.t,z˙2=−iτcω1+z2+λ¯1z2+h.o.t,z˙3=iτcω1−z3+λ2z3+h.o.t,z˙4=−iτcω1−z4+λ¯2z4+h.o.t,where(50)λ1=2D1q11∗ε2q12+q12∗−3rε2q11−3rcε2q12−τcε1+ε2Pccq14eiτcω2−+q13∗ε2q14+q14∗−3rε2q13−3rcε2q14−τcε1+ε2Pccq12eiτcω1+,(51)λ2=2D2q21∗ε2q32+q22∗−3rε2q31−3rcε2q32−τcε1+ε2Pccq34eiτcω2−+q23∗ε2q34+q24∗−3rε2q33−3rcε2q34−τcε1+ε2Pccq32eiτcω1+.To find the normal form up to third order, similarly, letM31 denote the operator defined in V36C4×kerπ, with(52)M31:V26C4↦V36C4,(53)M31Pz,ε=DzPz,εBz−BPz,ε,whereV36C4 denotes the linear space of the third order homogeneous polynomials in four variables z1,z2,z3,z4 with coefficients in C4.Decompose the spaceV36C4 as follows(54)V36C4=ImM31C4⊕ImM31C4c,where the complementary spaceImM31C4c spanned by the elements(55)ImM31C4c=spanz12z2000,z1z3z4000,0z22z100,0z2z3z400,00z1z2z30,00z32z40,000z1z2z4,000z42z3.Therefore, the normal form up to third-order terms is given by(56)z˙=Bz+12g21z,0,ε+13!g31z,0,0+h.o.t,where(57)13!g31z,0,0=13!I−P1,31f31z,0,0,andf31z,0,0 is the function giving the cubic terms in z,ε,v for V=0, and ε=0 is defined by the first equation of Equation (40).Finally, the normal form on the center manifold arising from Equation (40) becomes(58)z˙1=iτcω1+z1+λ1z1+a11z12z2+a12z1z3z4+h.o.t,z˙2=−iτcω1+z2+λ¯1z2+a¯11z22z1+a¯12z1z3z4+h.o.t,z˙3=iτcω2−z3+λ2z3+a21z32z4+a22z1z2z3+h.o.t,z˙4=−iτcω2−z4+λ¯2z4+a¯21z42z3+a¯22z1z2z4+h.o.t,where(59)a11=D1q12∗3B1q112q21+B2q112q22+2q21q12q11+3B3q122q22+B4q122q24eiτcω2−+2q12q22q14eiτcω2−+D1q14∗3B1q132q23+B2q132q24+2q23q14q13+3B3q142q24+B4q142q22eiτcω1++2q14q24q12eiτcω1+,(60)a12=D1q12∗6B1q11q31q41+B22q11q31q42+2q11q41q32+2q31q41q12+6B3q12q32q42+B42q12q32q44eiτcω2−+2q12q42q34eiτcω2−+2q32q42q14eiτcω2−+D1q14∗6B1q13q33q43+B22q13q33q44+2q13q43q34+2q33q43q14+6B3q14q34q42+B42q14q34q42eiτcω1++2q14q44q32eiτcω1++2q34q44q12eiτcω1+,(61)a21=D2q22∗3B1q312q41+B2q312q42+2q41q32q31+3B3q322q42+B4q322q44eiτcω2−+2q32q42q34eiτcω2−+D2q24∗3B1q332q43+B2q332q44+2q43q34q33+3B3q342q44+B4q342q42eiτcω1++2q34q44q32eiτcω1+,(62)a22=D2q22∗6B1q11q21q31+B22q11q31q22+2q11q21q12+2q31q21q12+6B3q12q32q22+B42q12q32q24eiτcω2−+2q12q22q34eiτcω2−+2q32q22q14eiτcω2−+D2q24∗6B1q13q33q23+B22q13q33q24+2q13q23q34+2q33q23q14+6B3q14q34q22+B42q14q34q22eiτcω1++2q14q24q32eiτcω1++2q34q24q12eiτcω1+.With the polar coordinates(63)z1=r1e−iρ1z2=r1eiρ1z3=r2e−iρ2z1=r2eiρ2then, the amplitude equation resulted from Equation (58) is(64)r˙1=α1+a¯r12+b¯r22r1r˙2=α2+c¯r12+d¯r22r1and the nature of bifurcation of Equation (64) can refer to [20]. We introduce new phase variables and rescale time according to R1=−a¯r1,R2=−d¯r1,τ=2t, and obtain(65)R˙1=α1−R1−θR2R1,R˙2=α2−δR1−R2R2,whereα1=Reλ1,α2=Reλ2,a¯=Rea11,b¯=Rea12,c¯=Rea21,d¯=Rea22,θ=b¯/d¯,δ=c¯/a¯.


## 4. Illustrations and Numerical Simulations

To give a more clear bifurcation picture, we chooser=4/3,c=94,Pc=30.21,τc=0.78 by referring [12], and yield ω1+=2.166,ω2−=1.847,ω1+:ω2−=52:6. By calculating, we can obtain α1=0.410ε2+0.605ε1, α2=−30.514ε2+24.351ε1, a¯=−4.729<0, b¯=−10.031, c¯=−.5.554, d¯=−2.003<0,θ=5.007, and δ=1.117. For the bifurcation behaviors of the original system Equation (20) in the neighborhood of the trivial equilibrium, the above critical bifurcation boundaries divide the parameter plane Pc,τc into six regions (see Figure 2(a)). In region A1, Equation (20) exists one trivial equilibrium, when the parameters are changed from region A1 to A6, the stable trivial equilibrium becomes a saddle and a stable periodic solution is bifurcated. In region A5, another stable periodic solution appears. When the parameters vary across line T2 from region A5 to A4, an unstable periodic solution bifurcates an unstable quasi-periodic solution which is a nobble, and there coexists with two stable periodic solutions and an unstable quasi-periodic solution. Further, in region A3, an unstable quasi-periodic solution disappears when parameters change. Finally, when the parameters are varied in the region A2, the periodic solution collides with the trivial solution and the trivial solution ranges from a source to a saddle. Here, in Figure 2, the bifurcation critical lines are, respectively(66)H1=α1,α2,α1=0:ε1=−0.698ε2,(67)H2=α1,α2,α2=0:ε1=1.253ε2,Figure 2
(a) critical bifurcation lines in theε1,ε1 parameter space near Pc,τc; and the corresponding phase portraits in the R1,R1 plane when Pc=30.21,τc=0.78, (b) to embed Figure 2(a) into Figure 1, and yielding complete bifurcation sets near Pc,τc.
(a)
(b)which stand for the Hopf bifuration lines,(68)T1=α1,α2,α2>0:ε1=1.262ε2(69)T2=α1,α2,α1>0:ε1=1.308ε2and which stand for the periodic solution of pitchfork bifurcation lines.We will present some numerical simulation results and choose the values of the parameters, respectively,ε1,ε2=0.4,−0.08∈A6 and ε1,ε2=−0.4,−0.08∈A1. Corresponding to a stable periodic solution shown in Figure 3, the equilibrium (0 0 0 0) is asymptotically stable, shown in Figure 4. Thus, the result shows that the numerical simulations are coincident with the analytical predictions.Figure 3
Simulated solution of Equation (20) for ε1,ε2=0.4,−0.08, showing a stable periodic solution.
(a)
(b)Figure 4
Simulated solution of Equation (20) for ε1,ε2=−0.4,−0.08, showing the equilibrium(0 0 0 0) is asymptotically stable.
(a)
(b)


## 5. Explanation for the Mechanical Background of the Model Using Mathematical Results

AsP=0, then Equation (2) can be written as(70)a˙−caa¨+32a˙2−a−3r+1−a˙3−3r−2a−3ra˙−2a˙=0,which shows the behavior of systems range from double bubbles coupling to single bubble. The zero solution of Equation (70) is asymptotically stable (see [12]).As noted in [12], τ=d/c shows that the distance of two bubbles d is proportional to the time delay τ when sound speed is constant. That coupling coefficient P decreases as the distance between bubbles increases (see [13]), and meaning the coupling coefficient P as a decreasing function of the delay τ. Clearly, if P is sufficiently small and τ is sufficiently big (τ is sufficiently small and P is sufficiently big), the behavior of Equation (2) at the equilibrium (0 0 0 0) will be close to Equation (70) (see Figures 5 and 6). And if neglecting relationship between P and τ and choosing them sufficiently big, Equation (2) will have a chaotic motion (see Figure 7). Therefore, the double bubble oscillators Equation (2) can show the possibility of more complicated dynamics.Figure 5
Simulated solution of Equation (2) for P=0.001,r=4/3,c=94,τ=1000 : the equilibrium (0 0 0 0) is asymptotically stable.Figure 6
Simulated solution of Equation (2) for P=1000,r=4/3,c=94,τ=0.001: the equilibrium (0 0 0 0) is asymptotically stable.Figure 7
Simulated solution of Equation (2) for P=1000,r=4/3,c=94,τ=100 : showing the chaotic motion.
(a)
(b)Bubble dancing [21] is an interesting phenomenon in bubble dynamics. It creates a pull or push under the sound pressure gradient between the bubbles, which makes them undergo a periodic motion in turn, and the periodic motion is called bubble dancing. As noted in Figure 3, it shows a stable periodic solution by taking value of parameters, which can depict periodic motion of bubble’ radius, and the change of bubble’ radius can describe the motion of bubble in the liquid. Thus, the stable periodic solution of system can be used to better explain this physical phenomenon (bubble dancing).


## 6. Conclusions and Discussion

In this work, we mainly have discussed the nonresonant double Hopf bifurcation in dynamics of microbubble with delay coupling and have used center manifold reduction methods to compute the normal form near the double Hopf bifurcation point. Moreover, with bifurcation analysis near the double Hopf critical point given, we confirm the existence of a stable periodic solution and an unstable quasi-periodic solution. We discover the system coexists with two stable periodic solutions and an unstable quasi-periodic solution, which is a new phenomenon with studying the double Hopf bifurcation.However, we only get the normal forms for nonresonant cases owing to the limitations of our knowledge, so analyzing the microbubble dynamics of the resonance case will be our next research.



---

*Source: 8059037-2020-01-03.xml*