Microbubbles coated by viscoelastic shells are important for ultrasound diagnosis using contrast agents, and the dynamics of single coated bubbles has been investigated in the literature. However, although a high number of contrast agents are used in practical situations, there has long been an absence of a nonlinear acoustic theory for multiple coated bubbles, except for our recent work by Kikuchi and Kanagawa [“Weakly nonlinear theory on ultrasound propagation in liquids containing many microbubbles encapsulated by viscoelastic shell,” Jpn. J. Appl. Phys. 60, SDDD14 (2021)], under several assumptions to be excluded. Aiming for generalization, in this study, we theoretically investigate weakly nonlinear propagation of ultrasound in liquid containing multiple bubbles coated by a viscoelastic shell with compressibility. Leveraging the method of multiple scales, both the Korteweg–de Vries–Burgers (KdVB) equation for a lowfrequency long wave and nonlinear Schrödinger (NLS) equation for a highfrequency short wave are derived from the volumetric averaged equations for bubbly liquids based on a twofluid model and the uptodate model for single coated bubbles with shell compressibility. Neglected factors in our previous paper, i.e., compressibility of the shell and liquid, drag force acting on bubbles, bubble translation, and thermal conduction, are incorporated in the present KdVB and NLS equations; the proposed model will be regarded as a generic physicomathematical model. The results show that shell compressibility attenuated ultrasound strongly and decreased nonlinearity of ultrasound. Finally, we compared the magnitudes of six dissipation factors (shell compressibility, shell viscosity, liquid compressibility, liquid viscosity, thermal effect, and drag force) for five typical ultrasound contrast agents, and a similar tendency between KdVB and NLS equations was revealed.
I. INTRODUCTION
Medical ultrasound has several important applications such as diagnosis via contrast images, tumor ablation therapy using high intensity focused ultrasound, and extracorporeal shock wave lithotripsy.^{1–8} Since the combination of micro (or nano) bubbles and ultrasound enhances an efficiency in medical practices,^{9,10} interaction of cavitating bubble oscillations and propagation of waves (i.e., shock, ultrasound, and pressure wave) has extensively been investigated (see, e.g., recent analytical and/or numerical progresses^{11–33}). Using microbubbles as ultrasound contrast agents (UCAs) has been shown to improve the contrast of the blood carrying tissue.^{34} UCAs were originally uncoated microbubbles used for echocardiography.^{35} However, in order to improve the stability and biocompatibility of the bubbles, protein coated bubbles such as Albumin was proposed, and then, lipidcovered microbubbles were developed as new UCAs and are now used to detect liver cancer.^{36} Further applications of coated microbubbles are expected to include therapeutic applications such as transportation in drug delivery systems.^{37–43}
Therefore, a precise understanding of the interaction between ultrasound propagation and oscillations of coated microbubbles as UCAs is desirable. Experiments on attenuation (or dissipation) and scattering were performed using various types of coatings to understand their properties,^{10,14–51} and molecular dynamics simulations were performed to clarify the properties of the shells.^{52,53} In addition, shell parameters were modeled based on intermolecular interactions.^{54,55} Particularly, Sojahrood et al. investigated the dynamics of coated (or uncoated) bubbles through bifurcation structure,^{56–59} sub and superharmonic behaviors,^{60–63} and the mechanics of nonlinear power dissipation.^{64–66}
Several physicomathematical models have been derived or proposed based on the equation of motion for coated bubbles. Equations of motion of single uncoated bubbles such as the Rayleigh–Plesset equation and Keller–Miksis equation^{67} have been modified previously.^{68–83} Furthermore, a pioneering work toward modeling single coated bubbles has been reported by de Jong^{45} and Church^{70} proposed a more rigorous theoretical model by treating the coating shell as an incompressible Kelvin–Voigt model. Hoff et al.^{75} made the assumption that the shell thickness is small enough compared with the bubble radius and derived a new model from the Church model.^{70} Previously, Sarkar et al. studied the viscous effect of the shell and developed a viscoelastic interface model for the thin coated bubble.^{80} Moreover, Marmottant et al.^{77} developed a model considering the physical properties of a lipid monolayer coating on a gas microbubble. Following this, other nonlinear models have been proposed with more complex rheological behaviors, such as strain softening and hardening^{78,83} or shear thinning.^{72} In addition, Liu et al. formulated the radial oscillations, translational motions, and deformations for double coated bubbles.^{84} Gubaidullin and Fedorov^{74} investigated heat transfer with coated bubbles. However, these studies^{10,45,58,69–88,112} were restricted to single coated bubbles or some (i.e., not a high number of) coated bubbles. Since multiple (i.e., a high number of) coated bubbles are used as contrast agents in practical applications, an understanding of multiple coated bubbles is desired. Although the basic theory for multiple coated bubbles in relation to ultrasound propagation should be established, our recent papers^{89,90} are among few that report on multiple coated bubbles, from the viewpoint of multiphase fluid dynamics (i.e., the formulation based on continuum model). Note that although Ma et al.^{91} obtained the parameters of acoustic nonlinearity for ultrasound in liquids containing multiple coated bubbles, they did not discuss the dissipation and various effects of propagation of nonlinear ultrasound.
The goal of the present study is (a) the consideration of the compressibility of the shell and (b) unified treatment of both lowfrequency long wave and highfrequency short wave bands. With regard to (a), we extend the previously established microscopic modeling into macroscopic dynamics of multiphase flows as a continuum approach (i.e., the extension of the dynamics of single coated bubbles to that of multiple coated bubbles by combining the equation of motion for single coated bubbles and the set of volumetric averaged equations for bubbly liquids). The most recent model^{68} incorporated the compressibility of coated shells and showed that neglecting the compressibility leads to an underestimation of the shear modulus. Based on the results of this paper, we found the nonnegligible contribution of shell compressibility to the ultrasound attenuation. We completed the theory by generalizing our previous work^{89} to incorporate other previously neglected factors, i.e., shell compressibility, liquid compressibility, influence of drag force on bubbles, translation of bubbles, and thermal conductivity. Equation for thermal conduction at the bubble–liquid interface^{92} is used. To the best of our knowledge, this is the first report considering both multiple coated bubbles and thermal conduction. With regard to (b), we focused on both lowfrequency long waves and highfrequency short waves, and derived the Korteweg–de Vries–Burgers (KdVB) and nonlinear Schrödinger (NLS) equations, respectively. Using both the low and high frequency modes is also important in a medical application. For the case of multiple uncoated bubbles, the KdVB and NLS equations have been derived independently (e.g., pioneering classical works^{93–96}). By changing the magnitude of the set of nondimensional ratios appropriately to low and highfrequency bands,^{97} we can systematically derive the KdVB and NLS equations (i.e., two types of nonlinear wave equations) from one set of basic equations in a unified manner. The resulting KdVB and NLS equations describe the weakly nonlinear propagation/evolution process of ultrasound in liquids containing multiple coated bubbles. The effect of the physical properties of the shell (e.g., compressiblity, shear modulus, viscosity) and thermal conductivity inside the bubbles on the propagation properties of ultrasound, i.e., dissipation, nonlinear, and advection effects, are clarified based on the KdVB and NLS equations. The coefficients in the dissipation, nonlinear, and advection terms of the KdVB and NLS equations are then compared for typical UCAs. The connection between the present study and our previous two studies on ultrasound in liquids containing multiple coated bubbles is as follows: (i) Although Ref. 89 originally derived only the KdVB equation without shell compressibility, the present study considers the unified derivation of the KdVB and NLS equations, making this the first derivation of an NLS equation for multiple coated bubbles. (ii) Reference 90 derived the KdVB equation incorporating the effect of shell compressibility but did not discuss the dependence of the initial void fraction as one of most important points in bubble dynamics. Here, we use twofluid model averaged equations to incorporate the dependence of the initial void fraction on each coefficient in both the KdVB and NLS equations.
The remainder of this paper is organized as follows. Section II formulates the problem by introducing the updated equation of motion for coated bubbles with shell compressibility^{68} and the volumetric averaged equations for bubbly liquids.^{98} Section III derives the KdVB equation for the lowfrequency long wave. Section IV derives the NLS equation for the highfrequency short wave. Section V discusses the dissipation, nonlinear, and advection effects in the resultant KdVB and NLS equations, and Sec. VI concludes the paper.
II. FORMULATION OF THE PROBLEM
A. Problem and assumptions
As shown in Fig. 1, uniform weakly nonlinear propagation of ultrasonic (or acoustic or pressure) waves in compressible liquids containing multiple spherical coated bubbles is studied theoretically. The number of coated bubbles (see Fig. 2) is significantly high, which can be formulated using the continuum model shown in (4) and (5) below. For simplicity, the following assumptions are used:

