Effects of Nonlinear Propagation of Focused Ultrasound on the Stable Cavitation of a Single Bubble

Many biomedical applications such as ultrasonic targeted drug delivery, gene therapy, and molecular imaging entail the problems of manipulating microbubbles by means of a high-intensity focused ultrasound (HIFU) pressure field; namely stable cavitation. In high-intensity acoustic field, bubbles demonstrate translational instability, the well-known erratic dancing motion, which is caused by shape oscillations of the bubbles that are excited by their volume oscillations. The literature of bubble dynamics in the HIFU field is mainly centered on experiments, lacking a systematic study to determine the threshold for shape oscillations and translational motion. In this work, we extend the existing multiphysics mathematical modeling platform on bubble dynamics for taking account of (1) the liquid compressibility which allows us to apply a high-intensity acoustic field; (2) the mutual interactions of volume pulsation, shape modes, and translational motion; as well as (3) the effects of nonlinearity, diffraction, and absorption of HIFU to incorporate the acoustic nonlinearity due to wave kinematics or medium—all in one model. The effects of acoustic nonlinearity on the radial pulsations, axisymmetric modes of shape oscillations, and translational motion of a bubble, subjected to resonance and off-resonance excitation and various acoustic pressure, are examined. The results reveal the importance of considering all the involved harmonics and wave distortion in the bubble dynamics, to accurately predict the oscillations, translational trajectories, and the threshold for inertial (unstable) cavitation. This result is of interest for understanding the bubble dynamical behaviors observed experimentally in the HIFU field.


Introduction
Acoustic cavitation is categorized into two main mechanisms: stable (non-inertial) cavitation and inertial (unstable) cavitation.Stable cavitation, the focus of this research, refers to the oscillations of bubbles in size and shape about their equilibrium radius with the insonation frequency, and are accompanied by higher harmonics, in addition to the linear driving frequency component.The inertial cavitation refers to sudden and rapid growth of bubbles (several times their initial radius), and is accompanied by a violent collapse [1][2][3][4].The importance of investigating stable cavitation is evident from the studies reporting that it reduces the risk of brain damage during the use of microbubbles in a focused ultrasound (FU)-induced blood-brain barrier (BBB) opening [5][6][7].Furthermore, the thermal effects as well as mechanical effects of stable cavitation, are considerably lower than that of inertial cavitation [8].In this research, we consider the stable cavitation caused by the oscillations of a gas bubble under an assumption that the bubble pre-existed in the FU field and is not generated by thermal effects.
The study of gas bubble dynamics during FU insonation has been the subject of interest in many applications such as medical fields in diagnostic and therapeutic procedures [1,9,10].The early observations of single-bubble sonoluminescence by Barber and Putterman [11] and Gaitan, et al. [12] motivated researchers to study the fundamentals of the coupled dynamics of ultrasonic cavitated gas bubbles spanning a wide range of nonlinear and multiphysics processes.A gas bubble excited by acoustic pressure moves either toward or away from the pressure antinode.When the amplitude of acoustic pressure approaches a threshold value (in high-intensity field), the excited bubble begins to translate erratically.This phenomenon is called dancing motion which is of great importance in many applications such as cell sorting [13][14][15].A brief review of the literature of the spherical bubble dynamics in a standing wave provides examples of proof-of-concept demonstrations.After initial investigations by Gaines [16] and Kornfeld and Suvorov [17], Eller and Crum [18] reported that the parametrically excited mode shape oscillations of a gas bubble caused by its radial pulsation in a standing sound wave, leads to its translational instability [19] and hence dancing motion.To the best of the authors' knowledge, Doinikov in 2002 [20], presented the comprehensive mathematical model for analyzing the translational motion of a spherical bubble in a high-intensity acoustic field.He modified the equation of radial motion to allow for the effects of acoustic medium compressibility.However, in this model [20], the interaction between the shape modes and translational motion was not investigated for a high-intensity field.In a later work in 2004, Doinikov [21] presented a mathematical model for the translational motion of a bubble undergoing shape oscillations while the effects of liquid compressibility are ignored.In addition, in none of the abovementioned works which employ high intensity ultrasound, the effects of finite-amplitude sound and nonlinear acoustics are incorporated.In HIFU, as a monochromatic finite-amplitude ultrasound wave is slanted during its propagation in time and space, higher harmonics of its fundamental frequency appear.Additional phenomena including shock formation, acoustic streaming (a large scale flow in a viscous fluid [22,23]), and acoustic saturation (the point at which the medium can no longer absorb any more energy from the wave) may take place.The wave distortion has local and cumulative effects due to localized pressure and temperature changes, that must be considered in the bubble dynamics simulations.
In this work, we extend Doinikov's works [20,21] for the purpose of taking account of (1) the liquid compressibility which allow us to apply high-intensity acoustic field; (2) the mutual interactions of volume pulsation, shape modes (up to four modes), and translation motion; and (3) the effects of nonlinearity, diffraction, and absorption of HIFU to incorporate the acoustic nonlinearity due to wave kinematics or medium-all in one model.The results aim to show that the radial response of the bubble and the onset of erratic dancing motion are not dependent only on the pressure amplitude but also on the harmonic components of the excitation.This feature differentiates the research in this paper from all previous studies where nonlinearity in excitation acoustic domain was not considered.
In all other previous modeling efforts to predict bubble dynamics [24][25][26][27], the wave kinematics and simulating the nonlinear focused acoustic waves were not taken into consideration.However, developing an approach based on the equations of continuum mechanics is of utmost importance to accurately incorporate the effects of acoustic nonlinearity, diffraction, and absorption; those strongly affect the acoustic wave propagation such as waveform distortion and change of pressure amplitude.The presence of higher harmonics due to transfer of energy from a fundamental component to higher harmonic components and distortion of waveform characterizes a nonlinear wave [28].The localization of the fundamental pressure field as well as the harmonics at the focal point is essential in evaluating the acoustic and thermal effects of nonlinear propagation of FU on the dynamics of stable cavitation; the amount of energy deposition is directly proportional to frequency-dependent absorption and to higher harmonics associated with the FU.Thus, stronger nonlinear effects lead to a stronger and higher number of harmonics creating enhanced localized energy deposition and therefore enhanced pressure amplitude and temperature rise in the medium surrounding the bubble.In diagnostic and therapeutic ultrasound applications, the precise prediction of acoustic wave propagation and pressure amplitude is of great importance.Moreover, to reduce the risk of tissue damage, for example, if a clinical practice is intended to avoid the violent collapse of bubbles, the resultant acoustic waves from FU transducers should be precisely modeled and analyzed.In this work, the results show the importance of incorporating all the involved harmonics and wave distortion in the bubble dynamics for accurately predicting the oscillations and translational trajectories, the onset of translational instabilities, and the threshold for inertial (unstable) cavitation.The results will be of interest for understanding the bubble dynamical behaviors observed experimentally in HIFU field.
The two common models used for simulating focused acoustic waves are Khokhlov-Zabolotskaya-Kuznetsov (KZK) and Westervelt equations [29][30][31].In both models, the approach is based on the equations of continuum mechanics.In this work, the KZK equation is chosen to model the acoustic field of FU to reduce the computational cost.We adopt the hybrid time-frequency approach using the operator-splitting technique [32] which accounts for the three effects of nonlinearity, diffraction, and absorption separately at every integration step (a linear analytical model for the pressure field due to focused transducers without including any nonlinear effects is also given by Bhargava et al. [28]).In our analysis, blood is chosen as an acoustic medium and different input parameters (such as input acoustic power and frequency) are selected to investigate the effects of nonlinearities imposed by the FU field.
Following the work done by Doinikov [21], bubble dynamics based on the Keller-Miksis model is adapted in which a second-order nonlinear mathematical model is employed to include high-intensity acoustic excitation of the bubble [33].Temporal waveform data obtained from the KZK equation is used to fit a Fourier series expansion and given as an input into the single-bubble dynamics equations.The nonlinear coupling between the radial pulsation, shape modes (up to fourth modes of oscillations), and translational motion of the oscillating bubble are analyzed and compared in linear and nonlinear pressure fields.In this research, we investigate a multiphysics model that incorporate the FU field as an input parameter into the bubble dynamics equations.The proposed model accounts for nonlinearities caused by the acoustic field as well as the bubble's oscillations.The scattered pressure in the far-field liquid domain originating from the bubble oscillations (acoustic streaming effect) was also studied, by the authors, in the FU field and various ranges of frequency [34].
This paper is organized as follows.In Section 2, we describe a multiphysics KZK-Keller-Miksis theoretical model that interconnects the nonlinear acoustic field from a FU transducer to the dynamics of a single microbubble at focal point.In addition, the numerical solution is presented.In Section 3, we present the numerical results which entail the linear and nonlinear input acoustic pressure waveforms along with the corresponding bubble volume oscillations, shape modes, and translational motion for varying frequencies of excitation (Section 3.1) and the input power to the focused transducer, which generates various input pressure at a fixed frequency (Section 3.2).Conclusions are summarized in Section 4.

