Weakly nonlinear focused ultrasound in viscoelastic media containing multiple bubbles

To facilitate practical medical applications such as cancer treatment utilizing focused ultrasound and bubbles, a mathematical model that can describe the soft viscoelasticity of human body, the nonlinear propagation of focused ultrasound, and the nonlinear oscillations of multiple bubbles is theoretically derived and numerically solved. The Zener viscoelastic model and Keller–Miksis bubble equation, which have been used for analyses of single or few bubbles in viscoelastic liquid, are used to model the liquid containing multiple bubbles. From the theoretical analysis based on the perturbation expansion with the multiple-scales method, the Khokhlov–Zabolotskaya–Kuznetsov (KZK) equation, which has been used as a mathematical model of weakly nonlinear propagation in single phase liquid, is extended to viscoelastic liquid containing multiple bubbles. The results show that liquid elasticity decreases the magnitudes of the nonlinearity, dissipation, and dispersion of ultrasound and increases the phase velocity of the ultrasound and linear natural frequency of the bubble oscillation. From the numerical calculation of resultant KZK equation, the spatial distribution of the liquid pressure fluctuation for the focused ultrasound is obtained for cases in which the liquid is water or liver tissue. In addition, frequency analysis is carried out using the fast Fourier transform, and the generation of higher harmonic components is compared for water and liver tissue. The elasticity supresses the generation of higher harmonic components and promotes the remnant of the fundamental frequency components. This indicates that the elasticity of liquid suppresses shock wave formation in practical applications.

