Hidden Attractors in a Model of a Bubble Contrast Agent Oscillating Near an Elastic Wall

A model describing the dynamics of a spherical gas bubble in a compressible viscous liquid is studied. The bubble is oscillating close to an elastic wall of finite thickness under the influence of an external pressure field which simulates a contrast agent oscillating close to a blood vessel wall. Here we investigate numerically the coexistence of chaotic and periodic attractors in this model. One of the tools applied for seeking coexisting attractors is the perpetual points method. This method can be helpful for localizing coexisting attractors, occurring in various physically realistic ranges of variation of the control parameters. We provide some examples of coexisting attractors to demonstrate the importance of the multistability problem for the applications.


Introduction
Ultrasound contrast agents are gas-filled bubbles of the order of 1-10 µm in diameter. Often they are encapsulated in a shell which can be a lipid, protein or polymer. Contrast agents are used in diagnostic ultrasound imaging because they can improve echogenecity. Bubble response to external pressure field, which is generated by the ultrasound, creates scattering sound waves. Some types of bubble response generate signals with recognizable signature. There are various works about the bifurcation structure of the contrast agent models, see e.g. [1,2] and references therein. However, the modified Rayleigh-Plesset equation for the description of a bubble oscillation near an elastic wall (see [3,4]) has not been studied numerically in detail. Moreover, the coexisting attractors governed by this equation have not been investigated at all. The choice of the initial data becomes a non-trivial problem when attractors can coexist. We apply the perpetual points method to overcome this problem [5,6]. It does not require calculations of the attractors basins and, thus, is much less CPU consuming. The system of equations determining the perpetual points consists of rational equations with respect to R,Ṙ and can be efficiently solved numerically.

Main equation and results
We investigate the following equation [3,4]: where R is the bubble radius, R 0 = 10 µm is the equilibrium bubble radius, P stat = 100 kPa is the static pressure, P v = 2.33 kPa is the vapor pressure, the notation P 0 = P stat − P v is used to simplify the equation view, P ac is the magnitude of the pressure of the external field, σ = 0.0725 N/m is the surface tension, ρ 1 = 1000 kg/m 3 is the density of the liquid inside the blood vessel, ρ 2 = 1060 kg/m 3 is the blood vessel wall density, ρ 3 = 1000 kg/m 3 is the density of the fluid surrounding the blood vessel, η = 0.001 Ns/m 3 is the viscosity of the liquid, c = 1500 m/s is the sound speed, γ = 4/3 is the polytropic exponent (the process is assumed to be adiabatic), ν is the Poisson's ratio, ν = 0.5 for an elastic wall and ν = 0 for a rigid wall, β = ρ 2 ν/(1 − ν) is a characteristic of the wall, h = 1 mm is the thickness of the wall, d denotes the distance between the wall and the center of the bubble and is close to R 0 . Let us numerically investigate the dynamics of the gas bubbles governed by equation (1). An attractor is called hidden if its basin of attraction does not intersect with small neighborhoods of equilibria. The dynamical system corresponding to (1) does not have any fixed points. Thus, all the attractors of system (1) are hidden according to this definition. However, here we focus on the scenarios of multistability, which are the most interesting. The perpetual points are defined by the system of equationsR = 0,Ṙ = 0. These points can be used as the initial conditions to achieve different attractors of the system [5,6]. Often, all attractors of a dynamical system can be obtained this way. However, sometimes we cannot use some of these points as the initial data, because the problem will be ill-posed. For example, one of the points can have ||Ṙ|| > c or R 1. Thus, the method does not always completely solve the problem of localizing all of the attractors.
The dynamics described by equation (1) can be quite complex. It can be shown that there can be a single periodic attractor, single chaotic attractor, coexisting periodic attractors and a chaotic attractor coexisting with a periodic attractor. It seems that the most interesting scenario is the one when periodic and chaotic attractors coexist. Examples of this behavior are given in Figs. 1 and 2. Points in the Poincare section are defined by the crossing of the phase trajectory with planes t = T k, where k ∈ N and T = 2π/ω. The original phase space is 3-dimensional, but here we demonstrate 2-dimensional projection on the (R/R 0 ,Ṙ) plane. We also place points of the Poincare section onto the phase portraits to demonstrate the period of the solution with respect to the period of the driving force (bold dots on the phase portraits). Another interesting scenario is happening in the following parameters range: P = 90 kPa, ω/2π ∈ [110, 125] kHz. Over the entire range there exists a 1-periodic attractor, similar to that shown in Fig. 2.
At ω/2π ≥ 120 kHz it coexists with a 2-periodic attractor (see Fig. 3, plate a). If the driving frequency is decreasing below 125 kHz, the 1-periodic attractor remains almost the same in the whole range, but the 2-periodic attractor goes through the period-doubling cascade and becomes a chaotic one, and then at approximately 110 kHz it becomes 2-periodic again (see Fig. 3, plate b). The corresponding bifurcation diagram is shown in Fig. 3, plate c. For any fixed control parameters in the set being discussed, there exists only one perpetual point that can be used as initial data. The diagram in Fig. 3, plate c, was obtained by applying this point. Now we consider other bifurcation diagrams for the same parameters region (see Fig. 4), but for another initial conditions. For these diagrams we use the state of the system at the last time step of the solution for the current frequency as the initial condition for the next step with respect to the frequency. For the diagram in Fig. 4, plate a, we use the increasing frequency sequence, and for plate b -the decreasing frequency sequence. For the first step with respect to the frequency, we use the perpetual point mentioned above as the initial data.
It can be seen that there are a lot of abrupt shifts from chaotic to periodic behavior and vice versa. This happens because the basins of the attractors deform under the change of the parameters of equation (1). If such a change is large enough, a point which was in a small neighborhood of one of the coexisting attractors can appear in the basin of another attractor at the next value of the parameter. The step size in the frequency axis, for which the diagrams are calculated, are about ∼1% of the whole range. From the physical point of view, these frequency changes can be treated  as finite disturbances in the system. Such finite perturbations always exist in real physical systems along with finite measurement errors in the control parameters. This fact highlights a problem for the applications: one can see that a disturbance of the order of ∼1% can cause unexpected shifts from a chaotic motion to a periodic one and vice versa. It is worth noting that signficant reduction of the step size will cause the system to stay on the same attractor. If we make the frequency steps ∼ 3-4 times smaller, we can achieve diagrams, completely according to the 1-periodic attractor, or to the another one (like in Fig. 3 (c)), if we start from the necessary attractor at the first step.

Conclusion
In this work we have studied an equation describing oscillations of a spherical bubble contrast agent near an elastic or rigid wall. We have numerically studied the considered model and have demonstrated that is exhibits complex dynamics. We have discussed some regions of parameters where chaotic and periodic attractors coexist. We have demonstrated that this problem can be important for applications, because the multistability can lead to unpredictable behavior in the bubble oscillations. These effects are undesirable for applications. Consequently, the corresponding areas of the control parameters should be avoided in practice.