Weakly nonlinear propagation of focused ultrasound in bubbly liquids with a thermal effect: Derivation of two cases of Khokolov–Zabolotskaya–Kuznetsoz equations

Highlights • Derivation of physico-mathematical simplified model for nonlinear focused ultrasound in liquids containing many bubbles.• Consistent incorporation of ultrasound propagation, bubble oscillation, and temperature fluctuation into single equation.• Spatially two- and three-dimensional cases.• Moderate temperature rise via a numerical analysis.


Introduction
Focused ultrasound is widely used in medical applications for both diagnosis and treatment of various conditions [1]. Ultrasound imaging [2] is a diagnosis method that can be performed in real time. Highintensity focused ultrasound (HIFU) treatment [3][4][5] is a lowinvasive treatment method that is utilized for tumor ablation therapy and shock wave lithotripsy (SWL). Both methods do not present a risk of radiation.
In recent reports, the use of bubbles has been shown to improve the efficiency of medical applications. For ultrasound imaging, microbubbles can be used as ultrasound contrast agents to significantly improve the resolution of images [43][44][45]. For tumor ablation therapy by HIFU, ciency of temperature rise [46][47][48][49][50][51]. Kaneko et al. [48] carried out an experiment using rabbit liver and reported that the use of bubbles approximately doubled the ablated volume and increased the temperature rise by approximately 7°C during 60 s of HIFU radiation. Chang et al. [51] carried out an experiment using a tissue-mimicking phantom and reported that the use of bubbles increased the ablated volume approximately sixfold during 60 s of HIFU radiation.
A physico-mathematical model for medical applications utilizing bubbles should describe the nonlinearity of both ultrasound propagation and bubble oscillation. Vanhille [30] focused on a medium containing a bubble cloud only in the pre-focal zone and proposed a method of solving the KZK equation for the area without bubbles and solving the linear wave equation and bubble dynamics equation (Rayleigh-Plesset equation) separately for the area near the bubble clouds. However, for the area near the bubble clouds, the nonlinearity of bubble oscillation was considered but that of ultrasound propagation was not. Hence, extending the KZK equation from pure liquid to liquid containing bubbles can reduce the computational cost and consistently describe the nonlinear effects of both ultrasound propagation and bubble oscillation. Khismatullin and Akhatov [52] derived a physico-mathematical model similar to the KZK equation for a liquid containing many bubbles, but dissipation effects were not considered (i.e., the Kadomtsev-Petviashvili equation). Our previous study [53,54] was the first to succeed in deriving a generalized KZK equation (i.e., the KZK equation including a dispersion term as third-order derivative due to bubble oscillations) with dissipation effects by the viscosity at bubble-liquid interface and the compressibility of the liquid phase, based on the basic equations for a liquid containing many spherical bubbles. However, temperature fluctuation and thermal conduction of the gas inside the bubbles was not accounted for. In some cases of tumor ablation therapy by HIFU, the temperature of the area near the focus can rise by over 80°C [5], making a description of temperature fluctuation and thermal effects necessary.
Recently, we proposed another generalized KZK equation describing the temperature fluctuation and thermal conduction of the gas inside bubbles [55,56] by introducing the energy equation for the gas inside bubbles [57]. However, this study [55,56] has some theoretical limitations: (i) Only the 2D KZK equation was derived because axisymmetric propagation was assumed. (ii) Volumetric averaged equations based on a mixture model [58] for bubbly liquids were used. However, it is reported [59] that a two-fluid model [60] can describe the dependence of the initial condition of the medium more accurately than the mixture model. (iii) The temperature gradient term of the energy equation [57] needs to be rewritten. Although many temperature gradient models have been proposed [61][62][63][64], only one model [62] was used in that study [55,56], and hence there is no knowledge of which model is the most accurate for modelling medical applications. In addition to these points, numerical solutions of the KZK equation for bubbly liquids have yet to be obtained.
The purpose of this study is the derivation of 2D and 3D generalized KZK equations for bubbly liquids using the volumetric averaged equations based on the two-fluid model [60], with a numerical example of the newly obtained equations provided. This paper is organized as follows: Section 2 introduces the basic equations; including the volumetric averaged equations based on a two-fluid model [60], the energy equation for the gas inside a bubble [57], and the temperature gradient models [61][62][63][64]. The perturbation expansions based on the multiple-scales method are demonstrated [65]. Section 3 presents the method for the derivation of the 2D and 3D KZK equations. The resultant 2D and 3D KZK equations are represented as the linear combination of nonlinear, dissipation, dispersion, and diffraction terms. The dissipation term is divided into three factors; the interfacial viscosity, liquid compressibility, and thermal conductivity of gas. Section 4 includes a numerical example of the newly obtained equations for the spatial distribution of temperature fluctuations. The finitedifference time-domain (FDTD) method developed by Lee and Hamil-ton [22,23] that is widely used as a numerical method for solving the KZK equation for pure liquids is extended to bubbly liquids in this study. The time development of the dissipation term is shown for the interfacial viscosity, liquid compressibility, and thermal conductivity of the gas. The calculation is carried out for three types of gas inside bubbles: Argon, air, and sulfur hexafluoride (SF 6 ). The results show that Argon and SF 6 gas are most effective for HIFU treatment and ultrasound imaging, respectively.