For these medical applications, the investigation of bubble dynamics in soft viscoelastic media like human body is essential. Newtonian fluid is the most famous model in the field of fluid dynamics, but only viscosity is considered. The Kelvin-Voigt [32][33][34][35][36]38,19,37,39,30,31] and Zener models [41,42,45,46,43,44,47] have been used as mathematical models of liquid with viscoelasticity. The Kelvin-Voigt model corresponds to the generalized Newtonian fluid with the elasticity introduced while the relaxation effect is omitted. Yang and Church [34] investigated a behavior of single spherical bubble based on the Kelvin-Voigt model and Keller-Miksis equation [34,52,41,35,42,36,38,19,43,37,44,39,40], and showed that liquid elasticity decreases the nonlinearity of bubble oscillation. Keller-Miksis equation is a model of bubble dynamics in compressible and viscoelastic liquid and often used in combination with Kelvin-Voigt model [34][35][36][38][39][40]. Murakami et al. [37] developed a model of the behavior of a non-spherical single bubble based on the Kelvin-Voigt model and Rayleigh-Plesset equation and investigated the shape stability of bubbles. Rayleigh-Plesset equation is the most famous model of bubble dynamics, but is the limiting case of Keller-Miksis equation assuming the incompressible liquid. Qin et al. [39] investigated the effects of the liquid viscoelasticity and the interaction of two bubbles on behavior of each bubble using the Kelvin-Voigt model and Keller-Miksis equation. The Maxwell viscoelastic model [48][49][50][51] is also the generalized Newtonian fluid with the relaxation effect of viecoelastic medium, however the elastic term is omitted. Fogler and Goddard [48] investigated the collapse of bububle in viscoelastic media based on the Maxwell model and Rayleigh-Plesset equation.
Further, the Zener model corresponds to the generalized Kelvin-Voigt and Maxwell model, including the elastic term and the relaxation effect of a viscoelastic medium. Warnez and Johnsen [42] developed a numerical method for bubble behavior based on the Zener model and Keller-Miksis equation, and they found that incorporating the relaxation time increases bubble growth. Zilonova et al. [46] investigated the interaction of two bubbles with the effects of drag force and translation based on the Zener model and Gilmore-Akulichev equation [53]. Gilmore-Akulichev equation has broader applications than the Keller-Miksis equation [46,47]. Filonets and Solovchuk [47] investigated the single bubble behavior excited by a dual-frequency ultrasound based on the Zener model and Gilmore-Akulichev equation, and they showed that the elasticity of a liquid increases the threshold pressure of bubble collapse.
The authors were involved in the first derivation of the KZK equation for a liquid containing multiple bubbles [88][89][90], using theoretical analysis based on volumetric averaged basic equations for liquids containing multiple bubbles. Kagami and Kanagawa [10] introduced the thermal effects of gas inside the bubbles into the KZK equation and obtained numerical solutions. However, these models [88][89][90]10] only considered viscosity and neglected the elasticity of the liquid. Recently, Hasegawa et al. [91] succeeded in introducing the elasticity of a liquid into the weakly nonlinear wave equation; however, this model was limited to the spatial one-dimensional form and could not be used for focused ultrasound. In addition, the viscoelasticity of the bubble-liquid interface was considered, but not the viscoelasticity of the entire liquid. In summary, a mathematical model that can describe the liquid viscoelasticity of both the bubble-liquid interface and the entire liquid, the nonlinear propagation of focused ultrasound, and the oscillations of multiple bubbles would be very useful.
In this study, the KZK equation describing the weakly nonlinear propagation of focused ultrasound in a viscoelastic liquid containing multiple bubbles is derived. In Section 2, volumetric averaged equations for a liquid containing multiple bubbles are introduced, based on the mixture model [96][97][98][99]. The viscoelasticities of the entire liquid and bubble-liquid interface are introduced, as shown in Fig. 1. The momentum conservation equation in the mixture model is combined with the Zener model to incorporate the viscoelasticity of the entire liquid, and the Keller-Miksis equation for bubble dynamics is used to incorporate the viscoelasticity of the bubble-liquid interface. Theoretical analysis based on the perturbation expansion and multiple-scales method [121] is carried out in Section 3, and the derivation of the KZK equation is demonstrated. This KZK equation is composed of terms representing the nonlinear, dissipation, dispersion, and diffraction effects of ultrasound propagation. The effects of the liquid elasticities on the nonlinearity, dissipation, dispersion, and phase velocity of focused ultrasound and the natural frequency of bubble oscillation are investigated by comparing cases in which the liquid is water (without elasticities) or liver tissue (with elasticities). By a comparison among the Zener model, the Maxwell model, and the Kelvin-Voigt model, the effects of the rigidity and the relaxation time are compared. In addition, a comparison of the liquid viscoelasticity of the entire liquid and bubble-liquid interface is conducted. In Section 4, the numerical solution of the newly obtained KZK equation is presented. The spatial distribution of the fluctuation of liquid pressure for the focused ultrasound is obtained. The dispersion effect on the numarical solution appeared. In addition, frequency analysis is conducted using the fast Fourier transform (FFT), and the generation of higher harmonic components is compared for water and liver tissue. As a result, the elasticity of liver tissue supresses the generation of higher harmonic components and promotes the remnant of the fundamental frequency components compare to the case of water, although the maximum value of the rise of liquid pressure are quite close. This result implies that the elasticity of liquid suppresses shock wave formation in practical applications.