Nonlinear Acoustic Pressure from FU Transducer
The KZK equation is used to predict the nonlinear acoustic pressure field generated by a FU transducer in thermoviscous fluid domains, solved by a numerical approach.This model incorporates diffraction, absorption, and nonlinearity in the parabolic wave equation.In cylindrical coordinates, the KZK equation for an axisymmetric sound beam propagating in the axial direction (positive z) is given as [35][36][37] ∂ 2 p ex ∂z∂t where p ex is an external acoustic pressure, z is the distance in axial direction, t = t − z/c is the retarded time, c is the small-signal sound speed in a fluid medium, and ρ is the ambient density of the fluid.The first term on the right-hand side of Equation (1) represents diffraction, the second term signifies absorption with δ as the diffusivity of sound and nonlinearity is given by the third term, where β = 1 + B/2A is the coefficient of nonlinearity and B/A is defined as the parameter of nonlinearity.The transverse Laplacian operator is given as ∇ 2 ⊥ = ∂ 2 /∂r 2 + (1/r)∂/∂r, where r is the transverse radial coordinate.
The focused source condition is given as p ex = p 0 f (r, t + r 2 /2cd) z=0,r≤a [38], where p 0 is the characteristic source pressure calculated as , d is the focal length, a and b are the characteristic inner and outer radius of sound source, respectively (shown in Figure 1), and P is the input acoustic power.The function f (r, t + r 2 /2cd) indicates dependence of source condition on time as a function of the radial coordinate.
A dimensionless form of the KZK equation is written as [39] ∂ 2 p ex ∂τ ∂z where p ex = p ex /p 0 is a dimensionless pressure and z = z/d is a dimensionless form of propagation distance.The retarded time is nondimensionalized with angular source frequency ω as τ = ω (t − z/c).Dimensionless parameter G, N, and A account for focusing gain, nonlinearity, and absorption, respectively.The parameter G = z 0 /d is the linear focusing gain, where z 0 = ωa 2 /2c is the Rayleigh distance.The term nonlinearity is N = d/z s , where z s = ρc 3 /βωp 0 is the plane-wave shock formation distance.Furthermore, the parameter A is defined as A = αd, where α = δω 2 /2c 3 is the thermoviscous attenuation coefficient.The dimensionless form of transverse Laplacian operator is given as ∇ 2 r = ∂ 2 /∂r 2 + (1/r)∂/∂r, where r = r/a.The dimensionless source condition becomes p ex = f (r, τ + Gr 2 ) z=0 , which is valid under the following three assumptions: sin 2 φ << 2, G << 8π/ sin 2 φ and 1 << k a << 16π/ sin 3 φ [39,40], where k = ω/c is the wave number.Here φ is the half aperture angle (Figure 1) and sin φ = a/d.While the first assumption is a relatively weak restriction which requires the half aperture angle to be less than 25 • , the other two are analogous and restrict the maximum focusing gain for a given aperture angle.Therefore, the given assumptions limit the domain of validity of the KZK equation for a focused transducer.Since, the transducers geometries used for our study satisfy the aforementioned conditions, the KZK equation is used to predict FU acoustic pressure field applied on the surface of the bubble.
In Equation ( 2), the diffraction is included by the equation Equation ( 3) is solved numerically using finite difference numerical approach to take into account the diffraction effects in the KZK equation.Absorption refers to the thermoviscous attenuation of the waves as they propagate through a medium leading to loss of energy of the wave.The KZK equation includes the absorption by the equation where A depends on the attenuation coefficient.The attenuation follows a power-law dependence on frequency given as [41,42] where α 0 and ν are empirical constants, with the value of ν typically given as 0 < ν ≤ 2. Nonlinearity results from variations of the pressure with variations of the density in a medium.Nonlinearity parameter / B A can also be expressed in terms of variations in the speed of sound.In nonlinear acoustics, as the wave propagates, waveform peaks may travel faster than troughs, thus distorting and steepening the profile of the wave.
The equation of state ( , ) ex p p s ρ = is expressed by Taylor series expansion as follows [30] 2 where p is the pressure, s is the specific entropy, 0 p and 0 ρ represent the pressure and density in the unperturbed state, respectively.The values of A and B , which are the coefficients of first and second order terms of Equation (6), are obtained as N integrates β into the KZK equation.Thus, the nonlinearity term in the KZK equation, while neglecting the diffraction and absorption effects, is expressed as Equation ( 7) has its origin in a simpler model known as the lossless Burger's equation [43].The KZK equation is an augmentation to the Burger's equation, in which it comprises the diffraction, absorption, and nonlinearity effects in a focused acoustic field.The numerical method applied in this work, to solve the KZK equation, is based on a split-step method [44] employed in MATLAB ® (2016a, Mathworks, Natick, MA, USA).In this method, the linear parts (Equations (3) and ( 4)) and nonlinear term (Equation ( 7)) are integrated separately.At each integration, the linear parts are calculated in the frequency-domain and then the solution is converted to the time-domain using a fast Fourier transform algorithm to solve the nonlinear part of the KZK equation.After integration, the result is again transformed into the frequency-domain.This method is repeated for each of the harmonics used in the solution.
To convert the KZK equation to frequency-domain representation, a time-harmonic representation of the pressure field is expressed as [45,46] *  Nonlinearity results from variations of the pressure with variations of the density in a medium.Nonlinearity parameter B/A can also be expressed in terms of variations in the speed of sound.In nonlinear acoustics, as the wave propagates, waveform peaks may travel faster than troughs, thus distorting and steepening the profile of the wave.
The equation of state p = p ex (ρ, s) is expressed by Taylor series expansion as follows [30] where p is the pressure, s is the specific entropy, p 0 and ρ 0 represent the pressure and density in the unperturbed state, respectively.The values of A and B, which are the coefficients of first and second order terms of Equation ( 6), are obtained as A = ρ 0 (∂p/∂ρ) ρ=ρ 0 ,s=s 0 ≡ ρ 0 c 2 and B = ρ 0 2 (∂ 2 p/∂ρ 2 ) ρ=ρ 0 ,s=s 0 , where s 0 indicates that Equation ( 6) is carried out in an isentropic process.
Consequently, the nonlinearity parameter B/A is written as B/A = ρ 0 /c 2 (∂ 2 p/∂ρ 2 ) ρ=ρ 0 ,s=s 0 and an alternative expression of B/A yields B/A = 2ρ 0 c(∂c/∂p)| ρ=ρ 0 ,s=s 0 .In Equation (2), the parameter N integrates β into the KZK equation.Thus, the nonlinearity term in the KZK equation, while neglecting the diffraction and absorption effects, is expressed as Equation ( 7) has its origin in a simpler model known as the lossless Burger's equation [43].The KZK equation is an augmentation to the Burger's equation, in which it comprises the diffraction, absorption, and nonlinearity effects in a focused acoustic field.The numerical method applied in this work, to solve the KZK equation, is based on a split-step method [44] employed in MATLAB ® (2016a, Mathworks, Natick, MA, USA).In this method, the linear parts (Equations ( 3) and ( 4)) and nonlinear term (Equation ( 7)) are integrated separately.At each integration, the linear parts are calculated in the frequency-domain and then the solution is converted to the time-domain using a fast Fourier transform algorithm to solve the nonlinear part of the KZK equation.After integration, the result is again transformed into the frequency-domain.This method is repeated for each of the harmonics used in the solution.
To convert the KZK equation to frequency-domain representation, a time-harmonic representation of the pressure field is expressed as [45,46] where u n is the dimensionless complex pressure corresponding to nth harmonic and u n * is the conjugate of dimensionless complex pressure.After substituting Equation (8) into Equation (2), the dimensionless KZK equation becomes where n u is the dimensionless complex pressure corresponding to nth harmonic and * n u is the conjugate of dimensionless complex pressure.After substituting Equation (8) into Equation (2), the dimensionless KZK equation becomes where n  is a complex number representing both absorption and sound velocity dispersion.Equation ( 9) accounts for frequency-domain representation of the dimensionless KZK equation.The real part of n  is proportional to the absorption parameter, which follows the frequency dependent power law (Equation ( 5)).The imaginary part of n  accounts for dispersion (phase velocity) effects.
The dispersion effects are necessary to maintain the causality of the system and are determined using Kramers-Kronig relations as [47] where 0 c is the sound velocity at reference frequency 0  and () c  represents the phase velocity.
Equation (10) takes into account the local generalized ultrasonic attenuation-dispersion relations.As a result, the linear part of Equation ( 9), written as frequency-domain while the remaining nonlinear term is solved in time-domain as mentioned earlier.
The boundary conditions, which satisfy parabolic approximations of a spherically converging harmonic wave having a uniform pressure distribution, are given as Soneson and Myers [45] proposed a criterion to determine the cases for which nonlinear effects can be neglected without significant errors to improve computing efficiency in time-domain.The criterion uses attenuation, a nonlinearity parameter, and linear focusing gain to determine a cutoff value and decide if the linear model is sufficient for the analysis.Thus, in this work, before the model computes for nonlinearity, a threshold cutoff amplitude of the solution is determined.If the amplitude is greater than the cutoff value, then the model accounts for nonlinearity and integrates Equation (7) using the upwind method.The number of harmonics () n included in the model is given by 2 k n  , where k is an integer.To accurately model nonlinear effects, 6 k  is used.To solve in frequency-domain, second-order, diagonally-implicit Runge-Kutta (DIRK), and the Crank-Nicolson (CN) method, which is based on trapezoidal rule [48], are used.The second-order DIRK method is used in the near-field in the boundary layer while CN method is used in the far-field region [44].