Problem statement
Weakly nonlinear propagation of focused ultrasound in a compressible liquid containing many spherical bubbles is considered in this study. The sound beam was radiated from a sound source in the bubbly liquid (Fig. 1). The surface shape of the sound source is not restricted to planar form, as it may be concave or convex; the former yields a focused beam and the latter a spreading beam. As shown in Fig. 1, the origin is set at the center of the sound source, the x Ã axis is normal to the surface of the sound source, the r Ã , y Ã , and z Ã axes are the distances from the x Ã axis. This provides the equation of The 2D and 3D KZK equations are derived by assuming the quasi-plane waves, i.e., weakly diffracted (weakly focused or weakly spreading) waves. As the assumption of quasi-plane waves for the derivation of 2D KZK equation, we assume that the typical wavelength L Ã is significantly longer than the initial bubble radius R Ã 0 , and the radius of the circular sound source a Ã r is significantly longer than L Ã , as in our previous work [53][54][55][56].
Although the diameter of the sound source was used in our previous work [53][54][55][56], the radius a Ã r is used in this study for applicability to the numerical calculation, but this does not significantly affect the result [53][54][55][56]. In this study, we extend the theory to the 3D KZK equation, where the profile of the sound source is not restricted to a circular form and may be elliptical or rectangular, hence, we introduce the assumption of where a Ã y and a Ã z are the typical lengths in the y Ã and z Ã directions, respectively.
In order to examine the thermal effects inside the bubble, the energy equation [57] is used (see (11) below) [55,56,66,67]. The pressure and temperature distributions of the gas inside the bubble [68]  are not examined directly. Instead, they are treated as the average values of the gas inside bubble. The temperature fluctuation of the gas inside the bubble is considered; however, that of the liquid phase is assumed constant. The bubble oscillation is spherically symmetric. The bubble-bubble interaction [69][70][71][72] is not considered. Bubbles do not coalesce, break up, appear, and disappear. Gas inside the bubbles is composed of only non-condensable gas, and the mass transfer induced by phase change at the bubble-liquid interface [77][78][79], which would be important in a scenario of high pressure and long pulses, is not considered. Further, heat transfer at the bubble-liquid interface [80] is not considered. Gas viscosity, viscosity of bulk liquid, and elasticity (from the viewpoint of medical applications) of the liquid phase [72,81] are not considered. Initially, the bubbly liquid is at rest and spatially uniform, except for the bubble distribution. The initial polydispersity of the bubbles [73][74][75][76] is not considered. Only the bubble distribution is spatially nonuniform in the initial state [53][54][55][56] because medical applications, such as HIFU treatment, often only use the bubbles at the focus.

