Nonlinear dynamics and bifurcation structure of ultrasonically excited lipid coated microbubbles

Highlights • Bifurcation structure of the lipid coated microbubbles (MBs) is studied• Effect of the coating is analyzed by comparing the dynamics of the lipid coated MBs to uncoated ones.• Buckling and rupture of the lipid coating enhances the nonlinear oscillations.• Period 2, Period 3 or chaos can be generated at pressures as low as 1 kPa.• In addition to subharmonic emissions, superharmonic emissions may also be enhanced.• With increasing pressure, period 2 is generated, then disappears and is regenerated again.

One of the first nonlinear phenomena detected with MBs in sound fields was through historical observations of Esche [21]. Esche reported the generation of a frequency peak at half the excitation frequency (f) in the power spectrum of the received signal [21]. In his investigation of MBs driven with 3Hz-3.3 MHz, he found the appearance of spectral lines at f/2 and in some cases f/3 for sufficiently high acoustic pressures. In a continuation of Esche's work, Bohn reported spectral lines down to f/4 [22]. In the chaotic (broadband noise) region of the sound emitted by the MBs, Holzfuss & Lauterborn [3] observed a surprisingly lowdimensional attractor with correlation dimension of about 2.5 which is the characteristic for driven damped nonlinear oscillators. Several other experimental studies investigated the nonlinear dynamics of ultrasonically excited MBs [23][24][25][26][27][28][29]; observing subharmonics, ultraharmonics and chaotic behavior. Numerical investigations have demonstrated the existence of multiple resonance peaks [2,30,31], period doubling route to chaos [32][33][34], strange attractors and chaotic behavior (e.g. [1][2][3]28,29]).
Within the last decade several studies have employed the methods of dynamical systems to study the behavior of MBs. There have been successful attempts in classification of some of the nonlinear dynamics of the MB oscillator [35][36][37][38]. Hegedüs [39] found numerical evidence for the existence of stable period 1 solutions beyond Blake's threshold [39]. Occurrence of higher order subharmonics (SHs) (f/3, f/4, f/5 etc) has been extensively investigated in [38,40] and for the case of ambient pressures slightly below the vapor pressure [39]. They are experimentally observed and numerically modelled in [41][42][43].
Hegedus [44] studied the topology of stable periodic solutions near Blake's threshold. The effect of high dissipation on the nonlinear evolution of the MB behavior is considered in [45,46] and it has been shown that MB becomes an over-damped oscillator suppressing collapse-like behavior. Moreover, they reported the existence of transient chaos [45]. Using two frequencies was proposed in [47] and extended in [48][49][50][51] to control the chaotic behavior of the MBs. The effect of multiple frequencies on the resonance behavior and nonlinear dynamics of the system was investigated in [52][53][54][55].
The influence of the pressure amplitude on the resonance frequency and bifurcation structure of the MBs which is driven by its resonance frequency is studied in [56]. It is shown that increasing the incident ultrasound pressure decreases the resonance frequency of the MB; when the MB is sonicated with its pressure dependent resonance frequency a saddle node bifurcation takes place at the corresponding pressure amplitude which enhances the nondestructive back-scattered pressure by the MBs. Non-spherical MB oscillations in a viscous liquid is studied in detail in [46] and its been shown that the increased rate of dissipation can significantly extend the stable domains in the acoustic excitation parameter planes. We have studied the ultraharmonic (UH) and super harmonic (SuH) behavior of the MB oscillator by introducing a more comprehensive method of construction of bifurcation diagrams [57]. Using this method, the bifurcation structure of the MBs undergoing period doubling and 1/2 order sub-harmonic emissions have been extensively studied [34]. It was found that sonication of MBs with twice their linear resonance frequency results in period doubling at a lower excitation and leads to non-destructive stable period 2 oscillations, however, sonication with resonance will most likely result in MB destruction before the appearance of period 2 oscillations. We showed in [58] that SH resonance frequency decreases with increasing pressure; and maximum SH strength is generated when the sonication frequency is 1.5-1.6 times the resonance frequency of the MBs.
In spite of numerous studies on the complex behavior of free (uncoated) MBs, the dynamics of the coated MBs have not been thoroughly studied. MBs stabilized by a coating in the form of phospho-lipid (e.g. Definity® [59]), or albumin (e.g. Optison [60]) or polymer (Point [61]) are designed to be used in clinical and pre-clinical medical ultrasound applications. Addition of the coating (more specifically in case of phospho-lipid coating) immensely increases the complexity of the MB oscillator. During MB oscillations phospho-lipid shell can undergo buckling and rupture [62] resulting in a dynamical system with varying stiffness. The dynamic stiffness of the nonlinear oscillator enhances the generation of nonlinear signatures in the oscillation of the coated MBs.
Buckling of the lipid shell has been shown to be one of the possible reasons for enhanced non-linearity [62][63][64][65][66][67][68][69][70][71]. Phospho-lipid shell MBs exhibit compression only behavior [67] during which MBs compress more than they expand. There exists a threshold behavior for the onset of oscillations [72]; the MB starts to oscillate only above a pressure threshold. It has been experimentally observed that phospholipid shell MBs can generate SH oscillations even at very low acoustic pressures (<30 KPa [64,66,73]). Such low threshold values not only contradict the predictions of the theoretical models for coated MBs [65,74,75], they are even below the threshold values expected for uncoated free MBs [65,76]. The low pressure thresholds are despite the increased damping due to the presence of the shell. Through experiments and numerical simulations it has been shown in [64] that the low pressure threshold for SH emissions is due to the compression only behavior of the MBs due to the buckling of the shell.
In [68] the lipid shell was found to enhance the nonlinear MB response at acoustic pressures as low as 10 kPa. The increase in acoustic pressure amplitude lead to a substantial decrease of the frequency of the maximum response even at very low acoustic pressures [68] resulting in a pronounced skewness of the resonance curve. Such shift in resonance has been postulated in [68] to be the origin of the 'thresholding' behavior [72]. Nonlinear resonance behavior of the lipid shell MBs was also observed in higher frequencies (5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15) in [77]. It is shown in [66] that the shell elasticity of the phospholipid shell varies with MB oscillation amplitude and the magnitude of 'compression only' behavior depends on the initial phospholipid concentration on the MB surface. Prosperetti [65] through theoretical analysis of the Marmottant model [62] attributed the lower SH threshold of the lipid MBs to the variation in the mechanical properties of the shell in the neighborhood of a certain MB radius (e.g. the occurrence of buckling).
In addition to the widely studied 1/2 order SHs, we have experimentally detected higher order SHs (1/3, 1/4 and 1/5) in the oscillations of lipid coated MBs at very low acoustic pressures and high frequencies (e.g. 25 MHz) [41][42][43]. Through analyzing bifurcation diagrams we concluded that buckling or rupture of the shell is responsible for the enhanced nonlinear behavior [41][42][43]. The closer the initial surface tension of the MB to the two limit values of the buckling and rupture of the shell, the lower the pressure threshold for nonlinear oscillations. Variation of the mechanical properties of the shell can also manifest itself in expansion dominated behavior in liposome-loaded lipid shells [69]. Expansion dominated oscillations occur for MBs with an initial surface tension near that of water [69,77]. Upon expansion, the stiffness of the coating weakens and the MB expands more than it compresses. Expansion-dominated behavior was used to explain the enhanced nonlinearity at higher frequencies (25 MHz) [78]. The Marmottant model effectively captures the behavior of the MB including expansiondominated behavior [59,43,77], compression only behavior [67], thresholding [72] and enhanced non-linear oscillations at low excitation pressures [41,43,63,64,66,68].
Previous studies (e.g. [2,34,36,[38][39][40]56,58]) investigated the bifurcation structure of the uncoated bubbles and bubbles coated with shells that exhibit linear viscoelastic behavior. However, the influence of the nonlinear viscoelastic behavior of the coating (e.g. buckling and rupture [59,62]) on the bifurcation structure of the bubble has not been investigated before. Due to the enhanced nonlinearity created by the behavior of the shell, it is important to rigorously investigate the impact of the exposure parameters on the MB oscillations. The current work addresses this problem for the first time. We perform a comprehensive analysis of the bifurcation structure of ultrasonically excited lipid coated MBs. Similar to our previous works in [34,38,56,58] we study the radius vs excitation pressure amplitude bifurcation structure of the lipid coated bubbles when the bubble is sonicated with multiples of its resonance frequency. Results are then compared to previous studies whereby we reveal the influence of the nonlinear shell viscoelasticity on the bubble behavior. We show that the buckling and/or rupture of the shell enhances the subharmonics (SHs), superharmonics (SuHs), ultraharmonics (UHs) and chaos at very low excitation pressures. The enhanced nonlinearity may disappear at moderate pressures. At higher pressures, nonlinear behavior may reappear in the bubble behavior exhibiting similar behavior to the uncoated bubbles and coated bubbles with linear viscoelastic behavior.
Knowledge of the effect of the shell behavior on the nonlinear response of the MB is essential to optimize the MBs response to an ultrasonic field. Moreover, the comprehensive knowledge that can be obtained through analyzing the bifurcation diagrams of the lipid coated MBs may help in revealing potential parameter spaces in which MB behavior can be beneficial to various applications. Last but not least, revealing the intricate behavior of the system and enhanced nonlinear effects is of potential interest in the field of nonlinear and chaotic dynamical systems.

