{'fraction_non_alphanumeric': 0.029000423805775867, 'fraction_numerical': 0.01749712417509233, 'mean_word_length': 4.599322033898305, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A comprehensive bifurcation method to analyze the super-harmonic and ultra-harmonic behavior of the acoustically excited bubble oscillator

Acoustically excited bubbles are involved in a wide range of phenomena and applications ranging from oceanography to sonoluminescence; they have applications in chemistry, medical imaging, and therapeutic ultrasound. The complexity of bubble dynamics and the limited understanding of their behavior restricts the exploration of their full potential. The bubble oscillator is a highly nonlinear system, which makes it difficult to generate a comprehensive understanding of its oscillatory behavior. One method used to investigate such complex dynamical systems is the bifurcation analysis. Numerous investigations have employed the method of bifurcation diagrams to study the effect of different control parameters on the bubble behavior. These studies, however, focused mainly on investigating the subharmonic (SH) and chaotic oscillations of the bubbles. Super-harmonic (SuH) and ultra-harmonic (UH) bubble oscillations remain under-investigated. One reason is that the conventional method used for generating bifurcation diagrams cannot reliably identify features that are responsible for the identification of SuH and UH oscillations. Additionally, the conventional method cannot distinguish between the UHs and SHs. We introduce a simple procedure for the generation of bifurcation diagrams to address this shortcoming. This method selects the maxima of the bubble oscillatory response and plots them alongside the traditional bifurcation points for the corresponding control parameter. Through applying this method, the oscillatory behavior of the bubble oscillator is analyzed, and stable SuH and UH bubble oscillations are investigated. Based on this new analysis, the conditions for the generation and amplification of UH and SuH regimes are discussed.

In this work, we introduce a more comprehensive and simple method to study the SuH and UH bubble oscillations. The bifurcation structure of the bubble oscillator is constructed by plotting the maxima of the stable oscillations of the bubble; which is plotted alongside the conventional bifurcation diagram. Using this method we were able to straightforwardly identify the SuH and UH oscillations and explore the conditions that are required to generate and amplify the SuH and UH oscillations.
This method establishes a framework that provides a more comprehensive understanding of the rich nonlinear behavior of bubbles; consequently, it may help in optimizing current applications and/or can be used to discover new nonlinear bubble behaviors that may result in new applications.
II. Methods

i. The Bubble model
The radial oscillations of the bubbles are numerically simulated by solving the well known Keller-Miksis equation [42]: In this equation, R is radius at time t, R 0 is the initial bubble radius, ρis the liquid density, c is the sound speed, p 0 is the atmospheric pressure, σ is the surface tension, µ is the liquid viscosity, P A and f are the amplitude and frequency of the applied acoustic pressure.
Oscillations of a bubble generates a backscattered pressure (P Sc ) which can be calculated by [43]: where d is the distance from the center of the bubble (and for simplicity is considered as 1m in this paper) [43]. Equation 1 is solved using the 4 th order Runge-Kutta technique; the control parameters of interest are R 0 , f and P A . The resulting radial bubble oscillations are visualized using the bifurcations diagrams. Bifurcation diagrams of the normalized bubble oscillations R R 0 are presented as a function of driving pressure. Detailed analysis is presented at select control parameters using a) the radius versus time curves, b) phase portrait analysis and c) the frequency spectrum of the backscattered pressure.
ii. Bifurcation diagrams For highly nonlinear systems like bubble oscillators, small changes in the initial conditions of the system or control parameters can result in large changes in the behavior of the system. However, many studies focus on analyzing the dynamics of the bubbles over discrete values of the control parameters. Due to the complexity and sensitivity of the bubble dynamics to the exposure parameters and initial conditions, the limited values of control parameters used in these studies do not provide a comprehensive understanding of the bubble dynamics. In addition, due to the discrete nature of the investigated parameters, many exposure parameter combinations have not been investigated. Bifurcation diagrams are valuable tools to analyze the dynamics of nonlinear systems where the qualitative and quantitative changes of the dynamics of the system can be investigated effectively over a wide range of the control parameters.

ii.1 Conventional bifurcation analysis
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 from a specific point in each driving period. The approach can be summarized in: Where Q denotes the points in the bifurcation diagram, R andṘ are the time dependentradius 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. The results are plotted for n = 100 − 150 to ensure a steady state solution has been reached.

ii.2 Method of peaks
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 R 0 ) 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.