Single-Bubble Dynamics in a FU Field
The surface of a spherical gas bubble, assuming axisymmetric deformations (with respect to xaxis, Figure 1) is defined as where r and  are the local spherical coordinates, originating at the center of the bubble, () Rt is the time-dependent radius of the bubble corresponding to the fundamental mode, () n at is the time- dependent amplitude of shape modes, and n P is the Legendre polynomial of order n.The dynamical equations describing the radial pulsation, translational motion, and shape oscillations of the bubble up to the fourth shape mode are derived and given as [21].
where EER REVIEW 6 imensionless complex pressure corresponding to nth harmonic and * n u is the sionless complex pressure.After substituting Equation ( 8) into Equation ( 2), the equation becomes mplex number representing both absorption and sound velocity dispersion.nts for frequency-domain representation of the dimensionless KZK equation.The roportional to the absorption parameter, which follows the frequency dependent on ( 5)).The imaginary part of n  accounts for dispersion (phase velocity) effects.cts are necessary to maintain the causality of the system and are determined using lations as [47]  into account the local generalized ultrasonic attenuation-dispersion relations.As part of Equation ( 9), written as while the remaining nonlinear term is solved in time-domain as mentioned earlier.
ditions, which satisfy parabolic approximations of a spherically converging ving a uniform pressure distribution, are given as z     yers [45] proposed a criterion to determine the cases for which nonlinear effects ithout significant errors to improve computing efficiency in time-domain.The uation, a nonlinearity parameter, and linear focusing gain to determine a cutoff the linear model is sufficient for the analysis.Thus, in this work, before the model linearity, a threshold cutoff amplitude of the solution is determined.If the er than the cutoff value, then the model accounts for nonlinearity and integrates the upwind method.The number of harmonics () n included in the model is here k is an integer.To accurately model nonlinear effects, 6 k  is used.To -domain, second-order, diagonally-implicit Runge-Kutta (DIRK), and the Crankthod, which is based on trapezoidal rule [48], are used.The second-order DIRK he near-field in the boundary layer while CN method is used in the far-field region ynamics in a FU Field a spherical gas bubble, assuming axisymmetric deformations (with respect to xfined as re the local spherical coordinates, originating at the center of the bubble, () Rt is t radius of the bubble corresponding to the fundamental mode, () n at is the time- de of shape modes, and n P is the Legendre polynomial of order n.The dynamical g the radial pulsation, translational motion, and shape oscillations of the bubble ape mode are derived and given as [21].
n is a complex number representing both absorption and sound velocity dispersion.Equation ( 9) accounts for frequency-domain representation of the dimensionless KZK equation.The real part of , 2 FOR PEER REVIEW 6 is the dimensionless complex pressure corresponding to nth harmonic and * n u is the f dimensionless complex pressure.After substituting Equation ( 8) into Equation ( 2), the ess KZK equation becomes is a complex number representing both absorption and sound velocity dispersion.
) accounts for frequency-domain representation of the dimensionless KZK equation.The n  is proportional to the absorption parameter, which follows the frequency dependent (Equation ( 5)).The imaginary part of n  accounts for dispersion (phase velocity) effects.sion effects are necessary to maintain the causality of the system and are determined using ronig relations as [47] is the sound velocity at reference frequency 0  and () c  represents the phase velocity.0) takes into account the local generalized ultrasonic attenuation-dispersion relations.As e linear part of Equation ( 9), written as is used.To quency-domain, second-order, diagonally-implicit Runge-Kutta (DIRK), and the Crank-N) method, which is based on trapezoidal rule [48], are used.The second-order DIRK sed in the near-field in the boundary layer while CN method is used in the far-field region ubble Dynamics in a FU Field rface of a spherical gas bubble, assuming axisymmetric deformations (with respect to x-1) is defined as nd  are the local spherical coordinates, originating at the center of the bubble, () Rt is pendent radius of the bubble corresponding to the fundamental mode, () n at is the time- amplitude of shape modes, and n P is the Legendre polynomial of order n.The dynamical escribing the radial pulsation, translational motion, and shape oscillations of the bubble urth shape mode are derived and given as [21].
n is proportional to the absorption parameter, which follows the frequency dependent power law (Equation ( 5)).The imaginary part of Acoustics 2019, 2 FOR PEER REVIEW 6 where n u is the dimensionless complex pressure corresponding to nth harmonic and * n u is the conjugate of dimensionless complex pressure.After substituting Equation (8) into Equation ( 2), the dimensionless KZK equation becomes where n  is a complex number representing both absorption and sound velocity dispersion.Equation ( 9) accounts for frequency-domain representation of the dimensionless KZK equation.The real part of n  is proportional to the absorption parameter, which follows the frequency dependent power law (Equation ( 5)).The imaginary part of n  accounts for dispersion (phase velocity) effects.
The dispersion effects are necessary to maintain the causality of the system and are determined using Kramers-Kronig relations as [47]  Equation (10) takes into account the local generalized ultrasonic attenuation-dispersion relations.As a result, the linear part of Equation ( 9), written as Soneson and Myers [45] proposed a criterion to determine the cases for which nonlinear effects can be neglected without significant errors to improve computing efficiency in time-domain.The criterion uses attenuation, a nonlinearity parameter, and linear focusing gain to determine a cutoff value and decide if the linear model is sufficient for the analysis.Thus, in this work, before the model computes for nonlinearity, a threshold cutoff amplitude of the solution is determined.If the amplitude is greater than the cutoff value, then the model accounts for nonlinearity and integrates Equation ( 7) using the upwind method.The number of harmonics () n included in the model is given by 2 k n  , where k is an integer.To accurately model nonlinear effects, 6 k  is used.To solve in frequency-domain, second-order, diagonally-implicit Runge-Kutta (DIRK), and the Crank-Nicolson (CN) method, which is based on trapezoidal rule [48], are used.The second-order DIRK method is used in the near-field in the boundary layer while CN method is used in the far-field region [44].

