Nonlinear ultrasound in liquid containing multiple coated microbubbles: effect of buckling and rupture of viscoelastic shell on ultrasound propagation

With promising applications in medical diagnosis and therapy, the behavior of shell-encapsula-ted ultrasound contrast agents (UCAs) has attracted considerable attention. Currently, second-generation contrast agents stabilized by a phospholipid membrane are widely used and studies have focused on the dynamics of single phospholipid shell-encapsulated microbubbles. To improve the safety and the efficiency of the methods using the propagation or targeted ultrasound, a better understanding of the propagation of ultrasound in liquids containing multiple encapsulated microbubbles is required. By incorporating the Marmottant–Gompertz model into the multiple scale analysis of two-phase model, this study derived a Korteweg–de Vries–Burgers equation as a weakly nonlinear wave equation for one-dimensional ultrasound in bubbly liquids. It was found that the wave propagation characteristics changed with the initial surface tension, highlighting two notable features of the phospholipid shell: buckling and rupture. These results may provide insights into the suitable state of microbubbles, and better control of ultrasound for medical applications, particularly those that require high precision.

For applications in medical treatment, a comprehensive understanding of the interaction of encapsulated bubbles with ultrasound is required. Several models have been proposed to explain the effect of shells on bubble oscillation, which encapsulate the gas core and provide important functions, such as improvement of stability against Laplace-pressure driven dissolution, gas diffusion, and coalescence. For example, de Jong proposed a pioneering model for viscoelastic shells [31][32][33][34], Church et al. [35] presented a model for a viscoelastic shell that also incorporates shell thickness, and Hoff et al. constructed a model by considering the limit of shell thickness as zero [35][36][37]. Furthermore, Chatterjee and Sarkar [38] created a model by assuming that shells behave as Newtonian viscous fluids. Marmottant et al. [39] proposed a model that incorporates the buckling and rupture phenomena, observed in phospholipid monolayer shells (it should be noted that the "rupture" phenomenon defined by Marmottant is not the irreversible collapse in solid mechanics). Recently, a nonlinear viscosity model was suggested by Doinikov et al. [40] and a model incorporating the effect of shell compressibility and anisotropy was developed by Chabouh et al. [41]. Extensive reviews on comparisons and discussions of the models have also been published [42][43][44][45]. Under simplifications and assumptions, the existing models show reasonable agreement with experimental results, mostly in the linear regime. However, nonlinear behavior of encapsulated bubbles is complicated, and attempts have been made to better understand these behaviors [46][47][48][49][50][51][52][53][54][55][56]. As Doinikov et al. [57] pointed out, available shell models do not possess required predictive capability for wide range of conditions. This problem may be solved with researches at the molecular scale, using theoretical [58,59] and simulation [60][61][62][63] approaches.
Research regarding single encapsulated bubbles in pressure waves is a foundation for studies on ultrasound propagation in a liquid with a large number of encapsulated bubbles. Accordingly, several models have been proposed as an extension of the analysis by van Wijngaarden and Caflisch, and multiple imaging techniques have been developed [64][65][66][67]. Ma et al. [68] derived a nonlinear evolution equation for ultrasound propagation in liquids containing multiple encapsulated bubbles, and Xia [69] theoretically studied the attenuation coefficient of ultrasound propagation. These studies, however, were limited to linear case and cannot holistically reflect nonlinear effects such as rapid attenuation change above the Blake threshold or the dependence on pressure of attenuation and sound speed. To address such limitations, in his pioneering paper, Louisnard [70,71] constructed a mechanical energy balance equation from fully nonlinear Caflisch equations. Here, energy loss was computed numerically by simu-lating bubble radial dynamics Rayleigh-Plesset equation. Consequently, a nonlinear Helmholtz equation for wave propagation in liquids with uncoated bubbles was derived. Accordingly, pressure dependent attenuation was derived from the imaginary part of wave number. This equation is relatively easier to solve than the fully nonlinear Caflisch equations and can predict more realistic attenuation and acoustic pressure values compared to the fully linearized models. Similarly, Jamshidi et al. [72] considered compressibility of the liquid using Keller-Miksis equation, and the additional attenuation effect of acoustic radiation was obtained in addition to small modifications to thermal and viscosity attenuation terms. Later, the Jamshidi model was shown to have non-physical values and a critical modification to the calculation of the damping terms was provided by Sojahrood et al. [73]. These values had better agreement with the linear model and resolved the non-physical values. Additionally, they extended previous studies and were the first to introduce a full nonlinear model capable of simultaneously calculating the pressure dependence of sound speed, which affected local acoustic pressure amplitude, and attenuation [74]. This model was verified numerically [75] and in the first controlled observation of the pressure dependence of sound speed and attenuation for coated bubbles [76]. To describe wave propagation in a liquid containing a high number of viscoelastic shellencapsulated microbubbles, we have adopted another approach using two-phase model. Consequently, the Korteweg-de Vries-Burgers (KdVB) equations were obtained [77][78][79]. Although our studies are limited to small pressure amplitude, the proposed models can account for the nonlinear propagation of waves, a characteristic neglected in models based on the Helmholtz model.
The buckling and rupture of the shell were first modeled by Marmottant et al. [39] for large-amplitude oscillations of coated bubbles. Researchers using this model have reported that it can predict important behaviors exhibited by phospholipid monolayer shells. De Jong et al. [80] observed a highly nonlinear response of phospholipid-coated UCA termed "compression-only" behavior, wherein microbubbles compressed but barely expanded beyond its initial radius. Later, Emmer et al. [81] linked the "compression-only" behavior to the enhanced second harmonic behavior, and Sijl et al. [48] theoretically showed that the subharmonic behavior and its threshold pressure can be explained through the value and rate of change of elasticity with radius, which are important features of the Marmottant model. The Marmottant model has been demonstrated to have predicted the intensification of the generation of 1/2 order subharmonics via simulation and experiment on bubbles in the buckling region [49]. Several researchers have reported the change in resonance curves and onset of vibrations [82,83]. These two prominent nonlinear effects were explained later using the Marmottant model by Overvelde et al. [47]. Recently, the generation of higher order subharmonics (e.g., 1/3, 1/4) has been demonstrated by analyzing bifurcation structure and verified experimentally by Sojahrood et al. [84]. Since many approved UCAs for human use such as Definity (2001), SonoVue (2001), and Sonazoid (2007) are coated by phospholipid shells, the incorporation of the two physical properties (buckling and rupture) is necessary to model the interaction between ultrasound and phospholipid shell encapsulated microbubbles.
The aim of the present paper is to extend previous studies on ultrasound propagation in liquids containing multiple microbubbles encapsulated with viscoelastic shell [77] or with compressible viscoelastic shell [78,79], by incorporating buckling and rupture phenomena. To the best of our knowledge, this is the first theoretical study that considers effect of buckling and rupture of phospholipid shell on ultrasound propagation in a liquid containing multiple encapsulated bubbles that includes nonlinearity parameters through twophase model. Instead of the original model proposed by Marmottant et al. [39], we consider the Marmottant-Gompertz model [85] to describe the behavior of the surface tension of the phospholipid shell, as the latter offers important features relevant to our study (as described in Sect. 2.2). Following the method described in our previous papers [86,87], we derived a KdVB equation as a physico-mathematical model for ultrasound propagation.
The remainder of this paper is organized as follows. Section 2 introduces basic equations for bubbly liquids based on a two-fluid model [88] and lipid-encapsulated bubble dynamics including shell buckling and rupture. Section 3 presents derivation of linear propagation for the first-order problem and a KdVB equation for the second-order problem. Section 4 presents a parametric analysis conducted to explore the effect of the initial void fraction and initial surface tension on the characteristics of the propagation (i.e., advec-tion, nonlinearity, attenuation and dispersion). Furthermore, the results obtained after the analysis proposed by Katiyar and Sarkar [50], with the upper and lower limits of the surface tension removed, are presented. The limitations of derived model are then discussed. Finally, Sect. 5 presents the conclusions of this paper.