Two propagation modes, that is, lowfrequency long wave (KdVB equation in Sec. III) and highfrequency short wave (NLS equation in Sec. IV), are investigated.

The shell coating the bubble is a viscoelastic body that obeys the Kelvin–Voigt model.

Although the shell is compressible,^{99} we regard the thickness of shell as a constant, for simplicity.

The propagation of onedimensional (plane) waves is treated.

The liquid is compressible.

The bubbles do not coalesce, break up,^{100} appear, or disappear, and bubble–bubble interaction^{101–105} is neglected.

The bubble oscillations are spherically symmetric.

Initially, the gas and liquid phases flow with a constant uniform velocity, and initial nonuniform velocity distributions of each phase^{106} are not considered.

The effect of liquid viscosity is considered only at the bubble–liquid interface.

The phase change^{99,107,108} is neglected.

The mass transport across the bubble–liquid interface^{109} is neglected.
B. Basic equations
To describe the dynamics as the expansion and contraction of coated bubbles, we used the updated model, i.e., Chabouh et al.'s equation,^{68} for spherically symmetric volumetric oscillations in a compressible liquid with coated bubbles. Furthermore, we incorporated the translation of bubbles [i.e., the last term in the righthand side in (1)] into Chabouh model.^{68} We then extended Chabouh equation^{68} for single coated bubbles without translation to the case of multiple coated bubbles with translation as follows:
Here, the pressures of surrounding liquid are described by the balance equation of normal stresses as follows:
Here, as an independent variable, $ t \u2217$ is the time; as dependent variables, $ R \u2217$ is the bubble radius, $ u G \u2217$ and $ u L \u2217$ are the gas and liquid velocities, respectively, $ p G \u2217$ and $ p L \u2217$ are the volumetricaveraged gas and liquid pressures, respectively, and $ P *$ is the surface averaged liquid pressure; as constants, $ G \u2217$ is the shear modulus of shell, $ K \u2217$ is the bulk modulus of shell, $ \mu S \u2217$ is the shear viscosity of the shell, $ \mu K \u2217$ is the viscosity that describes the friction losses due to volume change of the shell, $ \rho L 0 \u2217$ is the initial liquid density, $ p G 0 \u2217$ and $ p L 0 \u2217$ are the initial gas and liquid pressures, respectively [see also (12) below], $ \mu L \u2217$ is the liquid viscosity, $ c L 0 \u2217$ is the sound speed in the unperturbed liquid, $ d 0 \u2217$ is the initial shell thickness, $ R 0 \u2217$ is the initial bubble radius, and $ \sigma 1 \u2217$ and $ \sigma 2 \u2217$ are the surface tensions at the internal and external boundaries of the shell, respectively; the superscript * represents dimensional quantities, the subscript 0 denotes the initial unperturbed state, and the subscripts $G$ and $L$ represent the gas and liquid phases, respectively; $ D i / D t \u2217 ( i = G , L )$ are material derivatives for the gas and liquid phases, respectively,
where the independent variable $ x *$ is the spatial coordinate.
We defined uncoated bubbles by setting $ d 0 \u2217 = 0$. For the case of uncoated bubbles without translation, (1) coincides with the Keller–Miksis equation.^{67} We treated R^{*} not only as a time variable but also as a field variable [i.e., $ R \u2217 = R \u2217 ( x \u2217 , \u2009 t \u2217 )$], to extend (1) to the case of multiple bubbles. Since $ R \u2217$, $ p G \u2217$, $ p L \u2217$, $ P \u2217$, $ u G \u2217$, and $ u L \u2217$ in (1) are unknown variables, we should use conservation equations as a continuum model to close the set of equations. As volumetricaveraged equations based on a twofluid model,^{98} the conservation laws of mass and momentum for the gas and liquid phases are introduced as
where α is the void fraction (i.e., volumetric fraction of gas phase; $ 0 < \alpha < 1$), $ F \u2217$ is the interfacial momentum transport, and $ D *$ is the drag force. The following model of virtual mass force for $ F \u2217$ is used:^{110,111}
where for a spherical bubble, the values of the coefficients β_{1}, β_{2}, and β_{3} are 1/2. For the case of spherical bubbles, D^{*} is given as follows:^{112–114}
where $ C D$ is the drag coefficient for a single spherical bubble.
To describe thermal conduction at the bubble–liquid interface, the following equation for thermal effect inside bubbles^{92} is introduced:
where $ T G \u2217$ is the temperature of the gas phase, κ is the ratio of specific heats, $ r \u2217$ is the radial distance from the center of the bubble, and $ \lambda G \u2217$ is the thermal conductivity of the gas inside the bubble. For the temperaturegradient as the first term on the righthand side of (10), we use the model proposed by Sugiyama et al.,^{115}
where $ T 0 \u2217$ is the initial temperature. Furthermore, the resonant frequency of the linear spherical symmetric oscillation of a single bubble, $ \omega B \u2217$, is obtained from (1),
with the effective polytropic exponent $ \gamma e$, the initial effective viscosity $ \mu e 0 \u2217$, and the complex numbers $ \Gamma N$ and $ \alpha N$ are given by
where $i$ denotes the imaginary unit, and the complex number with the length dimension, $ L \u0303 P \u2217$ in (11), is given as
Finite amplitude propagation of pressure waves in liquids containing a large number of oscillating coated bubbles with a finite amplitude is the present target. For simplicity, we utilize the linear resonant frequency for the single bubble $ \omega B \u2217$ to distinguish low and highfrequency bands as shown in (33), we provide some comments for the linear resonant frequency for interacting bubbles. Guédra^{102} proposed an analytical form for the resonance frequency of the interacting bubbles by using $ \omega B \u2217$ and coupling strength. Yasui et al.^{88} proposed a simple parameter called the coupling strength describing an approximation of the pressure of other bubbles at the location of each bubble, which is calculated from the distance between the bubbles and the number of bubbles. In a forthcoming paper, we will incorporate these effects in $ \omega B *$.
Equations (1), (4)–(7), and (10) contain ten unknown variables, i.e., $ R \u2217$, $ p G \u2217$, $ p L \u2217$, $ P \u2217$, $ u G \u2217$, $ u L \u2217$, α, $ \rho G \u2217 , \u2009 \rho L \u2217$, and $ T G \u2217$, and should be closed by the Tait equation of state for the liquid phase, the equation of state for an ideal gas, and the conservation equation of mass inside the bubble, as follows:
where n is the material constant (e.g., n = 7.15 for water).
Note that our previous paper^{89} omitted the compressibility of the shell and liquid, drag force acting on bubbles, translation of bubbles, and thermal conductivity.
C. Multiple scales analysis
Independent variables are first nondimensionalized as
where $ T \u2217$ and $ L \u2217$ are the typical period and wavelength, respectively. Using t and x, the nearfield, farfield I, and farfield II are described by extended independent variables in the method of multiple scales as follows:^{116}
where ϵ is a nondimensional wave amplitude $ ( 0 < \u03f5 \u226a 1$, i.e., a weakly nonlinear problem^{116}). The differential operators are immediately expanded as
It is worth noting that farfield II is not used in Sec. III.
The dependent variables are now regarded as functions of (22)–(24) and expanded in a power series in ϵ as follows:
where $ U \u2217 \u2009 ( \u2261 L \u2217 / T \u2217 )$ is the typical propagation speed of a wave, and the nondimensional pressures for the gas and liquid phases of O(1) are given as $ p G 0 = p G 0 \u2217 / ( \rho L 0 \u2217 U \u2217 2 )$ and $ p L 0 = p L 0 \u2217 / ( \rho L 0 \u2217 U \u2217 2 )$. The ratio of the initial densities of the gas and liquid phases is $ \rho G 0 \u2217 / \rho L 0 \u2217 = O ( \u03f5 3 )$. The effect of initial polydispersity^{117–119} is neglected in this paper for simplicity. The expansion of liquid density is categorized into two cases.^{97}
The scaling relations of the nondimensional ratios among the physical parameters appropriate for the low frequency and long wave and the high frequency and short wave are defined as follows:^{97}
where $ \omega \u2217$ is a typical angular frequency of an incident wave, and Δ, Ω, and V are constants of O(1). Furthermore, the shell shear viscosity, viscosity describing the friction losses due to volume change, shell shear modulus, shell bulk modulus, and shell thickness are nondimensionalized as follows (see also Ref. 89):
where $ \mu S$, $ \mu K$, G, K, and $ d 0$ are constants of O(1).
The nondimensional liquid viscosity $ \mu L \u2009$ is defined as follows:^{97}
where $ \mu L$ is a constant of O(1). The nondimensionalization of initial effective viscosity and some parameters appearing in the energy equation are as follows:^{113,120,121}
where $ \mu L$, $ \zeta STM 1$, and $ \zeta STM 2$ are constants of O(1). The definition of the drag coefficient is as follows:^{114}
where Λ is a constant of O(1).
III. LOW FREQUENCY AND LONG WAVE
First, we focus on the low frequency and long wave that can be regarded as a weakly dispersive band in comparison with the analysis in Sec. IV and derive the KdVB equation. Our previous KdVB equation^{90} is successfully generalized by considering the additional effects, i.e., of drag force acting bubbles, thermal effect inside bubbles, and initial flow velocitiy; furthermore, the present KdVB equation is derived from twofluid model (not mixture) equations to incorporate α_{0} dependence.^{123} The derivation is performed by using the same mathematical procedure in our previous papers.^{89,90}
A. $ O ( \u03f5 )$: Linear propagation at nearfield
Substituting (25)–(39) into (1) and (4)–(7) and equating coefficients of like powers of ϵ in the set of resultant equations, we have the following set of linearized firstorder equations:
The following symbol is introduced as
and this also appears in Sec. IV, where
Eliminating α_{1}, $ u G 1$, $ u L 1$, and $ p L 1$ from (40)–(45), the single linear wave equation for firstorder perturbations of bubble radius, R_{1}, is obtained
with the phase velocity $ v p$,
with the constant coefficient,
The definition of the linear Lagrange derivative $ D / D t 0$ is
For simplicity, the initial velocities of both phases are assumed to be equal (i.e., $ u G 0 = u L 0 = u 0$). However, the perturbations of each velocity are not equal [i.e., $ u G i \u2260 u L i \u2009 ( i = 1 , 2 )$].
Setting $ v p = 1$ gives the explicit form of $ U \u2217$ as
with
Focusing on the rightrunning wave alone, a phase function $ \phi 0$ is defined as
All firstorder perturbations (i.e., α_{1}, $ u G 1 , \u2009 u L 1 , \u2009 p L 1$, and $ T G 1$) are expressed with regard to R_{1} as follows:
As a result, the dissipation and dispersion effects of waves and the effect of bubble translation and drag force do not appear in the leading order of approximation. The resultant physics of wave propagation are governed by the linear wave equation in (48) (including linear advection in $ D / D t 0$).
B. $ O ( \u03f5 2 )$: Nonlinear propagation at farfield
As in the case of Sec. III A, from the secondorder approximation, the following single inhomogeneous equation for R_{2} is derived:
where the explicit forms of the inhomogeneous terms $ K i$ ( $ i = 3 , 4 , 5 , 6$) are presented in Appendix A, and $ K 1$ and K_{2} are the same as the counterparts in Ref. 118. The solvability condition to the inhomogeneous equation in (60),^{97} K = 0 is imposed. By using (25), the independent variables $ \phi 0$, t_{1}, and x_{1} in (60) are restored into t and x
Finally, we have the KdVB equation with correction (i.e., fifth and sixth) terms
via the variable transformation
Equation (62) describes the competition of nonlinearity, dissipation due to five factors, dispersion, thermal dissipation, and dissipation due to drag force. Here, $ \Pi 0$ is the advection coefficient (i.e., the propagation speed in a moving coordinate), $ \Pi 1$ is the nonlinear coefficient for the second order nonlinear effect, $ \Pi 2$ is the dissipation coefficient, $ \Pi 3$ is the dispersion coefficient, $ \Pi 4$ is the dissipation coefficient due to thermal conduction,^{120} and $ \Pi 5$ is the nonlinear dissipation coefficient due to the drag force.^{114} The explicit form of the nonlinear coefficient $ \Pi 1$ is presented in Appendix B, and the other coefficients $ \Pi i \u2009 ( i = 0 , 2 , 3 , 4 , 5 )$ are given by
Here, $ \Pi 2$ and $ \Pi 3$ are positive, whereas $ \Pi 0$ and $ \Pi 1$ are negative. It should be noted that obtaining the analytical expression of the signs of $ \Pi 0$ and $ \Pi 1$ is challenging. In (64), the symbol C_{2} is introduced as
and this also appears in Sec. IV.
The three coefficients, the thermal conduction $ \Pi 4$, the drag force $ \Pi 5$, and other effects $ \Pi 2$, representing the dissipation are decomposed into
The first, second, third, fourth, and fifth terms in (65), $ \Pi 2 i$ ( $ i = 1 , 2 , 3 , 4 , 5$), describe the dissipation due to the liquid viscosity at the interface, liquid compressibility, shell viscosity, shell compressibility, and thermal effect (not thermal conduction), respectively. The detailed discussion is presented in Sec. V A.
Note that our previous KdVB equation^{89,90} neglected the thermal effect, drag force acting bubbles, bubble translation, and initial flow velocity.
IV. HIGH FREQUENCY AND SHORT WAVE
The most important difference between our previous^{90} and present papers is the focus on the high frequency and short wave that can be regarded as a strongly dispersive band in comparison with the analysis in Sec. III A. Weakly nonlinear modulation of the carrier wave is described by an envelope wave. Although some previous studies derived the NLS equation for the case of liquids containing uncoated bubbles (e.g., Refs. 97 and 113), here, for the first time, we successfully derive the NLS equation for liquids containing multiple coated bubbles. The focus on the wave frequency band near the resonance frequency of bubble oscillations will be a spark in medical applications.
A. $ O ( \u03f5 )$: Linear propagation of carrier wave at nearfield
As in the leading order of approximation in the low frequency case (Sec. III A), we first have the following set of linear equations:
Here, C_{1} is the same as that in (46). The difference between the present case (i.e., highfrequency short wave) and the previous case (lowfrequency long wave; Sec. III A) is the presence of the secondorder derivative of R_{1} on the lefthand side of (79), because the present case has a strong dispersive effect compared with the previous case [see (44)]. As in the previous case, although the initial velocities of both phases are assumed to be the same, the perturbations of each velocity are not equal [i.e., $ u G i \u2260 u L i$ ( $ i = 1 , 2 , 3$)].
Equations (75)–(80) are reduced to the following linear equation with a dispersion term:
The dispersion term (i.e., last term) on the lefthand side of (82) due to the dispersion term in the linearized Keller equation (79) did not appear in the previous case [see (48) in Sec. III A].
We assumed the solution of (81) as a quasimonochromatic wave train, which evolves into a slowly modulated wave packet^{97,116}
where A is the complex amplitude with a slow variation, which depends only on slow scales [i.e., t_{i} and x_{i} ( $ i = 1 , 2$)], and is obviously a constant in the nearfield characterized by fast scales (i.e., t_{0} and x_{0}); $ k ( \u2261 k \u2217 L \u2217 = k \u2217 U \u2217 / \omega B \u2217 )$ is the nondimensionalized wavenumber; and $ c . c .$ denotes the complex conjugate. Hereafter, two types of waves, $ e i \theta $ corresponding to the carrier wave and A to the envelope, are considered.^{97} The following linear dispersion relation connects Ω and k:
which is rewritten explicitly as
The positive Ω in (86) corresponds to the rightrunning carrier wave. The nondimensional phase velocity $ V p$ is readily obtained as
The nondimensional group velocity $ V g$ is also readily obtained as
We finally determine the typical propagation speed of a wave $ U \u2217$ by setting $ v p = 1$ under Ω = 1,^{97}
The typical length $ L * \u2261 U \u2217 T \u2217 = U \u2217 / \omega B \u2217$ is simultaneously determined.
with
As a result, the dispersion effect appeared in the nearfield, in contrast to Sec. III A.
B. $ O ( \u03f5 2 )$: Linear propagation of envelope wave at farfield I
From the secondorder of approximation, we have the following single inhomogeneous equation for R_{2}:
with
where M_{i} ( $ i = 1 , 2 , 3 , 4 , 5 , 6$) are the same as the counterparts in our previous case for uncoated bubbles without drag force.^{113} Substituting (83) and (91) into (93), the single inhomogeneous term M can be obtained as
with
where the real constant Γ is explicitly shown in Appendix B.
From the solvability condition for (96), the coefficient of $ e i \theta $ on the righthand side of (96) must be zero. Thus, we have a linear wave equation with a correction (i.e., third) term for the envelope wave A:
and its solution is obtained as follows:^{97}
with
Using (100) to set the secondorder equations, we have the other secondorder perturbations as follows:
with the new real constants e_{i} ( $ i = 1 , 2 , 3 , 4$) given by
Here, the explicit forms of the real constants are c_{i} ( $ i = 1 , 2 , 3 , 4 , 5$), d_{i} ( $ i = 1 , 2 , 3 , 4$), and $ c s$ are the same as those in our previous uncoated case.^{113}
Although the present result clarified the nonlinear propagation of carrier waves, the envelope wave was restricted to linear propagation. Furthermore, the wave propagation was not affected by the thermal effect, drag force, and bubble translation. We now proceed to the next order of approximation to elucidate the nonlinear propagation of envelope waves.
C. $ O ( \u03f5 3 )$: Nonlinear propagation of envelope wave and resultant NLS equation
The thirdorder of approximation clarifies the weakly nonlinear propagation process of the envelope wave. The single inhomogeneous equation for R_{3} is described as
where
with
Explicit forms of C_{5}, C_{6}, and C_{7} are not shown, because they do not affect the results of the present study. The inhomogeneous term N can be rewritten as
where $ \Lambda i \u2009 ( i = 1 , 2 , 3 )$ are the complex variables comprising A; $ \Lambda 1$, and $ \Lambda 2$ are not shown because they are not essential to deriving the NLS equation. Here, $ \Lambda 3$ is explicitly given as
and this should be zero from the solvability condition.
Finally, we can obtain the NLS equation with dissipation (i.e., fourth) and advection (i.e., second) terms
via the variable transformation
noting that complex amplitude A is assumed to be a sinusoidal solution [i.e., $ A = exp \u2009 ( i k \u0302 \xi \u2212 i \omega \u0302 \tau )$, where $ k \u0302$ and $ \omega \u0302$ are the frequency and wavenumber of the envelope wave, respectively].
The second, third, fourth, and fifth terms on the lefthand side of (114) represent the advection (i.e., wave propagation in moving coordinate), thirdorder nonlinear, dissipation, and dispersion effects of the envelope wave, respectively; ν_{02} represents the advection effect. The two advection coefficients, ν_{01} and ν_{02}, are given by
with the real constants
The explicit form of the nonlinear coefficient ν_{1} is presented in Appendix B. The dispersion coefficient ν_{3} is given by
The dissipation coefficient ν_{2} is represented by a linear combination of six factors as follows:
where
The detailed discussion is presented in Sec. V A.
V. RESULTS AND DISCUSSION
For the case of the KdVB equation, Figs. 3–5 show the advection, nonlinear, and dissipation coefficients $ \Pi 0 , \u2009 \Pi 1$, and $ \Pi 2$, respectively, vs the initial void fraction α_{0}. For the case of the NLS equation, Fig. 7 shows the two advection coefficients ν_{01} and ν_{02}, Fig. 8 shows the nonlinear coefficient ν_{1}, and Figs. 9 and 10 show the dissipation coefficient ν_{2} (since $ \nu 2 th$ is significantly large, the dissipations except for thermal dissipation, $ \nu 2 \u2212 \nu 2 th$ are shown in Fig. 9 and $ \nu 2 th$ is shown in Fig. 10). A quantitative discussion is presented using the physical properties of five wellknown UCAs, i.e., Albunex, Sonazoid, SonoVue, Levovist, and Optison, as summarized in Table I. In the following, uncoated bubbles (i.e., without shell) are defined by setting $ d 0 \u2217 = 0$.
UCA .  $ G S \u2217 \u2009 [ MPa ] $ .  $ \mu S \u2217 [ Pa \xb7 s ] $ .  $ d S 0 \u2217 \u2009 [ nm ] $ .  Reference . 