Single-Bubble Dynamics in a FU Field
The surface of a spherical gas bubble, assuming axisymmetric deformations (with respect to xaxis, Figure 1) is defined as where r and  are the local spherical coordinates, originating at the center of the bubble, () Rt is the time-dependent radius of the bubble corresponding to the fundamental mode, () n at is the time- dependent amplitude of shape modes, and n P is the Legendre polynomial of order n.The dynamical equations describing the radial pulsation, translational motion, and shape oscillations of the bubble up to the fourth shape mode are derived and given as [21].
n accounts for dispersion (phase velocity) effects.The dispersion effects are necessary to maintain the causality of the system and are determined using Kramers-Kronig relations as [47] where c 0 is the sound velocity at reference frequency ω 0 and c(ω) represents the phase velocity.Equation (10) takes into account the local generalized ultrasonic attenuation-dispersion relations.As a result, the linear part of Equation ( 9), written as where n u is the dimensionless complex pressure corresponding to nth harm conjugate of dimensionless complex pressure.After substituting Equation ( 8) in dimensionless KZK equation becomes where n  is a complex number representing both absorption and sound Equation ( 9) accounts for frequency-domain representation of the dimensionless real part of n  is proportional to the absorption parameter, which follows the fr power law (Equation ( 5)).The imaginary part of n  accounts for dispersion (ph The dispersion effects are necessary to maintain the causality of the system and a Kramers-Kronig relations as [47] frequency-domain while the remaining nonlinear term is solved in time-domain a The boundary conditions, which satisfy parabolic approximations of a sph harmonic wave having a uniform pressure distribution, are given as ( , 0, ) Soneson and Myers [45] proposed a criterion to determine the cases for wh can be neglected without significant errors to improve computing efficiency i criterion uses attenuation, a nonlinearity parameter, and linear focusing gain to value and decide if the linear model is sufficient for the analysis.Thus, in this wo computes for nonlinearity, a threshold cutoff amplitude of the solution is amplitude is greater than the cutoff value, then the model accounts for nonline Equation (7) using the upwind method.The number of harmonics () given by 2 k n  , where k is an integer.To accurately model nonlinear effect solve in frequency-domain, second-order, diagonally-implicit Runge-Kutta (DI Nicolson (CN) method, which is based on trapezoidal rule [48], are used.The method is used in the near-field in the boundary layer while CN method is used i [44].

Single-Bubble Dynamics in a FU Field
The surface of a spherical gas bubble, assuming axisymmetric deformation axis, Figure 1) is defined as  n u n = 0 is solved in frequency-domain while the remaining nonlinear term is solved in time-domain as mentioned earlier.The boundary conditions, which satisfy parabolic approximations of a spherically converging harmonic wave having a uniform pressure distribution, are given as p ex (r, 0, τ) = f (r, τ + Gr 2 ) and Soneson and Myers [45] proposed a criterion to determine the cases for which nonlinear effects can be neglected without significant errors to improve computing efficiency in time-domain.The criterion uses attenuation, a nonlinearity parameter, and linear focusing gain to determine a cutoff value and decide if the linear model is sufficient for the analysis.Thus, in this work, before the model computes for nonlinearity, a threshold cutoff amplitude of the solution is determined.If the amplitude is greater than the cutoff value, then the model accounts for nonlinearity and integrates Equation (7) using the upwind method.The number of harmonics (n) included in the model is given by n = 2 k , where k is an integer.To accurately model nonlinear effects, k ≥ 6 is used.To solve in frequency-domain, second-order, diagonally-implicit Runge-Kutta (DIRK), and the Crank-Nicolson (CN) method, which is based on trapezoidal rule [48], are used.The second-order DIRK method is used in the near-field in the boundary layer while CN method is used in the far-field region [44].

Single-Bubble Dynamics in a FU Field
The surface of a spherical gas bubble, assuming axisymmetric deformations (with respect to x-axis, Figure 1) is defined as where r and θ are the local spherical coordinates, originating at the center of the bubble, R(t) is the time-dependent radius of the bubble corresponding to the fundamental mode, a n (t) is the time-dependent amplitude of shape modes, and P n is the Legendre polynomial of order n.
The dynamical equations describing the radial pulsation, translational motion, and shape oscillations of the bubble up to the fourth shape mode are derived and given as [21].
where P sc accounts for scattered pressure in the surface of the bubble, which is expressed as where the coefficient A 0 represents the external acoustic pressure and the coefficients A n (n = 1, . . ., 4) indicate the acoustic excitations of the shape modes owing to radial pulsation of the bubble.Moreover, P 0 is the ambient (hydrostatic) pressure, σ is the surface tension of the infinite incompressible host liquid medium, µ is the dynamic viscosity of the liquid, γ is the specific heat ratio of the gas within the bubble, and .
x = v, which is the translational velocity of the bubble.It should be noted that the left-hand side of Equation ( 12) originates from the well-known dynamical model of bubble oscillation called Rayleigh-Plesset equation [49], which is written as follows: However, it is observed that a high-intensity acoustic field leads to nonlinear oscillations of a bubble with a higher amplitude.As a result, the radial velocity of the bubble can no longer be assumed small compared to the speed of sound in the liquid.Moreover, the effects of liquid compressibility in nonlinear and high-intensity oscillations of the bubble must be taken into consideration.To account for the aforementioned effects and acoustic radiation, the left-hand side of Equation ( 18) is replaced by the Keller-Miksis model [50] (Equation ( 12)).
It should be pointed out that the equations of motion are derived up to the second-order terms due to the fact that amplitudes of the surface modes are assumed to be small compared to its radius [51].For simplicity, the time dependency of R, a n , x and A n is excluded in the equations.The coefficients A n are calculated based on Equations ( 27)-(30) [21].
where P ac | S accounts for external acoustic pressure applied on the surface of the bubble.The external acoustic pressure originating from the plane sanding wave at any arbitrary point in the surrounding liquid is calculated from where and P a is the amplitude of driving acoustic pressure, ω is the angular driving frequency and k = 2π/λ is the wave number, where λ is the wavelength.The wavelength is expressed as λ = c/ f , where f is the excitation frequency.Here x(t) is the distance between the center of bubble to the nearest acoustic pressure antinode (x = 0), r defines the position vector originated at the center of the bubble and x 0 denotes initial position of the bubble centroid (Figure 1).Substituting Equation (28) into Equation ( 27) at r = R(t), one obtains the coefficients A n (t) as follows: A where j n is the spherical Bessel function of order n.
To incorporate p ex from FU transducer into Equations ( 28) and (30), temporal waveform data from the KZK equation is used to fit a Fourier series expansion.In our study, we assume that the bubble is located at the focal point.The corresponding Fourier series for p ex , which is a sum of sine and cosine functions and forms a periodic signal, can be described by where X 0 , X i and Y i are the Fourier constant coefficients.Using Curve Fitting toolbox in MATLAB ® , we obtain the Fourier coefficients in Equation ( 31) for the linear and nonlinear pressure fields.

Numerical Results
In this section, a comprehensive numerical investigation is carried out to analyze the effects of nonlinear ultrasonic pressure due to the change of frequency and input transducer power, at the focal point, on the bubble dynamics.We investigate the dependence of the bubble dynamics, including the bubble growth and its translational motion on the imposed acoustic pressure.Although there have been considerable theoretical and experimental studies concerning the radial oscillation and translational motion of a bubble in an acoustic field of high intensity [20,21,[25][26][27]52], a model that incorporates the effects of nonlinearity, diffraction, and absorption of FU field has not yet been investigated.In addition, fluid compressibility and shape modes (up to four modes) are considered in this work.
In the usage of FU-excited microbubbles, focused acoustic waves are strongly affected by the mutual influences of diffraction, absorption, and nonlinearity of the medium.Therefore, investigation of these effects in the focal area and the resultant FU pressure field is clearly of importance.In nonlinear acoustics, energy is transferred from a fundamental component to higher harmonic components.Generation of higher harmonics due to the increased energy transfer results in waveform steepening and distortion.The results in this section will show that the radial response of the bubble and the translational trajectories are dependent on both the pressure amplitude as well as the harmonic components of the excitation.
The parameters of the acoustic medium (blood) used in the KZK model, accounting for attenuation and nonlinearity effects, density, and sound speed are listed in Table 1 [53,54].The geometric properties of the focused transducer used in the simulations, including the inner and outer radius (a and b) and focal length, d, are given as 25, 10 and 80 mm, respectively.The transducer is based on the case study explained by Soneson [55].

Effects of Excitation Frequency on Bubble Dynamics
The numerical simulations are done assuming the bubble is excited by an external acoustic pressure at three frequencies.The excitation frequencies are selected to include fundamental resonance frequency, below and above fundamental resonance.The natural resonance frequencies of a bubble for the fundamental radial mode and axisymmetric modes of oscillations are given by [56,57] where n is the number of modes.Accordingly, the value of fundamental frequency in blood domain is estimated as f 0 = 0.5 MHz for a gas bubble with the initial radius R 0 = 7 µm and the given properties, mentioned later.
To analyze the effects of change of frequency of FU excitation on single-bubble dynamics, the external acoustic pressure waveform generated by the focused transducer is obtained in the blood domain as an input parameter to the bubble dynamics.Later, based on Equation ( 31), the equivalent harmonic excitation from Fourier series expansion is derived and included in the bubble dynamics.The input acoustic pressures, considering linear and nonlinear excitations at f = 0.3, 0.5 and 3 MHz (P = 2 W) are illustrated in Figure 2a.Furthermore, the corresponding frequency spectrum of nonlinear acoustic excitations, in which β = 4.1, are shown in Figure 2b.
( ) ( )( ) where n is the number of modes.Accordingly, the value of fundamental frequency in blood domain is estimated as 0 0.5 MHz f = for a gas bubble with the initial radius 0 7 μm R = and the given properties, mentioned later.
To analyze the effects of change of frequency of FU excitation on single-bubble dynamics, the external acoustic pressure waveform generated by the focused transducer is obtained in the blood domain as an input parameter to the bubble dynamics.Later, based on Equation ( 31), the equivalent harmonic excitation from Fourier series expansion is derived and included in the bubble dynamics.The input acoustic pressures, considering linear and nonlinear excitations at 0.3, 0.5 and 3 MHz f = ) are illustrated in Figure 2a.Furthermore, the corresponding frequency spectrum of nonlinear acoustic excitations, in which  It is observed that the nonlinearity at 3 MHz f = is more pronounced leading to a greater distorted pressure waveform, which in turn gives rise to the number and strength of harmonics, predominantly in the second and third harmonics.Thus, this study is focused on involving all the harmonics in dynamics of FU-excited microbubbles.Comparing the acoustic pressure waves in blood It is observed that the nonlinearity at f = 3 MHz is more pronounced leading to a greater distorted pressure waveform, which in turn gives rise to the number and strength of harmonics, predominantly in the second and third harmonics.Thus, this study is focused on involving all the harmonics in dynamics of FU-excited microbubbles.Comparing the acoustic pressure waves in blood domain with that obtained in water domain in our previous study [34], one can see the effect of attenuation on the pressure amplitude only at higher frequencies where smaller value of the attenuation parameter of water, leads to higher pressure amplitude at f = 3 MHz; however, at f = 0.3 and 0.5 MHz, the amplitudes of acoustic pressure are comparable.
The coupled equations of motion, which govern radial pulsation, translational motion, and shape oscillations of a single air bubble in a free-field blood domain, are solved numerically by the program package in MATHEMATICA.The theory is given for the general case of a gas bubble in a host liquid where specific heat ratio of air in Equation ( 17) and blood properties are given for the following results.The position of the bubble is set to x 0 = 0.001λ and the surface tension and dynamic viscosity of blood are defined as σ = 0.053 N/m [58] and µ = 0.0035 kg/(m•s) [59].Also, ambient pressures are given as P 0 = 101.3kPa and γ = 1.4.
Dynamics of the bubble under nonlinear acoustic excitation are analyzed using Equations ( 12) to (16).The equations are nonlinearly coupled owing to functions H n (t), where n = 0, 1, 2, 3, 4 and I n (t), where n = 2, 3, 4. The second-order interactions of two neighboring surface modes, namely a 2 and a 3 as well as a 3 and a 4 in the term H 1 (t), contribute to translational motion of the bubble (Equation ( 13)).The nonlinear interactions of two neighboring surface modes are also observed in the term H 3 (t), concluding that odd modes play important roles in deformations of all other modes.Similarly, the evidence of nonlinear interactions between translational motion and all other modes becomes apparent in the terms I n (t) and Equations ( 12) and ( 14) due to the velocity squared term, . x 2 .
In addition, there are strong two-way couplings between even modes (due to the terms a 2 and a 4 in H 2 (t) and H 4 (t)) which are also coupled with the radial mode as well.However, it is seen that the translational motion is mostly affected by the odd modes.Figures 3-5 show the time evolution of the amplitudes of radius and second to fourth modes of oscillations as well as translational motion of the bubble excited by different frequencies (at resonance and off-resonance frequency) for the input power of P = 2 W. The corresponding linear and nonlinear acoustic pressures are given in Figure 2a and as an example prolonged dynamical motions of the bubble excited at f = 3 MHz is shown in Appendix A, Figure A1.It is worthwhile to note that the bubble deformation dominates by the radial pulsation, which supports the assumption of small amplitude of surface modes compared to the radius of the bubble.In Figures 3-5, it is shown that the evolution of the third and fourth shape modes, along with the translational motion and velocity of the bubble are significantly larger or more unstable in the nonlinear excitation than that of the linear case and we observe more instability in the oscillations of the second mode.Moreover, the acoustic nonlinearity expedites the instability of the radial pulsation of the bubble in Figures 3a, 4a and 5a.The significant effects of acoustic nonlinearity on the translational motion of the bubble is seen at resonance frequency (see Figure 4e,f) where the bubble abruptly turns back toward the pressure antinode and shows an erratic translational (dancing) motion.This behavior of a bubble is favorable where in-line forward and backward movement of bubbles (or particles/cells) is needed in microfluidic devices for sorting applications.Among the three ranges of frequency investigated in this study, the results show that the changes in the dynamical behavior of the bubble and the translational instability due to acoustic nonlinear effects are more significant for higher frequencies due to the stronger and higher number of harmonics creating enhanced localized energy deposition surrounding the bubble.
Previous studies show that microbubbles exhibit nonlinear responses; however, the results in this study (Figures 3-5) show that the responses of the bubble are also significantly dependent on the excitation harmonics and hence they are different for linear and nonlinear acoustic excitations.The changes in the dynamical behavior of the bubble and the translational instability due to acoustic nonlinear effects are even more significant at higher frequencies (above resonance frequency).As shown in Figure 6, the Fourier spectra of the radial oscillation of the bubble (based on Figure A1a) verifies the nonlinear response of the bubble for both linear and nonlinear excitations.The spectral analysis reveals that in addition to the driving frequency component 0 ( ) f and second harmonic 0 (2 ) f , the bubble generates higher harmonics f [52,64].Although both of the excitations lead to the nonlinear response of the bubble, when the bubble is driven by the nonlinear excitation, the intensity of all the frequency components exceeds that of the linear excitation and the spectrum develops into a more complicated structure.Additionally, an approximate bandgap oscillation is found between 4 to 6 MHz frequencies for linear excitation whereas accounting for acoustic medium nonlinearity ( 0 β ≠ ) causes a broadband improvement of the oscillatory response of the bubble.When focused ultrasound is applied, the level of bubbles' oscillations in size and shape (stable cavitation) and accurately predicting translational motion are of utmost importance in biomedical applications such as blood-brain barrier (BBB) opening to enhance drug delivery to the brain.Hence, taking account of the nonlinearity parameter in ultrasound is important to predict bubbles' oscillatory behaviors more accurately.The results obtained by classical theory show that in a weak acoustic standing field, the bubbles congregate at pressure nodes or antinodes depending on the driving frequency [60].However, later studies confirm that, in the high-intensity acoustic field, bubbles show an erratic dancing motion caused by shape oscillations that are excited by the bubble volume pulsations [18].Indication of the translational instability is seen in the translational velocity and motion of the bubble (Figure 3e,f, Figure 4e,f and Figure 5e,f).The swift growth of the unstable third mode, in Figure 3c, Figure 4c, and Figure 5c, leads to instability of the second mode in Figure 3b, Figure 4b, and Figure 5b due to their nonlinear interaction.Moreover, the nonlinear interaction of the third and fourth modes contributes further in erratic translational motions of the bubble.The initiation of the erratic dancing motion of the bubble in the acoustic field in Figure 3e, Figure 4e, and Figure 5e,f is compatible with the onset of instability of other modes.Moreover, based on Figures 3-5, one can anticipate the inertial (unstable) cavitation threshold, in which the expansion ratio can be assumed around 2.3, i.e., R max /R 0 = 2.3 [61][62][63] and the bubble may collapse.Accordingly, the inertial cavitation may occur at f = 0.3 and 0.5 MHz, after few cycles of radial pulsation of the bubble, whereas the bubble exhibits stable cavitation at f = 3 MHz.
In Figures 3-5, it is shown that the evolution of the third and fourth shape modes, along with the translational motion and velocity of the bubble are significantly larger or more unstable in the nonlinear excitation than that of the linear case and we observe more instability in the oscillations of the second mode.Moreover, the acoustic nonlinearity expedites the instability of the radial pulsation of the bubble in Figures 3-5.The significant effects of acoustic nonlinearity on the translational motion of the bubble is seen at resonance frequency (see Figure 4e,f) where the bubble abruptly turns back toward the pressure antinode and shows an erratic translational (dancing) motion.This behavior of a bubble is favorable where in-line forward and backward movement of bubbles (or particles/cells) is needed in microfluidic devices for sorting applications.Among the three ranges of frequency investigated in this study, the results show that the changes in the dynamical behavior of the bubble and the translational instability due to acoustic nonlinear effects are more significant for higher frequencies due to the stronger and higher number of harmonics creating enhanced localized energy deposition surrounding the bubble.
Previous studies show that microbubbles exhibit nonlinear responses; however, the results in this study (Figures 3-5) show that the responses of the bubble are also significantly dependent on the excitation harmonics and hence they are different linear and nonlinear acoustic excitations.The changes in the dynamical behavior of the bubble and the translational instability due to acoustic nonlinear effects are even more significant at higher frequencies (above resonance frequency).As shown in Figure 6, the Fourier spectra of the radial oscillation of the bubble (based on Figure A1a) verifies the nonlinear response of the bubble for both linear and nonlinear excitations.The spectral analysis reveals that in addition to the driving frequency component ( f 0 ) and second harmonic (2 f 0 ), the bubble generates higher harmonics (3 f 0 , 4 f 0 , . ..), subharmonics (less than the driving frequency component such as 0.5 f 0 ) and ultraharmonics (such as 1.5 f 0 ) [52,64].Although both of the excitations lead to the nonlinear response of the bubble, when the bubble is driven by the nonlinear excitation, the intensity of all the frequency components exceeds that of the linear excitation and the spectrum develops into a more complicated structure.Additionally, an approximate bandgap oscillation is found between 4 to 6 MHz frequencies for linear excitation whereas accounting for acoustic medium nonlinearity (β = 0) causes a broadband improvement of the oscillatory response of the bubble.cases.In both cases, the pressure amplitude increases almost linearly by increasing the transducer input power, however at each power, the value of the pressure amplitude is greater in the nonlinear acoustic field along a steeper slope.As the input power is increased at a particular frequency ( f = 3 MHz), the amplitudes of all the harmonics increase while considering the nonlinearity parameter (β = 4.1); this leads to the wave distortion and adding more energy available to each harmonic, predominantly in the second and third harmonics (Figure 7b). Figure 7c compares the pressure amplitude at focal point for the two linear (β = 0) and nonlinear (β = 4.1) cases.In both cases, the pressure amplitude increases almost linearly by increasing the transducer input power, however at each power, the value of the pressure amplitude is greater in the nonlinear acoustic field along a steeper slope.