Marmottant model
The dynamics of the coated MBs undergoing buckling and rupture can be effectively modeled using the Marmottant equation [62]: In this equation, R is the bubble radius at time t, R 0 is the initial bubble radius, Ṙ is the wall velocity of the bubble, R is the wall acceleration, ρ is the liquid density (998 kg m 3 ), c is the sound speed (1481 m/s), P 0 is the atmospheric pressure, σ(R) is the surface tension at radius R, μ L is the liquid viscosity (0.001 Pa s), k s is the viscosity of the coating, k is the A.J. Sojahrood et al. polytropic index for the gas and P a (t) is the acoustic driving force P a (t) = P a sin(2πft) where P a and f are the amplitude and frequency of the applied acoustic pressure. The values in the parentheses are for pure water at 293 0 K. In this paper the gas inside the bubble is C3F8 and water is the host media.
The surface tension σ(R) is a function of radius and is given by: where σ 0 is the initial surface tension (at R = R 0 ), σ water is the water surface tension and χ is the shell elasticity. R r and R b are the rupture and the buckling radius respectively where In this work similar to [78], R breakup = R r . In this paper, simulations were run for different values of σ 0 . The initial surface tension σ 0 is a property of the lipid coated bubble and varies when using different manufacturing methods [79,80]. Moreover, σ 0 can be altered by varying the ambient pressure in the liquid [64,80]. Variations in σ 0 changes the R b and R r which in turn change the dynamical behavior of the bubble. In this paper, for simplicity, and similar to [78,80] we have assumed σ rupture = σ water . Fig. 1 shows a representation of buckling and rupture and the dependence of the effective surface tension (σ(R)) on microbubble radius.

Keller-Miksis model
Dynamics of the uncoated bubbles were also visualized alongside the lipid coated bubbles to highlight the effect of the lipid shell on the bubble dynamics. To model the uncoated bubble dynamics the Keller-Miksis model [79] is used: where G = . In both models we have neglected the effects of thermal damping. This is to decrease the problem complexity and to better highlight only the shell effects. Moreover, we have shown in [89] that in case of C3F8 gas cores thermal damping is significantly smaller compared to air. Moreover, in case of coated bubbles with C3F8 gas cores, thermal effects maybe be fully neglected. However, in case of the uncoated bubble effects of thermal damping at higher pressures should be considered using full ODEs [90] that account for the thermal damping. We have shown in [89] that the generally used linear assumptions [91] for thermal effects may lead to inaccuracies at pressures as low as ≈ 40 kPa. However, since the main focus of the paper is to highlight the coating effects and because the thermal effects of the C3F8 are weak [88], we have neglected the thermal effects in this paper.
It should be noted that the Keller-Miksis model (Eq. (3)) has some additional terms compared to the the Marmottant model (Eq. (1)). In this paper, the purpose of the qualitative comparison between the two models is to demonstrate the effect of shell on the MB dynamics. The behavior of the bubble in the absence of shell is used as a reference to reveal the effect of the enhanced non-linearities due to coating at low excitation pressures. The Marmottant model [62] is written a popular form [59] and addresses the inaccuracies of the Keller-Miksis model when |Ṙ| c ≈ 1. We show in Appendix A, that in the absence of the shell terms in the Marmottant model the radial oscillations of the bubbles as predicted by the two models are relatively in good agreement. Radial oscillation amplitude of the periodic behavior, the pressure threshold of the onset of various nonlinear regimes and the chaotic behavior are in good agreement. However, when the oscillations are chaotic the radial oscillation amplitudes as predicted by both models are not always equal.

Investigation tools
Bifurcation diagrams are valuable tools to analyze the dynamics of nonlinear systems since qualitative and quantitative changes of the dynamics of the system can be investigated effectively over a wide range of control parameters. In this paper, we employ a more comprehensive bifurcation analysis method introduced in [73,74].

Conventional bifurcation analysis (Poincaré cross section at each driving period)
When dealing with systems responding to a driving force, to generate the points in the bifurcation diagrams vs. the control parameter, one option is to sample the R(t) curves using a specific point in each driving period. The approach can be summarized in: where P denotes the points in the bifurcation diagram, R and Ṙ are the time dependent radius and wall velocity of the bubble at a given set of control parameters of (R 0 ,P 0 ,P A ,c,k,μ, σ,f) and Θ is given by n f . Points on the bifurcation diagram are constructed by plotting the solution of R(t) at time points that are multiples of the driving acoustic period. In this work, the results are plotted for n = 100 − 150 to ensure a steady state solution has been reached.

Method of maxima
As a more general method, bifurcation points can be constructed by setting one of the phase space coordinates to zero: In this method, the steady state solution of the radial oscillations for each control parameter is considered. The maxima of the radial peaks (Ṙ = 0) are identified (determined within 100-150 cycles of the stable oscillations) and are plotted versus the given control parameter in the bifurcation diagrams. The bifurcation diagrams of the normalized bubble oscillations ( R R0 ) are calculated using both methods a) and b). When the two results are plotted alongside each other, it is easier to uncover more important details about the SuH and UH oscillations, as well as the SH and chaotic oscillations.

Resonance curves
Compared to uncoated bubbles and coated bubbles with pure viscoelsatic behavior (e.g. Keller-Miksis model [81], Hoff model [82], Morgan model [83]), the resonance behavior of lipid coated bubbles are more complex. This is due to the buckling and rupture of the shell and dynamic variation of the effective surface tension of the bubble. As an example [68,77] have shown numerically and experimentally that a pressure increase leads to a significant displacement of the main resonance (frequency of maximum response) of the bubble leading to a significant shift of the resonance curve.  Similar to our previous work in [34,37,38,56], in this work we will attempt to classify the nonlinear dynamics of the lipid bubbles as a function of pressure amplitude when they are sonicated with fractions or multiples of their f r . However, the initial sharp decrease of the resonance frequency with pressure will make the classification difficult. Moreover, characterization of the coating parameters of the bubbles in experiments are generally through attenuation measurements of the bubble solution when there is an excitation pressure amplitude above 1 kPa is applied. As an instance negative peak pressure amplitude of 25 kPa, 12.5 kPa, 30 kPa, 10 kPa & 5 kPa were applied respectively in [85][86][87]80,84] and peak to peak pressures of 33 kPa were applied in [87]. Very low pressures can not be applied experimentally due to the signal to noise constraints of the measurements systems.
To simplify the classification method and to have a better comparison with published experimental data we have calculated the resonance frequency at P a = 10 kPa and used it for further study. Thus, in this paper for coated bubbles f r refers to the resonance frequency at P a = 10 kPa.