III. Results
We have considered a bubble with initial diameter of 4 microns. The linear resonance frequency (f r ) of the bubble is~2.04 MHz. We have studied the bifurcation structure of the bubble in the form of R R 0 as a function of the driving acoustic pressure for (0.2f r <f<3f r ) and (1kPa<P A <3 MPa). Results are shown in Figure 1. for f=0.7, 1.2 and 2.6 MHz respectively. To focus on more practical and stable oscillation regimes, we have omitted the parameter ranges that result in chaotic oscillations or bubble destruction ( R R 0 > 2 as discussed in detail in [16]). The bifurcation structures that are produced using the conventional method are presented in blue, and the ones produced by the method of peaks are shown in red. Fig. 1a shows the response of the bubble when f=0.7 MHz. The conventional method reveals a period 1 solution for P A <118 kPa and detects the generation of period 2 solution for P A >120 kPa. On the other hand, the peaks method reveals the generation of two maxima at 24 kPa<P A <56 kPa and three maxima for 56 kPa<P A <118 kPa. The three maxima undergo period doublings (PDs) resulting in a period 6 solution for PA > 118 kPa. Examination of the frequency content of the backscattered pressure reveals the underlying phenomenon which results in the discrepancy between the two methods. Figure 1d shows the amplitude of the harmonics and UHs of the backscattered signal. The occurrence of the maxima correlates with the resonance of the harmonic contents of the signal. After a pressure threshold (~25 kPa), the 3 rd SuH of the backscattered signal becomes stronger than the fundamental and other SuH harmonics, showing a 3 rd SuH resonance. The 3 rd SuH saturates for P A > 56 kPa concomitant with the occurrence of 3 maxima in the peaks methods. The SH and UH contents of the backscattered signal are shown in Fig. 1g. The simultaneous appearance of period doublings (PDs) in the blue curve and multiple PDs in the red curve are coincident with a sharp increase in the SH and UH content of the backscattered signal (Fig. 1g, arrow); the backscatter at the 7 2 (purple) and 9 2 UHs (green) are the strongest whit the 1 2 SH (blue) the weakest component. Fig. 1b shows that when sonication frequency is 1.2 MHz the conventional method depicts a similar behavior to f=0.7 MHz (fig. 1a); a linear response is observed for P A < 110kPa and radial oscillations undergo a PD for P A > 110 kPa. The method of peaks reveals a solution with one maximum for PA < 51 kPa which is similar to the conventional method; above this pressure, however, 2 maxima occur in the bifurcation diagram up until P A = 110kPa which are not detected in the conventional method. For PA > 110 kPa, the oscillations undergo two concomitant PDs resulting in a solution with 4 maxima. The conventional method predicted the same behavior for the two frequencies (0.7 and 1.2 MHz), however, the method of peaks revealed more intricate details of the bubble dynamics. Fig. 1e shows that the second harmonic of the backscattered signal has the strongest amplitude and saturates concomitant with the generation of the initial two maxima at~52 KPa. Fig. 1h illustrates a sharp increase in the amplitude of SH and UHs concomitant with the generation of PDs (arrow) in the blue and red curves as shown in Fig 1b. UH components of the signal are stronger than the SHs (~20 dB) with 5 2 and 3 2 UHs being the strongest while the 1 2 SH is the weakest component. Fig. 1c shows that when f= 2.6 MHz the conventional method (blue) and method of peaks (red) depict a similar behavior to f=0.7 MHz and f=1.2 MHz (Fig. 1a, 1b). The oscillations are of period 1 in both graphs until P A =243 kPa; above this pressure, PD occurs in both methods. The fundamental frequency is the strongest frequency component in the backscattered signal as is shown in Fig. 1f. Figure 1i indicates that the SH and UH components of the signal increase rapidly at pressures at which the PD occurs in Fig. 1c ; the SH component is stronger than all the UHs.
To gain a better insight of the oscillation characteristics that the conventional bifurcation analysis method was unable to reveal, a) radial oscillations vs. driving acoustic periods, b) phase portraits and c) backscattered frequency spectra are examined in detail. Exposure parameters were chosen for which both methods give similar predictions and exposure parameters for which the predictions diverged. Fig. 2a depicts the radial oscillations of the bubble for P A =125 kPa and f=2.6MHz. Inspection of Fig. 1c indicates that a period 1 (P1) oscillation regime is expected. The yellow stars represent the amplitude of radial oscillations after every period, and the red circle illustrates the maxima of the curve. There exists only one value for all red circles and yellow stars; therefore, the signal is a P1 signal with one maximum. The corresponding phase portrait in Fig. 2b is a semi circular orbit, and the fundamental component of the P Sc is stronger than the SuHs. Fig. 2d shows the R-T curve that corresponds to P A =78 kPa and f=1.2 MHz. The signal has two maxima, while the amplitude of the signal at each driving period remains the same. In other words, this suggests a P1 signal with two maxima. Fig. 2e shows that the phase portrait of the bubble undergoes an internal bend; depicting a similar behavior when SHs are present in the curve. However, the absence of SHs are evident in the frequency spectra of the corresponding P Sc (Fig. 2f) while the 2 nd harmonic has the highest value (2 nd harmonic resonance). In this case examination of the maxima provides more complete information about the oscillation characteristics indicating a 2 nd SuH resonance. Figure 3 compares the characteristics of two P2 oscillations; one has two maxima while the other has 4 maxima. Fig 3a shows the R-T curve of the bubble for f=2.6 MHz and P A =275kPa. The signal has two maxima (two red dots), and the signal repeats its pattern once every two acoustic driving periods (two yellow stars). The phase portrait (Fig 3b) consists of two circular orbits with one creating another and enclosing it within itself. Fig. 3c depicts the existence of 1 2 order  SH which is stronger than the UH components. Fig. 3d shows the R-T curve of the bubble for f=1.2 MHz and P A =145kPa; the signal is of P2 but with 4 distinct maxima. The radial oscillations repeat their pattern once every two acoustic periods (two yellow stars), and each pattern has 4 maxima (4 red dots). The phase portrait has two circular orbits similar to Fig 3b; however, each of these circular orbits underwent an internal bend. The frequency content of the P Sc in Fig. 3f has 1 2 order SHs as well as UHs; the 2 nd order SuH is the strongest signal, and 5 2 and 7 2 UHs are stronger than 1 2 order SH and other UHs.