Effects of Acoustic Pressure (Input Power to the Transducer) on Bubble Dynamics
Figure 8 shows the bubble dynamics at P = 6 W, insonified at focal point at the excitation frequency of 3 MHz.Referring to Figure 7c, the pressure amplitude for P = 2 W in the nonlinear case is very close to the pressure amplitude at P = 6 W in the linear field, however comparing the bubble dynamics in Figures 5 and 8 show the differences in the bubble growth, shape modes, and translational instability of the bubble which is evident in the velocity waveform and translational motion.This confirms that the radial growth and the onset of erratic motion of the bubble in a focused ultrasound field are reliant on both pressure amplitude and the harmonic components of the excitation.
input power levels ( 3 8 W P = − ), respectively.As the input power is increased at a particular frequency ( 3 MHz), f = the amplitudes of all the harmonics increase while considering the nonlinearity parameter ( 4.1); this leads to the wave distortion and adding more energy available to each harmonic, predominantly in the second and third harmonics (Figure 7b). Figure 7c compares the pressure amplitude at focal point for the two linear ( 0 β = ) and nonlinear ( β = ) cases.In both cases, the pressure amplitude increases almost linearly by increasing the transducer input power, however at each power, the value of the pressure amplitude is greater in the nonlinear acoustic field along a steeper slope.in the linear field, however comparing the bubble dynamics in Figures 5 and 8 show the differences in the bubble growth, shape modes, and translational instability of the bubble which is evident in the velocity waveform and translational motion.This confirms that the radial growth and the onset of erratic motion of the bubble in a focused ultrasound field are reliant on both pressure amplitude and the harmonic components of the excitation.
The radial growth of the bubble for varying input power (leading to varying excitation acoustic field) is given in Figure 9a,b, for linear and nonlinear acoustic excitation, respectively.In Figure 9a, it is seen that, for all values of power, the bubble radius reaches to the maximum value at the same time.However, in the nonlinear acoustic domain (Figure 9b), as the input power increases, the radial growth of the bubble to the maximum value is delayed.Moreover, based on Figure 9b, the bubble may collapse and show inertial (unstable) cavitation when the input power increases to 7 W and more.Hence, studying nonlinear acoustics in bubble dynamics is essential for detecting the threshold for inertial cavitation in biomedically enhanced delivery of drugs into tissue [65].In Figure 9c, the maximum radius of the bubble corresponds to each input power and acoustic pressure is given.The results show that ignoring nonlinearities of focused ultrasound can yield 56% error in the prediction of the bubble dynamics models.For low amplitude of pressure, we see a slight difference between the linear and nonlinear model which results in small error in the prediction of the maximum radius of the bubble.However, when the applied pressure increases considerably, the The radial growth of the bubble for varying input power (leading to varying excitation acoustic field) is given in Figure 9a,b, for linear and nonlinear acoustic excitation, respectively.In Figure 9a, it is seen that, for all values of power, the bubble radius reaches to the maximum value at the same time.However, in the nonlinear acoustic domain (Figure 9b), as the input power increases, the radial growth of the bubble to the maximum value is delayed.Moreover, based on Figure 9b, the bubble may collapse and inertial (unstable) cavitation when the input power increases to 7 W and more.Hence, studying nonlinear acoustics in bubble dynamics is essential for detecting the threshold for inertial cavitation in biomedically enhanced delivery of drugs into tissue [65].The effects of acoustic nonlinearity on the translational instability of the spherical gas bubble for various pressure fields are investigated in Figure 10, considering both linear and nonlinear excitations.To the best of our knowledge, most of the theoretical studies on the translational motion of a bubble in high-intensity sound wave fields are performed [20,21] while the effects of absorption, diffraction, and nonlinear distortion in the medium are not considered.The studies show that in a high-intensity field, the translational motion of the bubble occurs in more complicated trajectories such that the instability can develop very quickly.Figure 10b shows the evidence of the onset of In Figure 9c, the maximum radius of the bubble corresponds to each input power and acoustic pressure is given.The results show that ignoring nonlinearities of focused ultrasound can yield 56% error in the prediction of the bubble dynamics models.For low amplitude of pressure, we see a slight difference between the linear and nonlinear model which results in small error in the prediction of the maximum radius of the bubble.However, when the applied pressure increases considerably, the results show significant divergences; the linear model predicts a stable cavitation, whereas the nonlinear model reaches the threshold for unstable cavitation.
The effects of acoustic nonlinearity on the translational instability of the spherical gas bubble for various pressure fields are investigated in Figure 10, considering both linear and nonlinear excitations.To the best of our knowledge, most of the theoretical studies on the translational motion of a bubble in high-intensity sound wave fields are performed [20,21] while the effects of absorption, diffraction, and nonlinear distortion in the medium are not considered.The studies show that in a high-intensity field, the translational motion of the bubble occurs in more complicated trajectories such that the instability can develop very quickly.Figure 10b shows the evidence of the onset of unstable movement by the steep change in translational velocity of the bubble.Considering acoustic nonlinearity, one can see that the translational instability of the bubble accelerates and develops more intensely compared to that of linear excitation.Moreover, at higher pressure amplitudes, the shape oscillations develop strongly and reach high amplitudes, which leads to acceleration and amplification of the translational instability and dancing motion that can be seen in Figure 10a,b.The effects of acoustic nonlinearity on the translational instability of the spherical gas bubble for various pressure fields are investigated in Figure 10, considering both linear and nonlinear excitations.To the best of our knowledge, most of the theoretical studies on the translational motion of a bubble in high-intensity sound wave fields are performed [20,21] while the effects of absorption, diffraction, and nonlinear distortion in the medium are not considered.The studies show that in a high-intensity field, the translational motion of the bubble occurs in more complicated trajectories such that the instability can develop very quickly.Figure 10b shows the evidence of the onset of unstable movement by the steep change in translational velocity of the bubble.Considering acoustic nonlinearity, one can see that the translational instability of the bubble accelerates and develops more intensely compared to that of linear excitation.Moreover, at higher pressure amplitudes, the shape