Basic equations
Although simplified conservation equations based on a mixture model [58] were used in our previous works [55,56], the volumetric averaged conservation laws of mass and momentum for the gas and liquid phases, based on a two-fluid model [60], are used in this study: where t Ã is the time, α is the void fraction, ρ Ã is the density, u Ã is the velocity vector, p Ã is the pressure, and P Ã is the liquid pressure averaged at the bubble-liquid interface [60]. The subscripts G and L denote volume-averaged variables in the gas and liquid phases, respectively. The subscript * denotes a dimensional quantity. The virtual mass force model [60] is introduced as the interfacial momentum transport term F Ã where β 1 , β 2 , and β 3 are constants and may be set equal to 1/2 for the spherical bubble. For simplicity, the drag force [82][83][84], lift force, and gravitation are neglected. The Lagrange derivatives D G =Dt Ã and D L =Dt Ã are defined as The Keller equation for spherical bubble oscillation in a compressible liquid [85] is given as where R Ã is the bubble radius, c Ã is the sound speed, and the subscript 0 represents the physical quantity in the initial state. The energy equation [57] for thermal conduction at the bubble-liquid interface is used to express the thermal effect inside the bubble where κ is the ratio of specific heats, λ Ã G is the thermal conductivity of the gas phase. T Ã G is the temperature of the gas phase, and r Ã G is the radial distance from the center of the bubble. The temperature gradient term @T Ã G =@r Ã G j r Ã G ¼R Ã , in the first term on the right-hand side of (11), is rewritten using the following four models [66,67]: where T Ã 0 is the initial temperature of the liquid and gas phases, D Ã G is the thermal diffusivity of the gas phase, and Re and Im denote the real and imaginary parts, respectively. The main features of the temperature models (12)- (15) are summarized in our previous study [66]. Note that the STM model (15) incorporates the effect of a phase difference of time between the average temperature inside the bubble and the temperature gradient at the bubble-liquid interface. ω Ã B is the linear natural frequency of a single bubble oscillation [64], given as ð16Þ where γ e is the effective polytropic exponent, μ Ã e0 is the initial effective viscosity, μ Ã L is the liquid viscosity, Γ N and α N are complex numbers, and i denotes the imaginary unit. The explicit form of (16) is different from that used in our previous studies [53][54][55][56]59,86,87]. e L Ã P in (14) and (15) is the complex number with the length dimension, given as Note that the form of ω Ã B in (16) is valid only for the linear oscillations of symmetric waves with a sinusoidal shape. However, this study investigates the nonlinear propagation, wherein the wave form will develop an asymmetric form. Additionally, the natural frequency of bubble oscillations is affected by the pressure amplitude [88,89]. Furthermore, bubble-bubble interaction and dual-frequency ultrasound, which are not considered in this study, also affect the natural frequency of bubble oscillations [69][70][71][72][90][91][92], the latter of which can be employed to enhance the effects of bubble oscillations. Since incorporating these effects in the derivation process is quite difficult, the natural linear frequency of bubble oscillations ω Ã B is used in this study for simplicity. To close the set of equations, we use Tait's equation of state for the liquid phase, the equation of state for the ideal gas, the conservation equation of mass inside the bubble, and the balance of normal stresses across the bubble-liquid interface: where n is the material constant (e.g., n ¼ 7:15 for water). Although (10), (11), (24), and (25) were originally derived for a single bubble, this study deals with them as the volumetric averaged equations for liquids containing many bubbles. Hence, the bubble radius R Ã t Ã ð Þ, depending on only the time t Ã in the original equations, is treated as a volumetric averaged variable and extended as a multivariable function of the temporal and spatial variables, i.e., Þ , for the 2D and 3D KZK equations, respectively.
The liquid viscosity and initial effective viscosity [66] are nondimensionalized as Nondimensionalizations for the energy equations [66] are as follows: for (12) for (13) 3 for (14) 3 for (15) 3 where ζ SMK , ζ LSM , ζ PCB , ζ STM1 and ζ STM2 are the constants of O 1 ð Þ; ζ STM2 in (36) denotes the effect of the phase difference of the time between the average temperature inside the bubble and the temperature gradient at the bubble-liquid interface.