Problem statement
The weakly nonlinear propagation of focused ultrasound in a liquid containing multiple bubbles is considered. In this study, the viscoelasticity of the liquid phase is considered in two regions: the entire liquid and the bubble-liquid interface (see Fig. 1).
Focused ultrasound is radiated from a sound source in a viscoelastic medium containing multiple bubbles (Fig. 2). The center of the sound source is set as the origin. The x * -axis represents the distance from sound  source, and the r * -, y * -, and z * -axes denote the distances from the x * -axis. The following relation is provided: Note that the surface shape of the sound source is not only focusing, but can also be in a planar or spreading form, and the shape profile of the sound source is not only circular, but can also be rectangular or elliptical.
The KZK equation is derived by assuming quasi-plane propagation, i. e., weakly diffracted (focused or spreading) waves. In this study, two types of KZK equations are derived: two-dimensional (2D) and threedimensional (3D) spatial forms, as in our previous work [10]. Although the 2D KZK equation is assumed for axial symmetric ultrasound, the 3D KZK equation can handle asymmetric propagation, such as ultrasound radiated by an elliptical or rectangular sound source or propagated in a nonuniform medium, such as the human body. The following relation is provided as a mathematical formula of the quasiplane waves [88,90,93,10].
where R * 0 is the initial bubble radius, L * is the typical ultrasound wavelength, and a * i (i = r,y,z) are typical lengths in each direction. a * r is set as the radius of circular sound source, and a * y and a * z are set as the long and short diameters of elliptic sound sources, or the width and height of rectangular sound sources, etc.
In the theoretical analysis, the following assumptions are introduced for simplicity.
(a-2)The viscosity and elasticity of the gas phase are neglected because they are significantly smaller than those of the liquid phase. (a-3)The basic equations based on the mixture model [96][97][98][99] of bubbly liquids are used; [88,94,95,93] those based on the two-fluid model [99,100] were used in our previous work [88][89][90]10] because the two-fluid model needs many assumptions and experimental laws to incorporate the viscoelasticity of the liquid phase. Note that the mixture model requires some assumptions such as (a-4) and (a-6) described below. (a-4)The difference of velocity between the gas and liquid phases is sufficiently small, and then the drag force [102,103] is not considered.
(a-5)The initial void fraction is sufficiently small. (a-6)The viscosity and elasticity of the gas-liquid mixture are modeled by only those of the liquid phase, because the viscosity and elasticity of the gas phase are sufficiently smaller than those of the liquid phase according to Assumptions (a-2) and (a-5).
(a-7)The bubble-bubble interaction [104][105][106]46,44,39], which may be dominant under high void fraction, is not considered. (a-8)Bubble oscillation is spherically symmetric. (a-9)Bubbles do not collapse, appear, or coalesce. (a-10)The phase change [107][108][109] and the heat transfer [110] at the bubble-liquid interface is not considered. (a-11)The distributions of the pressure and temperature of the gas inside the bubbles [101] are not considered; they are treated as the averages inside each bubble. (a-12)The temperature of the liquid phase is assumed constant, whereas the temperature of the gas phase is assumed to fluctuate. (a-13)At initial state, the liquid is at rest and spatially uniform, except for the spatial distribution of bubbles. The nonuniformity of the bubble size [111][112][113][114][115] is not considered. The spatial nonuniformity of the bubble distribution at the initial state is introduced [88,90,93] because bubbles are only used in focus in medical applications.

Basic equations
The conservation laws of mass and momentum based on the mixture model [96][97][98][99] are used: where t * is time, ρ * is density, u * is the velocity vector, d * is the displacement vector, p * is pressure, λ * v,L is the second viscosity constant, λ * e,L is the Lamé constant, and τ * is the deviatoric stress tensor. The subscripts G and L mean the values of the gas and liquid phases, respectively. In Eq. (4), the pressure of the mixture is assumed to be that of the liquid phase [89,94,93]. The density of the gas phase is generally sufficiently small; hence, the density of the mixture is defined by the liquid density as follows: where α is the void fraction where R * is the bubble radius, n * is the bubble number density, and the following conservation law is imposed: The conservation of mass in each bubble is The second viscosity λ * v,L and Lamé's constant λ * e,L can be rewritten as follows: where K * v,L and K * e,L are the bulk viscosity and elasticity, respectively; μ * L is the viscosity; and G * L is the rigidity (shear modulus). The bulk viscosity K * v,L was set to zero based on the Stokes assumption [94,93]. The physical quantities in Eqs. (9) and (10) were originally that of the gas-liquid mixture; however, they have been substituted by that of the liquid phase from Assumption (a-6) in Section 2.1.
The subscript 0 represents initial values. Note that the Keller-Miksis equation can be used for cases in which the Mach number is less than one [46,47], and the Gilmore-Akulichev model [46,47] can be used with larger Mach numbers. The variable q * L is given by where r * ′ is the distance from the center of each bubble. Note that r * ′ is based on the spherical coordinate of each bubble and differs from r * , which is based on the macroscopic coordinate as shown in Fig. 2. The balance of the normal stresses across the bubble-liquid interface is To obtain the constitutive equation for q * L , the Zener Eq. (16) is integrated from the bubble-liquid interface to infinity: [41,42,47,44,38,39] From Eqs. (17)- (20), the effect of liquid viscoelasticity at the bubble-liquid interface is introduced. The energy equation describing the thermal conduction at the bubble-liquid interface [116] is used: [94,93,95,10] where κ is the ratio of the specific heat of the gas phase, λ * is the thermal conductivity of the gas phase, and T * is the temperature. The temperature gradient ∂T * G /∂r * ′ | R * is substituted into the following model [117], which is used for analysis.
Note that many temperature gradient models have been proposed, such as Refs. 117-120; however, the model in Eq. (22) [117] derived the most appropriate numerical solution of the models in our previous work [10] and is used again in this study. To close the set of equations, the Tait equation of state for the liquid phase and the equation of state for an ideal gas are used: where n is the material constant (e.g., n = 7.15 for water).