Conclusions
The studying of microbubble's oscillations in size and shape during FU insonation (ultrasound stable or non-inertial cavitation) is of immense significance in biomedical applications, such as enhancing drug delivery rate in the blood-brain barrier opening.Hence, accounting for the FU-induced acoustic nonlinearity in the mathematical models of bubble dynamics is essential to accurately predict the bubble oscillations and acoustic pressure due to streaming.Following the previously established mathematical model [21], this paper aims to investigate the effects of nonlinear propagation of a FU-induced pressure field into the dynamics of ultrasound microbubbles.A theoretical framework (based on the KZK equation) is used to predict the acoustic pressure field due to focused transducers considering the three effects of acoustic nonlinearity, diffraction, and absorption in the medium.Then the nonlinear acoustic model is coupled with the bubble dynamics (based on the Keller-Miksis model) to predict the oscillation, shape modes and translational motion of a single air bubble excited in a FU field at a focal point.The model is used to explore the effects of medium properties (nonlinearity and absorption) and input parameters (power and frequency) on the pressure waveform distortion and amplitude.The observations related to medium properties show that the coefficient of nonlinearity and absorption play significant roles in distorting the waveform and generating more harmonics, thus increasing the energy deposition at the focal point [34].
The KZK-Keller-Miksis model incorporates radial pulsations of a bubble along with its small shape oscillations and translational motion in a FU field.Nonlinear coupling between fundamental radial mode, second to fourth modes of oscillations, and translational motion of the bubble are analyzed to explain the translational instability of the bubble.It is concluded that second-order interactions of two neighboring surface modes impose instability to the system, which in turn leads to dancing motion of the bubble.Furthermore, the analysis show that in addition to resonance response of the bubble, the off-resonance excitation leads to interesting results.Particularly, the acoustic pressure fluctuations (due to acoustic streaming effect [34]) and translational instability of the bubble may increase at above or below fundamental resonance excitation.The translational instability has broad applications in cell sorting which resembles the bubble dancing motion in a microfluidic device [13,66,67].Therefore, studying the dynamics of off-resonant bubbles is essential in acoustic bubble manipulation.In addition, the results reveal the importance of studying nonlinear acoustics and consideration of all the involved harmonics in dynamics of FU-excited microbubbles.A summary of the observations is as follows: 1.
The results show the dynamical motions of the bubble in the FU field generated by the focused transducer at three different range of frequencies (including the resonance and off-resonance cases).Effects of the change of excitation frequency are investigated and it is shown that in all the frequency ranges, the nonlinearity expedites and amplifies the radial growth and translational motion of the bubble where the indication of the translational instability is observed by a steep rise in the translational velocity of the bubble.It is seen that the initiation of the erratic dancing motion of the bubble in the acoustic field is well-matched with the time of instability of other modes.Moreover, acoustic nonlinearity amplifies almost all the shape oscillations and the instability of the modes.Among the three ranges of frequency investigated in this study, the results show that the changes in the dynamical behavior of the bubble and the translational instability due to nonlinear effects are more significant for higher frequencies, which must be taken into account while investigating an erratic translational (dancing) motion of the bubble.However, a momentous effect of nonlinearity on the translational motion of the bubble is seen at resonance frequency where the bubble initially moves away from the pressure antinode and then makes the translation jumps around the antinode.It is seen that as the excitation frequency increases, the inertial (unstable) cavitation occurs earlier in the nonlinear acoustic field; the nonlinear excitation causes more energy deposition at higher harmonics at the focal point (where the bubble is located); therefore, the radius growth of the bubble is more pronounced at higher frequencies.Accordingly, studying nonlinear acoustics in bubble dynamics is essential for inertial cavitation prediction in biomedically enhanced delivery of drugs into tissue [65] as well as the prediction of the induced impact pressure as a shock loading to a structure, which might even cause mechanical destruction [68].