Multiple scale analysis
The independent variables t Ã , x Ã , r Ã , y Ã , and z Ã are nondimensionalized as where T Ã is the typical period of a wave and is related to the wave speed U Ã and wave length L Ã by U Ã ¼ L Ã =T Ã . Next, the nondimensionalized time t and the spatial coordinate x, for the direction of the ultrasound propagation, are expanded to the near field and far field using the nondimensional wave amplitude ε, for both 2D and 3D KZK equations [65] Owing to the assumption of weakly diffracted waves in (2), the spatial variations in the radial distance from the x axis appear only in the far field and not the near field [53][54][55][56]87]. Then, the spatial coordinates r, y, or z for the radial distance from x are defined only for the far field as Although the relationship in (41) for 2D KZK is the same as in our previous study [53][54][55][56]87], the relationship for 3D KZK is introduced in this study. For (41), the relationships are taken as All the dependent variables are regarded as functions of the extended independent variables of (37)- (42). The differential operators are thus expanded as [65] @ @t The dependent variables are expanded in powers of ε where u Ã is the velocity component in the x Ã direction. The velocity components in the r Ã , y Ã , and z Ã directions are v Ã r , v Ã y , and v Ã z , respectively. The velocity components are expanded for the gas and liquid phases as The magnitude of the variations in the velocity components in the r Ã , y Ã , and z Ã directions vertical to the x Ã axis are assumed to be smaller than those in the x Ã direction [53][54][55][56]87]. Then, the expansions of v Ã r , v Ã y , and v Ã z begin with a higher order than that of u Ã in (49)- (52). The density of the gas phase at the initial state is assumed to be significantly smaller than that of the liquid phase The nondimensional pressures of the liquid and gas phases at the initial state are defined as The expansion of the liquid density is given as which is determined using (22), (26), and (48) [86].
To consider the effect of the weak nonuniformity of the bubble number density at the initial state, the void fraction α is expanded as [53][54][55][56]87] where δ is the known variable and denotes the nonuniformity of the bubble number density at the initial state. The effect of the nonuniformity of the bubble number density is assumed to appear only at the far field from the sound source, thus δ is only a function of x 1 [53][54][55][56]87].

First order of approximation
From the approximations of O ε ð Þ, linear equations are obtained for (4)-(7), (10), and (11) as These equations are different from our previous study [55,56], which was based on the mixture model equations. (57)-(62) are combined into a single equation for R 1 as where the phase velocity v P is given as and s 4 is given as For simplicity, we set v P ≡ 1, and the typical wave speed U Ã is given as where s Ã 4 is given as Next, the following variable transformation is introduced: Then, we obtain the equation describing the propagation of rightrunning waves as By introducing the variable transform (68) into the approximated (57)-(62), all dependent variables can be written in terms of R 1 ¼ f as with

Approximation of radial direction
At the approximations of O ε 2=3 À Á , the equations obtained from the radial directions of the momentum conservation equations (6) and (7) are The independent variables r 1=2 , y 1=2 , and z 1=2 first appear here. The variable transformation (68) is introduced into (75) and (76), and then we obtain the relations (i = G, L) of where the constant coefficients C i (i ¼ G; L) are given as

Second order of approximation
From the approximations of O ε 2 ð Þ, equations are obtained for (4)-(7), (10), and (11) as @α 2 @t 0 À 3 where the explicit forms of the inhomogeneous terms K i i ¼ 1; 2; 3; 4; 5; 6 ð Þ are presented in Appendix A. Equations (80)-(85) are combined into the single equation The inhomogeneous term K is given as The differential operators (43)- (45), the relation in the near field (69), the relation in the radial direction (77)- (79) and (88) are combined as @ @t Finally, the 2D and 3D KZK equations are obtained: where Δ ? is the Laplacian with respect to the radial direction: The following variable transformations are used: where τ is the retarded time.

Coefficients of the KZK equation
In the KZK equation (90), Π 0 and Π 4 are the advection coefficients, Π 1 is the nonlinear coefficient, Π 21 and Π 22 are the dissipation coefficients, and Π 3 is the dispersion coefficient. The right hand side of the KZK equation (90) represents the diffraction (focusing) effect, and these terms are expressed differently for the 2D and 3D KZK equations. The nonuniformity of the bubble number density δ appears only in the variable transformation (92) with the advection coefficient Π 4 , thus, it only influences the advection effect [53][54][55][56]87]. Figure 2 shows the dependence of the (a) nonlinear, (b) and (c) dissipation, and (d) dispersion coefficients versus the initial void fraction α 0 for the normal condition of the air-water system. Note that the effect of bubble-bubble interaction [69][70][71][72] would be dominant when the void fraction increases; however, this effect is not considered in this study.
The nonlinear coefficient Π 1 is given as with In Fig. 2 (a), the nonlinear coefficient Π 1 in our previous study [55,56] derived from the mixture model, does not depend on α 0 . However, it does depend on α 0 in this study because the two-fluid model [60] is used. Furthermore, Π 21 , Π 3 , and Π 4 are given as The first term of the dissipation coefficient Π 21 (103) relates to the interfacial viscosity, the second and third terms relate to the liquid compressibility. Fig. 2 (b) shows the dependence of Π 21 divided into the effects of the interfacial viscosity and the liquid compressibility, showing that the interfacial viscosity is dominant. The viscous effects of bulk liquid were considered in our previous study [55,56]; however, these effects are omitted in this study because the basic equations based on the two-fluid model that incorporate the effect of the viscosity of the bulk liquid have yet to be developed. The advection coefficient Π 0 depends on which temperature gradient model from (12)-(15) is used (i) For (12)- (14), (ii) For (15), ◂ From (107), a phase difference between the average temperature inside the bubble and the temperature gradient at the bubble-liquid interface affects the advection effect. The dissipation coefficient Π 22 is given for each temperature gradient model in (12)-(15) (j = SMK, LSM, PCB, STM1) as Fig. 2 (c) shows the dependence of the dissipation coefficient Π 22 for cases where the gas inside the bubble is air, Argon (Ar) and sulfur hexafluoride (SF 6 ). Shimada et al. [61] (SMK) temperature gradient model in (12) is used. Although Π 21 is the dissipation coefficient owing to the interfacial viscosity and the liquid compressibility, Π 22 is the dissipation coefficient without differentiation owing to the thermal conduction of the gas inside the bubble. Then, Fig. 2 shows that the effects of thermal conduction decrease from Argon to air to SF 6 . This result agrees with the analysis for the single-bubble oscillation proposed by Matsumoto et al. [93].