Parameter scaling
The following scale parameters are introduced to consider the low frequency and long wave as in our previous studies [92,88,93,94,89,10], using the dimensionless amplitude ∊ (≪1): where ω * ,L * , and U * are the typical angular frequency, wavelength, and propagation speed, respectively, and Ω, Δ, and V are the dimensionless constants of O(1).
The assumption of weakly diffracted waves given by Eq.
(2) is rewritten as follows (i = r, y, z): where the dimensionless constant Γ i of O(1) represents the degree of diffraction for each direction.
The liquid viscosity μ * L , rigidity G * L , bulk elasticity K * v,L , and relaxation time λ * relax,L are nondimensionalized as where T * is the typical period of the wave, which can be related to the wavelength L * and wave speed U * as U * = L * /T * . The nondimensionalization of the energy Eq. (22) is as in Refs. 94,10: where λ G is the dimensionless constant of O(1).

Multiple scale analysis
The independent variables t * ,x * ,r * ,y * , and z * are nondimensionalized as Next, t and x are extended to the near and far fields using the dimensionless wave amplitude ∊: [121] By the assumption of weakly diffracted waves in Eq. (2), the dependence of unknown variables on the radial direction (r,y, z) is smaller than that on the propagative direction x [88,90,93,89]. Then, the spatial coordinates r, y, and z are defined in the far field only as All the unknown variables can be regarded as functions of the extended independent variables of Eqs. (35)- (38). The derivative operators are expanded as follows: [121] ∂ ∂t ∂ ∂x ∂ ∂r The unknown variables are expanded as a power series of ∊: where u * x ,u * r ,u * y , and u * z are the components of velocity vector u * in the x * , r * ,y * , and z * directions, respectively. The magnitudes of u * r ,u * y , and u * z are assumed to be smaller than that of the u * x direction [88,90,93,89,10]. Then, the expansions of u * r ,u * y , and u * z begin with a higher order than that of u * x in Eqs. (46)- (49). The components of displacement vector d * are expanded in the same manner as u * : The components of the deviatoric stress tensor τ * are expanded as where other components that do not affect the result are omitted. The gas density is assumed to be significantly smaller than the liquid density in the initial state: ) .
The dimensionless liquid pressure is defined as The liquid density is expanded as [92] ρ * where the first order of the expansion is determined using Eqs. (23), (27), and (44).
To incorporate the weak nonuniformity of the spatial distribution of bubbles in the initial state, void fraction α is expanded as [88,90,93 ,89,10] α α 0 where δ is the known variable that represents the spatial nonuniformity of the void fraction in the initial state. The effect of the spatial nonuniformity of the void fraction is assumed to appear only in the far field; then, δ is the function of x 1 only [88,90,93,89,10].

Results of theoretical analysis
The scale parameters of Eqs.