Sonazoid  52  0.99  4  80 
Albunex  15  0.05  15  70 
SonoVue  20  0.6  4  50 
Levovist  80  1.3  6  44 
Optison  20.7  1.7  7.5  69 
Note that the dependence of each coefficient on α_{0} exists qualitatively, as is clear from the explicit forms of coefficients in (64), (65), (70)–(74), (116), (117), (120)–(126), (B1), and (B10). However, as shown in Figs. 4–11, the dependence is quantitatively quite minimum for the case of low α_{0}.
A. Dissipation
1. Four factors: Importance of shell compressibility
First, we discuss the four factors of dissipation, i.e., the liquid viscosity, liquid compressibility, shell viscosity, and shell compressibility. Figures 6 and 11 represent each component of $ \Pi 2$ in the KdVB equation and ν_{2} in the NLS equation, respectively: (i) orange curves represent the liquid viscosity at the interface $ \Pi 21$ and ν_{21}, (ii) green curves represent the shell viscosity $ \Pi 23$ and ν_{23}, (iii) red curves represent the shell compressibility $ \Pi 24$ and ν_{24} and (iv) purple curves represent the liquid compressibility $ \Pi 22$ and ν_{22}.
The following tendency was observed for both the KdVB and NLS equations: The order of sizes of each dissipation factor differs depending on the UCAs. The dissipation due to the shell compressibility is always greater than that due to the shell viscosity. In contrast, the dissipation due to the liquid compressibility is always smaller than that due to the liquid viscosity. This reflects the importance of the consideration of shell compressibility.
The predominance of dissipation by the liquid phase (i.e., $ \Pi 21 + \Pi 22$ or $ \nu 21 + \nu 22$) or that by the shell (i.e., $ \Pi 23 + \Pi 24$ or $ \nu 23 + \nu 24$) is determined by the shell viscosity $ \mu S \u2217$, because $ \mu K \u2217 / \mu S \u2217$ is equal to $ K \u2217 / G \u2217$.^{68} Also, $ \mu K \u2217 = 0.7 \mu S \u2217$^{68} is used.
2. Thermal effect and drag force
The discussion of dissipation due to thermal effects and drag force is complex compared with the dissipation due to compressibility and viscosity (Sec. V A 1). Dissipation due to thermal effect is represented by both $ \Pi 25$ and $ \Pi 4$ in the KdVB equation and only $ \nu 2 th$ in the NLS equation. Dissipation due to drag force is represented by $ \Pi 5$ in the KdVB equation and ν_{25} in the NLS equation (pink curve in Fig. 11).
For the case of the NLS equation, the dissipation due to thermal conduction $ \nu 2 th$ is the largest, and that due to the drag force $ \nu 25$ is the smallest in ν_{2}. In the case of the KdVB equation, the dissipation owing to thermal effect $ \Pi 25$ is the smallest in $ \Pi 2$. Although ν_{2} is represented by a linear combination of every factor of dissipation, there are three different dissipation terms (coefficients $ \Pi 2 , \u2009 \Pi 4$, and $ \Pi 5$) in the KdVB equation in (62). It is worth noting that (i) $ \Pi 2$ is decomposed into five factors and (ii) although $ \Pi 2$ is the coefficient of the secondorder derivative, $ \Pi 4$ and $ \Pi 5$ are the coefficients without derivatives in (62).
The gas inside contrast agents can be regarded as isothermal behavior in clinical practice.^{40} In fact, the dissipation coefficient owing to thermal conduction, i.e., $ \Pi 4$ and $ \nu 2 th$, vanishes in the isothermal case. Furthermore, the difference of other components of dissipation coefficients between the isothermal and nonisothermal cases was negligible.
3. Comparison among UCAs
In the dissipation coefficients $ \Pi 2$ and ν_{2}, Optison with the largest shell viscosity $ \mu S *$ (see Table I) shows the largest dissipation and Albunex with the smallest $ \mu S *$ shows the smallest dissipation (Figs. 5 and 9, respectively). The order of size of each factor of $ \Pi 2$ and ν_{2} is the same in all the UCAs notably, although the effect of drag force is included in ν_{2} and not in $ \Pi 2$, this order is not affected.
In comparison with the case of uncoated bubbles, it turns out to be increased in all the UCAs, as summarized in Table II for the case of the KdVB equation and Table III for the case of NLS equation, including the advection and nonlinear coefficients.
UCA .  Albunex .  Sonazoid .  SonoVue .  Levovist .  Optison . 