2.
As the input power is increased, the amplitudes of all the harmonics increase.Comparing the bubble dynamics for two equal pressure amplitudes in linear and nonlinear cases show significant differences in the shape modes and translational instability of the bubble.This confirms that the radial growth and onset of the erratic motion of the bubble in a focused ultrasound field are reliant on both pressure amplitude and the harmonic components of the excitation.It is seen that when the applied pressure increases, the results of the linear model give up to 56% error for predicting the threshold for unstable cavitation (maximum radius of the bubble).In nonlinear acoustics, when the input power increases, the bubble may collapse and shows unstable cavitation as well as the amplified and accelerated translational instability.

Figure 1 .
Figure 1.Schematic representation of acoustic waves from a focused source to a bubble at the focal point.
s indicates that Equation (6) is carried out in an isentropic process.Consequently, the nonlinearity parameter / B A is written as

Figure 1 .
Figure 1.Schematic representation of acoustic waves from a focused source to a bubble at the focal point.

c
 represents the phase velocity.

where 0 c
is the sound velocity at reference frequency 0  and () c  represents the phase velocity.
) takes into account the local generalized ultrasonic attenuation-dis a result, the linear part of Equation (9), written as

Figure 2 .
Figure 2. Acoustic pressure generated by the focused transducer as the input parameter to bubble dynamics at 2 W P = (a) for various excitation frequency; temporal waveform and (b) the corresponding Fourier spectra for 4.1 β = .