Approximation of first order
The velocity and displacement are related by Eq. (14) as From the approximation of the Keller-Miksis Eq. (17), the linear natural frequency of bubble oscillation ω * B is obtained: In Eq. (69), the effect of liquid elasticity at the bubble-liquid interface [91] is introduced to our previous works [88,94,95,93]. Note that ω * B of Eq. (69) is obtained from the linear approximation, then the actual frequency of bubble oscillation will be shifted by accumulation of nonlinearity and wave amplitude [122,123]. In addition, bubble-bubble interaction [104][105][106]46,44,39] and dual frequency ultrasound [129][130][131] also affect the frequency of bubble oscillation, however these are not considered in this study.
Eqs. (62)-(67) are combined into a single equation of R 1 as where the phase velocity v P is given by For simplicity, v P ≡ 1 is imposed and the typical wave speed U * is given by In Eq. (72), the effect of the elasticity of the entire liquid is newly introduced to our previous works [88,94,95,93]. Next, the variable transformation is introduced as Then, the equation describing the right-running waves is obtained.
The variable transformation of Eq. (73) is introduced into approximated Eqs. (62)-(68)), and all unknown variables are written in terms of R 1 = f as follows. where

Approximation of radial direction
From the approximations of O(∊ 2/3 ), the radial components of the momentum conservation Eq. (13) are (1 − α 0 ) ∂u z1 ∂t 0 The components of the deviatoric stress tensor (16) related to the radial directions are The velocity and displacement are related by Eq. (14) as The variable transformation Eq. (73) is introduced into Eqs. (82)- (88): where the forms of (89)-(91) are the same as in our previous work [93].

Approximation of second order
From the approximations of O(∊ 2 ), the following equations are obtained for Eqs. (11)-(13) as ∂T G2 ∂t 0 The explicit forms of the inhomogeneous terms K i (i = 1, …, 6) are shown in Appendix A. Eqs. (92)-(97) are combined into a single equation: The inhomogeneous term K is given by As the solvable condition for Eq. (98), K = 0 is imposed [88,90,92,89], and the following relation is obtained for 2D and 3D spatial cases: Finally, the KZK equation is obtained: where Δ ⊥ is the Laplacian operator given for 2D and 3D spatial cases by The following variable transformations are used: where τ is the retarded time. In the KZK Eq. (102), the right-hand side with the Laplacian operator represents the diffraction (focusing) effect.

Coefficients of KZK equation
Π 01 , Π 02 , and Π 4 are the advection coefficients given by By the effect of the elasticity of the entire liquid, the new coefficient Π 02 is introduced in this study and terms are added to Π 01 and Π 4 from our previous work [10]. The spatial nonuniformity of the initial bubble distribution δ(x 1 ) only appears in the variable transformation of Eq.
(104), in which it only affects the advection term [88,90,93,89]. Π 1 is the nonlinear coefficient given by where Π 21 is the dissipation coefficient divided into five terms for each factor and region: Π 22 and Π 3 are the dissipation and dispersion coefficients, respectively, given by Whereas Π 21 depends on the compressibility, viscosity, and elasticity of the liquid phase, Π 22 depends on the thermal conductivity of the gas phase. Fig. 3 shows the coefficients for cases in which the liquid phase is water (black) and liver tissue (red). The properties of the water and liver tissue used in the calculations are shown in Table 1. Unless otherwise stated, these values are used in all calculations in Sections 3.5, 3.6, 3.7, and 4. As shown in Fig. 3, elasticity decreases the magnitudes of the nonlinear coefficient Π 1 , dissipation coefficients Π 21 and Π 22 , and dispersion coefficient Π 3 , and increases the phase velocity U * and linear natural frequency of bubble oscillation ω * B . Decrease of nonlinearity by the liquid elasticity is qualitatively consistent with the previous result [34] focusing of single bubble behavior.

Comparison of three viscoelastic models
In Fig. 4