Radial oscillations as a function of time and the corresponding changes in the σ(R)
In Fig. 2, we observed the generation of SuH as well as SH resonances at very low pressures in case of the coated bubbles. In this section, the enhanced nonlinear oscillations and their relationship with the bubble surface tension are briefly investigated to have a better insight on the mechanisms of enhanced nonlinearity. Fig. 5 shows the radial oscillations of the uncoated bubble as a function of 10 acoustic driving periods (100)(101)(102)(103)(104)(105)(106)(107)(108)(109)(110). The left column shows the radial oscillations when P a = 1 kPa and f = 0.3f r , 2f r and 3f r in Fig. 5a, c and e respectively. right column shows the radial oscillations when P a = 60 kPa and f = 0.3f r , 2f r and 3f r in Fig. 5b, d and f respectively. The red circles locate the amplitude of the radial oscillations at each period. This is the Poincaré cross section at each driving period which is used to generate the bifurcation diagram using the method introduced in 2.3.1. The bubble oscillations in Fig. 5 are period 1 (P1) and the red circles have the same value at each driving periods. This indicates the absence of any SHs. Only 3rd order SuHs are seen (P1 oscillations with 3 maxima) when pressure amplitude is 60 kPa in Fig. 5b.
Red circles correspond to the location of R(t) at each period.