Problem statement
We consider a weakly nonlinear (i.e., finite but small amplitude) propagation of plane and progressive pressure waves radiating from a source in a bubbly liquid.
In this study, the following assumptions were made for simplicity: Except for assumption (xii), which is the main topic of this study, the other assumptions are identical to those in our previous papers [86][87][88]. The buckling and rupture of the shell were modeled using the Marmottant-Gompertz model, first suggested by Gümmer et al. [85]. In this study, we use a two-fluid model [88,89] structured by a set of basic equations. This set includes the following equations: (i) Mass conservation law in the gas phase (ii) Mass conservation law in the liquid phase (iii) Momentum conservation law in the gas phase (iv) Momentum conservation law in the liquid phase where α is the void fraction (the gas fraction) (0 < α < 1), ρ * is the density, u * is the velocity, p * is the pressure, and the subscripts G and L indicate the volume-averaged variables in the gas and liquid phases, respectively. The righthand sides of Eqs. (3) and (4) show the interfacial momentum transport, denoted as F * , following the model of virtual mass force in a compressible liquid [89][90][91] DGρ * L Dt * , (5) For spherical bubble case, the coefficients β 1 , β 2 and β 3 can be set to 1/2. We refrained from explicitly using these values to present the contribution of each term to the result. The total derivatives are defined as follows: (v) Modified Rayleigh-Plesset equation for spherical oscillations of bubbles in a slightly compressible liquid [39] where R * is the bubble radius, ρ * L0 is the liquid density in unperturbed state, c * L0 is the initial sound velocity in the liquid phase and P * is the difference between volume-averaged liquid pressure and liquid pressure at bubble surface. Equations (1)-(5) and (7) are closed using the following equations: (vi) Tait equation of state for liquid (vii) Polytropic equation of state for gas (viii) Conservation law of mass inside a bubble where p * L0 and R * 0 are the liquid pressure and bubble radius, respectively, in the initial undisturbed state; p * G0 and ρ * G0 are the gas pressure and density inside the bubble in the initial state, respectively; γ is the polytropic exponent; and n is the material constant (e.g., n = 7.15 for water).
Here, μ * is the liquid viscosity, and κ * s is the shell dilatational viscosity derived for shells with small but finite and constant thickness. To incorporate the buckling and rupture phenomena in our study, instead of using the original model of surface tension proposed by Marmottant et al. [39], in which there are two discontinuities at R * = R * buck and R * = R * rupt , we use the continuous Marmottant-Gompertz model constructed by Gümmer et al. [85], the detailed form of this model is given in "Appendix A" (Fig. 2).
The Marmottant model can simulate certain nonlinear characteristics such as "compression-only" and subharmonic behaviors of phospholipid-coated microbubbles [43][44][45][47][48][49]. The Marmottant surface tension model has two important features: the first are an upper limit (ruptured state) and a zero lower limit (buckled state). The effects of these limits on the excitation threshold for subharmonic generation were studied by Katiyar and Sarkar [50]. The second is the rapid change of the elastic coefficient χ * between the elastic-buckled state and elastic-ruptured state with two singularities of dχ * /dR * at R * = R * buck and R * = R * rupt . These two singularities make the Marmottant model sensitive to the time step used in numerical methods [48,85]. To resolve this sensitivity, Silj et al. [48] proposed two additional quadratic crossover functions to smoothen the curve and control the two transition points. The Marmottant-Gompertz model proposed by Gümmer et al. [85] is based on the Gompertz function. This model can smoothen the curve and eliminate the singularities using the same set of initial parameters (i.e., R * 0 , χ * and σ * 0 ), while retaining the two significant features of the Marmottant model. The elimination of these two singularities allows for the calculation of the continuous first and second derivatives of the function. Accordingly, the linearization for theoretical analysis can be obtained. Furthermore, without the singularities, the sensitivity to the time step in numerical calculations is dismissed. As an inherent characteristic of the Gompertz function, there is a more rapid change of surface tension near buckled radius than near ruptured radius, providing good qualitative approximation to the behavior of the lipid monolayer observed in experiments. Additionally, the Marmottant-Gompertz equation is shown to have no significant difference to the original Marmottant model for acoustic emissions and overall dynam-ics of the bubble while demonstrating smoother bubble radius transition in and out of the elastic regime. The result also has good agreement to in vitro experiments under similar conditions. Readers are referred to Gümmer et al. [85] for further discussion of the Marmottant-Gompertz model.