IV. Summary and Conclusion
We have shown that the conventional method of generating bifurcation diagrams cannot reveal the hidden details of the oscillations that typically result in UH and SuH resonance. We have introduced a simple alternative method to generate the bifurcation diagrams; the method extracts the peaks of the oscillations and plots it as a function of the given control parameter. When this method is applied alongside the conventional method one can reveal hidden details of the bubble oscillations and identify the parameter ranges where SuH, UH or SH oscillations occur. We can briefly categorize the scenarios shown in this paper as follows: 1. The conventional method depicts a P1 oscillation regime, and the maxima method only reveals one maximum. In this case, the oscillation has a P1 resonance and the fundamental frequency component is the strongest in the backscattered signal. 1 2 order SH and UHs are generated concomitant with PD in both graphs (constructed by conventional method and method of maxima) and 1 2 order SHs, or 3 2 UHs are stronger than other UHs.
2. The conventional method depicts a P1 oscillation regime, but maxima method reveals n=2,3,. . . maxima. In this case, the n-th order SuH frequency component is the strongest in the backscattered signal. Generation of PD in the conventional method is concomitant with the generation of n-PDs in the diagram constructed by the peaks method; this correlates with an UH resonance behavior of (2n−1) 2 or (2n+1) 2 ; in other words, these UHs are stronger than the 1 2 order SHs and other UHs.
3. The conventional method depicts a P2 oscillation regime, but the maxima method reveals only 1 maximum; in this case, we have a P2 resonance; 1 2 order SH frequency component is the stronger than UHs.
4. UH and SHs only exist when the conventional method predicts a P2 oscillations; however, the method of maxima needs to be applied alongside of traditional method to determine whether SH or UH resonance are present, as well as the order of UHs.
The nonlinear behavior of the lipid shell enhances the generation of SHs ( 1 2 , 1 3 , 2 3 , 1 4 , 3 4 . . . ) at very low acoustic pressures and over a more extensive frequency range [45,46]. The behavior of lipid shell MBs are more complex due to the nonlinear response of the encapsulating shell. Implementation of the proposed method in this paper can shine a brighter light on the behavior of lipid coated MBs. Consequently, this approach can be used to optimize the wide range of applications that employ coated MBs. For example, the techniques presented can be used to optimize contrast enhanced techniques employing super harmonics [46,47] or ultraharmonics [48].
Detailed studies on the effect of initial conditions (ICs) (R(0) andṘ(0)) on the dynamical evolution of the bubble oscillations [28][29][30][31][32][33][34] have resulted in the discovery of new nonlinear features [29][30][31][32][33][34]. For example, depending on the ICs, the bubble has shown to exhibit period 1 (P1), P2 or P3 oscillations [31,33]. Application of the method proposed in this paper can help to better understand and categorize these nonlinear features, especially in the SuH and UH oscillation regimes. Results using this approach may be used to optimize applications by sending the proper pre-conditioning pulses to manipulate the ICs of the bubble to achieve the desirable behavior.