Comparison of elasticity of entire liquid and bubble-liquid interface
In Fig. 5, the effects of the liquid elasticity of the entire liquid and bubble-liquid interface are compared. As the effects of each factor become larger, the red (elasticity of entire liquid) and blue (elasticity of bubble-liquid interface) lines differ from the black line (both elasticities not considered). (See Fig. 6).
As shown in Fig. 5, the nonlinear coefficient Π 1 is decreased by both elasticities; the effect of the bubble-liquid interface is large at the wide range of the initial void fraction α 0 , while the effect of the entire liquid gradually increases as the value of α 0 increases. The magnitudes of the dissipation coefficients Π 21 and Π 22 and the dispersion coefficient Π 3 are affected by the elasticity of the bubble-liquid interface, whereas the effect of the elasticity of the entire liquid is only slight and increases with the initial void fraction α 0 . The phase velocity U * and linear natural frequency of the bubble oscillation ω * B increase when the elasticity of the bubble-liquid interface is considered, whereas the effects of the elasticity of the entire liquid are not observed.

Limitation of KZK Eq. (102)
The resultant KZK Eq. (102) includes the nonlinear term with coefficient Π 1 owing to nonlinearity of ultrasound propagation, while the dissipation and dispersion terms are limited to linear form. However, effects of nonlinear dissipation [124][125][126][127] will become large particularly in case of high amplitude. Two main methods to incorporate the nonlinear dissipation are (i) extending to over third order analysis from the second order of this study, (ii) nondimensionalization with lower order of ∊ (≪1) for liquid viscosity of Eq. (29). These extensions will significantly change the framework of theoretical analysis and the resultant equation will be no longer the form of original KZK equation, Table 1 Material properties of liquid phase.

Water
Liver tissue [47,128] Rigidity  but be very important for practial applications then be presented in our future work.

Method
The KZK equation of the spatial 2D form of Eq. (102) is numerically solved. As in our previous work [10], the finite-difference time-domain scheme [58,57] is used, which has been widely used to simulate focused ultrasound in single phase liquid [57][58][59][60][61][62][63][64][65][66][67][68]. To solve the fluctuation of liquid pressure p L1 , the KZK Eq. (102) is rewritten using the relation of the fluctuation of bubble radius f and liquid pressure p L1 in Eq. (75): where the focusing gain G r is given by The variable transformations (104) are rewritten again as The step sizes of each direction are Δτ = 2π/240, ΔX = 0.00025, and ΔR = 0.01, and the calculation regions are τ min ≦τ≦τ max , 0≦X≦X max , and 0≦R≦R max . To suppress the numerical oscillation, the calculation regions need to be sufficiently large for the direction of τ and R; then, τ min = − (G r + 24π), τ max = 34π, and R max = 400 are set.
Next, the KZK Eq. (125) is split into three parts and solved term by term: The diffraction term of Eq. (130)  In addition, the error caused by separately introducing the diffraction, dissipation, dispersion, and nonlinear effects is of order (ΔX) 2 . The total error is estimated as ΔX +(Δτ) 2 +ΔR [58,57].
The boundary condition at X = 0 for the focused sound source is given by where p * 0 is the amplitude of the source pressure. The wave is only given on the sound source, 0≦R≦1. As shown in Fig. 5, the wave at the sound source (x * = 0 mm and r * = 0 mm) is amplitude-modulated by the form of the exponential functuion as in Eq. (133). In Eq. (133), T d and m are the effective duration and the rise time of pulses, respecrively [57,58,60], and set as T d = 14π and m = 4 to radiate approximately 10 pulses.

Results
In Figs

Generation of higher harmonic
The dispersion effect is represented in the form of the third derivative with coefficient Π 3 in the resultant KZK Eq. (102) and obtained by introducing the bubble oscillation [10,93,25,26,28,29,27,115,88,114,9 5,89-92,94,102,103,113]. By the dispersion effect, the propagation speed of waves shifts depending on the frequency of each wave. Hence when the dispersion effect is indroduced, higher harmonic components generated by the nonlinear effect will propagate at a speed different from the fundamental frequency component. Fig. 9 shows the numerical solution that only the dispersion coefficient Π 3 is virtually set as zero to clarify the dispersion effect. In Figs. 7b and 8b with the dispersion effect, some peak points with low values newly appeared among the discontinuous points owing to the dispersion effect. In contrast to Figs. 7b and 8b, the time development of Fig. 9b  does not have the new peak points. In addition, tha spatial distribution of Fig. 9a without dispersion effect has the narrower peak region of pressure rise than the Figs. 7a and 8a.