Multiple-scale analysis
For weakly nonlinear problems, the nonlinear effect becomes apparent at a large distance from the sound source relative to the wavelength. In these problems, different time and length scales are considered. First, the independent variables are nondimensionalized as where T * and L * are the typical period and wavelength of the propagating wave, respectively. Then, the nearfield (i.e., the temporal and spatial scales of O (1)) is defined as and the far-field (i.e., the temporal and spatial scales of O(1/ )) is described as where is the nondimensional wave amplitude under the assumption (0 < 1). Using Eqs. (13) and (14) with the derivative-expansion method [92], the differential operators are expanded as follows: Further, the dependent variables are nondimensionalized and expanded in the power series of : where α 0 is the initial void fraction and U * is the typical wave propagation speed, which is related to the wavelength L * and period T * through U * ≡ L * /T * . The nondimensional pressures for the gas and liquid phases in the unperturbed states p G0 and p L0 are defined as: The liquid viscosity μ * , shell viscosity κ * s and initial surface tension σ * 0 are nondimensionalized as For simplicity, we define the following two parameters: Using Marmottant-Gompertz model, the surface tension σ * can be expanded as follows: where the explicit forms of N 1 , N 22 and N 21 are Substitute Eq. (28) into the modified Rayleigh-Plesset equation (Eq. (7)), and follow the analysis of Van der Meer et al. [82], the eigenfrequency of the bubble is obtained through linearization: We focus on the long-range propagation of nonlinear waves in low-frequency and long-wavelength bands. In this case, the appropriate scaling relations [86] are where V, and are of order unity, and ω * is the typical angular frequency.  4) and (7) and collecting the 1 terms, the following set of linearized equations for the first-order problem is obtained: (ii) Mass conservation law in the liquid phase (iii) Momentum conservation law in the gas phase (iv) Momentum conservation law in the liquid phase (v) Modified Rayleigh-Plesset equation The set of linearized equations, i.e., Eqs. (35)-(39) contains five dependent variables (i.e., α 1 , p L1 , u G1 , u L1 and R 1 ). By eliminating α 1 , p L1 , u G1 and u L1 , the linear wave equation for the first-order perturbation of the bubble radius R 1 is derived as where the phase velocity v p is given by Equation (40) indicates a linear and nondispersive wave motion, described in terms of t 0 and x 0 , while Eq. (41) shows a proportional relation between the phase velocity v p and 1/ √ α 0 (1 − α 0 ), a feature similar to the classical speed of sound in bubbly liquids [14,15]. The expressions in Eqs. (35)-(41) are identical to some of our previous results for the uncoatedbubble case [86,87] but are different from other results [78,[93][94][95] that incorporate effects such as those of polydispersity, thermal effect, and drag force as well as results based on other shell models (the Church-Hoff and Chabouh models). This similarity indicates that if the same set of nondimensionalized parameters, i.e., and is used, the derived solution of the first-order problem will be identical to that of the uncoated bubble. However, because the expression of ω * B changes with the contribution of σ * 0 and χ * , the behavior of the solution is different. By substituting the definitions of p G0 in Eq. (22), in Eq. (33), and in Eq. (34) into Eq. (41), we formulated a typical propagation speed U * as follows: Focusing on the right-running wave in the leadingorder of approximation, a phase function ϕ 0 for an arbitrary value of v p can be introduced as follows: By setting Rewriting Eqs. (35)-(39) using ϕ 0 and integrating them with respect to ϕ 0 , the other first-order perturbations α 1 , u G1 , u L1 and p L1 can be expressed in terms of f (ϕ 0 ): , Here, the integration constants are omitted because of the boundary conditions at x 0 → ∞ where the bubbly liquid is uniform and at rest. Furthermore, from the relations in Eq. (45), it is apparent that the firstorder perturbations in the near field characterized by t 0 and x 0 (i.e., α 1 , u G1 , u L1 , p L1 and R 1 ) are governed by Eq. (44).

Second order of approximation
Following the procedure presented in Sect. 3.1, but for the 2 terms, we derive a single inhomogeneous equation for R 2 : The explicit form of inhomogeneous terms K i (i = 1, 2, 3, 4, 5) is provided in "Appendix B". From the solvability condition of the inhomogeneous equation in Eq. (46), which is also the non-secular condition for the asymptotic expansions in Eqs. (16)- (21), we have K = 0, that is, Finally, through variable transformation, the KdVB equation for nonlinear propagation in the far-field can be obtained as follows: Here, 0 , 1 , 2 and 3 represent the advection, nonlinear, attenuation, and dispersion effects, respectively. The explicit forms of the coefficients i (i = 0, 1, 2, 3) can be expressed as where the explicit forms of k i (i = 1, 2, 3, 4, 5) can be written as: From this result, it is apparent that 0 is negative while 2 and 3 are positive. However, the sign of 1 is challenging to analyze; this case is evaluated with some examples in Sect. 4. A comparison with our previous results for the uncoated bubbles reveals certain discrepancies in the expressions of k 5 of 1 and 2 . The difference in k 5 of 1 originates from the complicated Marmottant-Gompertz surface tension model. Moreover, the difference in 2 arises from the incorporation of the shell dilatational viscosity κ * s , which leads to the additional nondimensional term κ s , and the use of the modified Rayleigh-Plesset equation in Eq. (7) results in the omission of term − 3 V /(6α 0 2 ). Notably, if we had followed the procedure presented in our previous paper using the Keller

Effect of shell
In this section, the effect of the initial surface tension on the coefficients of the derived KdVB equation is first investigated. For all the numerical calculations presented below, we set v p = 1 because our method cannot predict the typical speed of sound in a bubbly liquid U * . Figure 3 shows the coefficients versus the initial void fraction for the uncoated bubbles case (red) and three cases for lipid-shell encapsulated bubbles (i.e., near the buckled state, in elastic state, and near the ruptured state). The parameters used for each case are listed in Table 1. It should be noted that, in our calculation, the elasticity coefficient appears in denominators in certain relations and therefore zero value will make the program yields error. Here, since we would like to use the same model for every case (even for uncoated case) to avoid any unnecessary differences, the shell elasticity for uncoated bubble case is nonzero. Even though uncoated bubble case has shell elasticity, the contribution is neglegible and the result for uncoated bubble case is the same with our previous studies [86]. In this range of parameters, the sign of 1 is positive for all the cases. Furthermore, 0 diverges rapidly as gas fraction goes to zero. This behavior arises from the divergence of typical wave velocity U * chosen for v p = 1. For all terms except attenuation term 2 , the difference between the ruptured case (green) and the uncoated case (red) is relatively small. Moreover, the shell enhances all the effects of wave propagation, and for a larger void fraction, the absolute values of the coefficients decrease, although at different rates. It should be noted that, in medical applications and experimental studies, the void fraction is of the order of 10 −6 . Although a higher void fraction is difficult to achieve in experiments and applications, we draw graphs at very high void fractions to illustrate the behavior of the solutions.
To investigate the effect of the initial surface tension σ * 0 on wave propagation, graphs of σ * 0 versus i for different initial void fractions α 0 are plotted as solid lines in Fig. 4. The dashed line in Fig. 4 depicts the result of the coefficients i subtracted by its counterpart in the uncoated bubbles case (at α 0 = 0.05). Since for  Fig. 4 is the α 0 = 0.05 solid curve that shifts vertically and demonstrates the effect of the shell. From Fig. 4, we can deduce that the effect of the shell is relatively insignificant for the advection effect 0 near buckled and ruptured state, and dispersion effect 3 . However, for the attenuation effect 2 , the effect of the shell near buckled and ruptured state is significant. Moreover, for the nonlinear effect 1 , the effect of the shell is the main component. Another observation is that for the attenuation term 2 and dispersion term 3 , the coefficients vary insignificantly in the linear regime. In the linearbuckled and linear-ruptured transition regimes, σ * 0 N 1 and σ * 0 N 21 (i.e., the first-and second-order derivatives of the surface tension at R * /R * 0 = 1) change rapidly, which may explain this tendency. From the expression of the coefficients, it is apparent that only N 1 affects the values of 0 , 2 , 3 through the change in ω * B , as discussed in Sect. 3. Therefore, these values were determined by N 1 . The coefficient N 21 , however, implicitly defines the shape of the graphs because σ * 0 N 21 shows the rate of change of σ * 0 N 1 . Physically, σ * 0 N 1 shows shell stiffness. For regions with high σ * 0 N 1 , the shells are expected to be stiffer, making bubbles oscillation more difficult and impeding the rate of radius change of the bubbles. This tendency is reflected through the attenuation term 2 . As shown in Fig. 4c, the region of low absolute value of 2 corresponds to the region where shell stiffness is large. It should be noted that in our study, only the attenuation effect of shell viscosity and liquid near the bubbles are considered while other attenuation mechanisms such as thermal or radiation damping are neglected.
The behavior of the nonlinear term 1 is more complicated and is characterized by the absence of symmetry as in the other coefficients. Note that in the transition regime between linear-ruptured, an insignificant uphill region can be observed. These two features arise from the fact that the values of 0 , 2 and 3 are determined by N 1 and from the symmetry of N 1 , the graphs of these coefficients are symmetrical. However, the expression of component k 5 of 1 indicates that there is also a contribution of N 21 . In the transition regime from the buckled state to the elastic state, σ * 0 N 1 changes rapidly from zero to χ * (i.e., σ * 0 N 21 is positive), whereas in the transition regime from the elastic state to the ruptured state, σ * 0 N 1 changes rapidly back to zero (i.e., σ * 0 N 21 is negative). These rapid rates of change result in the dominance of σ * 0 N 21 in the nonlinear coefficient. For reference, we drawn the ratio of (σ * 0 N 21 )/(σ * 0 N 1 ) and present the result in Fig. 5.

Effect of upper and lower limits
To investigate the effect of the upper and lower limits of surface tension, each of the limit is removed by modifying the original Marmottant-Gompertz model. First, the upper limit is removed by setting asymptote value (i.e., free gas-water surface tension σ * c ) to a higher value while maintaining the lower limit. Second, the lower limit is removed while maintaining the upper limit through two steps: (i) The graph is translated down A unit (i.e., σ * = −A + σ * ); (ii) The values of the surface tension (i.e., σ * c and σ * 0 ) are set to a new value according to A (i.e., σ * c to σ * c = σ * c + A and σ * 0 to σ * 0 = σ * 0 + A). It should be noted that the change in the values of surface tension only occurs in the expansion terms of surface tension σ * and therefore, does not affect the zero-order term. In addition, the value of χ * is varied to obtain a good fit to the elastic region of the original model. The new parameters are listed in Table 2.
The modified models, after the removal of the limits, are shown in Fig. 6. Notably, since the Marmottant-Gompertz model is based on the Gompertz function, it uses three parameters: an asymptote parameter, a displacement along the horizontal-axis parameter, and a growth rate parameter. The Marmottant-Gompertz model varies these three parameters through (σ * 0 , σ * c and χ * ), and only two parameters can be set since we have already fixed σ * 0 . This results in high values of χ * used in this section. Moreover, we cannot find a set of parameters that can remove both the upper and lower limits with which the linear region is comparable with the original model. Therefore, a direct comparison for the case of the surface model without upper and lower limits cannot be achieved in this study, and readers are referred to our previous papers [77][78][79] for the results of the analysis using the Church-Hoff shell model.
In Figs. 7 and 8, the results for the upper-limitremoved case and original case are illustrated as dashed curves and solid curves, respectively. In Fig. 7, there is a shift in the ruptured case toward a higher absolute value for 0 , 1 , and 3 , whereas for the attenuation term 2 , the ruptured case moves toward a smaller absolute value. The relatively small differences in 0 , 2 , and 3 for elastic lines and buckled lines between two cases may arise from the small discrepancies between the upper-limit-removed model and the original model, specifically the value of σ * 0 N 1 (i.e., the slope of the graph). However, the considerable difference in the nonlinear coefficient 1 might be the result of two factors: the difference in σ * 0 N 1 and the difference in σ * 0 N 21 . Because the value and variation of σ * 0 N 1 are insignificant, the difference in σ * 0 N 21 for each case is dominant in determining the difference in 1 , as mentioned in Sect. 4.1. In Fig. 8, the case of α 0 = 10 −6 , 0.01 is omitted to highlight the results. The significant shape change for all coefficients can be explained as follows. Since the differences in σ * 0 N 1 and the differences in σ * 0 N 21 between the elastic and ruptured regimes are small for upper-limit-removed case, the curves of coefficients flatten in the region with a relatively large value of σ * 0 . This illustrates our discussion of the effects of σ * 0 N 1 and σ * 0 N 21 on the wave propagation coefficients.
In Figs. 9 and 10, the results for the lower-limitremoved case and original case are illustrated as dashed curves and solid curves, respectively. In Fig. 9, the buckling curves vary significantly while the elastic and rupture curves do not vary significantly for 0 , 2 and 3 . This is expected owing to the dependence of those curves on the slope of σ * . By removing the lower limit, the slope of σ * (σ * 0 N 1 ) near the buckled region becomes larger than the original case, explaining the lower dissipation absolute value for the buckling curve, since a stiffer shell makes bubbles harder to compress. The notable differences of the 1 curves can be explained similarly with the upper-limit-removed case. In Fig. 10, the i (σ * 0 ) curves are plotted. A similar trend with the upper-limit-removed case can be observed. However, in this case the curves flatten in the region with a relatively large σ * 0 value.   [85], and a same set of parameters was used to define and calculate the surface tension curve. Figure 5 shows that in Marmottant-Gompertz model, σ * 0 N 21 is an order of magnitude larger than σ * 0 N 1 . However, Sijl et al. [48] demonstrated that σ * 0 N 21 must be at least three orders of magnitude larger than σ * 0 N 1 to illustrate the abrupt elasticity change of the collapsing phospholipid monolayer. This characteristic is crucial for explaining the subharmonic behavior of phospholipid-encapsulated microbubbles [48] and corresponds with the experimentally determined value. From our analysis, it is noted that a change in σ * 0 N 21 directly changes the value of the nonlinear coefficient 1 and affects the tendency of i (σ * 0 ) graphs. The alternative model, which can control the slope of surface tension curve and its rate of change, was proposed by Sijl et al. [48]. They introduced a new parameter (i.e., ζ to manipulate σ * 0 N 21 ), and the model is shown as follows: In this study, we did not use the model proposed by Sijl et al. [48] since it requires the analysis to be conducted for each separated region; this is because, in this model, the surface tension is divided into different The dashed curves represent the case where the lower limit was removed, and the solid curves represent the original model. The other parameters are same with Fig. 3. The high void fractions are graphed to demonstrate the behavior of the solutions regions, which are expressed using different equations. Furthermore, in each separated region, there may not be sufficient information about the other regions. For example, in the region of R * elas2 , there is no information about the upper limit (i.e., σ * c ). Additionally, the KdVB equation derived in this study neglects the contributions of several factors. For instance, the effect of the bubble-bubble interaction, [96][97][98] that changes the resonance frequency of the bubbles, can alter the coefficients of the KdVB equation. Such interactions have been shown to affect the oscillation of the microbubbles, resulting in milder pulsation, and reduction in minimum pressure threshold at resonance [99,100]. Furthermore, Sojahrood et al. [76,101] proved theoretically and experimentally that the bubble-bubble interaction is crucial even at low concentrations and influences the pressure-dependent attenuation and sound speed of the bubbly liquid significantly. This effect explains the increase in scattering power until a plateau with increasing void fraction and then decrease, as predicted in [55] and observed experimentally in [102]. Therefore, the bubble-bubble inter-action will be addressed in future studies, particularly for the void fraction curves.
Certain assumptions mentioned in Sect. 2, such as the initial uniform bubble radius distribution, thermal effect, and drag force acting on the bubbles for uncoated bubbles, have been considered in our previous studies [94,103,104]. Since the inclusion of the aforementioned effects may obscure the effects of buckling and rupture, these effects are excluded from this study.

Conclusions and prospect for future research
Using multiple-scale analysis for the weakly nonlinear propagation of one-dimensional ultrasound, we derived a KdVB equation. This equation characterizes wave propagation based on advection, nonlinear, attenuation, and dispersion effects. The analysis showed that shell with buckling and rupture phenomena increases the absolute values of all the effects through the increase in eigenfrequency caused by the elastic coefficient and variation of surface tension value. Moreover, the behavior of the nonlinear term is significantly affected by the  Fig. 10 Coefficients of the KdVB equation versus initial surface tension i (σ * 0 ) curves of α 0 = 0.05 (blue) and the effects of the shell at α 0 = 0.05 (orange). The curves denoted as "Shell effect" represent difference of i 's value for coated bubble case and that for uncoated bubble case. The dashed curves and solid curves represent the case where the lower limit was removed and original model case, respectively. The other parameters are same with Fig. 3. The high void fractions are graphed to demonstrate the behavior of the solutions rate of change of the elastic coefficient. For the dispersion terms, the shell contribution is relatively small, whereas for the attenuation, the contribution is medium, more prominent near buckled and ruptured states, and for nonlinear term the contribution of the shell is dominant.
From these results, a basic understanding of the effect of the phospholipid shell on the nonlinear propagation of ultrasound was obtained. In the future, analysis of ultrasound propagation incorporating other shell models and other assumptions, such as bubble-bubble interaction and effect of non-Newtonian fluid [52], will be conducted. Numerical simulations, evaluation of the stability of solutions, and inclusion of elastic continuum mechanics, such as plastic deformation and anisotropy [41], will also be considered. Finally, a comprehensive model including the investigated assumptions will be proposed.
Funding This work was partially carried out with the aid of the JSPS KAKENHI (22K03898), based on results obtained from a project subsidized by the New Energy and Industrial Technology Development Organization (NEDO)(JPNP20004), and Ono Charitable Trust for Acoustics and supported by JKA and its promotion funds from KEIRIN RACE.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Declarations
Competing interests 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.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Marmottant-Gompertz model
The full form of Marmottant-Gompertz model [85] is given as follows: where σ * 0 is the initial surface tension at R * = R * 0 , σ * c is the surface tension of the clean gas-liquid interface, χ * is the surface elasticity of the lipid monolayer, and R * buck can be expressed as: Although not shown in Eq. (A.1), it is helpful to show the radius at which the bubble begins to rupture, that is, rupture radius: By using equation for balance of the normal stress in Eq. (11) and modified Rayleigh-Plesset equation in Eq. (7), we have the following equation: Finally, we can express the surface tension term in Eq. (A.4) by Eq. (A.1) and we obtain the full form of the equation for spherical oscillations of bubbles.