Comparison with original KZK equation
The original KZK equation [6,7] for pure liquid in the dimensional form is given as follows: where β and δ Ã are the nonlinear and dissipation coefficients, respectively. The most important difference between the present KZK equation (90) and the original (109) is that the former has the dispersion term as a third-order partial derivative. Further, in the present equation (90), the dissipation term without differentiation was discovered due to the bubble oscillations and the thermal conduction of gas inside bubbles.

Limitation of present model
The present KZK equation (90) is derived by perturbation analysis up to second-order approximation, in which nonlinear propagation is included as the term with the coefficient Π 1 . However, the viscosity, compressibility, and thermal conductivity dissipation terms are represented only in linear form. Although the nonlinear dissipation effects [94][95][96][97] will be important particularly when the pressure amplitude increases, they have been omitted in the present equation. Perturbation analysis with greater than third-order approximation will be necessary to incorporate the nonlinear dissipation effects and will be performed in a forthcoming study.

Method
We numerically solve the newly obtained 2D KZK equation (90) using the FDTD method developed by Lee and Hamilton [22,23]; this scheme has been widely utilized for the problem of focused ultrasound radiated by a circular sound source in pure water [9,[13][14][15]24,[26][27][28][29][30][31]. For consistency between theoretical studies [53][54][55][56]87] and the numerical studies of Lee and Hamilton [22,23], a scaling relation is assumed as where d Ã f is the focal length, and χ is the nondimensional constant of O 1 ð Þ. Hereafter, we set χ ¼ 1 for simplicity. For numerical calculations, the KZK equation (90) is rewritten by the definite integral regarding the retarded time τ, and the dependent variable is changed from the nondimensionalized bubble radius fluctuation f to the nondimensionalized temperature fluctuation of the gas inside the bubble T G1 through (70) as In previous studies of focused ultrasound in a pure liquid, the KZK equation is derived in the form describing the distribution of liquid pressure. However, in this study, we succeeded in describing the temperature distribution of the gas inside the bubble T G1 by introducing the energy equation (11). In (111), G r is the focusing gain, which is given as where (29) and (110) are used. The variable transformations (92) are also rewritten for consistency with Lee and Hamilton [22,23] as Next, the KZK equation (111) is solved term by term to incorporate the effects of diffraction, nonlinearity, dissipation, and dispersion separately in the same manner as Lee and Hamilton [22,23]. However, the dissipation term Π 22 T G1 and the dispersion term Π 3 @ 3 T G1 =@τ 3 in (111) are obtained by consideration of the effect of bubbles in this study. While Π 22 T G1 can be introduced easily, Π 3 @ 3 T G1 =@τ 3 is discretized for the finite difference approximation as where Δτ and i are the step size and the index in the τ direction, respectively. The boundary condition at X ¼ 0 for the focused sound source is given as where p Ã 0 is the effective pressure amplitude of the sound source, and E τ ð Þ is the amplitude modulation function. Note that the boundary condition in the algorithm proposed by Lee and Hamilton [22,23] was given for the liquid pressure. The example of the wave form at the sound source (i.e., the boundary condition) can be seen in Fig. 4 (a).
The step sizes for each direction are set as Δτ ¼ 2π=240, ΔX ¼ 0:001, and ΔR ¼ 0:01. In the calculations, there are three main sources of error [22,23]; (i) finite difference approximations of the diffraction and dissipation terms, (ii) incorporating nonlinear effects based on the implicit analytical solution, (iii) including the diffraction, dissipation, and nonlinear effects separately. Considering these three main sources of error, the total error is estimated as Δτ ð Þ 2 þ ΔX þ ΔR. Further details regarding the error are discussed by Lee and Hamilton [22,23]. Fig. 3 shows the spatial distribution of the first order dimensional temperature fluctuations of the gas inside the bubble T Ã G1 , where T Ã G1 ¼ εT G1 T Ã 0 and the initial temperature T Ã 0 is set as the normal temperature of the human body. The parameters used for calculations are; p Ã 0 ¼ 100 kPa, the frequency of the sound source f Ã ¼ 100 kHz, d Ã f ¼ 300 mm, a Ã r ¼ 150 mm, R Ã 0 ¼ 0:05 μm and α 0 ¼ 0:001. The Shimada et al. [61] (SMK) temperature gradient model in Eq. (12) is used. In Fig. 3, the gas inside the bubble is Argon in (a), air in (b), and SF 6 in (c). In the area near the focus at r Ã ¼ 0 mm and x Ã ¼ 300 mm, an intensive temperature rise is shown. The maximum temperature rise is approximately 28:1 K in (a), 17:7 K in (b), and 4:13 K in (c), with the Argon gas exhibiting the most effective temperature increase. Fig. 4 shows the time development of the temperature fluctuation T Ã G1 at (a) r Ã ¼ 0 mm, x Ã ¼ 0 mm (at the sound source) and (b) r Ã ¼ 0 mm, x Ã ¼ 300 mm (focus) of Fig. 3 (a), where the gas inside the bubble is Argon. In Fig. 4 (b) at the focus, wave distortion due to nonlinear effects is shown. In Fig. 4, the time development of two types of dissipation terms; Π 21 @ 2 T Ã G1 =@t Ã2 owing to the interfacial viscosity and liquid compressibility, Π 22 T Ã G1 owing to the thermal conductivity of gas inside bubbles; are also shown. Although in Fig. 4 (a) at the sound source, the values of the two dissipation terms are of the same order, in Fig. 4 (b) at the focus, Π 21 @ 2 T Ã G1 =@t Ã2 becomes dominant especially at the discontinuous point of T Ã G1 . Note that when we use the temperature gradient models, except for the Shimada et al. [61] (SMK) model of (12), the waves rapidly attenuate even near the sound source.  Fig. 3 (a), where the gas inside the bubble is Argon. The black, red, and blue curves represent the waveforms of the temperature fluctuation, the dissipation term owing to the interfacial viscosity and liquid compressibility, and the dissipation term owing to the thermal conductivity of gas.