Figure 2 .
Figure 2. Acoustic pressure generated by the focused transducer as the input parameter to bubble dynamics at P = 2 W (a) for various excitation frequency; temporal waveform and (b) the corresponding Fourier spectra for β = 4.1.

Figure 6 .
Figure 6.Fourier spectra of the radial oscillation of the bubble in blood domain shown in Figure A1a.The curve is on a log scale.

Figure 7a ,
Figure 7a,b show the acoustic pressure at focal point in time and frequency domains for various input power levels ( 3 8 W P = − ), respectively.As the input power is increased at a particular frequency ( 3 MHz), f = the amplitudes of all the harmonics increase while considering the nonlinearity parameter ( 4.1); β = Figure 7c compares the pressure amplitude at focal point for the two linear ( 0 β = ) and nonlinear (

Figure 6 .
Figure 6.Fourier spectra of the radial oscillation of the bubble in blood domain shown in Figure A1a.The curve is on a log scale.

Figure 7 .
Figure 7. Acoustic pressure generated by the focused transducer as the input parameter to bubble dynamics at 3MHz f = (a) for various input power; temporal waveform and (b) the corresponding Fourier spectra for 4.1 β = .(c) Pressure amplitude vs. transducer input power in linear and nonlinear

Figure 7 .
Figure 7. Acoustic pressure generated by the focused transducer as the input parameter to bubble dynamics at f = 3 MHz (a) for various input power; temporal waveform and (b) the corresponding Fourier spectra for β = 4.1.(c) Pressure amplitude vs. transducer input power in linear and nonlinear cases.

Figures 8
Figures 8 shows the bubble dynamics at 6 W P = , insonified at focal point at the excitation frequency of 3 MHz .Referring to Figure 7c, the pressure amplitude for 2 W P = in the nonlinear case is very close to the pressure amplitude at 6 W P = in the linear field, however comparing the

Acoustics 2019, 2 FOR PEER REVIEW 16 Figure 9 .
Figure 9. Radius versus time plot of a radius 7 μm − − bubble at 3 MHz f = for various input power in (a) linear and (b) nonlinear acoustic field; (c) A comparison of the linear (dashed line) and nonlinear (solid line) model predictions for varying input power (black line) and acoustic pressure (red line).

Figure 10 .
Figure 10.Dynamics of translation of a radius 7 μm − − bubble at 3 MHz f = for various input power: (a) Translational motion and (b) translational velocity (the corresponding acoustic pressures are given in Figure 7).

Figure 9 .
Figure 9. Radius versus time plot of a 7 − µm − radius bubble at f = 3 MHz for various input power in (a) linear and (b) nonlinear acoustic field; (c) A comparison of the linear (dashed line) and nonlinear (solid line) model predictions for varying input power (black line) and acoustic pressure (red line).

Figure 9 .
Figure 9. Radius versus time plot of a radius 7 μm − − bubble at 3 MHz f = for various input power in (a) linear and (b) nonlinear acoustic field; (c) A comparison of the linear (dashed line) and nonlinear (solid line) model predictions for varying input power (black line) and acoustic pressure (red line).

Figure 10 .
Figure 10.Dynamics of translation of a radius 7 μm − − bubble at 3 MHz f = for various input power: (a) Translational motion and (b) translational velocity (the corresponding acoustic pressures are given in Figure7).

Figure 10 .
Figure 10.Dynamics of translation of a 7 − µm − radius bubble at f = 3 MHz for various input power: (a) Translational motion and (b) translational velocity (the corresponding acoustic pressures are given in Figure 7).

Table 1 .
Parameters of the blood medium used in the Khokhlov-Zabolotskaya-Kuznetsov (KZK) model.