On the threshold of 1/2 order subharmonic emissions in the oscillations of ultrasonically excited bubbles
A.J. Sojahrood, N.R. Shirazi, H. Haghi, R. Karshafian, M.C. Kolios
OOn the threshold of 1/2 order subharmonicemissions in the oscillations of ultrasonicallyexcited bubbles
A.J. Sojahrood , a , b 1 , H. Haghi , a , b N.R. ShiraziR. Karshafian , a , b M.C. Kolios , a , b a Department of Physics, Ryerson University, Toronto, Canada b Institute for Biomedical Engineering, Science and Technology (IBEST) apartnership between Ryerson University and St. Michael’s Hospital, Toronto,Ontario, Canada
Abstract
The pressure threshold for 1/2 order subharmonic (SH) emissions and period dou-bling during the oscillations of ultrasonically excited bubbles is thought to be mini-mum when the bubble is sonicated with twice its resonance frequency ( f r ). Thisestimate is based on studies that simplified or neglected the effects of thermaldamping. In this work, the nonlinear dynamics of ultrasonically excited bubblesis investigated accounting for the thermal dissipation. Results are visualized usingbifurcation diagrams as a function of pressure. Here we show that, and depending onthe gas, the pressure threshold for 1/2 order SHs can be minimum at a frequencybetween 0 . f r ≤ f ≤ . f r . In this frequency range, the generation of 1/2 orderSHs are due to the occurrence of 5/2 order ultra-harmonic resonance. The stabilityof such oscillations are size dependent. For an air bubble immersed in water, onlybubbles bigger than 1 µm in diameter are able to emit non-destructive SHs in thesefrequency ranges. An ultrasonically excited bubble is a nonlinear oscillator that is present in awide range of applications and phenomena from underwater acoustics to ma-terial science and medicine [1]. Starting with the pioneering work of Parlitzet al [2], the bifurcation structure of the bubble oscillators exhibiting perioddoubling and chaos has been the topic of many recent investigations [3,4,5,6,7]. Email: [email protected]
Preprint submitted to Elsevier 17 September 2020 a r X i v : . [ phy s i c s . a pp - ph ] S e p eriod doubling (PD) results in the emission of 1/2 order SHs and (3/2,5/2,etc) UHs by bubbles. These nonlinear oscillators have many applications indiagnostic and therapeutic ultrasound [8,9,10,11,12,13,14] including the SHand UH emissions by bubbles is in monitoring the blood brain barrier openingby ultrasound [15,16].Cavitation and bubble oscillations are generally classified into two types, sta-ble cavitation, which results in emissions at subharmonics of the main ex-citation frequency [16,17] and inertial cavitation, which is characterized bybroadband noise emissions [17]. Understanding the generation and amplifica-tion of SHs as well as differentiating between stable and inertial regimes arecritical for several applications [17,18,19,20,21,22,23]. This is because gentlebubble oscillations may result in more control over the bubble dynamics andless unwanted tissue damage. In non-medical applications like material clean-ing, stable non-destructive bubble oscillations are an additional mechanismof surface cleaning thus eliminating erosive bubble collapse [24,25]. Potentialunderlying mechanisms of cleaning has been investigated in [25] and one of thecleaning mechanisms is through steadily fatiguing the surface contaminationby the non-destructive stable micro-streaming around bubbles.Among many causes, SH emissions can be generated due to nonlinear volumet-ric oscillations of the bubbles [26,27,28,29], from the onset of surface modes[30,31] or from energy transfer from surface to volume oscillations [32].Many studies have investigated the pressure threshold of SH emissions by ul-trasonically excited bubbles [29,33,34,35,36,37,38,39]. These studies indicatethat, the pressure threshold for SH emissions is the lowest when the bub-ble is sonicated with a frequency that is twice its resonance frequency ( f r )[29,33,34,35,36,37,38,39]. More detailed numerical investigations in [37,38] re-vealed that the pressure threshold for SH oscillations shift to frequencies nearresonance for bubbles smaller than 0 . µm in diameter. The mechanism behindthe shift in the minimum pressure threshold is attributed to the increaseddamping due to viscosity on smaller bubbles. Katiyar and Sarkar [38,39] con-clude that the increased dissipation due to viscosity damps the bubble re-sponse more at twice the resonance than at resonance, leading to a shift ofthe minimum from twice the resonance frequency toward the resonance fre-quency. The addition of shell has been shown analytically and numerically in[38,39] to decrease the pressure threshold for SH generation and possibly shiftthe frequency of minimum pressure threshold towards resonance. The impactof bubble-bubble interaction on the threshold of SH emissions has been in-vestigated by Guerda et al [40]. They show that, for a given bubble radius,increasing the concentration of bubbles shifts the subharmonic resonance fre-quency towards lower values and reduces the minimum pressure threshold forSH generation.The previous studies either simplified the thermal damping using the linearmodels [41,42] or ignored thermal effects. To provide a more accurate un-derstanding of the effect of dissipation on the SH emissions by bubbles, wehave considered the effect of thermal damping using the full set of ordinary2 hermal parameters of the gases at 1 atmGas type L( WmK ) c p ( kJkgK ) c v ( kJkgK ) ρ g ( kgm )Air [47] 0.01165+5 . × × T differential equations (ODEs) that describe the thermal effects [43,44]. Usingthe more accurate description of thermal effects, the bifurcation structure ofthe radial oscillations of the bubbles is studied as a function of pressure fordriving frequencies between 0 . f r < f < . f r . We used a more comprehen-sive method for bifurcation analysis introduced in [45]. We show that contraryto the common belief and depending on the gas, insonation with frequenciesbetween 0 . f r < f < . f r leads to the minimum pressure threshold for 1/2order SH emissions. A more careful analysis of the oscillations reveals that 3/2or 5/2 ultraharmonic (UH) resonance is the mechanism behind the generationof SHs in this frequency range. The dynamics of the uncoated bubble including the compressibility effectsto the first order of Mach number can be modeled using Keller-Miksis (KM)equation[46]: ρ [(1 − ˙ Rc ) R ¨ R + 32 ˙ R (1 − ˙ R c )] = (1+ ˙ Rc + Rc ddt )( P g − µ L ˙ RR − σR − P − P A sin (2 πf t ))(1)In this equation, R is radius at time t, R is the initial bubble radius, ˙ R is thewall velocity of the bubble, ¨ R is the wall acceleration, ρ is the liquid density(998 kgm ), c is the sound speed (1481 m/s), P g is the gas pressure.In the absence of thermal effects P g is given by P g = ( P + 2 σR )( R R ) γ , σ isthe surface tension (0.0725 Nm ), µ L is the liquid viscosity (0.001 Pa.s), and P A and f are the amplitude and frequency of the applied acoustic pressure. Thevalues in the parentheses are for pure water at 293 K. In this paper the gasinside the uncoated bubble is either air or C3F8, and water is the host media.Thermal properties for the gases are given in Table 1.3 .2 Thermal effects2.2.1 ODE thermal model If thermal effects are considered, P g is given by Eq. 2 [43,44]: P g = N g KT πR ( t ) − N g B (2)Where N g is the total number of the gas molecules, K is the Boltzman constantand B is the molecular co-volume of the gas inside the bubble. The averagetemperature inside the gas can be calculated using Eq. 3 [43,44]:˙ T = 4 πR ( t ) C v (cid:32) L ( T − T ) L th − ˙ RP g (cid:33) (3)where C v is the heat capacity at constant volume, T =293 K is the initial gastemperature, L th is the thickness of the thermal boundary layer. L th is givenby L th = min ( (cid:114) DR ( t ) | ˙ R ( t ) | , R ( t ) π ) where D is the thermal diffusivity of the gas whichcan be calculated using a = Lc p ρ g where L is the gas thermal conductivity and c p is specific heat capacity at constant pressure and ρ g is the gas density.Predictions of the ODE thermal model have been shown to be in good agree-ment with predictions of the models that incorporate the thermal effects usingthe PDEs [45] that incorporate the temperature gradients within the bubble.To calculate the radial oscillations of the uncoated bubble while including thethermal effects, Eqs. 1 is coupled with Eq. 2 and 3 and is then solved usingthe ode45 solver of Matlab. The linear thermal model [41,42] is a popular model that has been widelyused in studies related to oceanography [49,50] and the modeling and charac-terization of coated bubble oscillations [51,52,53,54,55]. In this model thermaldamping is approximated through linearization by adding an artificial viscos-ity term to the liquid viscosity. Furthermore, a variable isoentropic index isused instead of the polytropic exponent of the gas.In this model P g is given by: P g = P g (cid:18) R R (cid:19) k (4)Where the polytropic exponent γ is replaced by isoentropic indice ( k ): k = 13 (cid:60) ( φ ) (5)4here (cid:60) denotes the real part.The liquid viscosity is artificially increased by adding a thermal viscosity ( µ th )to the liquid viscosity. The thermal viscosity ( µ th ) is given by: µ th = P g (cid:61) ( φ ) ω (6)In the above equations (cid:61) denotes the imaginary part and the complex term φ is calculated from φ = 3 γ − γ − iχ (cid:20)(cid:16) iχ (cid:17) coth (cid:16) iχ (cid:17) − (cid:21) (7)where γ is the polytropic exponent and χ = DωR represents the thermal dif-fusion length. D is the thermal diffusivity of the gas. D = Lc p ρ g where c p , ρ g ,and L are the specific heat capacity at constant pressure, density and thermalconductivity of the gas inside the bubble.To calculate the radial oscillations of the bubble while including the linearthermal effects, Eq. 1 is coupled with Eq. 4 and liquid viscosity is increasedby µ th (Eq. 6). Bifurcation diagrams are valuable tools to analyze the dynamics of nonlinearsystems where the qualitative and quantitative changes of the dynamics ofthe system can be investigated effectively over a wide range of the control pa-rameters. In this paper, we employ a more comprehensive bifurcation analysismethod introduced in [45].
When dealing with systems responding to a driving force, one option is tosample the R(t) curves using a specific point for each driving period. Theapproach can be summarized in: P ≡ ( R (Θ)) { ( R ( t ) , ˙ R ( t )) : Θ = nf } (8)where P denotes the points in the bifurcation diagram, R and ˙ R are the timedependent radius and wall velocity of the bubble at a given set of controlparameters of ( R , P A and f ), Θ is given by nf and n=1,2,....440. Points onthe bifurcation diagram are constructed by plotting the solution of R ( t ) attime points that are multiples of the driving acoustic period. The results are5lotted for n = 400 −
440 to ensure a steady state solution has been reachedfor all bubbles.
Another way of constructing bifurcation points is by setting one of the phasespace coordinates to zero: Q ≡ max ( R ) { ( R, ˙ R ) : ˙ R = 0 } (9)In this method, the steady state solution of the radial oscillations for eachcontrol parameter is considered. The maxima of the radial peaks ( ˙ R = 0 and¨ R >
0) are identified (determined within 400-440 cycles of the stable oscilla-tions) and are plotted versus the given control parameter in the bifurcationdiagrams.The bifurcation diagrams of the normalized bubble oscillations (
R/R ) arecalculated using both methods a) and b). When the two results are plottedalongside each other, it is easier to uncover more important details about theSuH and UH oscillations, as well as the SH and chaotic oscillations. Figure 1 shows the bifurcation structure of an air bubble with R = 10 µm as a function of pressure. The frequency of sonication are in multiples of res-onance frequency ( f r ). For each model the resonance frequency is calculatednumerically by solving the corresponding equations at P A = 1 kP a and findingthe frequency of the maximum radial expansion. The first column in Fig. 1 iscalculated by solving Eq. 1 and neglecting the thermal effects (non-thermalmodel). The results in the second column take into account the linear thermaldissipation by coupling Eq. 1 with Eqs. 4-7 (linear thermal model) and thethird column considers the full thermal effects through coupling Eq.1 withEqs. 2 and 3 (full thermal model). The red graph is generated through thePoincare method (2.3.1) and the blue graph is generated via the method ofmaxima (2.3.2).Fig. 1 shows that in the case of the non-thermal model (NTM) and linearthermal model (LTM), the minimum pressure threshold for PD occurs when f = 2 f r (Figs. 1a and 1b). However, when full thermal effects are consideredthe minimum pressure threshold is when f = 0 . f r (Fig. 1i).When f = 2 f r (Figs. 1a-c) the bubble undergoes period doubling (PD) at6a) (b) (c) Pressure (Pa) R / R R0=10 m,f=1*fr (linear thermal model)
PD, P A = 185kPa Pressure (Pa) R / R R0=10 m,f=1*fr (ODE thermal model)
PD, P A = 171 kPa (d) (e) (f) Pressure (Pa) R / R R0=10 m,f=0.6*fr (linear thermal model) (g) (h) (i) Fig. 1. Bifurcation structure of the
R/R of an air bubble with R = 10 µm as afunction of pressure (The red curves are generated through the Poincar method andthe blue curves are generated via the method of maxima).Top row: f = 2 f r , middlerow: f = f r and bottom row: f = 0 . f r . P A = 17 kP a , 48 kP a and 96 kP a predicted by NTM, LTM and FTM respec-tively. Period doubling is highlighted in a subset in each graph where thered graph bifurcates in to two branches. The FTM model predicts the high-est pressure threshold for the generation of PD and SH emissions. We haverecently shown [56,57] that linear models for dissipation effects in bubble oscil-lations may loose accuracy at acoustic excitation pressures as low as 20 kPa .Predictions of the the LTM and the FTM deviate as pressure increases andthe FTM should be used to model bubble oscillations [57,58]. This is becausethe LTM is derived assuming very small (1-2 percent expansion ratio), andsymmetric bubble oscillation amplitude. However, as pressure increases, theradial oscillations grow above the size limit where linear assumptions are validand the bubble expansion and compression phases are no longer symmetric.Higher internal temperatures can be created during significant compression orcollapse (larger R max /R min ) and the average surface area available for thermalconduction becomes smaller or larger depending on the regime of oscillations[56,57]. This causes the predictions of the LTM and FTM to diverge. When7a) (b) (c)(d) (e) (f)(g) (h) (i) Fig. 2. Bifurcation structure of the
R/R of a C3F8 bubble with R = 10 µm asa function of pressure (The red curves are generated through the Poincar methodand the blue curves are generated via the method of maxima). Top row: f = 2 f r ,middle row: f = f r and bottom row: f = 0 . f r . f = 2 f r the LTM model under-predicts the PD threshold. This is becauseit treats the behavior of the thermal dissipation similar to viscous dissipa-tion through adding a thermal viscosity ( µ th ) term (Eq. 6). Thus, thermaldissipation behaves similarly to liquid viscous dissipation ( ∼ µ th ˙ RR ) which isproportional to the wall velocity amplitude. However, when f = 2 f r , for bub-bles with R > . µm , PD occurs at very gentle wall velocities ( ≈ < f = 2 f r , the bubble average surface area availablefor thermal conduction increases, and subsequently thermal dissipation (Td)increases. Thus, as shown in [57], compared to radiation damping and liquidviscosity damping, Td exhibits the largest increase when PD occurs at f = 2 f r (e.g. Fig. B.1.g in [56]). When PD occurs at f = 2 f r , the wall velocity is small,and the bubble rebounds slowly after collapse ([6]). After the occurence of thePD, the radial oscillations have one maximum (the bubble only collapses onceevery two acoustic cycles). Even when the second maxima reappears at higherpressures, its amplitude is very close to that at the first maxima. Due to the8arge surface area in each period and slow rebound after every collapse, theaverage surface area for thermal conduction is larger, and thus the Td as pre-dicted by the FTM model is higher than the LTM. Thus, in case of the bubblewith R = 10 µm (Fig. 1c) due to the increased Td, FTM predicts a higherthreshold for PD.When f = f r , PD occurs at 174 kPa, 185 kPa and 171 kPa as predictedby NTM, LTM and FTM models respectively. In this case PD occurs at ahigher pressure in the case of the LTM model. This is because when f = f r due to the very large amplitude bubble velocities at PD (for R > . µm , | ˙ R | > m/s ), thermal dissipation is over estimated by the LTM model as Td ∼ µ th ˙ RR . However, in reality the average surface area for temperature escapebecomes limited. When f = f r and at PD, the oscillations have two distinctmaxima. The bubble collapses strongly (leading to increased temperature),and it rebounds back to a radius that is smaller than the previous maxima.Thus, the bubble surface are becomes smaller. Moreover, the bubble reboundsquicker, thus the time duration available for the temperature escape will belimited. These lead to reduction of the average surface area for thermal conduc-tion and thus Td decreases. FTM model incorporate these effects. Therefore,it predicts a lower pressure threshold for PD. The FTM model also predictsa slightly lower pressure threshold (171 kPa in Fig. 1f) for the generation ofPD tcompared to the NTM model (174 kPa in Fig. 1d). Our detailed analysisof the pressure dependent mechanisms in [58] reveals that when f = f r , theradiation damping (Rd), thermal damping (Td) and damping due to liquidviscosity ( Ld ) increase with pressure increase with the order of strength of T d > Ld > Rd . However, Rd grows faster than the other dissipation factorswith increasing pressure. When PD occurs, due to the decrease in averagesurface area for thermal conduction, Td decreases and Rd becomes the majordissipation mechanism (
Rd > T d > Ld ). In case of the NTM model and dueto the absence of Td, the bubble attains higher wall velocities and accelera-tion and due to the dominant role of Rd at PD ( Rd ∼ ( R ¨ R + 2 ˙ R )) the totaldissipation increases slightly which lead to a 3 kPa higher pressure for thegeneration of PD.When f = 0 . f r , the P1 oscillation amplitude increases with increasing pres-sure and a second maximum (two blue lines) appears above a pressure thresh-old of ≈ kP a . This leads to the enhancement of the 2nd harmonic in thepressure scattered by the bubble [46]. Above a second pressure threshold, PDoccurs (red line bifurcates to two branches) at 78 kPa (Fig. 1g) and 79 kPa(Fig. 1i) respectively for the NTM and FTM model. The LTM model (Fig.1h) does not predict any PD due to overestimation of the Td. Overestimationof the Td in case of the LTM is due to the very large wall velocities thatthe bubble attain in the regime of 5/2 UH emissions [58]. The P2 oscillationhas 4 maxima (Fig. 1d and 1f) indicating that the 5/2 UH component of thescattered pressure is stronger than other SH and UH components [46].Fig. 2 shows the bifurcations structure of a C3F8 bubble with R = 10 µm as a9a) (b)(c) (d)(e) (f) Fig. 3. Normalized pressure thresholds of the P2 oscillations of an air bubble asa function of frequency. The left column is for bubbles with initial radii between1 µm ≤ R ≤ µm and the right column is for bubbles with initial radii between0 . µm ≤ R ≤ . µm : (a) & (b)- Non-thermal model, (c) & (d)-linear thermal (a)& (f)- ODE thermal model. Solid black circles mark the stable oscillation regimes. function of pressure. In case of C3F8 gas core, due to the significantly weakereffects of Td [56,59], the 3 models predict the same behavior for the bubble.Thus in agreement with previous analytical and numerical studies that negel-cet or simplify Td, the lowest pressure threshold for the generation of PD isat f = 2 f r .The results indicate that the minimum pressure for the generation of PDstrongly depends on the Td. For gases like air that have high Td, the fre-quency of minimum pressure threshold for PD can shift to frequencies belowresonance. 10 .2 Pressure threshold of P2 oscillations as a function of size and frequency In this section the pressure threshold for the generation of the PD is extractedfrom the bifurcation diagrams and plotted as a function of frequency. The PDpressure threshold is normalized to the pressure threshold at f = 2 f r for sim-plicity of comparison. Moreover, the regions of non-destructive oscillations areidentified. Stable domains are identified as oscillations that satisfy the relation-ship of R max /R ≤ µm ≤ R ≤ µm for which thermal effects are stronger. The second category is for bubbles with0 . µm ≤ R ≤ . µm for which the viscous effects are stronger. Fig. 3a showsthat for air bubbles with 1 µm ≤ R ≤ µm and in the absence of Td, theminimum pressure threshold for PD occurs at f = 2 f r . This is in agreementwith the results of analytical studies [34,35,36,37,38] and numerical simula-tions [39,40] that neglect the thermal dissipation. For the studied frequencyrange, the NTM model predicts stable SH emissions for 0 . f r < f < . f r and f ≥ . f r . For bubbles with 0 . µm ≤ R ≤ . µm however, the minimumpressure threshold occurs for 0 . f r ≤ f ≤ . f r . Previous numerical studies[39,40] showed that the minimum pressure threshold for SH generation shiftsfrom twice the resonance frequency toward the resonance frequency for bubbleswith R ≤ . µm . These studies conclude that, increased damping dampensthe bubble response more at twice the resonance than at the resonance whichleads to the shift of the frequency of minimum pressure threshold. The resultspresented in Figs. 3b, 3d and 3f are in agreement with these studies and pre-dict a minimum near the resonance frequency for bubbles with R ≤ . µm in the frequency range of 0 . f r < f < . f r . However, the previous studiesdid not cover frequencies below 0 . f r and thus were not able to reveal theabsolute minimum that exists in 0 . f r ≤ f ≤ . f r . The occurrence of PDin this frequency range is due to 5/2 UH resonance (e.g. Fig. 1g (P2 with4 maxima)). The NTM model predicts that the SH emissions are stable for f ≥ . f r and 0 . f r < f < . f r (solid black circles).When linear thermal dissipation is considered, bubbles with 1 µm ≤ R ≤ µm exhibit a minimum at 2 f r , except the bubble with R = 1 µm . Theminimum threshold for SH emissions for the bubble with R = 1 µm is at f = 0 . f r . Since the smaller bubbles are more affected by viscous damping,the addition of the thermal viscosity ( µ th ) to the liquid viscosity in the LTMmodel significantly increases the total damping. Thus, the bubble behavior11ith R = 1 µm becomes similar to the behavior of the smaller bubbles withvery high viscous damping (Fig. 3b). The shift in the minimum frequency inthis case, thus, can be described with a similar explanation given in [38,39].For bubbles with 0 . µm ≤ R ≤ . µm the minimum pressure threshold oc-curs for frequencies near resonance, except for the bubble with R = 0 . µm .The addition of the µ th to the liquid viscosity significantly increase the totaldissipation in smaller bubbles (as Td becomes inversly proportional to bubbleradius) and makes the bubble oscillator an over-damped oscillator. In this case,UH oscillations are suppressed and the PD requires much higher pressures for f < f r . Moreover the P2 oscillations are not in non-destructive regime for f ≤ . f r . For bubbles with 0 . µm ≤ R ≤ . µm , the LTM predicts non-destructive SH emissions only for 1 . f r ≤ f r ≤ . f r .When thermal effects are considered with their full non-linearity and for bub-bles with 1 µm ≤ R ≤ µm (Fig. 3e), minimum pressure threshold for SHemission is when 0 . f r < f ≤ . f r . SH emissions are stable for f ≥ . f r and 0 . f r ≤ f ≤ . f r . The bigger bubble ( R = 10 µm ) may sustain non-destructive SH emissions for 0 . f r ≤ f ≤ . f r as well. For bubbles with0 . µm ≤ R ≤ . µm (Fig. 3f) the minimum pressure threshold for SH emis-sions occurs at 0 . f r ≤ f ≤ . f r , however, only the bubble with R = 0 . f r may exhibit stable SH emissions.In case of FTM, all the bubbles in this size range exhibit stable SH emissionswhen 1 . f r ≤ f ≤ . f r .Fig. 4 shows the pressure threshold of SH emissions as a function of fre-quency for C3F8 bubbles. In this case thermal damping is not as significantas air bubbles. For bubbles with 1 µm ≤ R ≤ µm all three models of NTM(Fig. 4e), LTM (Fig. 4c) and FTM (Fig. 4e) predict the minimum at f = 2 f r .The only exception is in case of the bubble with R = 1 µm where the FTMmodel predicts the minimum at 0 . f r ≤ f ≤ . f r . For C3F8 bubbles with1 µm ≤ R ≤ µm , SH emissions are nondestructive for 0 . f r ≤ f ≤ . f r and f ≥ . f r . The stability region can be extended to lower frequencies forlarger bubbles.For bubbles with 0 . µm ≤ R ≤ . µm , the NTM (Fig. 4b) and FTM (Fig.4f) predict approximately the same behavior. The minimum is located at f = 0 . f r for the bubbles with 0 . µm ≤ R ≤ . µm . At this frequency FTMpredicts stable SH emissions only for the bubble with R = 0 . µm . All threemodels predict stable SH emissions for f ≥ . f r .Predictions of the LTM (Fig. 4d) deviate from the predictions of the FTM(Fig. 4f). This is because of the overestimation of the Td byb the LTM forsmaller bubbles. The large dissipation as predicted by the LTM suppresses theUH resonance for bubbles R ≤ . µm .12a) (b)(c) (d)(e) (f) Fig. 4. Normalized pressure thresholds of the P2 oscillations of a C3F8 bubble asa function of frequency. The left column is for bubbles with initial radii between1 µm ≤ R ≤ µm and the right column is for bubbles with initial radii between0 . µm ≤ R ≤ . µm : (a) & (b)- Non-thermal model, (c) & (d)-linear thermal (a)& (f)- ODE thermal model. Solid black circles mark the stable oscillation regimes. In this study for highlighting the thermal dissipation effects we only consid-ered spherical oscillations of a single bubble and one set of initial conditions( R ( t = 0) = R , ˙ R ( t = 0) = 0 m/s ). However, non-spherical oscillations,especially in larger bubbles can significantly influence the bubble behavior[4,63,64,65] and threshold of SH emissions [31,32]. Bubble-bubble interac-tion is another important factor that changes the threshold of SH emissionsin mono-disperse bubble clouds [40] as well as poly-disperse bubble clouds[66,67,68]. Interaction with a boundary also influences the pressure threshold13f SH emissions [68]. In real applications and especially for high bubble con-centrations, the effect of bubble-bubble interactions should be considered foraccurate predictions of the stable regions of SH emissions. The initial condi-tion of the bubble is another factor that can change the minimum pressurethreshold for the occurrence of nonlinear behavior and SH emissions [3,69,70].Sonication with a multi-frequency acoustic source [7,71,72,73,74,75] signifi-cantly influences the nonlinear behavior of the bubbles and the threshold forSH emissions. For sufficiently long exposures, mass transfer effects become animportant mechanism of the damping and can significantly affect the bubblebehavior [76,77,78,79,80]. Incorporation of these effects are beyond the scopeof present study. For accurate modeling of the cavitation phenomenon theseeffects must be included in the future studies. In this study we demonstrated that the thermal damping significantly influ-ences the pressure threshold of P2 oscillations and 1/2 order SH emissions.For gases with higher thermal damping like air, full thermal effects should beincorporated in the model to accurately predict the nonlinear oscillations. Pre-dictions of the the popular linear thermal model [41,42,49,50,81] significantlydeviate from the predictions of the full thermal model that more accuratelyincorporates all the thermal effects. When thermal damping is considered withits full non-linearity, we show that the minimum pressure for period doublingand 1/2 order SH emissions for bubbles with large thermal damping (e.g. air)occurs when 0 . f r ≤ f ≤ . f r ( f r is the linear resonance frequency). Themechanism of the period doubling is via the occurrence of 5/2 UH resonance.Results show that P2 oscillations and 1/2 order SH emissions are stable when f ≥ . f r for all the studied bubble sizes (0 . µm ≤ R ≤ µm ). However,only bubbles with R ≥ . µm ( ≥ µm in diameter) may exhibit stable SHemissions when 0 . f r ≤ f ≤ . f r . The work is supported by the Natural Sciences and Engineering ResearchCouncil of Canada (Discovery Grant RGPIN-2017-06496), NSERC and theCanadian Institutes of Health Research ( Collaborative Health Research Projects) and the Terry Fox New Frontiers Program Project Grant in Ultrasound andMRI for Cancer Therapy (project eferences [1] Lauterborn, Werner, and Thomas Kurz. ”Physics of bubble oscillations.”Reports on progress in physics 73.10 (2010): 106501.[2] Parlitz, U., et al. ”Bifurcation structure of bubble oscillators.” The Journal ofthe Acoustical Society of America 88.2 (1990): 1061-1077.[3] Varga, R., Klapcsik, K. and Hegeds, F., 2020. Route to shrimps: Dissipationdriven formation of shrimp-shaped domains. Chaos, Solitons & Fractals, 130,p.109424.[4] Klapcsik, K. and Hegeds, F., 2019. Study of non-spherical bubble oscillationsunder acoustic irradiation in viscous liquid. Ultrasonics sonochemistry, 54,pp.256-273.[5] Macdonald, C.A. and Gomatam, J., 2006. Chaotic dynamics of microbubbles inultrasonic fields. Proceedings of the Institution of Mechanical Engineers, PartC: Journal of Mechanical Engineering Science, 220(3), pp.333-343.[6] Sojahrood, A.J., Earl, R., Kolios, M.C. and Karshafian, R., 2020. Investigationof the 1/2 order subharmonic emissions of the period-2 oscillations of anultrasonically excited bubble. Physics Letters A, p.126446.[7] Behnia, S., Sojahrood, A.J., Soltanpoor, W. and Jahanbakhsh, O., 2009.Suppressing chaotic oscillations of a spherical cavitation bubble throughapplying a periodic perturbation. Ultrasonics sonochemistry, 16(4), pp.502-511.[8] Shankar, P.M., Krishna, P.D. and Newhouse, V.L., 1998. Advantages ofsubharmonic over second harmonic backscatter for contrast-to-tissue echoenhancement. Ultrasound in medicine & biology, 24(3), pp.395-399.[9] Forsberg, F., Shi, W.T. and Goldberg, B.B., 2000. Subharmonic imaging ofcontrast agents. Ultrasonics, 38(1-8), pp.93-98.[10] Nam, K., Liu, J.B., Eisenbrey, J.R., Stanczak, M., Machado, P., Li, J., Li,Z., Wei, Y. and Forsberg, F., 2019. ThreeDimensional Subharmonic AidedPressure Estimation for Assessing Arterial Plaques in a Rabbit Model. Journalof Ultrasound in Medicine, 38(7), pp.1865-1873.[11] Delaney, L.J., Machado, P., Torkzaban, M., Lyshchik, A., Wessner, C.E.,Kim, C., Rosenblum, N., Richard, S., Wallace, K. and Forsberg, F., 2020.Characterization of Adnexal Masses Using ContrastEnhanced SubharmonicImaging: A Pilot Study. Journal of Ultrasound in Medicine, 39(5), pp.977-985.[12] Haworth, K.J., Mast, T.D., Radhakrishnan, K., Burgess, M.T., Kopechek,J.A., Huang, S.L., McPherson, D.D. and Holland, C.K., 2012. Passive imagingwith pulsed ultrasound insonations. The Journal of the Acoustical Society ofAmerica, 132(1), pp.544-553.
13] Haworth, K.J., Radhakrishnan, K. and Mast, T.D., 2014. Frequency-sumpassive cavitation imaging. The Journal of the Acoustical Society of America,135(4), pp.2310-2310.[14] OReilly, M.A. and Hynynen, K., 2012. Blood-brain barrier: Real-time feedback-controlled focused ultrasound disruption by using an acoustic emissionsbasedcontroller. Radiology, 263(1), pp.96-106.[15] Tsai, C.H., Zhang, J.W., Liao, Y.Y. and Liu, H.L., 2016. Real-time monitoringof focused ultrasound blood-brain barrier opening via subharmonic acousticemission detection: implementation of confocal dual-frequency piezoelectrictransducers. Physics in Medicine & Biology, 61(7), p.2926.[16] Mestas, J.L., Lenz, P. and Cathignol, D., 2003. Long-lasting stable cavitation.The Journal of the Acoustical Society of America, 113(3), pp.1426-1430.[17] Datta, S., Coussios, C.C., McAdory, L.E., Tan, J., Porter, T., De Courten-Myers, G. and Holland, C.K., 2006. Correlation of cavitation with ultrasoundenhancement of thrombolysis. Ultrasound in medicine & biology, 32(8), pp.1257-1267.[18] Bader, K.B. and Holland, C.K., 2012. Gauging the likelihood of stable cavitationfrom ultrasound contrast agents. Physics in Medicine & Biology, 58(1), p.127.[19] Radhakrishnan, K., Bader, K.B., Haworth, K.J., Kopechek, J.A., Raymond,J.L., Huang, S.L., McPherson, D.D. and Holland, C.K., 2013. Relationshipbetween cavitation and loss of echogenicity from ultrasound contrast agents.Physics in Medicine & Biology, 58(18), p.6541.[20] Bader, K.B., Gruber, M.J. and Holland, C.K., 2015. Shaken and stirred:mechanisms of ultrasound-enhanced thrombolysis. Ultrasound in medicine &biology, 41(1), pp.187-196.[21] Gudra, Matthieu, Corentin Cornu, and Claude Inserra. ”A derivation ofthe stable cavitation threshold accounting for bubble-bubble interactions.”Ultrasonics Sonochemistry 38 (2017): 168-173.[22] Cornu, C., Gudra, M., Bra, J.C., Liu, H.L., Chen, W.S. and Inserra, C., 2018.Ultrafast monitoring and control of subharmonic emissions of an unseededbubble cloud during pulsed sonication. Ultrasonics sonochemistry, 42, pp.697-703.[23] Gudra, M., Cleve, S., Mauger, C. and Inserra, C., 2020. Subharmonic sphericalbubble oscillations induced by parametric surface modes. Physical Review E,101(1), p.011101.[24] Lohse, D. and Zijm, W.H.M., 2003. Bubble puzzles. Physics Today, 56(2), pp.36-41.[25] Chahine, G.L., Kapahi, A., Choi, J.K. and Hsiao, C.T., 2016. Modelingof surface cleaning by cavitation bubble dynamics and collapse. Ultrasonicssonochemistry, 29, pp.528-549.
26] Eller, A. and Flynn, H.G., 1969. Generation of subharmonics of order onehalfby bubbles in a sound field. The Journal of the Acoustical Society of America,46(3B), pp.722-727.[27] Lauterborn, Werner, and Joachim Holzfuss. ”Acoustic chaos.” InternationalJournal of bifurcation and Chaos 1.01 (1991): 13-26.[28] Holt, R.G., Gaitan, D.F., Atchley, A.A. and Holzfuss, J., 1994. Chaoticsonoluminescence. Physical review letters, 72(9), p.1376..[29] Prosperetti, Andrea. ”A general derivation of the subharmonic threshold fornon-linear bubble oscillations.” The Journal of the Acoustical Society ofAmerica 133.6 (2013): 3719-3726[30] Longuet-Higgins, M.S., 1989. Monopole emission of sound by asymmetricbubble oscillations. Part 1. Normal modes. Journal of Fluid Mechanics, 201,pp.525-541.[31] Holt, R. Glynn, and Lawrence A. Crum. ”Acoustically forced oscillations of airbubbles in water: Experimental results.” The Journal of the Acoustical Societyof America 91.4 (1992): 1924-1932.[32] Gudra, M., Cleve, S., Mauger, C. and Inserra, C., 2020. Subharmonic sphericalbubble oscillations induced by parametric surface modes. Physical Review E,101(1), p.011101.[33] Prosperetti, Andrea. ”Nonlinear oscillations of gas bubbles in liquids:steadystate solutions.” The Journal of the Acoustical Society of America 56.3(1974): 878-885.[34] Prosperetti, Andrea. ”Application of the subharmonic threshold to themeasurement of the damping of oscillating gas bubbles.” The Journal of theAcoustical Society of America61.1 (1977): 11-16.[35] Krishna, P. D., P. M. Shankar, and V. L. Newhouse. ”Subharmonic generationfrom ultrasonic contrast agents.” Physics in Medicine & Biology 44.3 (1999):681.[36] Shankar, P. M., P. D. Krishna, and V. L. Newhouse. ”Subharmonicbackscattering from ultrasound contrast agents.” The Journal of the AcousticalSociety of America106.4 (1999): 2104-2110.[37] Kimmel, Eitan, et al. ”Subharmonic response of encapsulated microbubbles:Conditions for existence and amplification.” Ultrasound in medicine & biology33.11 (2007): 1767-1776.[38] Katiyar, Amit, and Kausik Sarkar. ”Effects of encapsulation damping on theexcitation threshold for subharmonic generation from contrast microbubbles.”The Journal of the Acoustical Society of America 132.5 (2012): 3576-3585.[39] Katiyar, Amit, and Kausik Sarkar. ”Excitation threshold for subharmonicgeneration from contrast microbubbles.” The Journal of the Acoustical Societyof America 130.5 (2011): 3137-3147.
40] Gudra, M., Cornu, C. and Inserra, C., 2017. A derivation of the stable cavitationthreshold accounting for bubble-bubble interactions. Ultrasonics Sonochemistry,38, pp.168-173.[41] Prosperetti, A., Crum, L.A. and Commander, K.W., 1988. Nonlinear bubbledynamics. The Journal of the Acoustical Society of America, 83(2), pp.502-514.[42] Commander, Kerry W., and Andrea Prosperetti. ”Linear pressure waves inbubbly liquids: Comparison between theory and experiments.” The Journal ofthe Acoustical Society of America 85, no. 2 (1989): 732-746.[43] Toegel, R., Gompf, B., Pecha, R. and Lohse, D., 2000. Does water vapor preventupscaling sonoluminescence?. Physical review letters, 85(15), p.3165.[44] Stricker, L., Prosperetti, A. and Lohse, D., 2011. Validation of an approximatemodel for the thermal behavior in acoustically driven bubbles. The Journal ofthe Acoustical Society of America, 130(5), pp.3243-3251.[45] Sojahrood, A.J., Wegierak, D., Haghi, H., Karshfian, R. and Kolios, M.C., 2019.A simple method to analyze the super-harmonic and ultra-harmonic behavior ofthe acoustically excited bubble oscillator. Ultrasonics sonochemistry, 54, pp.99-109.[46] Keller, J.B. and Miksis, M., 1980. Bubble oscillations of large amplitude. TheJournal of the Acoustical Society of America, 68(2), pp.628-633.[47] Lide, D.R. and Kehiaian, H.V., 1994. CRC handbook of thermophysical andthermochemical data. Crc Press.[48] http://cern.ch/Detector-Cooling/data/C3F8 Properties.pdf[49] Dogan, H., & Popov, V. Numerical simulation of the nonlinear ultrasonicpressure wave propagation in a cavitating bubbly liquid inside a sonochemicalreactor.
Ultrasonics sonochemistry
30 (2016): 87-97.[50] Mantouka, Agni, Hakan Dogan, P. R. White, and T. G. Leighton. ”Modellingacoustic scattering, sound speed, and attenuation in gassy soft marinesediments.” The journal of the acoustical society of America 140, no. 1 (2016):274-282.[51] Hoff, L., Sontum, P. C., & Hovem, J. M. Oscillations of polymeric microbubbles:Effect of the encapsulating shell.
The Journal of the Acoustical Society ofAmerica
TheJournal of the Acoustical Society of America
54] Shekhar, Himanshu, Nathaniel J. Smith, Jason L. Raymond, and Christy K.Holland. ”Effect of temperature on the size distribution, shell properties, andstability of Definity.” Ultrasound in medicine & biology 44, no. 2 (2018): 434-446.[55] Raymond, Jason L., Kevin J. Haworth, Kenneth B. Bader, KirthiRadhakrishnan, Joseph K. Griffin, Shao-Ling Huang, David D. McPherson, andChristy K. Holland. ”Broadband attenuation measurements of phospholipid-shelled ultrasound contrast agents.” Ultrasound in medicine & biology 40, no.2 (2014): 410-421.[56] Sojahrood, A.J., Haghi, H., Li, Q., Porter, T.M., Karshfian, R. and Kolios,M.C., 2020. Nonlinear energy loss in the oscillations of coated and uncoatedbubbles: Role of thermal, radiation and encapsulating shell damping at variousexcitation pressures. Ultrasonics Sonochemistry, p.105070.[57] Sojahrood, A.J., Haghi, H., Karshfian, R. and Kolios, M.C., 2020. Criticalcorrections to models of nonlinear power dissipation of ultrasonically excitedbubbles. Ultrasonics Sonochemistry, p.105089.[58] Sojahrood A.J., Haghi H., Karshafian R. and Kolios M.C. 2020. Classification ofthe mechanisms of ultrasound energy dissipation from the nonlinear oscillationsof coated and uncoated bubbles. Arxiv 2020[59] Flynn, H.G., Church, C.C.: Transient pulsations of small gas bubbles in water.J. Acoust. Soc. Am. 84, 985998 (1988)[60] Church, C.C., 2005. Frequency, pulse length, and the mechanical index.Acoustics Research Letters Online, 6(3), pp.162-168.[61] Church, C.C. and Yang, X., 2005. The mechanical index and cavitation in tissue.The Journal of the Acoustical Society of America, 117(4), pp.2530-2530.[62] Sojahrood, A.J., Falou, O., Earl, R., Karshafian, R. and Kolios, M.C., 2015.Influence of the pressure-dependent resonance frequency on the bifurcationstructure and backscattered pressure of ultrasound contrast agents: a numericalinvestigation. Nonlinear Dynamics, 80(1-2), pp.889-904.[63] Versluis, M., Goertz, D.E., Palanchon, P., Heitman, I.L., van der Meer, S.M.,Dollet, B., de Jong, N. and Lohse, D., 2010. Microbubble shape oscillationsexcited through ultrasonic parametric driving. Physical review E, 82(2),p.026321.[64] Dollet, B., Van Der Meer, S.M., Garbin, V., de Jong, N., Lohse, D. and Versluis,M., 2008. Nonspherical oscillations of ultrasound contrast agent microbubbles.Ultrasound in medicine & biology, 34(9), pp.1465-1473.[65] Gudra, M., Cleve, S., Mauger, C., Blanc-Benon, P. and Inserra, C., 2017.Dynamics of nonspherical microbubble oscillations above instability threshold.Physical Review E, 96(6), p.063104.
66] Haghi, H., Sojahrood, A.J. and Kolios, M.C., 2019. Collective nonlinearbehavior of interacting polydisperse microbubble clusters. Ultrasonicssonochemistry, 58, p.104708.[67] Dzaharudin, F., Suslov, S.A., Manasseh, R. and Ooi, A., 2013. Effectsof coupling, bubble size, and spatial arrangement on chaotic dynamics ofmicrobubble cluster in ultrasonic fields. The Journal of the Acoustical Societyof America, 134(5), pp.3425-3434.[68] Dzaharudin, F., Ooi, A. and Manasseh, R., 2017. Effects of boundary proximityon monodispersed microbubbles in ultrasonic fields. Journal of Sound andVibration, 410, pp.330-343.[69] Hegeds, F., 2016. Topological analysis of the periodic structures in aharmonically driven bubble oscillator near Blake’s critical threshold: Infinitesequence of two-sided Farey ordering trees. Physics Letters A, 380(9-10),pp.1012-1022.[70] Hegeds, F., Hs, C. and Kullmann, L., 2013. Stable period 1, 2 and 3 structuresof the harmonically excited RayleighPlesset equation applying low ambientpressure. The IMA Journal of Applied Mathematics, 78(6), pp.1179-1195.[71] Gudra, M., Inserra, C. and Gilles, B., 2017. Accompanying the frequency shiftof the nonlinear resonance of a gas bubble using a dual-frequency excitation.Ultrasonics Sonochemistry, 38, pp.298-305.[72] Gudra, M., Inserra, C. and Gilles, B., 2017. Accompanying the frequency shiftof the nonlinear resonance of a gas bubble using a dual-frequency excitation.Ultrasonics Sonochemistry, 38, pp.298-305.[73] Klapcsik, K., Varga, R. and Hegeds, F., 2018. Bi-parametric topology ofsubharmonics of an asymmetric bubble oscillator at high dissipation rate.Nonlinear Dynamics, 94(4), pp.2373-2389.[74] Hegeds, F., Lauterborn, W., Parlitz, U. and Mettin, R., 2018. Non-feedbacktechnique to directly control multistability in nonlinear oscillators by dual-frequency driving. Nonlinear Dynamics, 94(1), pp.273-293.[75] Zhang, Y. and Li, S., 2017. Combination and simultaneous resonances ofgas bubbles oscillating in liquids under dual-frequency acoustic excitation.Ultrasonics sonochemistry, 35, pp.431-439.[76] Zhang, Y., Guo, Z. and Du, X., 2018. Wave propagation in liquids withoscillating vapor-gas bubbles. Applied Thermal Engineering, 133, pp.483-492.[77] Zhang, Y., Gao, Y., Guo, Z. and Du, X., 2018. Effects of mass transferon damping mechanisms of vapor bubbles oscillating in liquids. Ultrasonicssonochemistry, 40, pp.120-127.[78] Prosperetti, A., 2015. The speed of sound in a gasvapour bubbly liquid. Interfacefocus, 5(5), p.20150024.
79] Gaitan, D.F. and Holt, R.G., 1999. Experimental observations of bubbleresponse and light intensity near the threshold for single bubblesonoluminescence in an air-water system. Physical Review E, 59(5), p.5495.[80] Holt, R.G. and Gaitan, D.F., 1996. Observation of stability boundaries in theparameter space of single bubble sonoluminescence. Physical review letters,77(18), p.3791.[81] Zhang Y. and Li, S., 2014. Thermal effects on nonlinear radial oscillations ofgas bubbles in liquids under acoustic excitation. International Communicationsin Heat and Mass Transfer, 53, pp.43-49.79] Gaitan, D.F. and Holt, R.G., 1999. Experimental observations of bubbleresponse and light intensity near the threshold for single bubblesonoluminescence in an air-water system. Physical Review E, 59(5), p.5495.[80] Holt, R.G. and Gaitan, D.F., 1996. Observation of stability boundaries in theparameter space of single bubble sonoluminescence. Physical review letters,77(18), p.3791.[81] Zhang Y. and Li, S., 2014. Thermal effects on nonlinear radial oscillations ofgas bubbles in liquids under acoustic excitation. International Communicationsin Heat and Mass Transfer, 53, pp.43-49.