Result
Hence, other temperature gradient models probably overestimate the dissipation effect of the thermal conduction of the gas inside bubbles. Table 1 shows the values of the calculation of Fig. 3. When the gas inside the bubble is SF 6 , the value of the temperature rise and the dissipation coefficient owing to thermal conductivity Π 22 are the lowest of the three types. This result agrees with the analysis for the singlebubble oscillation [93]. In addition, the value of the nonlinear coefficient Π 1 =s 5 is the highest. For the applications of bubbles as contrast agent for ultrasound imaging, a temperature rise is not required, but higher-harmonic generation induced by nonlinear effects is necessary. Hence, achieving a result where the temperature rise is lowest, and the nonlinear coefficient is highest for SF 6 is consistent with the applications for ultrasound imaging.

Conclusions
We theoretically and numerically examine a weakly nonlinear propagation of focused ultrasound in liquids containing many spherical microbubbles with a nonuniformity of bubble number density. The KZK equation [6,7], which has long been used  as the physico-mathematical model for focused ultrasound in pure liquids, is extended to bubbly liquids using the volumetric averaged equations for bubbly liquids based on the two-fluid model [60]. For the description of the temperature rise and the thermal effects of gas inside bubbles, the energy equation [57], (11), for the gas inside a bubble is used [55,56,66,67]. The temperature gradient term in (11) is rewritten by the four major models of (12)- (15). When we use the temperature gradient models, the numerical results show that the waves rapidly attenuate even near the sound source except for the Shimada et al. [61] (SMK) temperature gradient model of (12). Hence, it is likely that the other temperature gradient models of (13)-(15) overestimate the dissipation effect of the thermal conduction of the gas inside bubbles.
The resultant generalized KZK equation (90) for bubbly liquids is represented as the linear combination of nonlinear, dissipation, dispersion, and diffraction (focusing) terms. The dispersion term is owing to bubble oscillations. The diffraction terms are expressed differently for the 2D and 3D KZK equations. By using volumetric averaged equations based on a two-fluid model [60], the dependence of the nonlinear coefficient Π 1 (95) on the initial void fraction α 0 is obtained, as shown in Fig. 2 (a). The dissipation effect is divided into the two terms; the term with the second-order derivative is due to the interfacial viscosity and liquid compressibility, and the term without differentiation is due to the thermal conductivity of the gas inside bubbles. Hence, we succeed in consistently describing ultrasound propagation, bubble oscillations, and temperature increase in a single equation.
The KZK equation is then numerically solved by the FDTD method [22,23]. For comparison of the value of the two dissipation terms, Π 21 @ 2 T Ã G1 =@t Ã2 being based on the interfacial viscosity and the liquid compressibility and Π 22 T Ã G1 being based on the thermal conductivity of the gas inside the bubbles are of the same order at the sound source; however, the former becomes dominant at the focus. Finally, the spatial distribution of temperature fluctuations is obtained for Argon, air, and SF 6 as the gas inside the bubble. As shown in Table 1, Argon gas shows the highest temperature rise, making it effective for applications in tumor ablation therapy by HIFU. SF 6 gas exhibits the lowest temperature rise and the highest nonlinearity, making it effective for applications in ultrasound imaging.
In a future work, theoretical extensions of the KZK equation incorporating the viscous effects of bulk liquid, the elasticity of body tissues, mass transfer, heat transfer, and higher-order approximation to describe nonlinear dissipation will be carried out. Particularly, the description of blood vessels [98][99][100], effect of initial nonuniform distributions of velocities [101], and effect of a shell encapsulating microbubble toward ultrasound diagnosis and therapy [102]

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. 1 @t 0 þ 3 @R 1 u G1 @x 0 À @α 1 u G1 @x 0 þ δ x 1 ð Þ 3 @R 1 @t 0 À @u G1 @x 0 ; ðA3Þ K 2 ¼ 1 À α 0 ð Þ @u L1 @x 1 À α 0 @α 1 @t 1 À α 0 @α 1 u L1 @x 0 þ 1 À α 0 ð Þ @ρ L1 @t 0 À α 0 δ x 1 ð Þ @u L1 @x 0 : The remaining inhomogeneous terms K 3 , K 4 , K 5 , and K 6 are common to both the 2D and 3D KZK equations K 3 ¼ 3p G0 @R 1 @x 1 À p G0 @T G1 @x 1 À β 1 @u G1 @t 1 À @u L1 @t 1 Àβ 1 α 1 @u G1 @t 0 À @u L1 @t 0 À β 2 u G1 À u L1 ð Þ @α 1 @t 0 þ3p G0 α 1 @R 1 @x 0 À 6p G0 @R 2 1 @x 0 À β 1 u G1 @u G1 @x 0 The explicit form of K 6 depends on which temperature gradient model is used, i.e., for (12)-(15): [66] (i) For (12)- (14), (ii) For (15), wherê Appendix B. Derivation in our previous study We briefly review the derivation in our previous study [55,56]. The volumetric averaged conservation laws of mass and momentum based on a mixture model [58] was used: where ρ Ã is the density of mixture and defined as follows: where the density of the gas phase is sufficiently smaller than that of the liquid phase, thus it is omitted. The conservation equation of the bubble number density n Ã and the relationship between the bubble number density and the void fraction are used: The set of equations were closed by (10)- (25). The form of the linear natural frequency of a single bubble ω Ã B was different from (16) of this paper: As the same procedure of this paper, the same form of the resultant KZK equation (90) is derived. However, the coefficients of the resultant KZK equation are greatly different from this paper and given as ; ðB7Þ where the symbol ∼ means the coefficients of our previous study [55,56]. Note that in our previous study [55,56], only the 2D KZK equation was derived, and only the LSM temperature gradient model of (13) was used.