Bifurcation structure of the uncoated bubble
In this section, we briefly highlight the main nonlinear regimes of the dynamics of the uncoated bubble as a function of pressure amplitude at different frequencies. This data will be useful when analyzing the behavior of the lipid coated bubble by highlighting the shell effects on the coated bubble dynamics. The red curve is constructed using the method of maxima (Section 3) and the blue curve is constructed using the Poincaré cross section at each driving period. Fig. 11a shows the bifurcation structure of the uncoated bubble with R 0 = 2μm sonicated with f = 0.3f r . Pressure increase above ≈ 50 kPa leads to the generation of 3 maxima in the bubble oscillations (3 blue lines) for a period 1 (P1) oscillation regime. Thus 3rd order SuH regime [57] is generated. In the regime of 3rd order SuH oscillations, the frequency component at 3f is stronger than the rest of the frequency components of the scattered pressure. Oscillations undergo period doubling (PD) at about 124 kPa. The blue curve with 3 maxima undergoes PD concomitant with the 1 PD in the red curve; thus oscillations become P2 with 6 maxima and 7/2 order UH oscillations are generated (124kPa < P a < 178 kPa). When 7/2 order UHs occur the frequency component at 3.5f in the scattered pressure by the bubble is stronger than the frequency components at 0.5f, 1.5f, 4.5f, etc. The 3rd order region and the 7/2 order UH region are highlighted as an inset in Fig. 11a. Further pressure increase leads to SN bifurcation to 2nd order SuH oscillations of higher amplitude, followed by 5/2 UHs, and a small chaotic window. Finally a giant P1 resonance emerges out of the chaotic window undergoing further PDs at higher pressures.
When f = 0.5f r (Fig. 11b), as pressure amplitude increases above 14 kPa, 2 maxima are generated in the P1 oscillation regime (2nd order SuH). Further pressure increase results in a PD in both the blue and red graphs leading to a P2 oscillation with 4 maxima (5/2 UH oscillations). This region is highlighted as an inset in Fig. 11b. Chaos occurs in a small window above 160 kPa with a tiny window of periodic (P3 with 5 maxima) behavior within. Afterwards, a giant P1 resonance emerges out of the chaotic window. The P1 oscillations undergo a multiple cascades When f = 0.6f r (Fig. 11c) 5/2 UH oscillations (P2 with 4 maxima) are developed and then transition to P1 oscillations through a bubble in the pressure window of 116-150 kPa (highlighted in an inset). P1 oscillations then undergo a saddle node bifurcation to a P1 oscillation with higher amplitude at P a ≈ 166 kPa. This is due to the pressure dependent resonance behavior that has been discussed in detail in [56]. Further pressure increase leads to a PD to P2 oscillations (at 406 kPa) which is followed by a cascade of PDs to chaos at ≈ 614 kPa.
The dynamics of the bubble sonicated with f = 0.7f r (Fig. 11d) is similar to the case of f = 0.6f r ; however, 5/2 UH oscillations are not generated and SN bifurcation occurs at a slightly lower pressure amplitude (117 kPa). At this pressure amplitude the red curve meets the blue curve. This is the pressure dependent resonance and the wall velocity becomes in phase with the driving signal. This is discussed in detail with numerical and experimental observations in [86]. PD occurs at 326 kPa which is lower than the PD threshold in Fig. 11c. Chaos settles through a cascade of PDs at 504 kPa.
When f = f r (Fig. 11e) oscillations are P1 and the blue line and the red line have the same value (highlighted in the inset) which indicates that the wall velocity is in phase with the acoustic driving force due to the resonance (page 290 in [92]). The two curves start diverging as soon as pressure increases above 18 kPa and at 215 kPa the oscillations undergo PD. Oscillations become chaotic above 400 kPa with a small window of periodic behavior (P3 with 3 maxima).
When f = 1.2f r (Fig. 11f), we witness the similar behavior as the case of f = f r ; however, P2 oscillations are developed for R max /R 0 < 2, thus P2 oscillations are more likely stable [93].
When f = 1.5f r (Fig. 12a), P1 oscillations undergo PD with 2 maxima at 236 kPa. P2 oscillations undergo a SN bifurcation to P2 oscillations of higher amplitude at 347 kPa. The SN bifurcation is coincident with the pressure dependent SH resonance (Pdf sh ) [58]. This results in the oversaturation and enhancement of the SH signal from the pressure scattered by bubbles [58]. P2 oscillations undergo successive PDs to chaos at ≈ 494 kPa.
When f = 1.8f r (Fig. 12b) P1 oscillations undergo a SN bifurcation to P2 oscillations of higher amplitude at 155 kPa. The P2 oscillations amplitude Rmax R0 < 2 thus bubbles may have higher stability compared to Fig. 12a. Further pressure amplitude increase leads to chaos through successive PDs. At 931 kPa a giant P3 resonance emerges out of the chaotic window.
When f = 2f r (linear SH resonance frequency), PD occurs at the lowest pressure threshold of 77 kPa (highlighted in an inset) [34]. P2 oscillations undergo successive PDs and chaos appears at 400 kPa and extends to ≈ 600 kPa where giant P3 resonance emerges out of the chaotic window. Oscillations later become chaotic again through successive PDs. When f = 2.2f r (Fig. 12d), PD occurs at 189 kPa which is higher than the PD threshold when f = 2f r . P2 oscillations undergo PD to P4-2 at 445 kPa and then are followed by chaos through consecutive PDs at 482 kPa.
The case of f = 2.8f r is depicted in Fig. 12e. P1 oscillations undergo a SN to P3 oscillations at 390 kPa. P3 oscillations undergo PD to P6 at 489 kPa and a small chaotic window appears at 587 kPa. Chaos disappears and low amplitude P1 emerges out of the chaotic window at 588 kPa which later undergo a PD similar to Fig. 12d at 661 kPa. Further pressure increase results in the occurrence of P4 through a SN at 819 kPa. P4 oscillations undergo PD to P8 at about 900 kPa.
When f = 3f r (Fig. 12f), P3 occurs at 353 kPa through SN bifurcation. P3 extends to 567 kPa where P6 oscillations are generated through a PD. A small chaotic window appears before the low amplitude P1 which then undergoes a SN to P8 oscillations. Finally chaos is generated at ≈ 800 kPa. initially in the buckled and ruptured states receptively. The bubbles start oscillation in a P1 with two maxima (2nd order SuH) right from P a = 1 kPa. The following evolution 2nd order SuH → 3rd order SuH (P1 with 3 maxima) → 4th order SuH (P1 with 4 maxima) takes place as pressure amplitude increases (these are highlighted as insets in Fig. 13a-b). Compared to the uncoated bubble case, the 2nd order SuH appears at a very small pressure threshold (P a = 1 kPa). Wall velocity is in phase with the driving acoustic pressure for most of the pressures below 200 kPa. Further pressure amplitude increases results in the gradual disappearance of the maxima, and above 210 kPa, only two maxima remain in When f = 0.5f r ( Fig. 13c-d), oscillations start with 2nd order SuH oscillations (P1 with 2 maxima) right from the start at 1 kPa and this stretches to ≈ 50 kPa in both cases at which point 2nd maxima disappears (highlighted as insets in Fig. 13c-d). Compared to the uncoated bubble case, the coating at its ruptured or buckled state reduces the pressure threshold for SuH oscillations. UH oscillations, however, are suppressed and only occur at higher pressures and for a much shorter range of excitation pressures. The pressure threshold for the giant PD increases and chaotic oscillations are suppressed within the excitation pressure amplitude range that is examined here. This can be due to the increased damping due to the coating.
When f = 0.6f r (Fig. 13e-f     When f = 0.7f r (Fig. 13g-h), oscillations start in a similar manner to the case of f = 0.6f r . The growth rate of the P1 oscillation amplitude increases abruptly above a pressure threshold which is lower than the case of the f = 0.6f r (90 kPa for σ 0 = 0 N/m & 29 kPa σ 0 = 0.072 N/m).
Consequently, wall velocity becomes in phase with the driving pressure (highlighted as an inset in Fig. 13)g and further pressure amplitude increases result in the divergence of the blue and the red curve. PD occurs at P a = 495 kPa for σ 0 = 0 N/m & 295 kPa for σ 0 = 0.072 N/m. Chaotic oscillations are finally generated through successive PDs with some periodic windows within.
The cases of the coated bubbles in Fig. 13e-h are similar to the case of the uncoated bubble sonicated with f = 0.6f r & f = 0.7f r (Pdf r [56]). However, the pressure threshold for the SN bifurcation or the increase in the growth rate of the oscillations (inflection point) is much lower in case of the coated bubble with σ 0 = 0.072 N/m despite being excited with lower frequencies. Moreover, the pressure threshold for PD and chaotic oscillations are higher for the coated bubbles with PD occurring at a higher The Rmax R0 . This can be due to the increased damping in the bubble oscillations.
When f = f r (Fig. 14a-b) (note that in this paper in case of the lipid coated bubbles f r was considered the frequency of maximum response at 10 kPa) the red and blue curve have the same value for P a < 20 kPa. The P1 oscillations amplitude grows as pressure amplitude increases and the two curves diverge with amplitude increase. PD occurs at at P a = 267 kPa for σ 0 = 0 N/m & 317 kPa for σ 0 = 0.072 N/m which are higher than the PD pressure amplitude for the uncoated bubble (P a = 215 kPa (Fig. 11)e. Rmax R0 of the P2 oscillations of the coated bubble however, are below 2 while the oscillation amplitude of the P2 oscillations in uncoated bubble are above 2. In case of the bubble with σ 0 = 0 N/m a further pressure amplitude increase leads to P4 oscillations through another PD. P4 oscillations become P8 and then again P4 through a bubbling bifurcation; P4 oscillations later undergo a PD cascade to chaos. At P a ≈ 915 kPa a P4 oscillation emerges out of the chaotic window through reverse PD bifurcation. P4 becomes P2 through another SB. For the bubble with σ 0 = 0.072 N/m, At P a = 600 kPa the P2 oscillations undergo a SN bifurcation to P2 oscillations of higher amplitude. This is similar to the behavior of the uncoated bubble sonicated by its Pdf sh =≈ 1.5 − 1.9f r [58] & Fig. 12a (f = 1.5f r ). Thus, in case of the lipid coated bubble the buckling and rupture of the coating significantly decreases the Pdf sh .
When f = 1.2f r (Fig. 14c-d), the P1 oscillation amplitude increases with increasing pressure amplitude and PD occurs at P a = 314 kPa for σ 0 = 0 N/m & 238 kPa for σ 0 = 0.072 N/m. Pressure thresholds for PD are higher than the pressure threshold of PD (218 kPa) in the uncoated bubble case in Fig. 11f. In both cases, with increasing pressure amplitude a SN bifurcation takes place from P2 to another P2 with higher amplitude (P a = 796 kPa for σ 0 = 0 N/m & 314 kPa for σ 0 = 0.072 N/m). This is similar to the dynamics of the uncoated bubble sonicated by its Pdf sh [58] & Fig. 12a-b (f = 1.5f r & 1.8f r ). This shows that the dynamic variations of the effective surface tension including buckling and rupture decreases the Pdf sh . In the case of σ 0 = 0.072 N/m chaos appears through successive PDs, however, the bubble with σ 0 = 0 N/m does not exhibit chaotic oscillations in this pressure amplitude range. Additionally at a given pressure, Rmax R0 is higher for the bubble with σ 0 = 0.072 N/m because of the expansion dominated behavior of the bubble. This can be one of the reasons for lower pressure threshold of P2 and chaotic oscillations in case of the bubble in ruptured state.
When f = 1.5f r (Fig. 14e-f), the bubble behavior is similar to the uncoated bubble sonicated with its Pdf sh . The pressure threshold for P2 oscillations are P a = 338 kPa for σ 0 = 0 N/m & P a = 98 kPa for σ 0 = 0.072 N/m. In case of the bubble with σ 0 = 0.072 N/m pressure threshold for PD is lower than the case of the uncoated bubble (Fig. 12a). Increasing the pressure amplitude results in a SN bifurcation from a P2 regime to a higher amplitude P2 regime. In case of the uncoated bubble the SN bifurcation results in Rmax R0 > 2, however, here P2 oscillations remain below 2 when SN occurs. The P2 oscillations undergo successive PDs to P8 in both bubbles (Fig. 14e-f). However, only the bubble with σ 0 = 0.072 N/m, exhibits chaotic oscillations. Similar to the previous cases, Rmax R0 is higher for the bubble in the ruptured state due to expansion dominated behavior.
When f = 1.8f r a very interesting phenomenon is observed (Fig. 14gh). In both cases, the bubble starts oscillating in the P2 regime at the very low pressure threshold of 1 kPa. To our best knowledge, such a low excitation threshold for P2 oscillations in nonlinear oscillators is first reported here. The dynamic of the bubble exhibits three interesting stages. The generation of P2 oscillations (at very low pressure), the disappearance of P2 oscillations and regeneration of P2 oscillations. Such behavior has been observed experimentally in [73,94]. In [73], the disappearance of SH oscillations is referred to as an "unexpected standstill" of SHs. This will be discussed further in discussion. Within the initial P2 window, a very small P4 window occurs for both bubbles.  (Fig. 12b); however, the SN occurs at a higher pressure. Similar to the uncoated bubble, after the SN occurrence, the bubble with σ 0 = 0 N/m undergoes chaotic oscillations through successive PDs.
When f = 2f r (Fig. 15a-b), the dynamics are similar to Fig. 14g-h. P2 oscillations are generated at 1 kPa, and they disappear above 200 kPa.
For the bubble with σ 0 = 0 N/m, P2 oscillations re-emerge at ≈ 600 kPa and through a PD bifurcation. Similar to the coated bubble sonicated with its Pdf sh in [58], P2 oscillations undergo a SN bifurcation to P2 oscillations with higher amplitude. Further pressure amplitude increase results in chaotic oscillations through successive PDs. In case of the bubble with σ 0 = 0.072 N/m, soon after the disappearance of the P2 oscillations, a rather small window (293-310 kPa) of P2 oscillations is generated through a SN. P2 oscillations disappear and P1 oscillations undergo a SN to P3 at 707 kPa.This dynamical feature is similar to Fig. 12f where the uncoated bubble is sonicated with f = 3f r .
The dynamics of the bubble sonicated with f = 2.2f r (Fig. 15c-d) is similar to f = 2f r and the general dynamical features of the system stays the same.
The dynamics of the bubbles with σ 0 = 0 N/m & 0.072 N/m sonicated with 1.8f r ⩽f ⩽2.2f r exhibits three main stages. In stage one the bubble shows enhanced non-linearity by which P2 oscillations are generated at very low pressure thresholds. The P2 oscillations disappear by increasing the pressure amplitude however, they re-emerge as P2 or P3 oscillations above a pressure threshold higher than the uncoated counterpart, and in a similar fashion to the uncoated bubble sonicated by its Pdf sh or f = 2.8 − 3f r .
The bifurcation structure of the bubbles when f = 2.8f r is shown in  Fig. 11h. Thus, here we decided to present f = 3.1f r to highlight the generation of P3 at P a = 1 kPa). P3     (Fig. 13a-d). However, there are two differences: 1) Rmax R0 is generally lower than the initially buckled or the ruptured bubble over all pressures studied and, 2) The threshold for the start of SuH oscillations is ≈ 11 kPa which was 1 kPa in (Fig. 13a-d).
The pressure threshold for SuH oscillations is still lower than the case of uncoated bubble in Fig. 11a-b.
Sonication with f = 0.5f r is depicted in Fig. 16c- The case of sonication with f = 0.7f r is shown in Fig. 16g-h. There are two SN bifurcations with pressure amplitude increase. The initial SN occurs at ≈ 15 kPa and results in P1 oscillations of higher amplitude. After the first SN oscillation amplitude grows with increasing pressure amplitude and PD occurs in both cases. A small P4 window is generated within the P2 window. For the case of σ 0 = 0.01 N/m and at P a = 710 kPa P4 oscillations are regenerated and then transition to P2 via reverse PD at 982 kPa. For the bubble with σ 0 = 0.062 N/m at P a = 479 kPa P2 oscillations undergo a SN bifurcation to P2 oscillations with higher amplitudes. This is similar to the dynamics of the uncoated bubble sonicated with its Pdf sh (Fig. 12a-b).
Case of the f = f r is shown in Fig. 17a-b. At P a = 10 kPa a SN bifurcation takes place and oscillation amplitudes increase slightly (Pdf r at 10 kPa). Oscillation amplitude increases slowly with pressure amplitude and PD occurs at P a = 326 kPa & 148 kPa respectively for σ 0 = 0.01 N/ m & σ 0 = 0.062 N/m. After the SN, the dynamics of the bubble with σ 0 = 0.01 N/m (Fig. 17a) & σ 0 = 0.062 N/m (Fig. 17b) sonicated with f = f r are respectively similar to the dynamics of the bubble with σ 0 = 0 N/m (Fig. 14e) & σ 0 = 0.072 N/m (Fig. 14f) sonicated with f = 1.5f r . For the bubble with σ 0 = 0.01 N/m increasing pressure amplitude results in a SN bifurcation from P2 oscillations to a higher amplitude P2 oscillations at 378 kPa. This is similar to the dynamics of the uncoated bubble sonicated with its Pdf sh (Fig. 12a) [58]. P2 oscillations then grow in amplitude with pressure amplitude increase and oscillations become P4-2 through a PD at 624 kPa. Bubble continues with P4 oscillations with a P8 window within, which is created and annihilated through a bubbling bifurcation. The dynamics of the bubble with σ 0 = 0.062 N/m resembles the case of the uncoated bubble sonicated with f = 2f r (Fig. 12c) [34]. P2 oscillations spread between 148 − 555 kPa. At 555 kPa, P4-2 oscillations are generated via a PD and later undergo successive PDs to chaotic oscillations at 638 kPa.
The case for f = 1.2f r is shown in Fig. 17c-d. In both cases we witness the generation of the P2 oscillations, their disappearance and regeneration which is similar to the dynamics of the initially buckled and ruptured bubble in Fig. 14g-h. For bubbles with σ 0 = 0.01 N/m (Fig. 17c) & σ 0 = 0.062 N/m (Fig. 17d) P2 oscillations occur via a PD at P a = 47 & 48 kPa respectively. With pressure amplitude increase P2 oscillations transition to P1 (disappearance of 1/2 order SHs) at P a = 140 & 269 kPa respectively. P2 oscillations are then re-appear at P a = 418 & 545 kPa respectively. Dynamics of the coated bubbles in this pressure amplitude region is similar to the dynamics of the uncoated bubble sonicated with its Pdf sh (Fig. 12b). In case of σ 0 = 0.01 N/m, further pressure amplitude increase results in a cascade of PDs to chaos.
In case of the σ 0 = 0.062 N/m further pressure amplitude increase results in the appearance of P3 oscillations which later undergo PD to P6 oscillations.
The dynamic variation of the effective surface tension due to the lipid coating decreased the frequency of Pdf sh to frequencies close to resonance. Moreover, P3 oscillations are unexpectedly generated.  (Fig. 12e-f).
For f = 1.8f r (Fig. 17g-h A.J. Sojahrood et al. another SN. The bubble oscillates with P1 for the rest of the studied pressure amplitude range. In this case the dynamic variation of the effective surface tension of the lipid coating enhances the generation of P2 oscillations at very low pressures. The coating lowers the pressure threshold for the Pdf sh ; however, at higher pressures suppresses the nonlinear oscillations. For f = 2f r (Fig. 18a-b), PD is initiated at P a = 17 kPa for both bubbles. A SN bifurcations transition the P2 oscillations to P2 oscillations of higher amplitude at 26 kPa. P4 oscillations are generated and transition back to P2 oscillations through bubbling bifurcation. The P2 oscillations undergo SN bifurcation to P1 oscillations at P a = 268 kPa & 230 kPa respectively for the bubbles with σ 0 = 0.01 N/m & σ 0 = 0.062 N/m. Oscillations of the bubble stay at P1 for the rest of the pressure amplitude ranges that is studied here. Compared to the uncoated bubble in Fig. 12c the lipid coating enhances the P2 oscillations at low acoustic pressures; however, the P2 oscillations of the bubble is suppressed at higher pressures.
When f = 2.2f r (Fig. 18c-d), P2 oscillations are generated through a PD at P a = 17 kPa and then at P a = 20 kPa undergo a SN bifurcation to P2 oscillations of higher amplitude. Oscillations undergo a cascade of At 213 kPa, P2 becomes P3 via a SN which is then followed by a SN from P3 to P2 and reverse PD to P1 for the rest of the studied pressure amplitude range. For the bubble with σ 0 = 0.062 N/m within the pressure amplitude range of 140-376 kPa, there are intermittent transitions between P2 and P3 via SNs. At 401 kPa, P1 oscillations give birth to a P4 oscillations which then transition to P3 via a SN at 411 kPa. P1 oscillations emerge out of the P3 window via a SN at 425 kPa. Compared to the uncoated bubble, lipid coating enhances the P2, P3 and chaotic oscillations at very low acoustic pressures. Moreover, P4 oscillations appear at 3f r . In case of the uncoated bubble and for the same initial conditions however, P4 is expected to appear at frequencies near 4f r [38,40].   19a shows the case of sonication with f = 0.3f r . P1 oscillations undergo a SN at P a = 44 kPa and the bubble oscillations amplitude increases abruptly (This is similar to the dynamics of the bubble sonicated with its Pdf r in Fig. 11c-d). Wall velocities are in phase (blue curve meets the red curve) with the driving acoustic pressure for a range of acoustic excitation pressures and with increasing pressure amplitude the two curves diverge. PD occurs at 371 kPa followed by a cascade of PDs leading to chaos at ≈ 595 kPa. Further pressure amplitude increase results in the intermittent transition between chaos and periodic behavior. This behavior is similar to the dynamics of the uncoated bubble sonicated with its pressure dependent resonance frequency (Pdf r ) in Fig. 11a-d. The presence of the coating thus lowers the pressure threshold for the SN bifurcation. However, the pressure threshold for PD is higher and the bubble oscillation amplitude is generally smaller than the uncoated bubble.
When f = 0.5f r (Fig. 19b) the bubble undergoes two bifurcations that leads to two abrupt increases in the bubble oscillation amplitude. The first is a SN which takes place at 34 kPa transitioning the P1 oscillations to a P1 with higher amplitude. The second one is an inflection point at 460 kPa transitioning the P2 oscillations to P2 oscillations with slightly higher amplitude. Here, the system exhibits dynamics that are similar to two different regimes of the oscillations in the uncoated bubble. The low pressure amplitude transition is similar to the low pressure amplitude transition of the uncoated bubble sonicated with Pdf r (Fig. 11c-d). The second transition that occurs at a higher pressure amplitude resembles the dynamics of the bubble sonicated with its Pdf sh in Fig. 12b. When compared to the uncoated counterpart, for the coated bubble the first transition occurs at a lower pressure amplitude while the second transition occurs at a higher pressure.
When f = 0.6f r (Fig. 19c), we witness the same two pressure thresholds as the previous case. Two SN occur, one at P a = 29 kPa and the second one at 327 kPa. The first SN transition P1 to a P1 oscillation of higher amplitude (Pdf r ) while the second SN transition the P2 oscillations to P2 oscillations of higher amplitude (Pdf sh ). Further pressure amplitude increases leads to P4 with bubbles of P8 inside. Right after the bubble 4 small windows of chaos appear which transition to P4 and then again to chaos.
When f = 0.7f r (Fig. 19d) two SN takes place; the first SN transitions a P1 oscillation to a P1 oscillation of higher amplitude at 25 kPa (Pdf r ) and the second SN transition the P1 oscillation to P2 oscillations of higher amplitude at 277 kPa (Pdf sh ). Pressure amplitude increase leads to P4 through a PD at ≈ 600 kPa and later chaos at 671 kPa.
Looking at Fig. 19a-d, the dynamic variation of the surface tension of the lipid coating significantly decreases the frequencies of pressure dependent resonance (Pdf r ) & and pressure dependent SH resonance frequency (Pdf sh ). As an instance, Pdf sh typically occurs for 1/5f r < f < 2f r for the uncoated bubble ( [58] and Fig. 12a-b) while here Pdf sh occurred at frequencies as low as 0.5f r for the coated bubble with σ 0 = 0.036 N/m. When f = f r (Fig. 19e), an unexpected behavior is observed. P1 oscillations undergo a SN to P3 oscillations at 833 kPa. In case of the uncoated bubble (Fig. 12e-f) or bubbles with pure viscoelastic coating [38], this behavior only occurs for frequencies close to 3f r . Thus the lipid coating here, decreased the P3 resonance frequency by 200 %. The pressure threshold for P3 oscillations, however is higher for the coated bubble when compared to the uncoated counterpart.
When f = 1.2f r (Fig. 19f), nonlinear oscillations are suppressed to only a P1 oscillation for the studied pressure amplitude range.
Compared to the P3 oscillations in case of the uncoated bubble ( Fig. 12e-f), P3 occurs at lower pressure thresholds. For instance at f = 2.2f r P3 is generated at 157 kPa. This is however, higher than the pressure threshold for P3 oscillations in case of the coated bubbles with σ 0 = 0, 0.01, 0.62 & 0.072 N/m.
When f = 2.8f r (Fig. 20e), P3 is generated through a SN at 134 kPa and later transition to P1 via another SN at 269 kPa. P5 oscillations are generated at 373 kPa through a SN and transition to P1 at 433 kPa. P5 oscillations re-appear again for a short pressure window through the same mechanism at 647 kPa.
When f = 3f r (Fig. 20f), P3 oscillations start at 128 kPa and stretch up to 279 kPa with a short window of P1 oscillations within. P5 oscillations are generated at 351 and 592 kPa for two short pressure windows. Compared to the uncoated bubble sonicated with 3f r (Fig. 12f), the pressure threshold of P3 oscillations is lowered by about 276 %.
Coating with σ 0 = 0.036 N/m significantly reduced the frequency

Investigation of the mechanism of the disappearance (standstill) and regeneration of P2
In subSection 3.2 we showed that the enhancement of P2 oscillations at lower pressures can be caused by the asymmetric variations of the effective surface tension due to buckling or rupture. Here, we look into the possible reasons of the disappearance of the P2 oscillations when increasing pressure. Fig. 21a shows the radial oscillations as a function of the driving periods of the coated bubble in Fig. 17d (R 0 = 2μm and σ 0 = 0.062 N/m) at P a = 400 kPa. At this pressure amplitude the P2 oscillation regime disappeared. Radial oscillations are P1, and the red circles return only one value. The corresponding σ(R) curve, depicts a rather symmetrical variations in the buckling and rupture, the bubble spends the same approximate time in the buckled stage as the ruptured stage. For 10 driving periods, the bubble buckles 10 times and ruptures 10 times.
As pressure amplitude increases, P2 is regenerated (Fig. 17d). At 650 kPa the radial oscillations vs period curves have two maxima (Fig. 21)c and the red circles have 2 distinct values. σ(R) as a function of the driving periods (Fig. 21d) exhibits an asymmetrical behavior between the buckled and the ruptured state. The bubble spends a longer time duration at the ruptured stage than the buckled stage. As a result, the bubble buckles 5 times and ruptures 5 times within 10 driving periods. Thus oscillations become P2 again.

Sonication with f < f r
First the findings related to the sonications with frequencies smaller than resonance are presented. 1-SuH regimes are generated at lower excitation thresholds compared to the uncoated bubbles. The bubbles initially at the buckled or the ruptured stages exhibit SuH regime of oscillations at the lowest pressure threshold of 1 kPa. Thus applications of coated bubbles with initial surface tension close to 0 N/m or 0.072 N/m have the potential to increase the contrast in super harmonic imaging (e.g. [95]). Due to the lower threshold of the SuH generation, the amplitude of the generated harmonics in tissue will be smaller. Therefore, the contrast to tissue ratio may be higher.
2-The sudden increase in the bubble oscillation amplitude (SN bifurcation or the inflection point in bifurcation diagrams) occurs at lower excitation amplitudes when compared to the uncoated bubble and coated bubbles with linear viscoelastic behavior [56]. The SN bifurcation is more pronounced in case of the bubbles with σ 0 = 0.01 N/m and σ 0 = 0.062 N/m. The wall velocity stays in phase with the driving sound field for a larger pressure amplitude range. The reason for the lower pressure threshold for the SN and lower frequencies of Pdf r is the fast decrease of the resonance frequency with increasing pressure. Overvelde et al. [68] has experimentally and numerically shown that for coated microbubbles undergoing buckling, the nonlinear resonance behavior is enhanced at pressures as low as 10 kPa. Helfield and Goertz [77] experimentally observed the enhanced nonlinear resonance behavior of the lipid coated microbubbles at pressures of 13-25 kPa. The SN bifurcation can have applications in amplitude modulation techniques [96]. occur at frequencies as low as 0.6f r . In such cases two SN bifurcations are observed with increasing pressure. The first SN occurs at a lower pressure threshold and transfers a P1 oscillation to a P1 oscillation of higher amplitude. The second SN occurs at a higher pressure amplitude and transfers a P2 oscillation to a P2 oscillation of higher amplitude. In [58] we have shown that Pdf sh typically occurs at 1.5f r < f < 2f r and can lead to oversaturation of the 1/2 order SH component of the scattered signal. The enhanced nonlinear resonance behavior of the coating thus shifts the Pdf sh to frequencies below resonance. The occurrence of the two SNs may have potential applications in increasing the contrast in multi-pulse amplitude modulation techniques.

4.1.2.
Case of the coated bubble with σ 0 = 0.036 N/m 1-Compared to the uncoated bubble and coated bubbles with linear viscoelastic behavior (Pdf r is within 0.5f r < f < f r [56]), the frequency of Pdf r is much lower (as low as 0.3f r ).2-Pdf sh can occur even at f = 0.5f r . In case of the uncoated bubble Pdf sh occurs at 1.5f r < f < 2f r [58].3-For 0.5f r ⩽f ≤ 0.7f r and with increasing pressure, two SN occur; the first one transfers a P1 oscillation regime to a higher amplitude P1 and the second one which is at a higher pressure transfers a P1 or a P2 oscillation regime to a higher amplitude P2 regime.

f⩾f r
In this section findings of the sonications with frequencies above resonance are summarized. Such a frequency range is typically used in SH imaging of microbubbles in contrast enhanced ultrasound [20,68,78,96] 2-For the coated bubbles with σ 0 = 0 N/m & σ 0 = 0.072 N/m that are insonated with 1.8f r ⩽f⩽2.2f r , with increasing pressure, P2 oscillations are generated through a PD (at a pressure threshold of 1 kPa), disappear and then are regenerated at a higher pressure amplitude as a higher amplitude P2 or P3. The second P2 is similar to the dynamics of the uncoated bubble undergoing a SN to P2 when f = Pdf sh . The second P3 is similar to the dynamics of the uncoated bubble undergoing a SN to P3 when f ≈ 3f r . In [73], the disappearance of SH oscillations is referred to as an "unexpected standstill" of SHs. This means that, in the case of a bubble able to generate a stable subharmonic oscillation, the subharmonic emission disappears if the acoustic pressure is raised above a second pressure threshold. The subharmonic standstill however, is a reversible [73]; that is, if the acoustic pressure amplitude is decreased again, the bubbles start generating subharmonics one more time [73]. Thus, disappearance is not due to the bubble destruction. Prior works on subharmonics performed on a population of microbubbles did not report this kind of behavior because it was probably 'masked' by the overall response of the several other bubbles within the same sample volume that experience different pressure amplitudes [73]. The standstill of subharmonic emission also was not explained by the numerical studies of the nonlinear models of the bubble dynamics. Here, we show that the disappearance of the SHs is due to the symmetric buckling and rupture of the shell at moderate pressures. At higher pressures, similar to the lower pressures the buckling and rupture of the shell becomes asymmetric. This manifests itself in the re-generation of P2 signals. Above the second pressure threshold, the bubble spends more time in the ruptured stage than the buckling stage. This exposes the bare gas to water for a longer duration and thus can explain the reduced stability of SH oscillations when they were re-generated [73].
In sensitive therapeutic applications like blood-brain barrier opening, the SH components of the scattered pressure by microbubbles are commonly used as a signature for quantifying the nonlinear oscillations of the bubble cloud and treatment efficacy [97,98]. Due to the strong interplay between stable and inertial cavitation regimes, understanding the origin and stability of P2 oscillation regimes is crucial. Thus, the information on the generation, disappearance, amplification and stability of the P2 oscillations that is obtained here, provides a framework for the analysis of the optimization of SH oscillations in applications.
3-For the coated bubbles with σ 0 = 0.01 N/m & σ 0 = 0.062 N/m sonicated with f = 1.5f r , with increasing pressure, P2 oscillations are generated through a PD and then disappear via a SN. Above a second pressure threshold, a P3 oscillation regime occurs via a SN from a P1 regime. This is similar to the dynamics of the uncoated bubble undergoing a SN to P3 when f ≈ 3f r . The pressure threshold for PD is smaller than the uncoated bubbles [40] and the coated bubbles with linear viscoelastic behavior [38]. icated with 2.8f r ⩽f⩽3f r which is followed by the emergence of P3 out of the chaotic window. To our best knowledge, such low pressure thresholds for chaotic oscillations has not been observed in a bubble oscillator. The pressure threshold for P3 is approximately 5 times smaller than the uncoated counterpart.
Here we identified several different types of SN that occur with increasing pressure amplitude in the oscillations of the lipid coated bubbles. This information, can provide the fundamental framework for the optimization of amplitude modulation techniques and SH imaging procedures. Moreover, the enhanced P3 and higher order oscillations may find potential in mixing applications and drug delivery.
In the cases analyzed in this paper, Rmax R0 was higher for the bubbles with a higher σ 0 because of the expansion dominated behavior of the bubble. This can be one of the reasons for the lower pressure threshold of P2 and chaotic oscillations in case of the bubble in the ruptured state.

4.2.2.
Case of the coated bubble with σ 0 = 0.036 N/m 1-For 1.5f r ⩽f⩽3f r with increasing pressure amplitude a P3 occurs via a SN through a P1 oscillation regime. The pressure threshold for P3 is about half of the uncoated counterpart. P3 disappears via a SN. A second or 3rd SN may occur with pressure increase that can lead to the regeneration of P3 or the generation of P5 or P7 oscillations. Due to the wide range of the pressure amplitude and frequency of P3 behavior for the bubbles with σ 0 = 0.036 N/m, engineering of the coatings with such initial surface tensions may find potential in higher order SH imaging with potential higher resolution and contrast. In [59] we have shown that the 2/3 or 3/4 order SHs are stronger than 1/2 order SHs and due to their close proximity to the transducer central frequency they may be detected with superior sensitivity.

Limitations and future work
The goal of this paper was to study the influence of the lipid coating on the nonlinear dynamics of the MBs. Thus for simplicity we only analyzed the radial oscillations of the bubble. Future work can be extended by analyzing the scattered pressure of the bubbles to find the regions of SH power enhancement. Bubble-bubble interaction should also be considered as in applications bubbles exist in poly-disperse clouds. The bifurcation structure of the interacting bubbles has been studied extensively in [99][100][101][102][103][104][105]. These studies have shown that the bubble-bubble interaction significantly influences the dynamics of each bubble. Effects of coupling, bubble size, and spatial arrangement have been studied in [102] and effects of boundary proximity on the dynamics of a bubble cluster is investigated in [103]. We have shown in [105] that the bubble cluster may exhibit collective behavior dominated by the response of the larger bubbles. Future studies need to look into the effect of the interaction of lipid coated MBs and potential collective behavior of the lipid coated bubbles at lower excitation pressure amplitudes.
This study was limited to the case of sonication with a single frequency acoustic excitation. Dynamics of MBs that are sonicated with dual-frequency acoustic excitation have been investigated in several studies [47,49,50,52,51,55,[106][107][108]. Sonication with two frequency forcing can be used to suppress chaos [47,49,50] or enhance the nonlinearity of the system [47,49,50,52,51,55]. Dual frequency forcing may also be used to enhance the bubble expansion to achieve a higher chemical yield [106,107]. The dynamics of lipid coated MBs excited by multi-frequency acoustic excitations can be a subject of future studies. The enhanced non-linearity of such a system may find potential new applications and improvements in contrast enhanced imaging and therapy.
In this study non-spherical oscillations of the MBs were also neglected. Pioneering work of Holt and Crum [109] investigated the subharmonic behavior of larger bubbles (≈ 100 − 200 μm in size) and have experimentally observed the shape oscillations concomitant with subharmonic oscillations. They showed that the appearance of shape oscillations could be phenomenologically mistaken for a simple perioddoubling of the radial mode. Versluis et al. [110] through using high speed optical observations were able to identify shape oscillations of mode n = 2 to 6 in the behavior of single air bubbles with radii between 10 μm and 45 μm. Their study [111] showed that the close to resonance bubbles were found to be most susceptible toward shape instabilities. In A.J. Sojahrood et al. case of coated MBs, non-spherical bubble oscillations were investigated in [111] through high speed optical observations. It was shown that nonspherical bubble oscillations are significantly present in medically relevant ranges of MB radii and applied acoustic pressures. Nonspherical oscillations develop preferentially at the resonance radius and may exist during SH oscillations [111]. Recently Klapcsik and Hegedüs [46] through GPU accelerated large parameter investigations and 2D bifurcation diagrams, have shown that non-spherical oscillations can affect the subharmonic threshold and nonlinear behavior of bubbles. Most recently Guédra et al. [112] through optical observations have shown that subharmonic oscillations can be triggered by energy transfer from surface to volume oscillations and thus can change the pressure threshold for SH emissions. Thus, for a more accurate modeling of the nonlinear behavior of lipid coated MBs, deeper theoretical modeling of MB coating, accounting for membrane shear and bending is required [111].
One of the limitations of the current study is neglecting the effects of lipid shedding and mass transfer. These effects may become important at higher pressures and they should be considered for accurate modeling. The lipid coating undergoes buckling and rupture when the bubble oscillates. The lipid coating reseals quickly when the bubble contracts [113,114]. For a long enough pulse and depending on the applied pressure amplitude, the coating may shed some lipids while it undergoes buckling, rupture and reseal [115,116]. This leads to mass transfer and shrinkage of the bubble which reaches a stable size after a few cycles [114,115]. The incorporation of these effects can be the subject of future studies.
In this work we studied the bifurcations structure of the bubble oscillator using the standard methods of bifurcation analysis. The buckling and rupture of the shell, however, makes the lipid coated bubble a non-smooth system [117][118][119]. Similar features of the bifurcation structure of the lipid coated bubble may be seen in the behavior of the pressure relief valve model which is a non-smooth system. Future studies that focus on the nonlinear properties of the lipid coated bubble, can reveal more detailed information about the system behavior using the tools of non-smooth dynamics [118,119].

Conclusion
In this work, the bifurcation structure of the lipid coated bubbles undergoing buckling and rupture was studied extensively. Our results further confirmed that the rapid variation of the effective surface tension and buckling and rupture of the coating enhances the generation of nonlinear behavior including higher order SHs, SuHs and chaos. We showed for the first time that P2 and P3 can occur at pressures as low as 1 kPa (≈ 1% of the ambient pressure). Existence of chaos was confirmed at pressures as low as 10 kPa. The closer the initial surface tension of the bubble to the buckling stage or the ruptured stage, the lower the pressure threshold for the nonlinear behavior. We showed that rapid variations of the surface tension on the bubble may not be enough for enhanced non-linearity. In case of asymmetrical variations of the surface tension between buckling and rupture, nonlinear behavior is enhanced. However, symmetrical behavior of the effective surface tension may suppress the non-linearity.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The Marmottant model [62] in the absence of the coating reduces to: In this equation some of the terms of the Keller-Miksis (K-M) equation [81] are dropped. This is because K-M model is no longer valid when |Ṙ| c ≈ 1. [59]. The Marmottant model is in a popular form that is developed for sonoluminescence and ultrasound contrast agent (UCA) modeling [120,121]. The results of the Marmottant model were compared to the K-M uncoated bubble model to investigate the shell effects. However, as the forms of the K-M model and the Marmottant model in the absence of coating are not the same, the comparison may not seem justified. In this section we show that the K-M model in Eq. (3) and Eq. (A1) indeed predict the same overall behavior over the exposure parameter ranges that were studied in our paper. Thus the comparison between the models is justified. Figs. A.figurekk and A.figurekk show that the predictions of the bifurcation structure of R/R 0 vs pressure of the Marmottant model with no shell (Eq. (A.1)) are in good agreement with the Keller-Miksis model (Eq. (3)) within the exposure parameter ranges that were studied in this work. It should be noted that the comparison between the coated bubble model and the free bubble model are qualitative. In other words, the uncoated bubble behavior was used as a reference to reveal the nonlinear effects induced by the presence of the shell. However, we see here that despite the difference in some terms of the two equations, the quantitative behavior of (R/R 0 ) of the uncoated bubble as predicted by Eq. (A.1) and Eq. (3) are in agreement. The radial oscillation amplitude and the pressure threshold for the onset of nonlinear behavior (e.g. PD, chaos) are in agreement. However, although both models predict the same pressures for the onset and disappearance of chaos, when the oscillations are chaotic, the radial oscillation amplitudes are not always equal. For better visualizing the predictions of the two models at lower pressure amplitudes, the figure insets concentrate on the lower acoustic pressure regions P A ⩽250 kPa.  Appendix B. Influence of χ and σ rupture on the buckling and rupture radius Analysis in this paper showed that the nonlinear behavior of the lipid coated MBs is enhanced at very low pressures. This occurs due to the rapid variation of the effective surface tension when the bubble radius exceeds above the rupture radius or contracts below the buckling radius. Fig. B. figurekk.a shows that for a constant σ 0 and σ rupture as χ increases, R r decreases and R b increases (R r and R b converge to R 0 ); thus, a lipid MB with a higher χ which is not initially at buckled or ruptured state exhibits nonlinear behavior at a lower pressure threshold.
MBs with resistant shells that can withstand higher tensions before rupture (e.g. with σ rupture > 0.1 N/m) has been observed by Marmottant et al. [62]. Fig. B.figurekk.b shows that R r increases with increasing σ rupture . For a lipid bubble with initial surface tension of more than 0.036 N/m, the pressure threshold for enhanced non-linearity will increase as σ rupture increases. This is because R r increases linearly with σ rupture , thus higher pressure is needed to enhance the nonlinear oscillations. A detailed understanding of the influence of χ and σ rupture over a wide pressure and frequency range can be a subject of future investigations.  figurekk shows the pressure threshold for the onset of nonlinear oscillations for a C3F8 lipid coated MB with χ = 3.5 N/m and σ rupture = 0.072 N/m as a function of the shell viscosity k s and σ 0 . The frequency is 2f r . The pressure threshold is extracted from the radius vs pressure amplitude bifurcation diagrams for P2 or P3 oscillations. Increasing k s does not change the pressure threshold (P a = 1 kPa) for nonlinear oscillations for the MBs initially at buckled or ruptured stage. For MBs which are not initially at buckled or ruptured stage, the pressure threshold increases with increasing k s .
As an example for a MB with σ 0 = 0.062 N/m, the pressure threshold increases from 16 kPa to 71 kPa when k s is increased from 1 * 10 − 9 to 8 * 10 − 9 kg/s. Bubbles which are initially at a buckled or ruptured state start nonlinear oscillations at the lowest pressure threshold (P A = 1 kPa), regardless of the value of k s .