Advection $ \Pi 0$  0.41  0.69  0.71  0.49  0.63 
Nonlinearity $ \Pi 1$  0.70  0.88  0.90  0.73  0.42 
Dissipation $ \Pi 2$  1.4  5.6  3.9  6.9  14 
UCA .  Albunex .  Sonazoid .  SonoVue .  Levovist .  Optison . 

Advection $ \Pi 0$  0.41  0.69  0.71  0.49  0.63 
Nonlinearity $ \Pi 1$  0.70  0.88  0.90  0.73  0.42 
Dissipation $ \Pi 2$  1.4  5.6  3.9  6.9  14 
UCA .  Albunex .  Sonazoid .  SonoVue .  Levovist .  Optison . 

Advection $ \nu 01$  9.3  5.7  8.6  0.64  17 
Advection $ \nu 02$  66  62  66  38  91 
Nonlinearity $ \nu 1$  0.34  0.71  0.76  0.39  0.59 
Dissipation $ \nu 2 \u2212 \nu 2 th$  1.2  4.7  3.4  5.9  12 
Dissipation $ \nu 2$  0.086  0.69  0.79  0.16  0.45 
UCA .  Albunex .  Sonazoid .  SonoVue .  Levovist .  Optison . 

Advection $ \nu 01$  9.3  5.7  8.6  0.64  17 
Advection $ \nu 02$  66  62  66  38  91 
Nonlinearity $ \nu 1$  0.34  0.71  0.76  0.39  0.59 
Dissipation $ \nu 2 \u2212 \nu 2 th$  1.2  4.7  3.4  5.9  12 
Dissipation $ \nu 2$  0.086  0.69  0.79  0.16  0.45 
B. Nonlinearity
To discuss the size of nonlinearity, we observed the absolute value of the nonlinear coefficients, i.e., $  \Pi 1 $ in Fig. 4 and $  \nu 1 $ in Fig. 8, and compared the size of all the UCAs (except for the uncoated bubbles). Then, the same order was observed in both $  \Pi 1 $ and $  \nu 1 $. Furthermore, Albunex, having the highest shell thickness $ d S 0 *$, has the smallest $  \Pi 1 $ and $  \nu 1 $.
C. Advection
We discuss the advection coefficients, $  \Pi 0 $ in the KdVB equation (Fig. 3) and ν_{01} and ν_{02} in the NLS equation (Fig. 7), among the five UCAs. In the range of $ 10 \u2212 8 \u2264 \alpha 0 \u2264 10 \u2212 5$, Albunex has the smallest $  \Pi 0 $; however, Levovist has the smallest ν_{01} and ν_{02}. Conversely, SonoVue has the largest $  \Pi 0 $, although Optison has the largest ν_{01} and ν_{02}.
D. Case of high void fractions
Figure 12 provides an example of dependence of coefficients $ \Pi 0$ and ν_{1} on the initial void fraction α_{0}. Note that although Fig. 12 shows only a visible range of α_{0} dependence (i.e., high void fraction case), α_{0} dependence qualitatively exists regardless of the size of α_{0}.
Unlike the case of $ \Pi 0$, although the dependence of the other coefficients on α_{0} is quite small, noting that slight dependence of α_{0} exists. In fact, ν_{1} for large α_{0} slightly depends on α_{0} [Fig. 12(b)].
E. Dependence of physical properties
Figure 13 shows the dependence of physical properties (i.e., shell rigidity $ G S *$, surface tension $ \sigma 1 *$, and shell thickness $ d S 0 *$) on the nonlinear coefficients $ \Pi 1$ and ν_{1}. Notably, the range of physical properties (i.e., horizontal axis) does not correspond to a specific contrast agent; however, it corresponds to five typical contrast agents shown in Table I.
Clearly, the nonlinear coefficient strongly depends on these physical properties. As values, the degree of change is larger for ν_{1} than for $ \Pi 1$. Nevertheless, as is clear from explicit forms in (64)–(74) and (115)–(126), the remaining coefficients are also dependent on these physical properties.
VI. CONCLUSIONS
Apart from our previous research,^{89,90} no theoretical studies have been reported on the dynamics of multiple coated bubbles. In this paper, we theoretically investigated the weakly nonlinear propagation process of onedimensional ultrasound in a compressible flowing liquid uniformly containing multiple microbubbles coated by a viscoelastic body that obeys the Kelvin–Voigt model, especially focusing on (i) both lowfrequency long wave (Sec. III) and highfrequency short wave (Sec. IV) bands and (ii) the incorporation of the compressibility of the shell.^{68} As a result, the KdVB equation for the low frequency case and NLS equation for the high frequency case, which describe nonlinear evolutions of ultrasound, were successfully derived in a unified manner. Importantly, for the first time, we successfully derived the NLS equation for the case of multiple coated bubbles. Moreover, the compressibility of the liquid and shell, thermal conductivity, drag force acting on bubbles, and translation of bubbles, which were neglected in our previous study,^{89} were incorporated into the present theory. Especially, in the incorporation of shell compressibility, the updated equation of motion for a single coated bubble with shell compressibility^{68} was successfully extended to the case of multiple coated bubbles.
Focusing on the coefficients of advection, nonlinear, and dissipation terms in the resultant KdVB and NLS equations for five typical UCAs (i.e., Albunex, Sonazoid, SonoVue, Levovist, and Optison) and uncoated case, the following results were obtained: (i) The shell compressibility increased the dissipation effect but decreased the nonlinear effect; (ii) The dissipation due to thermal conduction was the largest of all the dissipation factors for the NLS equation; (iii) By comparing each dissipation factor for typical UCAs, both similarity and difference were observed; (iv) Nonnegligible contribution of shell compressibility to ultrasound dissipation was identified although it depends on the types of contrast agents, e.g., the dissipation due to shell compressibility is approximately ten times larger than that of liquid compressibility for the case of Optison [Fig. 6(e)]; (v) Focusing on the dependence of nonlinear and advection effects on the initial void fraction, a similar tendency was observed for both the KdVB and NLS equations. From these results, we elucidated a basic understanding of the effect of the shell of coated bubbles on the dissipation, nonlinear, and advection effects of weakly nonlinear ultrasound, for both low and high frequencies. Particularly, we discovered that the effect of shell compressibility, which has long been neglected, is crucial.
Finally, as with any study, the limitations of the present study should be acknowledged and addressed in future work: (i) the buckling of the shell;^{77} (ii) anisotropy of the shell;^{68} (iii) asymmetric oscillations of nonspherical coated bubbles; (iv) the effect of viscoelasticity of liquids such as tissue;^{125} (v) bubble–bubble interaction and the resultant change in eigenfrequency;^{58,102} (vi) pressure dependence of the eigenfrequency of bubble oscillations;^{81,82} and (vii) discussion of the nonlinear dissipation effect^{63,64,66,124} for a high pressure amplitude case. Also, the KdVB and NLS equations derived here will be numerically solved under practical conditions conducive to diagnostic (and/or therapeutic) applications of ultrasound.
ACKNOWLEDGMENTS
This work was partially carried out with the aid of the JSPS KAKENHI (Nos. 18K03942 and 22K03898), based on the results obtained from a project subsidized by the New Energy and Industrial Technology Development Organization (NEDO) (No. JPNP20004), and Ono Charitable Trust for Acoustics. We would like to thank referees for their valuable comments and suggestions, and Editage (www.editage.com) for English language editing.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Tetsuya Kanagawa: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (lead); Project administration (lead); Supervision (lead); Validation (equal); Writing – original draft (lead); Writing – review & editing (lead). Mitsuhiro Honda: Data curation (equal); Formal analysis (equal); Investigation (equal); Validation (equal); Writing – original draft (supporting); Writing – review and editing (supporting). Yusei Kikuchi: Data curation (supporting); Formal analysis (equal); Investigation (equal); Methodology (supporting); Validation (supporting).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX A: INHOMOGENEOUS TERMS
The inhomogeneous terms K_{i} ( $ i = 3 , 4 , 5 , 6$) in (60) are given by
where C_{2} was defined in (69).
APPENDIX B: EXPLICIT FORM OF SOME COEFFICIENTS
The nonlinear $ \Pi 1$ coefficient of the KdVB equation in (62) is given as follows:
where explicit forms of $ k i$ ( $ i = 1 , 2 , 3 , 4 , 5 , 6$) are the same as counterparts in Ref. 46.
where
Finally, the nonlinear coefficient of the NLS equation ν_{1} in (114) is given as follows:
where