Effect of elasticity on higher harmonic generation
In Fig. 7, the initial void fraction α 0 = 0.00025 and the liquid is water. The frequency domains of Figs. 7f show that the fundamental frequency component gradually shift into higher harmonic components. In x * = 29 mm and x * = 32 mm of Figs. 7e and 7f, the third harmonic components ecxeed the second harmonic components.
In Fig. 8, the initial void fraction α 0 = 0.0005 and the liquid is liver tissue. The other conditions including the frequency, the pressure amplitude and the radius of sound source, are same as Fig. 7. The maximum values of pressure fluctuations are 506 kPa and 505 kPa in Figs. 7 and 8, respectively, these values are quite close. As Fig. 7, the frequency domains of Figs. 8f show that the fundamental frequency component gradually shift into higher harmonic components, however, third harmonic components do not exceed the second harmonic components in contrast to Fig. 7. In addition, Figs. 8f retain more fundamental frequency components than Figs. 7f. These results imply that the elasticity of the liquid phase supresses the generation of higher harmonic components and promotes the remnant of the fundamental frequency components; therefore, the elasticity of the liquid phase suppresses shock wave formation in practical applications.

Conclusion
Weakly nonlinear propagation of focused ultrasound in viscoelastic liquid containing multiple bubbles was investigated using the volumetric averaged equations of liquid containing multiple bubbles based on the mixture model [96][97][98][99]. The viscoelasticity of the entire liquid was introduced to the momentum conservative equation in Eq. (4) and modeled with the Zener model in Eq. (16), whereas that of the bubble-liquid interface was considered in the Keller-Miksis equation in Eq. (17).
As a result of the theoretical analysis based on perturbation expansion with the multiple-scales method [121], the KZK equation in Eq. (102) describing the weakly nonlinear propagation in viscoelastic liquids containing multiple bubbles was derived. The resultant KZK Eq. (102) is composed of terms representing nonlinear, dissipation, dispersion, and diffraction effects of ultrasound propagation. As shown in Fig. 3, the liquid elasticity decreases the magnitudes of nonlinear coefficient Π 1 , dissipation coefficients Π 21 and Π 22 , and dispersion coefficient Π 3 , and increases the phase velocity U * and linear natural frequency of bubble oscillation ω * B . As shown in Fig. 4, the comparison among the viscoelastic models; Zener model, Maxwell model, and Kelvin-Voigt model, is conducted. Then, the effects of the rigidity is quite larger than the relaxation time. In addition, a comparison of the liquid viscoelasticity of the entire liquid and bubble-liquid interface was conducted. As shown in Fig. 5, the nonlinear coefficient Π 1 was

Table 2
Parameters and resultant values for five different initial void fractions α 0 , including the cases in Figs. 7 and 8. The absolute values of Π 21 and Π 3 are shown because they are originally negative. Maximum value of p * L1 is at τ * = 10 μs. α0 = 0.0001 (Fig. 7) α0 = 0.00025 (Fig. 8 decreased by both elasticities. However, the magnitudes of the dissipation coefficients Π 21 and Π 22 , dispersion coefficient Π 3 , phase velocity U * , and linear natural frequency of the bubble oscillation ω * B were strongly affected by the elasticity of the bubble-liquid interface, while the effects of the elasticity of the entire liquid were considerably small.
In Section 4, the numerical solution of the newly obtained KZK equation was shown for different cases of the initial void fraction α 0 . The dispersion effect introduced by the bubble oscillation are shown, by virtually setting the dispersion coefficient Π 3 = 0 of Fig. 9. In addition, a frequency analysis was carried out using FFT and the generation of higher harmonic components was compared for water and liver tissue. As shown in the frequency analysis of Figs. 7-8f, the elasticity of liver tissue supresses the generation of higher harmonic components and promotes the remnant of the fundamental frequency components compare to the case of water, although the maximum values of p * L1 are quite close. Therefore, this result implies that the elasticity of the liquid phase suppresses shock wave formation in practical applications.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  Fig. 8.