Bubbling transition as a mechanism of destruction of synchronous oscillations of identical microbubble contrast agents
BBubbling transition as a mechanism of destruction of synchronous oscillations of identical microbubble contrast agents
Ivan R. Garashchuk and Dmitry I. Sinelshchikov National Research University Higher School of Economics, Moscow, RussiaFebruary 5, 2021
Abstract
We study the process of destruction of synchronous oscillations in a model of twointeracting microbubble contrast agents exposed to an external ultrasound field. Com-pletely synchronous oscillations in this model are possible in case of identical bubbles,when the governing system of equations possess a symmetry leading to the existence of asynchronization manifold. Such synchronous oscillations can be destructed without break-ing the corresponding symmetry of the governing dynamical system. Here we describethe phenomenological mechanism responsible for such destruction of synchronization anddemonstrate its implementation in the studied model. We show that the appearance andexpansion of transversally unstable areas in the synchronization manifold leads to transfor-mation of a synchronous chaotic attractor into a hyperchaotic one. We also demonstratethat this bifurcation sequence is stable with respect to symmetry breaking perturbations.
Here we consider a model describing oscillations of two interacting microbubble contrast agents.Microbubble contrast agents are micrometer size encapsulated gas bubbles which are used forenhancing ultrasound visualisation [1–3]. There are also other applications of such bubbles suchas targeted drug delivery and noninvasive therapy [4, 5].It is known that the type of dynamics of contrast agents may be important for particularapplications (see, e.g. [3, 6]). For instance, chaotic dynamics may be beneficial for ultrasoundvisualization, while regular dynamic may be required for targeted drug delivery. Moreover,large amplitude hyperchaotic oscillations may be utilized for the purpose of quick destructionof the bubbles, e.g. in some drug delivery related applications [4]. Therefore, it is necessary todetermine both possible type of bubbles oscillations and their bifurcations. On the other hand,in a number of works it was shown that the contrast agents dynamics may be complicatedand strongly depend on both control parameters and initial conditions (see, e.g. [6–12]). Inworks [6–9] the case of one microbubble was studied and it was demonstrated that the dynamicsmay be either regular or chaotic and multistable in both cases, while, e.g. in [10–12], interactionsbetween bubbles were taken into account.In particular, in [11] it was found that the dynamics of two interacting identical con-trast agents may be regular, quasiperiodic, chaotic and hyperchaotic. Furthermore, authorsof Ref. [11] demonstrated that oscillations of bubbles may be both synchronous and asyn-chronous. It is worth noting that a synchronization manifold exist only in the case of two1 a r X i v : . [ n li n . C D ] F e b dentical bubble, when the governing system of equations possesses a symmetry [11, 12]. Thus,complete synchronization is possible only in the symmetrical case. In [12] authors studiedhow symmetry breaking affects various synchronous and asynchronous attractors. However,the destructions of synchronous oscillations without symmetry breaking in the model of twointeracting contrast agents has not been studied previously.Therefore, the main aim of this work is to study the process of destruction of synchronousoscillations in a system of two identical interacting contrast agents. We demonstrate that thetypical mechanism responsible for this is the bubbling transition [13–18] that was observedin various dynamical systems [19, 20]. The key point in this scenario is the appearance oftransversally unstable saddle orbits inside the synchronization manifold leading to the emer-gence of transversally unstable areas there and, eventually, to the disappearance of the syn-chronous attractor. At first, a trajectory can escape the synchronization manifold for a shorttime. Then, the more transversally unstable orbits emerge in the synchronization manifold, themore frequent and long these escapes become. Finally, the trajectory spends the majority oftime outside the synchronization manifold and the attractor becomes hyperchaotic. Since thisscenario leads to a hyperchaotic attractor, the bubbling transition can be considered as anotherroute to hyperchaos that is alternative to the scenario involving the emergence of the discreteShilnikov attractor [11] (see also [21–24]).It is worth noting that in this work we consider complete synchronization of two chaoticsystems (see, e.g. [17, 25–29]), each of which represents forced oscillations of a microbubblecontrast agent. Single forced model without coupling would demonstrate chaotic behaviour inthe considered regions of the control parameters.The rest of this work is organized as follows. In the next section we introduce the governingdynamical system and discuss some of its properties. In Section 3 we propose a phenomenologi-cal scenario of the destruction of synchronous oscillations in the symmetrical case of consideredmodel. Section 4 is devoted to numerical demonstration of the implementation of the proposedscenario. In the last Section we briefly summarize and discuss our results. We consider the following system of equations [10, 30–34] (cid:32) − ˙ R c (cid:33) R ¨ R + 32 (cid:32) − ˙ R c (cid:33) ˙ R = 1 ρ (cid:34) R c + R c ddt (cid:35) P − ddt (cid:32) R ˙ R d (cid:33) , (cid:32) − ˙ R c (cid:33) R ¨ R + 32 (cid:32) − ˙ R c (cid:33) ˙ R = 1 ρ (cid:34) R c + R c ddt (cid:35) P − ddt (cid:32) R ˙ R d (cid:33) , (1)where P i = (cid:18) P + 2 σR i (cid:19) (cid:18) R i R i (cid:19) γ − η L ˙ R i R i − σR i − P − χ (cid:18) R i − R i (cid:19) − κ S ˙ R i R i − P ac sin( ωt ) , i = 1 , . System (1) describes oscillations of two encapsulated gas bubbles, which interact throughthe pressure field, in a compressible viscous liquid. Here R and R are radii of the bubbles, d is the distance between the centers of bubbles, P stat = 100 kPa is the static pressure, P v = 2 . kPa is the vapor pressure, P = P stat − P v , P ac is the magnitude of the pressure of the externalfield, σ = 0 . N/m is the surface tension, ρ = 1000 kg/m is the density of the liquid, η L = 0 . Ns/m is the viscosity of the liquid, c = 1500 m/s is the sound speed, γ = 4 / isthe polytropic exponent (the process is assumed to be adiabatic), χ and κ s denote the shell2lasticity and shell surface viscosity respectively. Dynamical system (1) takes into accountliquid’s compressibility according to the Keller–Miksis [35] model and bubbles’ shell accordingto the de-Jong model [36–38].We use the following values of the parameters χ = 0 . N/m and κ S = 2 . · − kg/s thatcorrespond to the SonoVue contrast agent with equilibrium radius R i = 1 . µ m [39].In this work we study the case of bubbles with the same equilibrium radii, i.e. R = R = R . Under this condition the system (1) is symmetrical with respect to transformation R ←→ R . This leads to the existence of the synchronization manifold, defined by equations R = R , ˙ R = ˙ R , i.e. it is a three-dimensional plane in a five-dimensional phase space. Belowwe show that a synchronous attractor lying inside this manifold can be destructed through thebubbling transition mechanism. We also demonstrate that our results are stable with respectto the symmetry breaking perturbation.Parameters P ac and ω can be naturally considered as control ones, since they define proper-ties of the external ultrasound field. The parameter d is defined by the density of the injectedbubbles cluster and is also considered as control one. In this work, we focus on the influ-ence of the parameter d on the dynamics and fix P ac and ω as follows P ac = 1 . MPa and ω = 2 . s − . These values of P ac and ω correspond to biomedically relevant regions ofthese parameters [3].We perform all calculations in the following non–dimensional variables R i = R r i , t = ω − τ , where ω = 3 κP / ( ρR ) + 2(3 κ − σ/R + 4 χ/R is the natural frequency of bubbleoscillations. We also denote the non-dimensional radial speeds of the bubbles’ surfaces by u and u . We use the fourth-fifth order Runge–Kutta method [40] for finding numerical solutionsof the Cauchy problem for (1). For the calculations of the Lyapunov spectra we use the standardalgorithm by Bennetin et al [41], modified according to the specifics of our problem, see Sec.4 for details. Poincar´e maps are constructed as a cross-section at every period of the externalpressure field, i.e. we use stroboscopic section. In this section we propose a phenomenological scenario of the destruction of a synchronousattractor in model (1). This scenario is based on the bubbling transition mechanism [13–17].We consider a synchronous chaotic attractor in system (1), belonging to a set
F ix ( S ) ,obtained via a cascade of period-doubling bifurcations of a synchronous limit cycle. Supposealso we have a non-attractive hyperchaotic hyperbolic set SH (2 , near F ix ( S ) , see Fig. 1a.Notice that for better clarity we demonstrate the Poincar´e map instead of the full phase space.Thus, for instance, points in Fig. 1 correspond to periodic trajectories of the original system.Let us also remark that it does not matter, how the set SH (2 , is created (e.g. it is initiallynon-attractive, or appears after a crisis of a hyperchaotic attractor).There is an infinite number of saddle limit cycles of (3 , -type inside the synchronouschaotic attractor, which are transversally stable and have one-dimensional unstable manifoldembedded inside the synchronization manifold (see white points on F ix ( S ) in Fig. 1a). In thiscase we observe strong synchronization [17], since all periodic orbits are transversally stable.On the other hand, set SH (2 , outside of the synchronization manifold contains infinitelymany periodic saddle orbits of (2 , -type (see red points in Fig. 1a). Due to the symmetryof system (1) ( R ←→ R ), set SH (2 , consists of two symmetrical parts above and below F ix ( S ) . The set SH (2 , is hyperchaotic because two-dimensional unstable manifolds of (2 , saddle periodic orbits result in two directions of hyperbolic instability when a trajectory passesclose to them. This will be numerically demonstrated in Section 4.Then, if we increase a control parameter, set SH (2 , can approach F ix ( S ) and unstablemanifolds of its (2 , -type saddle orbits can form heteroclinic manifolds with transversally sta-3igure 1: Scheme of the main steps of the bubbling transition scenario. (a) Strongly syn-chronous state, (b) weakly synchronous state, (c) weakly asynchronous state.ble (3 , -type orbits lying inside F ix ( S ) and forming the skeleton of the synchronous attractor.Eventually, at a certain value of the control parameter these transversally stable (3 , saddleorbits from F ix ( S ) collide with (2 , -type saddles from SH (2 , merging via subcritical pitch-fork bifurcations and, thus, creating transversally unstable (2 , -type orbits O i inside F ix ( S ) (see Fig. 1b). One unstable invariant manifold of O i is embedded inside F ix ( S ) , while anotherone is transversal to it. Since the most of the orbits belonging to F ix ( S ) are still transversallystable, the ejection of a trajectory from F ix ( S ) is a rare event, which can be considered as anextreme event [20]. Moreover, after escaping from F ix ( S ) through the transversally unstableinvariant manifold of O i , the trajectory goes through SH (2 , and returns back to F ix ( S ) through a transversally stable area of F ix ( S ) . At this point the synchronous attractor be-comes a Milnor attractor with riddled basin [17]. Below, by computing transversal Lyapunovexponents [17, 29], we demonstrate that it is stable on average. Therefore, we have a weaklysynchronous state at these values of the control parameter.Consequently, an orbit naturally splits into two phases: synchronous and asynchronous.There are two logically distinct sets inside the attractive area of the phase space, throughwhich a trajectory goes: chaotic set inside F ix ( S ) and hyperchaotic saddle set SH (2 , , andattractor as a whole consists of synchronous and asynchronous components. If we further in-crease the value of the control parameter, the subcritical pitchfork bifurcations continue tohappen, and more and more transversally unstable synchronous orbits of (2 , -type appear4n F ix ( S ) . As the joint contribution of such transversally unstable orbits becomes more sig-nificant, substantial areas of F ix ( S ) become transversally unstable. Therefore, a trajectoryspends more and more time outside the synchronization manifold. At a certain point a balancecan be reached: trajectory spends roughly equal time inside and outside F ix ( S ) , and bothcomponents have similar impact on the behavior of a trajectory as the whole. At the pointwhen the largest transversal Lyapunov exponent becomes positive, the synchronous set F ix ( S ) turns transversally unstable on average, while there still exist transversally stable orbits insidethe synchronization manifold. Here we observe a weakly asynchronous state. Since the asyn-chronous component is hyperchaotic, this weakly asynchronous state can also be hyperchaotic.The implementation of this scenario in system (1) will be demonstrated numerically in Sec. 4.Increasing the control parameter further could lead to the blow-out bifurcation [14, 16, 20]:almost all of the orbits that used to form the skeleton of the synchronous attractor becometransversally unstable. The manifold F ix ( S ) becomes completely transversally unstable, anda trajectory never comes back to it. The set SH (2 , becomes attractive, and the attractorbecomes fully asynchronous and hyperchaotic. However, in the case of system (1) along theone-parametric route studied in this paper, the blow out bifurcation does not happen since thehyperchaoic attractor disappears earlier due to a crisis. In this section we demonstrate the implementation of the above proposed scenario of destruc-tion of a synchronous attractor in model (1). As we show below, this bifurcations sequence leadsto a transformation from a synchronous periodic attractor to a hyperchaotic one, consisting ofsynchronous and asynchronous parts, in accordance with the hypothesis in Sec. 3. In numericalcalculations we split the points of a trajectory into synchronous and asynchronous ones accord-ing to the following criterion: if | r − r | < − , | u − u | < − , then the corresponding pointbelongs to the synchronization manifold, otherwise it belongs to the asynchronous part of atrajectory. We also modify the scheme of calculating the Lyapunov spectrum in the same way:exponents for synchronous and asynchronous components are calculated separately. Lyapunovexponents of an attractor as a whole can be calculated by averaging along the whole trajectoryor as a weighted sum of both components with weights N sync /N and N async /N respectively,where N is total number of points used in calculations and, N sync and N async are numbers ofpoints belonging to synchronous and asynchronous parts of a trajectory, respectively.We perform calculations of the Lyapunov spectrum with the following time step T l = 0 . T between orthonormalizations, where T is the non-dimensional period of the external pressurefield. For each value of the control parameter we use . · steps of length T l , which cor-responds to the non-dimensional time . · T . We consider our estimations of Lyapunovexponents statistically correct only if at least · steps belong to each component of the tra-jectory (or at least of points belong to each phase). Thus, synchronous and asynchronousparts of the Lyapunov spectrum are demonstrated only for intervals of the control parameterfor which the number of points belonging to each component exceeds this threshold.As a quantitative estimation of transversal stability of the synchronization manifold weuse the largest transversal Lyapunov exponent [17, 27, 29]. It shows whether an attractor istransversally stable or unstable on average. If this exponent is negative, then the correspond-ing attractor is transversally stable on average, and otherwise if it is positive. We do notprovide the values of the second transversal Lyapunov exponent, because it is always negativewith a comparatively large absolute value. With the help of the combination of the calculatedquantities (number of synchronous/asynchronous points, Lyapunov and transversal Lyapunovexponents) we can separate different qualitative types of behavior along the route from syn-chronous to asynchronous state. If there are no asynchronous points along a trajectory and5 d/R -0.0500.050.1 f21 d/R -0.100.10.20.3 fsync1sync2async1async2 ab Figure 2: The dependence of the three largest Lyapunov exponents on d/R , ∈ [12 . , . :(a) calculated along the full trajectories, (b) calculated separately for synchronous and asyn-chronous phases of the trajectories (plotted only for intervals with enough points in bothphases).the largest transversal Lyapunov exponent is negative, we observe strongly synchronous state.Suppose that the asynchronous phase of the trajectory is present, but the largest transversalLyapunov exponent is still negative. Consequently, the synchronous attractor is still transver-sally stable on average and we have a weakly synchronous state. If there are asynchronous andsynchronous phases of a trajectory and the largest transversal Lyapunov exponent is positive,we deal with a weakly asynchronous state. If all the orbits inside the synchronization manifoldbecame unstable, one could observe strongly asynchronous state. However, in the consideredone-parametric route the hyperchaoic attractor disappears earlier due to a crisis.In what follows we fix ω and P ac and consider d/R as the control parameter and provideone-parametric maps of quantities characterizing the dynamics of the system with respect tothis parameter. In Fig. 2a we present the dependence of Lyapunov exponents calculated alongthe whole trajectory on the control parameter. To show possible types of behavior in differentphases we demonstrate Lyapunov exponents that are calculated separately for synchronous andasynchronous components of this trajectory in Fig. 2b. We also show how ratios N async /N and N sync /N change with the parameter d/R in Fig. 3a. They approximately represent fractionsof time that a trajectory spends inside synchronization manifold and outside of it, respectively.The graph of the largest transversal Lyapunov exponent is presented in Fig. 3b. Once can6 d/R sync /NN async /N12.1 12.2 12.3 12.4 d/R -0.0100.010.02 t r ab Figure 3: The dependence of (a) relative numbers of points in synchronous and asynchronouscomponents N async /N and N sync /N on d/R and (b) the maximal transversal Lyapunov expo-nent on d/R .clearly notice strong correlation between graphs in Fig.-s 3a,b.In order to calculate the Lyapunov spectrum for (1) we convert it into five-dimensionalautonomous dynamical system. Any attractor of such system, which is not a stationary point,has a zero Lyapunov exponent, that corresponds to the translations along this attractor. Wedenote this Lyapunov exponent by λ f and plot it on the graphs as a reference to convenientlyseparate positive and negative values of the exponents. Apart from λ f , we present only the twolargest exponents, because they help to distinguish between different types of attractors. Onthe other hand, the other two exponents are always negative (since system (1) is dissipative)and their values do not carry any substantial information. If we presented these Lyapunovexponents in the same plots with the two largest ones, we would significantly decrease theclarity due to the following reason. It would be necessary to change the scale of the graphsgreatly in order to fit those large negative values. As a consequence, it would be very difficultto distinguishing between various dynamical regimes since the two largest Lyapunov exponentsare often close to zero. To avoid ambiguity, below, we exclude λ f from the consideration, whenwe discuss numbers of positive and negative exponents. In Tables 1, 2 we present the valuesof Lyapunov exponents, the proportion of asynchronous points and the largest transversalLyapunov exponent for several particular values of parameters, that correspond to certainsteps in the phenomenological scenario proposed in the previous section. In these Tables we7lso exclude λ f from consideration.One can see that the increase in the parameter d/R leads to the emerging of asynchronouspoints even if one starts with synchronous initial conditions. This means that transversallyunstable saddle orbits indeed start appearing inside the synchronization manifold and a trajec-tory may leave the synchronization manifold when it passes close to one of these saddle orbits.Therefore, the suggested splitting of a trajectory on synchronous and asynchronous parts isquite natural. Notice also that the asynchronous phase is virtually always characterized bytwo positive Lyapunov exponents, which means that the dynamics in the asynchronous phaseis hyperchaotic, i.e. two directions of hyperbolic instability are present outside of the synchro-nization manifold. This supports our hypothesis of the existence of a saddle set outside thesynchronization manifold, containing saddle orbits of (2 , -type.Figure 4: Projections on the ( r , u , r ) space of (a) phase curve and Poincar´e map of thesynchronous 2-periodic limit cycle at d/R = 7 . , (b) synchronous chaotic attractor at d/R =12 . , emerging as a result of a cascade of period-doubling bifurcations of the limit cycle inpanel (a)In Fig.-s 4a-8b we present projections of the Poincar´e maps of the attractors at givenparameter values on ( r , u , r ) space. We use the Poincar´e sections instead of phase curves,because high orbits density for the chaotic attractors makes it very hard to interpret the results.The corresponding values of Lyapunov exponents are presented in Tables 1, 2. Fig. 4a isthe only one where we present the phase curve along with the Poincar´e section, because itshows a synchronous limit cycle at d/R = 7 . . If one increases the parameter d/R , thislimit cycle undergoes the cascade of period-doubling bifurcations and becomes a synchronouschaotic attractor. Poincar´e section of this attractor at d/R = 12 . is presented in Fig.8b. One can notice that all the points lie on the plane r = r , which is the projection ofthe synchronization plane onto the three-dimensional space ( r , u , r ) . Just as expected froma chaotic attractor embedded in the three-dimensional synchronization manifold, the largestLyapunov exponent of this attractor is positive and the second one is negative, see Table1. Taking into account negative values of the corresponding largest transversal Lyapunovexponents for these synchronous attractors, we conclude that they are transversally stable onaverage and represent strongly synchronous states.Table 1: Lyapunov exponents and N async /N ratio for certain values of parameters.Fig. d/R N async /N λ λ λ sync λ sync λ async λ async λ tr . − . − . − . − . – – − . . . , − . . − . – – − . .
05 0 .
716 0 . − . − . − . . . . .
12 0 .
206 0 . − . . − . . . − . If we further increase the value of d/R , we observe the start of the bubbling transition se-quence (see the initial increase of N async /N and λ tr in Fig.-s 3a,b inside . (cid:46) d/R (cid:46) . interval). The left border of this quick growth marks the starting point at which subcriticalpitchfork bifurcations start happening and the end of the existence of the synchronous attractor.We demonstrate in Fig. 5a an attractor that can be obtained at d/R = 12 . and representsa typical attractor in this region of values of the control parameter. The whole orbit is charac-terized by one positive Lyapunov exponent, i.e. behavior on average is chaotic. However, theasynchronous component has two positive Lyapunov exponents implying hyperchaotic dynam-ics in this phase, while in the synchronous component all exponents are negative (see Table1). Notice also the λ tr > for this attractor, which means that it is transversally unstable onaverage and, consequently, corresponds to a weakly asynchronous state. While there still existareas inside the synchronization manifold, where the Lyapunov sums are growing, the dynamicsin the synchronous component on average become non-chaotic. This probably happens becausethe trajectory inside the synchronization manifold does not pass through the regions where theone-dimensional divergence of the orbits used to be the strongest in the synchronous attractor.One can see that the transition from the synchronous state to weakly asynchronous one isnot monotonous: after initial sharp increase of the asynchronous component, the number ofasynchronous points on the attractor sharply drops almost to zero, and stays low for some timebefore beginning to grow again (see Fig. 3a). In Fig. 5b we show an attractor at d/R = 12 . that corresponds to the interval in d/R for which the presence of the asynchronous phase islow, but still visible. In this case there are two positive Lyapunov exponents in the asynchronousphase and one positive Lyapunov exponent in the synchronous phase. On average, across bothphases one exponent is positive and the others are negative. However, it is worth noting thatfor this attractor N async /N = 0 . , which means that it belongs to the local spike in the areaof low N async /N in Fig. 3a. Negative value of λ tr for this attractor implies that it is weaklysynchronous. Another attractor from the interval of low N async /N , that is closer to the localminimum of asynchronous points, is presented in Fig. 6a and corresponds to d/R = 12 . .In this case N async /N = 0 . and only a few points belong the asynchronous phase. Thus,we do not present the values of Lyapunov exponents for the asynchronous component, becausethere are not enough points to make statistically correct conclusions. The values of Lyapunovexponents averaged across both phases and for synchronous phase are given in Table 2. Thebehavior of an orbit is mostly determined by the synchronous phase and is chaotic on average.This attractor is also transversally stable on average.9igure 5: Projections of Poincar´e maps of attractors on the ( r , u , r ) space at (a) d/R =12 . , N async /N = 0 . high presence of asynchronous component, (b) d/R = 12 . , N async /N = 0 . substantial drop in the amount of points in the asynchronous componentIn Fig.-s 6a–8b we show a sequence of attractors with steadily growing amount of pointsbelonging to the asynchronous phase. The values of Lyapunov exponents corresponding to eachattractor are presented in Table 2. Combined with the graphs from Fig.-s 2a-3b, the generaltrend can be seen. In Fig. 6a the asynchronous component is barely seen, because there areonly a few asynchronous points and the vast majority of them are located very close to thesynchronization manifold. Thus, asynchronous points can be hardly visually distinguished fromsynchronous points. In Fig. 6b the asynchronous component becomes clearly noticeable, butstill the asynchronous points are mostly concentrated in a narrow band that is close to thesynchronization manifold. In these region of d/R the values of λ tr oscillate close to zero, witha steady trend towards positive direction. This means that the synchronous attractor is on theverge of being transversally unstable.Increasing d/R further, we observe the next step of the proposed scenario. More andmore transversally unstable saddle orbits of (2 , -type start appearing inside the synchro-nization manifold, which contributes into the formation of transversally unstable areas on thesynchronization manifold. This can be observed in Fig.-s 7a,b. The number of asynchronouspoints in Fig. 7a is only two times larger than the number of those presented in Fig. 6b,however, they spread much further away from the synchronization manifold, which makes themsignificantly more apparent. The largest transversal Lyapunov exponent is positive for the lasttwo attractors characterizing them as weakly asynchronous and keeps growing further with d/R .In Fig. 7b the asynchronous phase is presented by even more far-off array of points. Not10able 2: Lyapunov exponents and N async /N ratio for certain values of parameters.Fig. d/R N async /N λ λ λ sync λ sync λ async λ async λ tr .
18 0 .
037 0 . − . . − . – – − . .
20 0 .
150 0 . − . . − . . . . .
27 0 .
338 0 . − . . − . . . . .
275 0 . . . . − . . . . .
30 0 .
636 0 . . − . − . . . . .
36 0 .
900 0 . . − . − . . . . Figure 6: Projections of Poincar´e maps of attractors on the ( r , u , r ) space at (a) d/R =12 . , N async /N = 0 . the asynchronous component is almost absent and barely visible, (b) d/R = 12 . , N async /N = 0 . asynchronous component becomes a little more noticableonly they reach further away from the synchronization manifold, but also their density closerto the synchronization manifold increases significantly. At this point, the synchronous andasynchronous phases are roughly balanced in terms of time a trajectory spends in them, andby their impact on the overall dynamics on the attractor. The character of motion is very closeto becoming hyperchaotic.In the right side of the studied interval in Fig. 3a, one can see that the attractor becomesalmost asynchronous, i.e. N async /N becomes significantly larger than N sync /N and the trajec-tory rarely comes back to the synchronization manifold, spending most of the time outside of it.Thus, after a certain point, the dynamics in the asynchronous phase almost completely deter-mines the character of motion in the whole attracting area. Since the regimes of dynamics in theasynchronous phase are always hyperchaotic, the motion averaged across the whole attractor11igure 7: Projections of Poincar´e maps of attractors on the ( r , u , r ) space at (a) d/R =12 . , N async /N = 0 . there is substantial presence of the asynchronous component, theasynchronous points occur a little bit further away from the synchronization manifold, (b) d/R = 12 . , N async /N = 0 . , asynchronous points can go further away from the synchro-nization manifold, and their density has increasedbecomes hyperchaotic as well. This step of the proposed scenario corresponds to the last stagebefore the blowout bifurcation and corresponds to the scheme in Fig. 1c. Such attractors arepresented in Fig.-s 8a,b. They correspond to d/R = 12 . and d/R = 12 . respectively, thevalues of the largest Lyapunov exponents and the largest transversal Lyapunov exponent canbe found in Table 2. From Fig.-s 8a,b we see that asynchronous points spread out significantlyfurther away from the synchronization plane than before and the synchronous points are barelyvisible behind all the asynchronous ones. At d/R = 12 . , N async /N is equal to . (seeFig. 8a), which means that there is still significant presence of the synchronous component, butnonetheless the motion is already hyperchaotic. From Fig. 8b one can observe that althoughvisually asynchronous points do not extend further away for the synchronization manifold, incomparison to the asynchronous points in Fig. 8a, their density has become much higher.However, the trajectory still occasionally returns to the synchronization manifold after orbitingoutside of it for a long time. This means that we do not observe the blowout bifurcation in ourscenario. Since the largest transversal Lyapunov exponents are positive, these attractors corre-spond to weakly asynchronous states. Let us also remark that in the right side of the interval inFig. 2b the two largest Lyapunov exponents of the synchronous component start decreasing andboth λ sync and λ sync become negative for attractors in Fig.-s 8a,b. This happens because themajority of the (3 , orbits previously forming the skeleton of the synchronous attractor havebecome transversally unstable (2 , orbits (making most parts of the synchronization manifold12igure 8: Projections of Poincar´e maps of attractors on the ( r , u , r ) space at (a) d/R =12 . ( N async /N = 0 . , (b) d/R = 12 . , N async /N = 0 . , the asynchronous component isdominating, the synchronous points can be barely seen behind the scattered asynchronous onestransversally unstable). Thus, a lot of these unstable orbits now contribute into the ejection ofa trajectory outside of the synchronization manifold or prevent a trajectory from coming backinto it, instead of contributing into the one-dimensional instability inside the synchronizationmanifold. In Fig. 3a it is marked by the last growth of N async /N (with simultaneous drop of N sync ), after which N async /N ∼ . making N sync very small.One can also notice the drop of the values of Lyapunov exponents at the very right end ofthe studied interval in Fig. 2a. This corresponds to the crisis of the hyperchaotic attractor anda leap to an asynchronous limit cycle of period 4.We would also like to note that the attractor we are studying here is relatively stable tosymmetry breaking in model (1), corresponding to the transition from two identical bubbles toa pair of slightly different. Let us introduce a small parameter that characterizes this transitionas follows [12] R = εR , | (cid:15) − | (cid:28) , ε > . If ε (cid:54) = 1 , equations (1) are no longer symmetricalwith respect to R ←→ R transformation.Now we study the influence of symmetry breaking on several attractors representing typicalstages of the bubbling transition process. We numerically continue these attractors with respectto ε and calculate the corresponding Lyapunov spectra. We present the results of this contin-uations for four attractors along the way of the bubbling transition scenario: the synchronouschaotic attractor from Fig. 4b, d/R = 12 . , (see Fig. 9 a); the one with low presence ofthe asynchronous component from Fig. 6b, d/R = 12 . (see Fig. 9b); the attractor withmoderate presence of the asynchronous component from Fig. 7a, d/R = 12 . (see Fig. 9c);the hyperchaotic attractor for which the asynchronous component is dominant from Fig. 8b,13 .98 0.99 1 1.01 1.02ɛ-0.2-0.100.05 f 2 1 d -0.2-0.100.1 f 2 1 ɛ f 2 1 ɛ0.98 0.99 1 1.01 1.02-0.2-0.100.1 f 2 1 ɛ b a c Figure 9: The dependence of Lyapunov spectra on ε for continuation of (a) the completely syn-chronous attractor attractor at d/R = 12 . , (b) the attractor the in middle of the bubblingtransition process at d/R = 12 . , (c) the attractor in the middle of the bubbling transi-tion at d/R = 12 . , (d) the hyperchaotic attractor in the end of the bubbling transition at d/R = 12 . . d/R = 12 . (see Fig. 9d). The results for all these attractors are quite similar, except thatdynamics is hyperchaotic in some neighborhood of ε = 1 on the last two graphs.In [12] attractors were considered unstable with respect to symmetry breaking if they disap-pear through a crisis or underwent other sharp qualitative changes of dynamics almost instantlyafter the symmetry breaking is introduced (i.e. | ε − | (cid:28) . ). This corresponds to deviationsof the equilibrium radii from the symmetrical state approximately on the order of . or less.On the other hand, if an attractor is stable with respect to perturbation in ε on the order of or more (i.e. ε ∼ ± . ), it was considered stable to symmetry breaking because this isa sensible margin of error from a physical point of view. Here we follow the same approach.Graphs in Fig.-s 9a-d are very similar, and the sharp qualitative change in type of behaviorcan be observed at the values of ε slightly smaller then . , which corresponds to deviations14rom symmetrical case on the order of . − . . These regions of ε , though relatively small,can be considered as areas of stability with respect to symmetry breaking in comparison withthe cases when attractors disappear with arbitrary small perturbations in ε . Therefore, weconsider that all the discussed above attractors are stable with respect to symmetry breaking.As a consequence, dynamical regimes emerging in the entire bubbling transition scenario mightbe observed in a system of two slightly nonidentical contrast agents and the bubbling transitionprocess may be observed experimentally.Finally, we would like to note that the bubbling transition provides a new bifurcationscenario for the appearance of hyperchaos in system (1). This scenario complements the one thatis based on the appearance of a homoclinic chaotic attractor containing a saddle-focus periodicorbit with its two-dimensional unstable manifold [11] (see, also [21, 22] for other applications ofthe latter scenario). In this work we have studied the process of the destruction of a synchronous chaotic attractorin a model of oscillations of two interacting identical ultrasound contrast agents. We havephenomenologically described the main steps and driving mechanisms of the correspondingscenario. We have numerically shown that the implementation of such bifurcation scenariotakes place in the studied model.We have found that trajectories can be ejected from the synchronization manifold, whichconfirms the existence of transversal instabilities inside it. Moreover, we have demonstrated thatthe asynchronous phase is characterised by two positive Lyapunov exponents, which implies thepresence of (2 , -type saddle orbits in the saddle set outside of the synchronization manifold.Thus, this asynchronous saddle set is hyperchaotic. We have provided Poincar´e maps of indi-vidual attractors that highlight certain steps of the scenario and support the general pictureshown in Fig.-s 2a-3b.Implementation of the bubbling transition scenario discussed here is similar to the onepresented in work [20]. However, in [20] the asynchronous saddle set was simply chaotic.In our case this set is hyperchaotic and, as a result, motion in the asynchronous phase isalso hyperchaotic. Consequently, proposed scenario can be considered as another route tohyperchaos in systems of coupled oscillators.We have also studied the stability of the proposed scenario with respect to the symmetrybreaking perturbation. We have numerically shown that several attractors at the main steps ofthe scenario are stable with respect to symmetry breaking and their behaviour after symmetrybreaking is quite similar. This means that the entire scenario could be observed at values of ε close but not equal to 1, that correspond to the case of two bubbles with slightly differentequilibrium radii. This suggestion is a little counterintuitive, since for ε (cid:54) = 1 the symmetry ofsystem (1) no longer exists together with the synchronous attractors and the F ix ( S ) manifold,that play an important role in the scenario described in Section 3. Nevertheless, it seemsthat small shifts and deformations of the orbits in Fig. 1 do not eliminate the possibility ofimplementation of the scenario. Therefore, we believe that the dynamical regimes and scenariosdiscussed here may be observed in physically realistic systems.We would also like to note that hyperchaotic oscillations of contrast agents may be beneficialfor various applications, for example, for applications related to targeted drug delivery, whena quick rupture of the bubbles is desirable [3, 4, 6]. Moreorver, hyperchaotic attractors aremore robust to various perturbations including symmetry breaking. Therefore, we believe thatfinding new routes to hyperchaos may be important for certain applications of the contrastagents. 15 cknowledgments Authors are grateful to Alexey Kazakov and anonymous referees for useful discussions andremarks that helped to significantly improve the manuscript. This work (except for Section 3)was supported by Russian Science Foundation grant no. 19-71-10048. The results in Section3 were supported by Laboratory of Dynamical Systems and Applications NRU HSE, of theMinistry of science and higher education of the RF grant ag. N. 075-15-2019-1931
References [1] T. Szabo,
Diagnostic Ultrasound Imaging: Inside Out (Academic Press, Cambridge, 2004).[2] F. Goldberg, B.B., Raichlen, J.S., Forsberg,
Ultrasound Contrast Agents: Basic Principlesand Clinical Applications (Martin Dunitz, London, 2001).[3] L. Hoff,
Acoustic characterization of contrast agents for medical ultrasound imaging (Springer, Berlin, 2001).[4] A.L. Klibanov, Invest. Radiol. 41, (2006).[5] C.C. Coussios and R.A. Roy, Annu. Rev. Fluid Mech. 40, 395 (2008).[6] J.M. Carroll, M.L. Calvisi, and L.K. Lauderbaugh, J. Acoust. Soc. Am. 133, 2641 (2013).[7] U. Parlitz, V. Englisch, C. Scheffczyk, and W. Lauterborn, J. Acoust. Soc. Am. 88, 1061(1990).[8] S. Behnia, A. Jafari, W. Soltanpoor, and O. Jahanbakhsh, Chaos Soliton. Fract. 41, 818(2009).[9] I.R. Garashchuk, D.I. Sinelshchikov, and N.A. Kudryashov, Regul. Chaotic Dyn. 23, 257(2018).[10] C.A. Macdonald and J. Gomatam, Proc. Inst. Mech. Eng. Part C-Journal Mech. Eng. Sci.220, 333 (2006).[11] I.R. Garashchuk, D.I. Sinelshchikov, A.O. Kazakov, and N.A. Kudryashov, Chaos AnInterdiscip. J. Nonlinear Sci. 29, 063131 (2019).[12] I.R. Garashchuk, A.O. Kazakov, and D.I. Sinelshchikov, Nonlinear Dyn. 101, 1199 (2020).[13] P. Ashwin, J. Buescu, and I. Stewart, Phys. Lett. A 193, 126 (1994).[14] E. Ott and J.C. Sommerer, Phys. Lett. A 188, 39 (1994).[15] P. Ashwin, J. Buescu, and I. Stewart, Nonlinearity 9, 703 (1996).[16] S.C. Venkataramani, B.R. Hunt, and E. Ott, Phys. Rev. E 54, 1346 (1996).[17] A. Pikovsky, M. Rosenblum, and J. Kurths,