# Nonspherical oscillations of an encapsulated microbubble with interface
energy under the acoustic field

Nehal Dash  Ganesh Tamadapu<EMAIL_ADDRESS>
###### Abstract

The practical applications of gas-filled encapsulated microbubbles involve
inherent nonspherical oscillations under acoustic fields. The gas-
encapsulation and encapsulation-liquid interfaces significantly affect the
mechanics of the bubbles, especially of smaller radii, and their consideration
is vital for mimicking the experimental setting. In this paper, we apply the
interface energy model [N. Dash and G. Tamadapu, J. Fluid Mech. 932, A$26$
$(2022)$] to examine the nonspherical oscillations of an encapsulated
microbubble with a radius of $2\,\mu$m and $5\,\mu$m under an acoustic field.
Using the Lagrangian energy formulation, the coupled dynamical governing
equations for spherical and nonspherical modes are derived, incorporating the
effects of interface energy at the interfaces, shell elasticity, and
viscosity. Through a perturbation analysis based on the Krylov-Bogoliubov
method of averaging, a set of first-order differential (slow-time) equations
is obtained to conduct steady-state and conditional-stability analysis. The
stability analysis helped in determining the excitation pressure and frequency
of the acoustic field required for smaller radii bubbles to exhibit finite
amplitude shape oscillations. Direct numerical simulations of the governing
equations revealed that the parametrically forced even mode ($n=2$) excites
even modes, while the odd modes ($n=3$) excite both even and odd modes. For
smaller radii bubbles, we observe shape mode oscillations of finite non-zero
amplitudes only in the presence of interface parameters. The initial size-
dependent interface parameter and shell viscoelastic parameters are identified
as the key parameters that play a critical role in exhibiting finite shape
mode oscillations of the bubble.

##  I Introduction

Encapsulated microbubbles (EBs) have emerged as contrast agents in clinical
ultrasound imaging. (Hoff, 2001; Stride and Saffari, 2003; Postema  _et al._ ,
2004; Lindner, 2004; Klibanov, 2006) In addition to their established
applications, research has unveiled their potential in various other fields
such as microvascular imaging and blood flow, (Errico  _et al._ , 2015)
targeted drug delivery, (Tachibana and Tachibana, 1999; Tsutsui, Xie, and
Porter, 2004; Hernot and Klibanov, 2008; Hynynen, 2008; Kooiman  _et al._ ,
2014) sonochemistry, (Suslick, 1990; Blake, 1999) for biofilm removal,
(Agarwal  _et al._ , 2012; Seo  _et al._ , 2018) cavitation cleaning, (Ohl
_et al._ , 2006; Reuter  _et al._ , 2017) among others.

In practical applications, maintaining the stability of EBs throughout the
process is vital to ensure their efficacy. Different techniques have been
developed to fabricate microbubbles, and enhance their controlled stability.
(Lee  _et al._ , 2015) However, when exposed to a sufficiently intense
acoustic field, the oscillating surface of the EB suspended in a fluid becomes
unstable. These instabilities result in the initiation of nonspherical (shape)
oscillations in the EB when the driving acoustic pressure surpasses a critical
threshold. While the behavior of EBs inherently involves nonspherical
oscillations in their applications, for the purpose of simplification in
analysis, most of the studies on EBs focused on developing bubble shell models
undergoing radial symmetric oscillations. (de Jong  _et al._ , 1992; de Jong
and Hoff, 1993; Church, 1995) Subsequently researchers reported the presence
of nonspherical or surface mode oscillations in shell-free bubble models
(Brenner, Lohse, and Dupont, 1995; Hao and Prosperetti, 1999) demonstrating
their relevant applications in the field of sonoluminescence. However, there
has always been a question of whether such observations would also be observed
in shelled bubble models, given the effects of the bubble shell on the damped
oscillations. Versluis  _et al._ (2004) and van der Meer  _et al._ (2006)
described the existence of surface modes in shelled bubbles used as contrast
agents. They also conveyed that the surface modes evolve as parametric
instabilities. Later, researchers investigated the parametric shape
instabilities in the coated microbubbles. (Dollet  _et al._ , 2008; Vos  _et
al._ , 2011) Numerical models have also been developed to study the shape
oscillations of the encapsulated microbubbles. (Liu  _et al._ , 2011;
Tsiglifis and Pelekasis, 2011) Liu  _et al._ (2011) used a boundary-fitted
finite-volume method to model the movement of boundary of the bubble, where
the shell was considered as a neo-Hookean membrane with an energy dissipation
equation to capture the shell deformation. Tsiglifis and Pelekasis (2011)
developed the numerical model to study the parametric stability and dynamic
buckling of the encapsulated microbubble. Tamadapu, Grishenkov, and Eriksson
(2016) investigated the resonance characteristics of spherical and
nonspherical modes in a thick encapsulated bubble filled with air and
suspended in water, where the shell material was considered to linear
viscoelastic and quasi-incompressible.

A noticeable aspect of understanding these nonspherical oscillations is to
have a better understanding of the physics or behavior of a single spherical
EB suspended in a viscous medium. Another important aspect is establishing an
explanation for the instabilities observed in EBs, (Dollet  _et al._ , 2008;
Vos  _et al._ , 2011) which may lead to collapse/rupture/buckling of these
EBs. (Hilgenfeldt, Lohse, and Brenner, 1996; Brenner, Hilgenfeldt, and Lohse,
2002) With these motivations, Loughran, Eckersley, and Tang (2012) proposed a
reduced analytical model for bubble surface mode oscillations that takes into
account the effect of shell properties on the oscillations. Doinikov (2004)
conducted a study on the nonlinear coupling between spherical modes,
translational motion, and shape modes of an oscillating bubble, with a
particular focus on translational instability. They derived a set of coupled
equations that described spherical oscillations of a bubble, its translational
motion and shape oscillations evolving on the bubble surface. Taking this a
step further, Shaw (2006) studied the nonlinear interactions between the
axisymmetric shape distortions, the axial translational motion, and the volume
oscillations of a gas bubble in an inviscid, incompressible liquid. They
represented the surface deformation by a complete set of Legendre polynomials
and assumed that the deformation and translational motion are small. By using
a Lagrangian energy formulation, they derived a system of equations that
remained valid up to third order approximation in these interaction terms.
Shaw (2009) also studied the stability and nonspherical oscillations of a
bubble assuming that the bubble translation and deformation is small.
Additionally, Shaw (2009) discussed the stability and nonspherical
oscillations of a bubble, employing the assumption that both bubble
translation and deformation are small. They used a combination of Rayleigh
dissipation function and perturbation analysis to account for the effects of
viscosity. In the subsequent study Shaw (2017) considered the role on
nonlinear shape mode interactions in bubble dynamics. They identified the
parametrically identified excited shape mode together with nonlinear
excitement of other shape modes, modification of the spherical mode and
induced translation of the bubble. Introducing the damping effect in
compressibility and viscosity, they showcased that a nonlinear coupling
between parametrically forced shape mode and other modes is essential for the
bubble to achieve stable oscillatory shape deformation. Recently Guédra and
Inserra (2018) followed the mathematical formulation of Shaw (2006) and
derived analytical solutions for weakly nonlinear shape oscillations of finite
amplitude. Using perturbation methods they analysed the shape oscillations of
bubbles at small amplitudes via spherical harmonics. Under the steady-state
conditions, the equations yield analytical expressions of the modal
amplitudes, conditionally stable and absolutely stable thresholds for shape
oscillations are derived and analysed.

Lately, Dash and Tamadapu (2022a) developed a mathematical model based on
interface energy within the framework of surface continuum mechanics,
(Steigmann and Ogden, 1999; Gao  _et al._ , 2014) proposing its relevance for
radial oscillations of an EB. This model highlighted the significance of
interface energy at the gas-encapsulation and the encapsulation-fluid
interfaces in understanding and analyzing radial oscillations of EBs. The
interface energy model explained the influence of each of the interface
parameters, which exerted their effects through area strain, curvature,
initial size dependence, and the coupled effects of area strain and curvature.
This model, for the first time, also resolved the spurious dependency of shell
viscoelastic parameters on the initial size of the bubble, which remained
hitherto unexplained. Very recently this was experimentally clarified by
Cattaneo and Supponen (2023) where they characterized the shell dilatational
viscosity analysing the bubble’s time domain response using ultra-high-speed
microscopic imaging and optical trapping, instead of conventional bubble
spectroscopy approach. Expanding on the interface energy model, Dash and
Tamadapu (2022b) provided a description of the radial dynamics of an EB with a
nonlinear viscoelastic shell. The interface energy model has presented
numerous possibilities for studying the behavior of EBs, especially for
smaller radii bubbles. While several models have been developed to investigate
the nonlinear shape oscillations of an EB in viscous medium, the investigation
of nonspherical oscillations of a smaller radii EB under acoustic field,
considering the influence of interface energy at the two interfaces remains
unexplored.

With this motivation, we study the spherical (volumetric) and nonspherical
(shape) oscillations of a gas-filled microbubble using the interface energy
model. The bubble is encapsulated with a thin shell membrane and suspended in
an infinite incompressible medium. The influence of interface energy at both
the gas-encapsulation and encapsulation-liquid interfaces is taken into
consideration. Furthermore, the effects of elasticity and viscosity of the
thin shell membrane are incorporated into the analysis. The resulting coupled
governing equations for spherical and nonspherical modes of bubble
oscillations, which captures the effects of interface energy and shell
elasticity, are derived. The behavior of the spherical bubble in both the
spherical and nonspherical modes is analyzed using direct numerical
simulations, considering a reasonable set of physical and interface
parameters. In order to explore the conditional stability of the EB, a
perturbation analysis using Krylov-Bogoliubov averaging method is carried out.
While the model developed by Guédra and Inserra (2018) focused on analyzing
nonspherical oscillations in uncoated bubbles of larger radius, the present
study is directed towards encapsulated bubbles with smaller radii, where the
influence of interfaces becomes particularly significant. The novelty of the
present study unfolds in three key aspects: (i) it introduces an interface
energy-based mathematical model aimed at analyzing the nonspherical
oscillations of an EB with smaller radii, (ii) the model considers the shell
material as a thin viscoelastic membrane, introducing elasticity and viscosity
parameters into the mathematical model, (iii) through numerical simulations
and stability analysis the model highlights the importance of interface
parameters in existence of shape mode oscillations in smaller radii EBs with
interface effects.

This paper is organized as follows. Section II introduces the kinematics of
deformation of the EB. Section III presents the mathematical model formulation
and the resulting coupled governing equations for the present model. The
temporal evolution of the spherical and shape mode oscillations are discussed
in section IV. Section V provides the analysis of applying the Krylov-
Bogoliubov perturbation technique to the governing equation. This is followed
by the steady-state solutions of the slow-time equations and the conditional
stability analysis in section VI. The conclusions and an outlook of further
work are provided in section VII.

##  II Kinematics of deformation

Figure 1: The deformed configuration of a thin encapsulated microbubble filled
with gas and suspended in viscous medium.

Consider a gas-filled encapsulated spherical bubble suspended in an infinite
incompressible viscous medium. In the present analysis, the encapsulation is
treated as a thin mathematical surface, and the gas-encapsulation and
encapsulation-liquid interfaces are assumed to have zero thickness, following
the nonlinear continuum framework of Steigmann and Ogden (1999). These
interfaces are also referred to as Steigmann-Ogden interfaces (SOIs). (Dash
and Tamadapu, 2022a) The equilibrium radius and the encapsulation thickness of
the bubble are denoted by $R_{0}$ and $q$, respectively. A schematic of the
bubble’s deformed geometry suspended in a viscous medium with the frame of
reference is shown in figure 1.

###  II.1 Reference configuration

When no acoustic fields are applied, the bubble is considered to be in its
reference configuration $\mathfrak{D}_{0}$. Assume that $\\{X^{\rm 1},X^{\rm
2},X{\rm{}^{3}}\\}$ be the three-dimensional Cartesian coordinate system with
unit basis vector triad as $\\{\bm{e}_{1},\bm{e}_{2},\bm{e}_{3}\\}$, and that
$\\{r,\phi,\theta\\}$ be the principal spherical polar coordinates with origin
at the center and unit basis vector triad
$\\{\bm{e}_{r},\bm{e}_{\phi},\bm{e}_{\theta}\\}$. Here, $\\{r,\phi,\theta\\}$
corresponds to the radial distance from the centre of the bubble, polar angle
formed with $X^{3}-$axis, and azimuthal angle measured about $X^{3}-$axis in
the counter clockwise direction from $X^{1}-$axis, respectively. The
coordinates of a material particle at a point in the undeformed reference
configuration is given by the following relation between the Cartesian and the
spherical polar coordinates as

$\displaystyle\begin{split}X^{1}&=r\sin\phi\cos\theta,\\\
X^{2}&=r\sin\phi\sin\theta,\\\ X^{3}&=r\cos\phi.\end{split}$ (1)

###  II.2 Deformed configuration

Under the action of acoustic field, the bubble deforms and attains a deformed
configuration $\mathfrak{D}$. The coordinates of the same typical point at a
given time instant $t$ in an arbitrary deformed configuration, with the
Cartesian coordinates $\\{x^{\rm 1},x^{\rm 2},x{\rm{}^{3}}\\}$ and the
spherical polar coordinates $\\{\rho(\phi,t),\phi,\theta\\}$ is given by

$\displaystyle\begin{split}x^{1}&=\rho(\phi,t)\sin\phi\cos\theta,\\\
x^{2}&=\rho(\phi,t)\sin\phi\sin\theta,\\\
x^{3}&=\rho(\phi,t)\cos\phi.\end{split}$ (2)

###  II.3 Interface energy model

Steigmann and Ogden (1999) proposed a nonlinear continuum framework to
describe the kinematics of the elastic surface-substrate interactions. Later,
Gao  _et al._ (2014) developed a more general interface theory considering the
effects of curvature-dependent interface energy and the residual elastic field
in the bulk induced by this interface energy. Following these works, Dash and
Tamadapu (2022b) and Dash and Tamadapu (2022a) proposed an interface energy
model within the framework of surface continuum mechanics to study the radial
dynamics of an encapsulated microbubble. In this section, we will provide some
preliminary notations and definitions to describe the kinematics of
deformation of these interfaces. For detailed explanations and derivations the
reader may refer toSteigmann and Ogden (1999); Gao  _et al._ (2014); Dash and
Tamadapu (2022b, a).

Let $\theta^{1}=\phi$ and $\theta^{2}=\theta$ be the surface coordinates of
the spherical interface with radius $r$. Consider
$\bm{Z}(\theta^{1},\theta^{2})$ and $\bm{z}(\theta^{1},\theta^{2})$ as the
position vectors of the same point on the undeformed and deformed interfaces.
Assuming that the interface is convected by the deformation of the bulk of the
bubble, the same material point before and after the deformation can be
related using the deformation mapping $\bm{\chi}$ such that

$\displaystyle\bm{z}(\theta^{1},\theta^{2})=\bm{\chi}\left(\bm{Z}(\theta^{1},\theta^{2})\right).$
(3)

The respective tangent vectors $\bm{G}_{\alpha}$ and $\bm{g}_{\alpha}$ on the
undeformed and deformed interfaces induced by these coordinates, are given by

$\displaystyle\begin{split}\bm{G}_{\alpha}=\bm{Z}_{,\alpha},\quad\bm{g}_{\alpha}=\bm{z}_{,\alpha},\quad\alpha\in\\{1,2\\},\quad(\cdot)_{,\alpha}=\frac{\partial(\cdot)}{\partial\theta^{\alpha}}.\end{split}$
(4)

The components of the covariant metric tensor in the undeformed and deformed
interfaces, respectively, are given by

$\displaystyle G_{\alpha\beta}$
$\displaystyle=\text{diag}\left[r^{2},r^{2}\sin^{2}{\phi}\right],$ (5)
$\displaystyle g_{\alpha\beta}$
$\displaystyle=\text{diag}\left[{\rho(\phi,t)}^{2}+{\rho^{\prime}(\phi,t)}^{2},{\rho(\phi,t)}^{2}\sin^{2}\phi\right],$
(6)

where the prime ${(\cdot)}^{\prime}$ denotes the derivative with respect to
$\phi$. The components of the mixed right Cauchy-Green deformation tensor can
be obtained as

$\displaystyle{C^{\alpha}}_{\beta}=G^{\alpha\delta}\,g_{\delta\beta}=\text{diag}\left[\frac{{\rho(\phi,t)}^{2}+{\rho^{\prime}(\phi,t)}^{2}}{r^{2}},\frac{{\rho(\phi,t)}^{2}}{r^{2}}\right].$
(7)

The second fundamental forms which represent the normal curvatures for the
undeformed and deformed interfaces are denoted by $Q_{\alpha\beta}$ and
$q_{\alpha\beta}$, respectively, given by

$\displaystyle Q_{\alpha\beta}=\bm{N}\cdot\bm{G}_{\alpha,\beta},\quad
q_{\alpha\beta}=\bm{n}\cdot\bm{g}_{\alpha,\beta}.$ (8)

where $\bm{N}$ and $\bm{n}$ are the oriented unit normals to the undeformed
and deformed interfaces, respectively. By definition, the relative curvature
tensor $\kappa_{\alpha\beta}=-q_{\alpha\beta}$ and its co-variant components
are given by

$\displaystyle\kappa_{\alpha\beta}=\text{diag}\left[\kappa_{11},\kappa_{22}\right],$
(9)

where

$\displaystyle\kappa_{11}$
$\displaystyle=-\frac{{\rho(\phi,r)}^{2}+2{\rho^{\prime}(\phi,t)}^{2}-\rho(\phi,t)\rho^{\prime\prime}(\phi,t)}{\sqrt{{\rho(\phi,t)}^{2}+{\rho^{\prime}(\phi,t)}^{2}}},$
$\displaystyle\kappa_{22}$
$\displaystyle=-\frac{\rho(\phi,t)\sin(\phi)\left[\rho(\phi,t)\sin{\phi}-\rho^{\prime}(\phi,t)\cos{\phi}\right]}{\sqrt{{\rho(\phi,t)}^{2}+{\rho^{\prime}(\phi,t)}^{2}}}.$

The contra-variant components of the relative curvature tensor
$(\kappa^{\alpha\beta})$ are obtained using the contraction operation

$\displaystyle\kappa^{\alpha\beta}=G^{\alpha\gamma}G^{\beta\delta}\kappa_{\gamma\delta}.$
(10)

The adjugate $(\tilde{\cdot})$ of the symmetric right Cauchy-Green deformation
tensor $\tilde{C}^{\alpha\beta}$, and the relative curvature tensor
$\tilde{\kappa}^{\alpha\beta}$ are defined by Steigmann and Ogden (1999)

$\displaystyle\tilde{C}^{\alpha\beta}$
$\displaystyle=\left(g/G\right)g^{\alpha\beta},$ (11)
$\displaystyle\tilde{\kappa}^{\alpha\beta}$
$\displaystyle=\mu^{\alpha\beta}\mu^{\beta\lambda}\kappa_{\gamma\lambda},$
(12)

where $g={\rm{det}}\left(g_{\alpha\beta}\right)$,
$G={\rm{det}}\left(G_{\alpha\beta}\right)$,
$\mu^{\alpha\beta}=e^{\alpha\beta}/\sqrt{G}$, and
$e^{\alpha\beta}=e_{\alpha\beta}$ is the alternator symbol.

The gas-encapsulation and the encapsulation-liquid interfaces are considered
as hemitropic interfaces, and the energy density can be expressed as a
function of right Cauchy-Green interface deformation tensor $\bm{C}$ and the
relative curvature tensor $\bm{\kappa}$ (Fung, 1977; Steigmann and Ogden,
1999; Gao  _et al._ , 2014; Dash and Tamadapu, 2022a) as

$\displaystyle\gamma=\gamma(\bm{C},{\bm{\kappa}}),$ (13)

and satisfies the relation

$\displaystyle\gamma(\bm{C},{\bm{\kappa}})=\gamma(\bm{Q}\bm{C}\bm{Q}^{\rm
T},\bm{Q}{\bm{\kappa}}\bm{Q}^{\rm T}),$ (14)

where $\bm{Q}$ is a proper-orthogonal second-order tensor. The Cauchy
interface stress $(\bm{\sigma})$ and the moment $(\bm{m})$ tensors are
calculated using the relations

$\displaystyle
J{\bm{\sigma}}=2\,\frac{\partial\gamma}{\partial\bm{C}},\quad{\bm{m}}=\frac{\partial\gamma}{\partial\bm{\kappa}}.$
(15)

Further, the expressions for the components $(T^{\alpha\beta})$ of the
interface stress tensor $\left(\bm{T}\right)$ and the components
$(M^{\alpha\beta})$ of the bending moment tensor $\left(\bm{M}\right)$ are
obtained in the form of the following constitutive equations (for detailed
derivation see Steigmann and Ogden (1999); Gao  _et al._ (2014); Dash and
Tamadapu (2022a))

$\displaystyle\begin{split}T^{\alpha\beta}&=\frac{1}{2}J\sigma^{\alpha\beta}=\frac{\partial\gamma}{\partial
I_{1}}G^{\alpha\beta}+\frac{\partial\gamma}{\partial
I_{2}}\tilde{C}^{\alpha\beta}+\frac{\partial\gamma}{\partial
I_{5}}\kappa^{\alpha\beta}+\frac{1}{2}\frac{\partial\gamma}{\partial
I_{6}}\left(D^{\alpha\beta}+D^{\beta\alpha}\right),\\\
M^{\alpha\beta}&=Jm^{\alpha\beta}=\frac{\partial\gamma}{\partial
I_{3}}G^{\alpha\beta}+\frac{\partial\gamma}{\partial
I_{4}}\tilde{\kappa}^{\alpha\beta}+\frac{\partial\gamma}{\partial
I_{5}}C^{\alpha\beta}+\frac{1}{2}\frac{\partial\gamma}{\partial
I_{6}}\left(E^{\alpha\beta}+E^{\beta\alpha}\right),\end{split}$ (16)

where

$\displaystyle J=\sqrt{g/G},\quad
D^{\alpha\beta}=G_{\gamma\delta}\mu^{\alpha\gamma}\kappa^{\beta\delta},\quad
E^{\alpha\beta}=G_{\gamma\delta}\mu^{\alpha\gamma}C^{\beta\delta},$ (17)

and $(I_{1},I_{2},I_{3},I_{4},I_{5},I_{6})$ are the six basis invariants
(defined in appendix A) of the right Cauchy-Green interface deformation tensor
$\bm{C}$, the relative curvature tensor $\bm{\kappa}$, and the permutation
tensor-density $\bm{\mu}$ on the undeformed interface.

##  III Mathematical model

The following section presents the mathematical model formulation for the
present problem. Section III.1 outlines the energy density functions for the
shell material. The first variational formulations and the resulting governing
equations are discussed in subsections III.2 and III.3, respectively.

###  III.1 Material energy density function for the thin shell membrane

The shell material and surrounding fluid are considered to be homogeneous,
isotropic and incompressible. The bubble shell material assumed to be a
hyperelastic membrane following the constitutive relation of the neo-Hookean
material model. The neo-Hookean material is chosen because it offers a simple
and realistic model for a rubbber-elastic type material. For the
incompressible bulk neo-Hookean material model, the strain energy density
function (per unit undeformed volume) is given by

$\displaystyle\Psi_{\rm s}=C_{1}(\tilde{I}_{1}-3),$ (18)

where $C_{1}$ is the material elastic constant related to the shear modulus
$\mu=2\,C_{1}$, and $\tilde{I}_{1}$ is the first invariant of the right-Cauchy
Green deformation tensor $\bm{C}$. Thus, $\tilde{I}_{1}$ is given by

$\displaystyle\tilde{I}_{1}=\rm{tr}\,\bm{C}=\Lambda_{1}^{2}+\Lambda_{2}^{2}+\Lambda_{3}^{2},$
(19)

where $(\Lambda_{1},\Lambda_{2},\Lambda_{3})$ are the principal stretch ratios
in the radial (across the thickness), meridional, and azimuthal directions,
respectively. The expressions for the principal stretch ratios are given by

$\displaystyle\Lambda_{1}=\frac{r^{2}}{\rho(\phi,t)\sqrt{{\rho(\phi,t)}^{2}+{\rho^{\prime}(\phi,t)}^{2}}},\quad\Lambda_{2}=\frac{\sqrt{{\rho(\phi,t)}^{2}+{\rho^{\prime}(\phi,t)}^{2}}}{r},\quad\Lambda_{3}=\frac{{\rho(\phi,t)}}{r}.$
(20)

Here the stretch $(\Lambda_{1})$ across the radial (thickness) direction has
been calculated using the incompressibility constraint
$\Lambda_{1}\Lambda_{2}\Lambda_{3}=1$.

###  III.2 Variational formulation

In addition to those in Guédra and Inserra (2018), the new terms in the
governing equations related to the interface energy and viscoelastic shell are
obtained using a variational formulation of the energy functional. The energy
functional $(\Pi)$ for the present problem constitutes of three parts: the
interface energy $(\Pi_{\rm{ie}})$, and energy of the viscoelastic shell
$(\Pi_{\rm{s}})$, such that

$\displaystyle\Pi=\Pi_{\rm{ie}}+\Pi_{\rm{s}},$ (21)

and the expressions for the respective energy functional are given by

$\displaystyle\Pi_{\rm{ie}}=\int_{0}^{\pi}\gamma\left(\bm{C},\bm{\kappa}\right)\sqrt{G}\,\text{d}\phi,\quad\Pi_{\rm{s}}=\int_{0}^{\pi}\Psi_{\rm{s}}\;q\sqrt{G}\,\text{d}\phi,$
(22)

where $\sqrt{g}=r^{2}\sin\phi$. Using the variational principle, the total
energy functional is minimized $(\delta\Pi=0)$ to obtain respective terms in
the governing equilibrium equations. Since the interface energy density is a
function of $\bm{C}$ and $\bm{\kappa}$, its variation can further be written
as (Gao  _et al._ , 2014)

$\displaystyle\delta\left(\gamma\left(\bm{C},\bm{\kappa}\right)\right)=\frac{\partial\gamma}{\partial\bm{C}}:\delta\bm{C}+\frac{\partial\gamma}{\partial\bm{\kappa}}:\delta\bm{\kappa}=\bm{T}:\delta\bm{C}+\bm{M}:\delta\bm{\kappa}.$
(23)

Therefore, the first variation of the energy functional is given by

$\displaystyle\delta\Pi=\int_{0}^{\pi}\left(\bm{T}:\delta\bm{C}+\bm{M}:\delta\bm{\kappa}\right)\sqrt{G}\,\text{d}\phi+\int_{0}^{\pi}\delta\left(\Psi_{\rm{s}}\right)q\sqrt{G}\,\text{d}\phi.$
(24)

The expression for $\delta\bm{C}$ and $\delta\bm{\kappa}$ can be further
expressed in terms of the $\delta R$ and $\delta a_{n}$ to obtain additional
terms in the governing differential equation.

###  III.3 Equations governing the spherical mode and the shape mode

In this study, the wavelength of the acoustic field is considered to be
sufficiently large such that it can be assumed to behave uniformly on the
bubble surface. For simplicity, we neglect the translational motion of the
bubble. Assuming that the motion of the spherical bubble is dominated in the
radial direction, the perturbation to the spherical bubble can be expanded in
terms of spherical harmonics. Since we restrict attention to small
axisymmetric shape deformation of the bubble, the spherical harmonics reduces
to Legendre polynomials. At time $t$, the surface of the bubble
$\rho\left(\phi,t\right)$ is given by

$\displaystyle\rho\left(\phi,t\right)=R\left(t\right)+\sum_{n=2}^{\infty}\epsilon
a_{n}(t)P_{n}\left(\cos{\phi}\right),$ (25)

where $R(t)$ is the radius of the spherical bubble (spherical or volume mode),
$P_{n}(\cos\phi)$ represents the Legendre polynomial of order $n$, and
$\epsilon\,a_{n}(t)$ is the amplitude of the $n$th Legendre mode, commonly
referred to as shape mode amplitude. Here, the surface/shape distortion terms
of the bubble radius $R(t)$ are assumed to be small, therefore their
amplitudes $a_{n}(t)$ are scaled by the small parameter $\epsilon$, whereas no
restrictions were imposed on the spherical/volume oscillations.

Shaw (2006) developed the Lagrangian formulation for a gas bubble in an
incompressible liquid of infinite extent. By integrating the Lagrangian
density across the problem domain, the resulting Lagrangian $\mathcal{L}$ is
obtained and can be written as

$\displaystyle\mathcal{L}=\mathcal{C}-\mathcal{T}-\mathcal{K}.$ (26)

where $\left(\mathcal{C},\mathcal{T},\mathcal{K}\right)$ stands for the
constraint, potential energy and kinetic energy terms, respectively. Their
respective expressions are given by Shaw (2006)

$\displaystyle\mathcal{C}$ $\displaystyle=-2\pi\rho^{\rm
L}\int_{0}^{2\pi}\Phi\left(\rho\frac{\partial\rho}{\partial
t}\right)\rho\sin\phi\,\rm{d}\phi,$ (27) $\displaystyle\mathcal{T}$
$\displaystyle=\sigma_{(\cdot)}\,\rvert
S\rvert+\Pi-V\left[\frac{p_{g_{0}}}{1-k}{\left(\frac{V_{0}}{V}\right)}^{k}-p_{\infty}\right],$
(28) $\displaystyle\mathcal{K}$ $\displaystyle=-\pi\rho^{\rm
L}\int_{-1}^{1}\Phi f_{1}\,\rm{d}(\cos\phi),$ (29)

where

$\displaystyle\rvert S\rvert$
$\displaystyle=2\pi\int_{0}^{\pi}\rho^{2}\left[1+\frac{1}{\rho^{2}}{\left(\frac{\partial\rho}{\partial\phi}\right)}^{2}\right]^{1/2}\,\sin\phi\,\rm{d}\phi,$
(30) $\displaystyle f_{1}$
$\displaystyle=b_{0}(t)-\sum_{n=1}^{\infty}\frac{\epsilon
b_{n}(t)(n+1)}{\rho^{n}}P_{n}(\cos\phi)-(1-\cos^{2}\phi)$
$\displaystyle\qquad\times\left(\sum_{n=2}^{\infty}\epsilon
a_{n}(t)\frac{\text{d}P_{n}}{\text{d}(\cos\phi)}\right)\left(\sum_{n=2}^{\infty}\frac{\epsilon
b_{n}(t)}{\rho^{n+1}}\frac{\text{d}P_{n}}{\text{d}(\cos\phi)}\right),$ (31)

$\rho^{\rm L}$ is the density of liquid, $V$ is the volume at any instant of
time $t$, $V_{0}$ is the initial volume, $\Phi$ is obtained by evaluating
$\varphi$, defined in (32), on the bubble surface, and $\sigma_{(\cdot)}$ is
the respective effective surface tension parameter for spherical and
nonspherical modes.

Consistent with much of other works, they defined the bubble surface as in
(25) and the general fluid velocity potential $\varphi$ as

$\displaystyle\varphi=-\frac{b_{0}(t)}{r}+\sum_{n=1}^{\infty}\epsilon\frac{b_{n}(t)}{r^{n+1}}P_{n}(\cos\phi).$
(32)

The expansions for $\rho\left(\phi,t\right)$ and potential $\varphi$ evaluated
on the bubble surface are substituted into the respective expressions of
constraint, potential, and kinetic energy. The resulting expressions are
truncated to $\mathcal{O}(\epsilon^{3})$ to obtain the expression for
Lagrangian $\mathcal{L}$ in terms of parameters
$R(t),\dot{R}(t),a_{n}(t),\dot{a}_{n}(t),$ and $b_{n}(t)$ to order
$\epsilon^{3}$. The coefficients of $b_{n}(t)$ are then eliminated from the
Lagrangian $\mathcal{L}$ by applying the kinematic boundary conditions
imposing the set of following conditions

$\displaystyle\frac{\partial\mathcal{L}}{\partial b_{n}}=0,\qquad n=0,1,2,...$
(33)

From the above conditions, we get

$\displaystyle b_{0}(t)=R^{2}\dot{R},\qquad
b_{n}(t)=-\frac{R^{n+1}}{n+1}\left(2a_{n}\dot{R}+\dot{a}_{n}R\right)+\mathcal{O}(\epsilon).$
(34)

The resultant Lagrangian is then used in the Euler-Lagrange equations to
obtain the equations governing the spherical, and shape oscillations. A
detailed description of the Lagrangian formulation can be found in Shaw
(2006). By following the mathematical formulation presented in Shaw (2006),
and using the orthogonality of Legendre polynomials, the integrals are deduced
further to obtain additional terms in the governing equations. Incorporating
the effects of interface energy (through the interface tension), and viscosity
and elasticity of the shell membrane, the present model introduces additional
terms in (35) and (36) beyond those considered by Guédra and Inserra (2018).
Thus, the resulting equations governing the spherical mode $R(t)$ and the
nonspherical (shape) modes $a_{n}(t)$ accurate to second order approximation
in $\epsilon$ with $n\geqslant 2$ are given, respectively, by

$\displaystyle\ddot{R}+\frac{3}{2}\frac{\dot{R}^{2}}{R}-\frac{1}{\rho^{\rm
L}R}\left[p_{g_{0}}{\left(\frac{R_{0}}{R}\right)}^{3k}-p_{\infty}\right]+\frac{2\sigma_{\rm
sp}}{\rho^{\rm L}R^{2}}+4\nu^{\rm
S}\frac{\dot{R}}{R^{2}}+\frac{4C_{1}q}{\rho^{\rm
L}R^{2}}\left[1-\frac{R_{0}^{6}}{R^{6}}\right]$
$\displaystyle+\frac{\epsilon^{2}}{R}\sum_{n=2}^{\infty}\Bigg{\\{}\frac{1}{\mathcal{H}_{n}}\left[\left(n+\frac{3}{2}\right)\dot{a}_{n}^{2}+\left(n+3\right)a_{n}\ddot{a}_{n}-\left(n-3\right)\left(\frac{\ddot{R}}{R}a_{n}^{2}+\frac{1}{2}\frac{\dot{R}^{2}}{R^{2}}\,a_{n}^{2}+2\frac{\dot{R}}{R}a_{n}\dot{a}_{n}\right)\right]\Bigg{\\}}$
$\displaystyle+\frac{\epsilon^{2}}{\rho
R^{3}}\left[p_{\infty}-p_{g_{0}}\left(1-3k\right)\left(\frac{R_{0}}{R}\right)^{3k}\right]\sum_{n=2}^{\infty}\,\frac{a_{n}^{2}}{\left(2n+1\right)}+\epsilon^{2}\frac{6C_{1}qR_{0}^{6}}{\rho^{\rm
L}R^{10}}\left[\sum_{n=2}^{\infty}\,\frac{(n(n+1)-10)}{2n+1}a_{n}^{2}\right]=0,$
(35)

and

$\displaystyle\epsilon\ddot{a}_{n}+\epsilon\left(\frac{3\dot{R}}{R}+F_{\nu}\right)\dot{a}_{n}+\epsilon(n+1)\frac{\sigma_{\rm
nsp}^{0}}{\rho^{\rm
L}R^{3}}a_{n}+\epsilon\left(n-1\right)\left(G_{\nu}-\frac{\ddot{R}}{R}\right)a_{n}$
$\displaystyle+\epsilon\,(n+1)a_{n}\frac{2C_{1}q}{\rho^{\rm
L}R^{3}}\left[\big{\\{}(n-1)(n+2)\big{\\}}+\big{\\{}14-n(n+1)\big{\\}}\frac{R_{0}^{6}}{R^{6}}\right]+\epsilon^{2}\frac{\mathcal{H}_{n}}{2\rho^{\rm
L}R^{4}}$
$\displaystyle\times\left[\sum_{i=2}^{\infty}\sum_{j=2}^{\infty}\sigma_{\rm
nsp}^{1}a_{i}a_{j}\right]+\epsilon^{2}\frac{\mathcal{H}_{n}}{4R}\sum_{i=2}^{\infty}\,\sum_{j=2}^{\infty}\Bigg{\\{}\frac{\ddot{R}}{R}a_{i}a_{j}G_{d_{ijn}}+\frac{\dot{R}^{2}}{R^{2}}a_{i}a_{j}M_{a_{nij}}+\frac{\dot{R}}{R}a_{i}\dot{a}_{j}M_{b_{nij}}$
$\displaystyle+a_{i}\ddot{a}_{j}M_{c_{nij}}+\dot{a}_{i}\dot{a}_{j}M_{d_{nij}}\Bigg{\\}}+\epsilon^{2}\frac{\mathcal{H}_{n}}{2\rho^{\rm
L}R^{3}}\left[p_{\infty}-p_{g_{0}}\left(\frac{R_{0}}{R}\right)^{3k}\right]\sum_{i=2}^{\infty}\sum_{j=2}^{\infty}a_{i}a_{j}I_{a_{nij}}$
$\displaystyle-\epsilon^{2}\frac{6C_{1}q\mathcal{H}_{n}R_{0}^{6}}{\rho^{\rm
L}R^{10}}\sum_{i=2}^{\infty}\sum_{j=2}^{\infty}\left(5I_{a_{nij}}+\frac{1}{2}I_{c_{nij}}+I_{c_{ijn}}\right)a_{i}a_{j}=0.$
(36)

Here, $\mathcal{H}_{n}=(2n+1)(n+1)$, $R_{0}$ is the radius of bubble at static
equilibrium configuration, $\nu^{\rm L}$ and $\nu^{\rm S}$ denote the
kinematic viscosity of surrounding liquid and bubble shell material,
respectively, and $k$ is the polytropic expansion index. The gas pressure
inside the bubble $\left(p_{g_{0}}\right)$ is represented in terms of static
liquid pressure $(p_{0})$ and the surface tension parameter at static
equilibrium $(\sigma_{\rm eq})$ given by

$\displaystyle p_{g_{0}}=p_{0}+\frac{2\sigma_{\rm eq}}{R_{0}},$ (37)

with

$\displaystyle\sigma_{\rm eq}$
$\displaystyle=\gamma_{0}+2\gamma_{1}+2\gamma_{2}+\gamma_{3}\frac{1}{R_{0}}+\gamma_{4}\frac{1}{R_{0}^{2}}+3\gamma_{5}\frac{1}{R_{0}},$
(38)

where $\gamma_{k}=\gamma_{1k}+\gamma_{2k}$ for $k=\\{0,1,2,4\\}$ and
$\gamma_{l}=-\gamma_{1l}+\gamma_{2l}$ for $l=\\{3,5\\}$ are the interface
constants with $\gamma_{ij},i=1,2,\,j=1\,\text{to}\,5$ as interface parameters
at the gas–encapsulation and encapsulation–liquid interfaces, respectively. In
Dash and Tamadapu (2022a) and Dash and Tamadapu (2022b), the interface
parameters for the gas-encapsulation and encapsulation-liquid interfaces are
associated with the inner and outer radii of the bubble, respectively. In the
present work, mathematical surface treatment of the bubble oscillations in
terms of single radial parameter $R(t)$ leads to the net interface parameters
$\gamma_{i}$ in governing equations. It is important to note that the
orientation of the normal vector plays an important role for the interface
parameters connected to the curvature tensor. (Dash and Tamadapu, 2022a, b)

The liquid pressure $p_{\infty}\left(t\right)$ at any instant of time $t$ is

$\displaystyle p_{\infty}\left(t\right)=p_{0}-p_{a}\cos\left({2\pi
ft}\right),$ (39)

where $p_{a}$ is acoustic pressure excitation. The terms $(\sigma_{\rm
sp},\sigma_{\rm nsp}^{0},\sigma_{\rm nsp}^{1})$ represent the effective
surface tension in terms of interface parameters $\left(\gamma_{ij}\right)$
for the spherical $(\cdot)_{\rm sp}$ and nonspherical $(\cdot)_{\rm nsp}$
modes, respectively, expressed as follows

$\displaystyle\sigma_{\rm sp}$
$\displaystyle=\gamma_{0}+2\left[\gamma_{1}+\gamma_{2}\frac{R^{2}}{R_{0}^{2}}+\gamma_{5}\frac{R}{R_{0}^{2}}\right]+\frac{1}{R}\bigg{[}\gamma_{3}+\gamma_{4}\frac{R}{R_{0}^{2}}+\gamma_{5}\frac{R^{2}}{R_{0}^{2}}\bigg{]}$
$\displaystyle+\frac{\epsilon^{2}}{R^{2}}\sum_{n=2}^{\infty}\frac{a_{n}^{2}}{2n+1}\bigg{[}\big{\\{}n(n+1)+6\big{\\}}\gamma_{2}\frac{R^{2}}{R_{0}^{2}}-\frac{n(n+1)}{2R}\gamma_{3}+\big{\\{}2n(n+1)+3\big{\\}}\gamma_{5}\frac{R}{R_{0}^{2}}\bigg{]},$
(40a) $\displaystyle\sigma^{0}_{\rm nsp}$
$\displaystyle=(n-1)(n+2)\Bigg{[}\gamma_{0}+2\gamma_{1}+2\gamma_{2}\frac{R^{2}}{R_{0}^{2}}+2\gamma_{3}\frac{1}{R}+5\gamma_{4}\frac{1}{R_{0}^{2}}+8\gamma_{5}\frac{R}{R_{0}^{2}}\Bigg{]}$
$\displaystyle+4\gamma_{2}\frac{R^{2}}{R_{0}^{2}}+8\gamma_{4}\frac{1}{R_{0}^{2}}+16\gamma_{5}\frac{R}{R_{0}^{2}},$
(40b) $\displaystyle\sigma_{\rm nsp}^{1}$
$\displaystyle=N^{2}_{nij}\gamma_{2}\frac{R^{2}}{R_{0}^{2}}+N^{3}_{nij}\gamma_{3}\frac{1}{R}+N^{4}_{nij}\gamma_{4}\frac{1}{R_{0}^{2}}+N^{5}_{nij}\gamma_{5}\frac{R}{R_{0}^{2}}.$
(40c)

In the present study, the viscous dissipation in the surrounding fluid and the
shell membrane has been introduced in a ad hoc manner into equations (35) and
(36). Following the discussions presented in Guédra and Inserra (2018), these
terms are introduced as classical viscous terms to lower orders. In the
equation governing the spherical mode of oscillation, the influence of
viscosity arising from both the bubble shell and the surrounding fluid is
introduced using the classical term, similar to the equation governing the
radial oscillations of the bubble. The viscous effects associated with the
surrounding fluid in the shape mode equation are incorporated using
established functions within the boundary layer approximation such that

$\displaystyle F_{\nu}$ $\displaystyle=(n+2)\left[(2n+1)\nu^{\rm
S}-2n(n+2)\frac{\delta}{R}\nu^{\rm L}\right]\frac{2}{R^{2}},$ (41a)
$\displaystyle G_{\nu}$ $\displaystyle=(n+2)\left[\nu^{\rm
S}+2n\frac{\delta}{R}\nu^{\rm L}\right]\frac{2}{R^{2}}\frac{\dot{R}}{R},$
(41b)

where $\delta=\sqrt{\nu^{\rm L}/\omega}$ represents the viscous boundary layer
thickness. (Brenner, Lohse, and Dupont, 1995) A comprehensive mathematical
analysis has been conducted by Shaw (2009), taking into account the combined
impact of the Rayleigh dissipation function and perturbation analysis, in
order to address the effects of viscosity. This would result in nonlinear
viscous terms at higher orders in (35) and (36). However, Guédra and Inserra
(2018) considered fluid viscosity as a small perturbation, resulting in the
neglect of such higher-order terms in the subsequent asymptotic expansion of
both the spherical and nonspherical equations anyway. Following a similar
approach, we consider that both the fluid and shell viscosity act as small
perturbations. This simplifies the analysis and drops the nonlinear viscous
terms appearing at the higher-orders expansions. The set of integrals (see
appendix B), and the other coefficients appearing the subsequent calculations
are listed in appendix C.

Table 1:  The numerical values of inner and outer interface parameters such as $(\gamma_{11},\gamma_{21},\gamma_{12},\gamma_{22})\,$N/m, $(\gamma_{14},\gamma_{24})\,$Nm and $(\gamma_{13},\gamma_{23},\gamma_{15},\gamma_{25})\,$N of the bubble chosen for numerical simulations. $\gamma_{11}$ | $\gamma_{21}$ | $\gamma_{12}$ | $\gamma_{22}$ | $\gamma_{13}$ | $\gamma_{23}$ | $\gamma_{14}$ | $\gamma_{24}$ | $\gamma_{15}$ | $\gamma_{25}$
---|---|---|---|---|---|---|---|---|---
0.01 | 0.01 | 0.01 | 0.01 | 0.04 | 0.035 | 0.03 | 0.01 | 0.04 | 0.03

The primary focus of this work is to understand the nonspherical oscillations
of smaller radii microbubbles, where the influence of interface parameters
becomes particularly important. Therefore, we consider the equilibrium radius
of the bubble $R_{0}=2\,\mu$m with shell membrane thickness of $q=20\,$nm. The
density of the surrounding liquid is considered to be $\rho^{\rm
L}=1000\,$kg/${\rm{m^{3}}}$, viscosity of the fluid $(\nu^{\rm L})$ and shell
$(\nu^{\rm S})$ are considered to be $(10^{-6})\,\rm{m^{2}/s}$, the shell
material elastic constant $C_{1}=0.1\,$MPa, and static liquid pressure
$p_{0}=0.1\,$MPa. The values of these physical and material parameters remain
consistent throughout the numerical simulations, unless otherwise specified.
Given our interest in studying an EB with a radius of
$\mathcal{O}\left(10^{-6}\right)$m and a thickness of
$\mathcal{O}\left(10^{-9}\right)$m, we assume the order of interface
parameters
$\gamma_{10},\gamma_{20},\gamma_{11},\gamma_{21},\gamma_{12},\gamma_{22}$ as
$\mathcal{O}(1)\,$N/m, $\gamma_{13},\gamma_{23},\gamma_{15},\gamma_{25}$ as
$\mathcal{O}(10^{-6})\,$N and $\gamma_{14},\gamma_{24}$ as
$\mathcal{O}(10^{-12})\,$N m. The order of these interface parameters are such
that the effective interface tension parameters possesses reasonable values of
$\mathcal{O}(1)$.(Dash and Tamadapu, 2022b, a) The values of these interface
parameters for the present study are tabulated in table 1. The coupled
governing equations (35) and (36) indicate the complex interactions between
spherical and shape deformation of an EB. It is also apparent that numerical
studies using such mathematical model may offer many possibilities that can be
explored. Nevertheless, the present model introduces the interface parameters
as a preliminary, yet comprehensive, study that highlights their substantial
influence on the oscillations of smaller radii bubbles.

##  IV Temporal evolution of spherical and shape mode amplitudes

In this section, the variation of the spherical and shape mode amplitudes for
an EB are discussed in detail. The coupled governing equations are solved for
the interface parameters tabulated in table 1. The two quantities $x$ and
$\epsilon s_{n}$, respectively, given by

$\displaystyle x=\frac{R(t)}{R_{0}}-1,\quad\epsilon
s_{n}=\epsilon\frac{a_{n}(t)}{R_{0}},$ (42)

are plotted against time $(t)$. Eight modes are retained in the computations
and the initial conditions for all the shape modes are set to
$s_{n}(0)=10^{-2}$ and $\dot{s}_{n}(0)=0$. Prior to these calculations the
steady state simulation of the Rayleigh-Plesset radial equation has been run
to set the initial conditions for $R(0)$ and $\dot{R}(0)$. From the direct
numerical simulations of the coupled governing equations it is observed that
the parametrically forced even mode $(n=2)$ only excites the even modes whilst
the odd modes $(n=3)$ can excite both even and odd modes, as shown in figures
2 and 3, respectively. This aligns with that of the observations highlighted
by Shaw (2006) and Guédra and Inserra (2018). However, it is important to
emphasize that the interface parameters play a crucial role in this.

Figure 2: Steady state oscillations of the spherical and the seven first shape
modes for the EB with $R_{0}=2\,\mu$m driven above the first parametric
resonance of the $n=2$ mode with $p_{a}=0.37\,$MPa, $f=4.96\,$MHz, and
$\Omega_{0,n}\sim 1$ (corresponding to the point $\bullet$ in figure 5) with
the set of interface parameters in table 1.

For the case of $n=2$, specific observations can be made regarding the
consequences of setting certain interface parameters to zero. This enables us
to interpret the effect of the specific interface parameter in the numerical
simulations while keeping all other working parameters unchanged. When the
interface parameter $\gamma_{3}$, which captures the curvature effects, is set
to zero, the even shape modes still get excited. However, the amplitude of
shape mode oscillations is slightly higher in this case. When the interface
parameter $\gamma_{4}$, which accounts for the initial size-dependent effects,
is set to zero, it leads to all the shape modes being zero. The dominant
nature of the initial size of the bubble in finite shape mode oscillations can
be attributed to the fact that the radius of the EB directly influences its
natural frequency of oscillation. Smaller EBs have higher natural frequencies,
making them more responsive to external perturbations, resulting in more
noticeable shape mode oscillations. When interface parameter $\gamma_{5}$ is
set to zero, the even shape mode exhibit a finite amplitude, though with
slightly smaller amplitudes. It is also observed that the interface parameters
$(\gamma_{1},\gamma_{2},\gamma_{4})$ play a more significant role compared to
$(\gamma_{3},\gamma_{5})$. This aspect can also be understood by looking at
their net contributions in the expressions of effective interface tension
parameters in (40). The net effect of $(\gamma_{1},\gamma_{2},\gamma_{4})$ at
both the inner and outer interfaces gets added up since
$\gamma_{i}=(\gamma_{1i}+\gamma_{2i}),\,i=\\{1,\,2,\,4\\}$, resulting in the
increase of their net contributions. Whereas in the case of
$(\gamma_{3},\gamma_{5})$, the outer interface possesses a negative sign due
to the orientation of the normal vectors, such that
$\gamma_{j}=(\gamma_{1j}-\gamma_{2j}),\,j=\\{3,\,5\\}$, resulting in the
decrease of their net contributions.

Figure 3: Steady state oscillations of the spherical and the seven first shape
modes for the EB with $R_{0}=2\,\mu$m driven above the first parametric
resonance of the $n=3$ mode with $p_{a}=0.92\,$MPa, $f=8.52\,$MHz, and
$\Omega_{0,n}\sim 1$ (corresponding to the point $\bullet$ in figure 5) with
the set of interface parameters in table 1.

Moreover, we extend this analysis for a relatively larger EB radii to enable
us to discuss the influence and significance of interface parameters. The
steady state oscillation of spherical and first seven shape modes for an EB
with radius $R_{0}=5\,\mu$m is shown in figure 4. Also as the bubble’s radius
increases, the importance of interface parameters diminishes and the analysis
is predominantly governed by the initial size of the bubble and other physical
parameters (like viscosity and inside gas pressure), making it easier to
simulate the equations without encountering numerical difficulties. However,
for smaller radius bubbles, where a small length scale is involved, the
interface parameters exhibit a dominant behavior in conjunction with other
physical parameters. Therefore, it is essential to consider their effects,
particularly while analysing small radii bubbles. Among other physical
parameters, the viscosity and elastic material constant of the shell have
consistently played a crucial role in the behavior of EBs, influencing both
their radial and nonspherical oscillations. In the present model, their
importance remains prominent.

Figure 4: Steady state oscillations of the spherical and the seven first shape
modes for the EB with $R_{0}=5\,\mu$m driven above the first parametric
resonance of the $n=2$ mode with $p_{a}=0.07\,$MPa, $f=1.17\,$MHz, and
$\Omega_{0,n}\sim 1$ (corresponding to the point $\bullet$ in figure 7) with
the set of interface parameters in table 1.

##  V Perturbation analysis using Krylov-Bogoliubov averaging method

In this section, the Krylov-Bogoliubov asymptotic perturbation technique is
used to analyze the system further. By canceling out the secular terms in the
coupled governing equations, a set of first-order differential equations is
derived. These equations capture the essential dynamics of the system and
provide a foundation for further analysis. In the subsequent section, these
equations are utilized to perform a steady-state analysis and investigate the
conditional stability of the present model.

Let us represent the spherical mode $R(t)$ and the shape mode $a_{n}(t)$ in
terms of non-dimensional coordinates $x(t)$ and $s_{n}(t)$, respectively, such
that

$R(t)=R_{0}[1+x(t)],\quad a_{n}(t)=R_{0}s_{n}(t).$ (43)

Additionally, the new time scale $\tau$ is defined as

$\tau=\frac{\omega t}{2}.$ (44)

In this analysis, the non-dimensional coordinate $x$, pressure amplitude
$p_{a}$, viscous damping measuring parameter of surrounding liquid
$(\delta/R_{0})^{2}$, and shell $\nu^{\rm S}$ are all considered as small
quantities of $\epsilon$ order. By performing a power series expansion of $x$
and considering terms up to the second orders in $\epsilon$, the equations
(35) and (36) can expressed as

$\displaystyle\ddot{x}$
$\displaystyle+\Omega_{0}^{2}x=A\cos{2\tau}+\epsilon\mathcal{F}_{1}$ (45)
$\displaystyle\ddot{s}_{n}$
$\displaystyle+\Omega_{0,n}^{2}s_{n}=\epsilon\mathcal{F}_{2}.$ (46)

where $\Omega_{0}^{2}=\omega_{0}^{2}+4B$,
$\Omega_{0,n}^{2}=\omega_{0,n}^{2}+4B(n+1)$, $\mathcal{F}_{1}$ and
$\mathcal{F}_{2}$ are, respectively, given by

$\displaystyle\mathcal{F}_{1}$
$\displaystyle=\bigg{[}\left(Q+22B\right)x^{2}-\frac{3}{2}\dot{x}^{2}-xA\cos{2\tau}-F_{\nu,0}^{(0)}\dot{x}\bigg{]}-\frac{1}{\mathcal{H}_{n}}\sum_{n=2}^{\infty}{\left(n+\frac{3}{2}\right)\dot{s}_{n}^{2}}$
$\displaystyle+\frac{1}{\mathcal{H}_{n}}\sum_{n=2}^{\infty}{\left[(n+3)\Omega_{0,n}^{2}-(n+1)U-(n+1)(n(n+1)-10)B\right]s_{n}^{2}},$
$\displaystyle\mathcal{F}_{2}$
$\displaystyle=\Bigg{[}(n-1)A\cos{2\tau}+\bigg{\\{}3\Omega_{0,n}^{2}-(n-1)\Omega_{0}^{2}+2(n+1)\big{\\{}14-n(n+1)\big{\\}}B$
$\displaystyle-(n+1)I_{\rm
p}\bigg{\\}}xs_{n}-3\dot{x}\dot{s}_{n}-F_{\nu,0}^{(n)}\dot{s}_{n}\Bigg{]}+\mathcal{H}_{n}\sum_{i=2}^{\infty}\sum_{j=2}^{\infty}\Bigg{[}\Bigg{\\{}\Omega_{0,j}^{2}\frac{M_{c_{nij}}}{4}+W_{0}I_{a_{nij}}$
$\displaystyle+B\left(5I_{a_{nij}}+I_{c_{ijn}}+\frac{I_{c_{nij}}}{2}\right)+N_{nij}\Bigg{\\}}s_{i}s_{j}-\frac{M_{d_{nij}}}{4}\dot{s}_{i}\dot{s}_{j}\Bigg{]}.$

Here $\dot{(\cdot)}$ represents the time derivatives with respect to $\tau$.
The other parameters in (45) and (46) are given by

$\displaystyle A$
$\displaystyle=\left(\frac{2}{\omega}\right)^{2}\frac{p_{a}}{\rho^{\rm
L}R_{0}^{2}},\quad B=\left(\frac{2}{\omega}\right)^{2}\frac{6C_{1}q}{\rho^{\rm
L}R_{0}^{3}},\quad W_{0}=\left(\frac{2}{\omega}\right)^{2}\frac{\sigma_{\rm
eq}}{\rho^{\rm L}R_{0}^{3}},$ (47a) $\displaystyle W$
$\displaystyle=\left(\frac{2}{\omega}\right)^{2}\frac{1}{\rho^{\rm
L}R_{0}^{3}}\Bigg{[}(n-1)(n+2)\bigg{\\{}\gamma_{0}+2\gamma_{1}+2\gamma_{2}+5\gamma_{4}\frac{1}{R_{0}^{2}}+2\gamma_{3}\frac{1}{R_{0}}+8\gamma_{5}\frac{1}{R_{0}}\bigg{\\}}$
$\displaystyle\qquad\qquad+4\gamma_{2}+8\gamma_{4}\frac{1}{R_{0}^{2}}+16\gamma_{5}\frac{1}{R_{0}}\Bigg{]},$
(47b) $\displaystyle Q$
$\displaystyle=\left(\frac{2}{\omega}\right)^{2}\frac{1}{\rho^{\rm
L}R_{0}^{2}}\Bigg{[}\frac{9}{2}k(1+k)p_{g_{0}}-\frac{2}{R_{0}}\bigg{\\{}2\gamma_{0}+4\gamma_{1}-2\gamma_{2}+5\gamma_{3}\frac{1}{R_{0}}+2\gamma_{4}\frac{1}{R_{0}^{2}}\bigg{\\}}\Bigg{]},$
(47c) $\displaystyle U$
$\displaystyle=\left(\frac{2}{\omega}\right)^{2}\frac{1}{\rho^{\rm
L}R_{0}^{2}}\Bigg{[}3kp_{g_{0}}-\frac{2}{R_{0}}\bigg{\\{}\gamma_{0}+2\gamma_{1}-(n(n+1)+4)\gamma_{2}+\left(\frac{n(n+1)+2}{2}\right)\gamma_{3}\frac{1}{R_{0}}$
$\displaystyle\qquad\qquad+\gamma_{4}\frac{1}{R_{0}^{2}}-2n(n+1)\gamma_{5}\frac{1}{R_{0}}\bigg{\\}}\Bigg{]},$
(47d) $\displaystyle I_{\rm p}$
$\displaystyle=\left(\frac{2}{\omega}\right)^{2}\frac{1}{\rho^{\rm
L}R_{0}^{3}}\Bigg{[}(n-1)(n+2)\bigg{\\{}4\gamma_{2}-2\gamma_{3}\frac{1}{R_{0}}+8\gamma_{5}\frac{1}{R_{0}}\bigg{\\}}+8\gamma_{2}+16\gamma_{5}\frac{1}{R_{0}}\Bigg{]},$
(47e) $\displaystyle\omega_{0}^{2}$
$\displaystyle=\left(\frac{2}{\omega}\right)^{2}\frac{1}{\rho^{\rm
L}R_{0}^{2}}\Bigg{[}3kp_{g_{0}}-\frac{2}{R_{0}}\bigg{\\{}\gamma_{0}+2\gamma_{1}-2\gamma_{2}+2\gamma_{3}\frac{1}{R_{0}}+\gamma_{4}\frac{1}{R_{0}^{2}}\bigg{\\}}\Bigg{]},$
(47f) $\displaystyle\omega_{0,n}^{2}$ $\displaystyle=(n+1)W,$ (47g)
$\displaystyle N_{nij}$
$\displaystyle=\left(\frac{2}{\omega}\right)^{2}\frac{1}{2\rho^{\rm
L}R_{0}^{3}}\Bigg{[}N^{2}_{nij}\gamma_{2}+N^{3}_{nij}\gamma_{3}\frac{1}{R_{0}}+N^{4}_{nij}\gamma_{4}\frac{1}{R_{0}^{2}}+N^{5}_{nij}\gamma_{5}\frac{1}{R_{0}}\Bigg{]}.$
(47h)

The functions representing the viscous dissipation in the spherical
$F_{\nu,0}^{(0)}$ and nonspherical $F_{\nu,0}^{(n)}$ mode equations, at the
lowest order, are given by

$\displaystyle F_{\nu,0}^{(0)}$
$\displaystyle=\left(\frac{2}{\omega}\right)\frac{4\nu^{\rm S}}{R_{0}^{2}},$
(48a) $\displaystyle F_{\nu,0}^{(n)}$
$\displaystyle=\left(\frac{2}{\omega}\right)\frac{2}{R_{0}^{2}}\left[(n+2)\left\\{(2n+1)\nu^{\rm
S}-2n(n+2)\frac{\delta}{R_{0}}\nu^{\rm L}\right\\}\right].$ (48b)

Consistent with the discussions in Guédra and Inserra (2018), Krylov-
Bogoliubov averaging method (Nayfeh, 2008) is followed to find the solutions
to (45) and (46). In order to obtain the solutions of (45) and (46), the
following form of asymptotic expansions for $x(\tau)$ and $s_{n}(\tau)$,
respectively, are assumed

$\displaystyle x(\tau)$ $\displaystyle=x_{0}+\epsilon
x_{1}+\epsilon^{2}x_{2}+\cdot\cdot\cdot\;\;,$ (49) $\displaystyle s_{n}(\tau)$
$\displaystyle=s_{n,0}+\epsilon
s_{n,1}+\epsilon^{2}s_{n,2}+\cdot\cdot\cdot\;\;.$ (50)

For a shape mode $n$ close to a parametric resonance, such that
$\Omega_{0,n}\sim a$ ($a$ is an integer), the solution to
$\mathcal{O}(\epsilon^{0})$ equation of (46) can be assumed to be of the form

$\displaystyle s_{n,0}=S_{n}\cos{\theta_{n}},\quad\theta_{n}=a\tau+\phi_{n}.$
(51)

In this analysis, the first parametric resonance of the shape mode $n$ is
considered exclusively, hence $a=1$. In the case of spherical oscillations of
the bubble, it is assumed that the oscillations are significantly far from the
harmonic resonances, ensuring that $\Omega_{0}\neq 2a$. Considering that the
steady state of the spherical oscillations has already been achieved at the
initial time, the solution to (45) in the non-resonant region can be assumed
to be of the form

$\displaystyle x_{0}=X\cos{\theta_{x}},\quad\theta_{x}=2\tau+\phi_{x}.$ (52)

To account for the dynamics of other modes where $m\neq n$, it is assumed that
these secondary modes are sufficiently far from their parametric resonances
$\Omega_{0,m}\neq a)$ and that these $m$ modes oscillate at the frequency of
the excitation forcing at the lowest order, such that

$\displaystyle
s_{m,0}=S_{m}\cos{\theta_{m}},\quad\theta_{m}=2\tau+\phi_{m}\;\;(m\neq n).$
(53)

In order to apply Krylov-Bogoliubov method, the amplitudes $(X,S_{n},S_{m})$
and phases $(\phi_{x},\phi_{n},\phi_{m})$ are assumed to exhibit slow
variations over time $\tau$. Also, assuming that

$\dot{x}_{0}=-X\sin{\theta_{x}},\quad\dot{s}_{n,0}=-S_{n}\sin{\theta_{n}},\quad\dot{s}_{m,0}={S}_{m}\cos{\theta_{m}}$

leads to the following three set of equations

$\displaystyle\dot{X}\cos\theta_{x}-X\sin\theta_{x}\dot{\phi}_{x}=0,$
$\displaystyle\dot{S}_{n}\cos\theta_{n}-S_{n}\sin\theta_{n}\dot{\phi}_{n}=0,$
(54)
$\displaystyle\dot{S}_{m}\cos\theta_{m}-S_{m}\sin\theta_{m}\dot{\phi}_{m}=0.$

Substituting (51), (52) and (53) back into (45) and (46) leads to the
following equations

$\displaystyle-2\dot{X}\sin\theta_{x}-2X\cos\theta_{x}(2+\dot{\phi}_{x})+\Omega_{0}^{2}X\cos\theta_{x}=A\cos
2\tau+\epsilon\mathcal{F}_{1}(Z),$ (55a)
$\displaystyle-\dot{S}_{n}\sin\theta_{n}-S_{n}\cos\theta_{n}(1+\dot{\phi}_{n})+\Omega_{0,n}^{2}S_{n}\cos\theta_{n}=\epsilon\mathcal{F}_{2}(Z),$
(55b)
$\displaystyle-2\dot{S}_{m}\sin\theta_{m}-2S_{m}\cos\theta_{m}(2+\dot{\phi}_{m})+\Omega_{0,m}^{2}S_{m}\cos\theta_{m}=\epsilon\mathcal{F}_{2}(Z).$
(55c)

where
$Z=\\{X\cos\theta_{x},S_{n}\cos\theta_{n},S_{m}\cos\theta_{m},-X\sin\theta_{x},-S_{n}\sin\theta_{n},-S_{m}\sin\theta_{m}\\}$.
Solving the system of equations (54) and (55) leads to six set of first order
ordinary differential equations in terms of
$(X,\phi_{x},S_{n},\phi_{n},S_{m},\phi_{m})$ as

$\displaystyle-2\dot{X}$
$\displaystyle=\sin\theta_{x}\Big{[}(4-\Omega_{0}^{2})X\cos\theta_{x}+A\cos
2\tau+\epsilon\mathcal{F}_{1}(Z)\Big{]},$ (56a)
$\displaystyle-2X\dot{\phi}_{x}$
$\displaystyle=\cos\theta_{x}\Big{[}(4-\Omega_{0}^{2})X\cos\theta_{x}+A\cos
2\tau+\epsilon\mathcal{F}_{1}(Z)\Big{]},$ (56b) $\displaystyle-\dot{S}_{n}$
$\displaystyle=\sin\theta_{n}\Big{[}\big{(}1-\Omega_{0,n}^{2}\big{)}S_{n}\cos\theta_{n}+\epsilon\mathcal{F}_{2}(Z)\Big{]},$
(56c) $\displaystyle-S_{n}\dot{\phi}_{n}$
$\displaystyle=\cos\theta_{n}\Big{[}\big{(}1-\Omega_{0,n}^{2}\big{)}S_{n}\cos\theta_{n}+\epsilon\mathcal{F}_{2}(Z)\Big{]},$
(56d) $\displaystyle-2\dot{S}_{m}$
$\displaystyle=\sin\theta_{m}\Big{[}\big{(}4-\Omega_{0,m}^{2}\big{)}S_{m}\cos\theta_{m}+\epsilon\mathcal{F}_{2}(Z)\Big{]},$
(56e) $\displaystyle-2S_{m}\dot{\phi}_{m}$
$\displaystyle=\cos\theta_{m}\Big{[}\big{(}4-\Omega_{0,m}^{2}\big{)}S_{n}\cos\theta_{m}+\epsilon\mathcal{F}_{2}(Z)\Big{]}.$
(56f)

The first approximation of the Krylov-Bogoliubov begins with the Fourier
expansion of the right hand side of (56) in terms of $\theta_{x},\theta_{n}$
and $\theta_{m}$ and identifying the leading order term. These terms are also
called as averaged equations or slow-time equations given by

$\displaystyle-4\dot{X}$
$\displaystyle=2F_{\nu,0}^{(0)}+A\sin{\phi_{x}}+S_{n}^{2}\frac{\beta_{n}}{2}\sin{(\phi_{x}-2\phi_{n})},$
(57a) $\displaystyle-4X\dot{\phi}_{x}$
$\displaystyle=(4-\Omega_{0}^{2})X+A\cos{\phi_{x}}+S_{n}^{2}\frac{\beta_{n}}{2}\cos{(\phi_{x}-2\phi_{n})},$
(57b) $\displaystyle-2\dot{S}_{n}$
$\displaystyle=F_{\nu,0}^{(n)}S_{n}+\frac{S_{n}}{2}\Bigg{[}(n-1)A\sin{2\phi_{n}}+\zeta_{n}X\sin{(2\phi_{n}-\phi_{x})}$
$\displaystyle\quad+\sum_{m\neq
n}W_{nm}S_{m}\sin{(2\phi_{n}-\phi_{m})}\Bigg{]}$ (57c)
$\displaystyle-2S_{n}\dot{\phi}_{n}$
$\displaystyle=(1-\Omega_{0,n}^{2})S_{n}+\frac{S_{n}}{2}\Bigg{[}(n-1)A\cos{2\phi_{n}}+\zeta_{n}X\cos{(2\phi_{n}-\phi_{x})}$
$\displaystyle\quad+\sum_{m\neq
n}W_{nm}S_{m}\cos{(2\phi_{n}-\phi_{m})}\Bigg{]},$ (57d)
$\displaystyle-4\dot{S}_{m}$
$\displaystyle=2F_{\nu,0}^{(m)}S_{m}+O_{mn}\frac{S_{n}^{2}}{2}\sin{(\phi_{m}-2\phi_{n})},$
(57e) $\displaystyle-4S_{m}\dot{\phi}_{m}$
$\displaystyle=(4-\Omega_{0,m}^{2})S_{m}+O_{mn}\frac{S_{n}^{2}}{2}\cos{(\phi_{m}-2\phi_{n})},$
(57f)

where the coefficients $\beta_{n},\zeta_{n},W_{nm},O_{mn}$ are listed in
Appendix C. In the subsequent section, the first-order differential equations
in (57) are used to provide the steady-state solutions and derive analytical
expressions for the amplitude of shape oscillations.

##  VI Steady-state solutions and stability analysis

The set of differential equations in (57) is written in terms of a complex
form by introducing the complex variables
$\bar{X}=X\text{e}^{\text{i}\phi_{x}}$,
$\bar{S}_{n}=S_{n}\text{e}^{\text{i}\phi_{n}}$ and
$\bar{S}_{m}=X\text{e}^{\text{i}\phi_{m}}$ to get (Guédra and Inserra, 2018)

$\displaystyle 4\text{i}\dot{\bar{X}}$
$\displaystyle=\bar{\Delta}_{0}\bar{X}+A+\frac{\beta_{n}}{2}\bar{S}_{n}^{2},$
(58a) $\displaystyle 2\text{i}\dot{\bar{S}}_{n}$
$\displaystyle=\bar{\Delta}_{n}\bar{S}_{n}+\left[\zeta_{n}\bar{X}+(n-1)A+\sum_{m\neq
n}W_{nm}\bar{S}_{m}\right]\frac{\bar{S}_{n}^{*}}{2},$ (58b) $\displaystyle
4\text{i}\dot{\bar{S}}_{m}$
$\displaystyle=\bar{\Delta}_{m}\bar{S}_{m}+O_{mn}\frac{\bar{S}_{n}^{2}}{2},$
(58c)

where $*$ denotes the conjugate, and the complex quantities
$\bar{\Delta}_{(\cdot)}$ are given by

$\displaystyle\bar{\Delta}_{0}$
$\displaystyle=4-\Omega_{0}^{2}-2\text{i}F_{v,0}^{(0)},$ (59a)
$\displaystyle\bar{\Delta}_{n}$
$\displaystyle=1-\Omega_{0,n}^{2}-\text{i}F_{v,0}^{(n)},$ (59b)
$\displaystyle\bar{\Delta}_{m}$
$\displaystyle=4-\Omega_{0,m}^{2}-2\text{i}F_{v,0}^{(m)}.$ (59c)

In the steady-state regime, the first temporal derivatives in equations (58a)
and (58c) can be cancelled and the complex amplitudes of the spherical and
secondary shape modes can directly be determined as

$\displaystyle\bar{X}$
$\displaystyle=-\bar{\Delta}_{0}^{-1}\left(A+\frac{\beta_{n}}{2}\bar{S}_{n}^{2}\right),$
(60) $\displaystyle\bar{S}_{m}$
$\displaystyle=-\bar{\Delta}_{m}^{-1}O_{mn}\frac{\bar{S}_{n}^{2}}{2}.$ (61)

Further, the above relations in (60) and (61) are substituted in (58b) to get

$\displaystyle\left[2\bar{\Delta}_{n}-\frac{S_{n}^{2}}{2}\left(\frac{\zeta_{n}\beta_{n}}{\bar{\Delta}_{0}}+\sum_{m\neq
n}\frac{W_{nm}O_{mn}}{\bar{\Delta}_{m}}\right)\right]\bar{S}_{n}=-\left[(n-1)-\frac{\zeta_{n}}{\bar{\Delta}_{0}}\right]A\bar{S}_{n}^{*}.$
(62)

The trivial solution for the above equation is $S_{n}=0$, such that

$\displaystyle\bigg{\rvert}\bar{y}_{1}-\frac{S_{n}^{2}}{2}\bar{y}_{2}\bigg{\rvert}=\rvert\bar{y}_{3}\rvert
A,$ (63)
$\displaystyle\text{e}^{2\text{i}\phi_{n}}=-\frac{\bar{y}_{3}A}{\bar{y}_{1}-\cfrac{S_{n}^{2}}{2}\bar{y}_{2}}.$
(64)

Here the complex quantities $\bar{y}_{(\cdot)}$ are defined as

$\displaystyle\bar{y}_{1}$ $\displaystyle=2\bar{\Delta}_{n},$ (65)
$\displaystyle\bar{y}_{2}$
$\displaystyle=\bar{\Delta}_{0}^{-1}\zeta_{n}\beta_{n}+\sum_{m\neq
n}\bar{\Delta}_{m}^{-1}W_{nm}O_{mn},$ (66) $\displaystyle\bar{y}_{3}$
$\displaystyle=(n-1)-\bar{\Delta}_{0}^{-1}\zeta_{n}.$ (67)

It is worth mentioning that the solutions obtained in this analysis include
the onset threshold for the parametric excitation of the shape mode $n$ in the
vicinity of the first resonance. Following a similar approach as discussed in
Guédra and Inserra (2018), an analytical expression for the absolute stability
threshold (in terms of driving amplitude pressure) is derived by neglecting
the quadratic terms in (60) and (63), which can be written as

$\displaystyle
A_{\text{th}}=\Bigg{\rvert}\frac{\bar{y}_{1}}{\bar{y}_{3}}\Bigg{\rvert}=\cfrac{2\rvert\bar{\Delta}_{n}\rvert\rvert\bar{\Delta}_{0}\rvert}{\rvert\zeta_{n}-(n-1)\bar{\Delta}_{0}\rvert}.$
(68)

Similarly, an expression can also be deduced in terms of the spherical mode
driving amplitude

$\displaystyle
X_{\text{th}}=\cfrac{2\rvert\bar{\Delta}_{n}\rvert}{\rvert\zeta_{n}-(n-1)\bar{\Delta}_{0}\rvert}.$
(69)

Figure 5: The plot shows the driving amplitude $A$ as a function of the
frequency $\Omega_{0,n}$, and $n=2$ for the EB with $R_{0}=2\,\mu$m. Curves
corresponds to a specific stability threshold and fixed values of $S_{n}$ in
the steady state with the set of interface parameters in table 1.

Upon comparing the relations given in (68) and (69) with those in Guédra and
Inserra (2018), it is evident that the structure of these relations remains
unchanged. However, the definitions of the involved quantities differ due to
the inclusion of interface parameters and the shell elasticity constant in the
present model. The steady-state solutions of (63) can be written as

$\displaystyle
S_{n}=\frac{\sqrt{2h_{0}}}{y_{2}}\left[1\pm\sqrt{1+\cfrac{y_{2}^{2}(y_{3}A-y_{1})(y_{3}A+y_{1})}{h_{0}^{2}}}\right]^{1/2}.$
(70)

Here $y_{i}=\rvert\bar{y}_{i}\rvert$, and

$\displaystyle
h_{0}=\Re(\bar{y}_{1})\Re(\bar{y}_{2})+\Im(\bar{y}_{1})\Im(\bar{y}_{2}).$ (71)

The conditional stability threshold $(A^{\prime}_{\text{th}})$ is the
condition for $A$ that ensures the solutions are real $(A>0)$, then

$\displaystyle
A>A^{{}^{\prime}}_{\text{th}}=\frac{1}{y_{3}}\sqrt{\left(y_{1}-\frac{h_{0}}{y_{2}}\right)\left(y_{1}+\frac{h_{0}}{y_{2}}\right)}.$
(72)

The driving amplitude $(A)$ corresponding to isovalues of shape mode amplitude
$(S_{n})$ can directly be calculated using the relation in (63). The absolute
stability threshold $(A_{\text{th}})$ and the conditional stability threshold
$(A^{\prime}_{\text{th}})$ can be calculated from (68) and (72), respectively.
The driving amplitudes for various values of $S_{n}$ (isolines), along with
the absolute and conditional stability thresholds, are plotted against
$\Omega_{0,n}$ for $n=2$ and $n=3$ in figures 5 and 6, respectively.

Figure 6: The plot shows the driving amplitude $A$ as a function of the
frequency $\Omega_{0,n}$, and $n=3$ for the EB with $R_{0}=2\,\mu$m. Curves
corresponds to a specific stability threshold and fixed values of $S_{n}$ in
the steady state with the set of interface parameters in table 1.

In figure 5, the point marked with $(\bullet)$ on the isoline $S_{n}=0.2$
represents the value of $A$ corresponding to $\Omega_{0,n}\sim 1$. By using
the definitions of $A$ and $\omega_{0,n}^{2}$ as given in (47), the excitation
pressure and frequency values are determined for this specific point.
Subsequently, these calculated pressure and frequency values are considered in
the direct numerical simulations, which lead to the determination of the
finite shape mode amplitudes illustrated in figure 2. It is also apparent that
the amplitude of oscillation of shape mode $\epsilon a_{2}=0.2$ corresponds to
the isoline $S_{n}=0.2$.

A similar analysis is conducted for $n=3$, at the point denoted by $(\bullet)$
on the isoline $S_{n}=0.3$ in figure 6, and the calculated values of pressures
and frequency is used in the direct numerical simulations depicted in figure
3. It is observed that a slight change in the excitation pressure would still
yield in finite amplitude oscillations in the time-series analysis, while a
significant change would cause the finite non-zero amplitudes to drop to zero.
To study the temporal variation, it is crucial to ensure that the simulations
are conducted within and around the conditional stable zone close to the
corresponding $\Omega_{0,n}\sim 1$. Hence, for a specific value of $n$,
conditional stability plots can be generated, and based on those, appropriate
working values of excitation pressure and frequency are calculated at
$\Omega_{0,n}\sim 1$ to conduct further temporal evolution analysis.

As an extended analysis, we consider slightly larger radii bubbles with an
initial radius of $R_{0}=5\,\mu$m and a shell thickness of $q=20\,$nm. The
stability curves for $n=2$ using the same interface material parameters as
before are illustrated in figure 7. Subsequently, direct numerical simulations
are conducted for this case, using the calculated values of excitation
pressure and frequency, as depicted in figure 4. Although the analysis remains
unchanged, the larger bubble radius leads to a decrease in the effects of
interface parameters.

Figure 7: The plot shows the driving amplitude $A$ as a function of the
frequency $\Omega_{0,n}$, and $n=2$ for the EB with $R_{0}=5\,\mu$m. Curves
corresponds to a specific stability threshold and fixed values of $S_{n}$ in
the steady state with the set of interface parameters in table 1.

It is essential to emphasize that the direct numerical simulations of the
governing equations are influenced by the choice of numerical values of
interface parameters and the viscosity of the shell. In certain cases, these
simulations may exhibit a blow-up of the solution, indicating the presence of
a finite time singularity in the problem. One such possibility is when there
are no interface parameters. The blow-up of the solution can also be
interpreted as an unstable configuration of the EB, because it represents a
sudden and uncontrollable escalation of the bubble’s oscillations, which can
lead to unpredictable behavior and potential collapse of the bubble. The
introduction of interface effects through the surface tension parameter plays
an important role in averting this blow-up. However, it’s crucial to carefully
balance the benefits of these interface parameters against the potential
complexity introduced by them in the system. Hence, the behavior of the bubble
is strongly affected by these interface parameters, a finding that has also
been reported in other interface energy models studying the radial dynamics of
EBs.(Dash and Tamadapu, 2022b, a) Properly chosen interface parameters
effectively address and prevent the occurrence of finite time singularities in
the problem, as demonstrated in the current study. The interface parameters
modify the EB’s behavior, contributing to stability and preventing the
undesirable blow-up phenomenon. Additionally, it is worth noting that higher
shell viscosity values can lead to the emergence of such finite time
singularities and special numerical techniques are required to handle these
singularities.

##  VII Conclusion

This work is focused on developing a mathematical model to investigate the
nonspherical oscillations of smaller radii EB suspended in fluid. This model
takes into account the essential interfacial mechanics at the gas-
encapsulation and encapsulation-fluid interfaces, which contributes
significantly to the mechanics of smaller radii bubbles. (Dash and Tamadapu,
2022a, b) The shell material is treated as a thin membrane with both elastic
and viscous effects. The coupled dynamical equations governing the spherical
and shape mode oscillations of the EB are derived using Lagrangian energy
formulation. These governing equations are then analyzed using direct
numerical simulations.

The study reveals that the parametrically forced even mode $(n=2)$ excites
only the even modes, while the odd modes $(n=3)$ excites both even and odd
modes. This has also been observed in the larger radii bubbles. (Guédra and
Inserra, 2018) But in the case of smaller radii EBs, the interface parameters
play a crucial role in the analysis, unlike in larger bubbles. For instance,
Guédra and Inserra (2018) reported the finite amplitude oscillations of larger
radii bubbles with conventional surface tension parameter without considering
the effect of interface energy. But in the present analysis with smaller radii
bubbles, we observe that the interface parameters are essential in identifying
the stable finite amplitude shape mode oscillations. In the absence of
interface parameters, the EB demonstrates an unstable behavior, possibly
causing solutions to blow up. However, the EB shows finite amplitude
oscillations when interface parameters are introduced. The interface energy
models describing the radial dynamics of EBs have proven the significant
influence of interface parameters.(Dash and Tamadapu, 2022a, b) The present
analysis further highlights this notion, illustrating that interface
parameters play an even more substantial and influential role when it comes to
nonspherical oscillations. This is substantiated by the fact that finite
amplitude shape mode oscillations are observed in smaller radii EBs only when
all interface parameters are considered; otherwise, they are not present.
Alongside the interface parameters, the elasticity and viscosity of the shell
material are also essential factors in the numerical simulations.
Additionally, perturbation analysis has been used to derive the equations and
analyze the system in the vicinity of the first parametric resonance of shape
mode $n$. The perturbation method based on the Krylov-Bogoliubov method of
averaging and steady-state analysis proves invaluable in calculating the
required excitation pressure and frequency for the existence of shape mode
oscillations. It helps reveal insights into how the resonance and the
introduced perturbations affect the EB’s stability and amplitude of
oscillations. This information is particularly valuable in determining optimal
excitation pressures and frequencies for small radii bubbles, which hold
significant promise for various medical applications.

Future work includes several directions. The present study can be extended to
investigate the nonspherical oscillations of EBs suspended in any biological
fluids, which can be assumed as linear or nonlinear viscoelastic fluids.
(Allen and Roy, 2000a, b) In the context of biomedical applications, when EBs
approach blood vessels, their dynamic behavior alters. (Doinikov, Zhao, and
Dayton, 2009; Suslov, Ooi, and Manasseh, 2012) Another interesting extension
is exploring the nonspherical oscillations of EBs near blood vessel walls
using the interface energy model. Interesting recent developments in bubble
dynamics involve strategies like coating bubbles with magnetic nanoparticles
to enhance their targeting efficiency in drug delivery.(Stride  _et al._ ,
2009; Jamburidze  _et al._ , 2019; Malvar, Gontijo, and Cunha, 2018) Moreover,
by suspending EBs in magnetic fluids, researchers gain insights into how they
respond under the combined influence of acoustic and magnetic fields. The
interface energy model can also be used to study the behavior of these
magnetic microbubbles.

## Acknowledgements

We are grateful to P. G. Senapathy Center for Computing Resource, IIT Madras,
for providing the computational resources. This work was supported by the
Science and Engineering Research Board (SERB), Department of Science and
Technology, Government of India, under the project number CRG/2022/005775.

##  Appendix A Basis invariants for the interface

The interface energy density can be expressed in terms of six basis invariants
of the right Cauchy-Green interface deformation tensor $\bm{C}$, the relative
curvature tensor $\bm{\kappa}$, and the permutation tensor-density $\bm{\mu}$
on the undeformed interface as

$\displaystyle\gamma=\gamma(I_{1},I_{2},I_{3},I_{4},I_{5},I_{6}),$ (73)

and the invariants are given by

$\displaystyle\begin{split}I_{1}&={\rm
tr}\,\bm{C}=G^{\alpha\beta}C_{\alpha\beta},\\\ I_{2}&={\rm
det}\,\bm{C}=J^{2}=g/G,\\\ I_{3}&={\rm
tr}{\,\bm{\kappa}}=G^{\alpha\beta}\kappa_{\alpha\beta},\\\ I_{4}&={\rm
det}{\,\bm{\kappa}}=\frac{1}{2}\mu^{\alpha\beta}\mu^{\gamma\delta}\kappa_{\alpha\gamma}\kappa_{\beta\delta},\\\
I_{5}&={\rm
tr}(\bm{C}\bm{\kappa})=C_{\alpha\beta}\kappa^{\alpha\beta}=C^{\alpha\beta}\kappa_{\alpha\beta},\\\
I_{6}&={\rm
tr}(\bm{C}\bm{\kappa}\bm{\mu})=G_{\alpha\beta}C_{\gamma\delta}\kappa^{\alpha\gamma}\mu^{\beta\delta}=G_{\alpha\beta}\kappa_{\gamma\delta}C^{\alpha\gamma}\mu^{\beta\delta}.\end{split}$
(74)

##  Appendix B Nonlinear interaction coefficients

The integrals and the nonlinear interaction coefficients are defined through
the following where $\mu=\cos{\phi}$:

$\displaystyle I_{a_{nij}}$
$\displaystyle=\int_{-1}^{1}P_{n}P_{i}P_{j}\text{d}\mu,$ (75) $\displaystyle
I_{c_{nij}}$
$\displaystyle=-\int_{-1}^{1}(1-\mu^{2})P_{n}\frac{\text{d}P_{i}}{\text{d}\mu}\frac{\text{d}P_{j}}{\text{d}\mu}\text{d}\mu,$
(76) $\displaystyle I_{g_{nij}}$
$\displaystyle=-\int_{-1}^{1}\mu(1-\mu^{2})\frac{\text{d}P_{n}}{\text{d}\mu}\frac{\text{d}P_{i}}{\text{d}\mu}\frac{\text{d}P_{j}}{\text{d}\mu}\text{d}\mu,$
(77) $\displaystyle L_{a_{nij}}$
$\displaystyle=(n+1)(n+i+1)I_{a_{nij}}+I_{c_{nij}},$ (78) $\displaystyle
L_{b_{nij}}$
$\displaystyle=\frac{1}{2}(n+1)(n^{2}+4n+4)I_{a_{nij}}+(n+2)I_{c_{inj}},$ (79)
$\displaystyle L_{c_{nij}}$
$\displaystyle=(n+i+3)(n+i+2)I_{a_{nij}}+I_{c_{nij}}+I_{c_{inj}},$ (80)
$\displaystyle G_{c_{nij}}$
$\displaystyle=\frac{2n^{2}-n+1}{n+1}I_{a_{nij}}-\frac{2}{n+1}\left(L_{b_{nij}}-\frac{2L_{a_{nij}}}{i+1}\right),$
(81) $\displaystyle G_{d_{nij}}$
$\displaystyle=\frac{2L_{c_{nji}}}{(n+1)(j+1)}-\frac{L_{b_{jin}}-j(j-1)I_{a_{nij}}}{j+1}-\frac{4(n-1)}{n+1}I_{a_{nij}},$
(82) $\displaystyle G_{e_{nij}}$
$\displaystyle=\frac{\left[n(n-j)+3(j+1)\right]I_{a_{nij}}+I_{c_{nji}}}{(n+1)(j+1)},$
(83) $\displaystyle M_{a_{nij}}$
$\displaystyle=G_{d_{ijn}}-G_{c_{nij}}-G_{c_{inj}}-G_{c_{ijn}},$ (84)
$\displaystyle M_{b_{nij}}$
$\displaystyle=G_{d_{jin}}+G_{d_{ijn}}+2\left(G_{e_{nij}}+G_{e_{jin}}\right)-G_{d_{nij}}-G_{d_{inj}},$
(85) $\displaystyle M_{c_{nij}}$ $\displaystyle=G_{e_{nij}}+G_{e_{jin}},$ (86)
$\displaystyle M_{d_{nij}}$
$\displaystyle=G_{e_{nij}}+G_{e_{jin}}-G_{e_{inj}}$ (87) $\displaystyle
N^{2}_{nij}$ $\displaystyle=12I_{a_{nij}}-2I_{c_{nij}}-4I_{c_{ijn}},$ (88)
$\displaystyle N^{3}_{nij}$
$\displaystyle=\left[\left(\frac{n(n+1)}{2}+1\right)I_{c_{nij}}+(i(i+1)+2)I_{c_{ijn}}\right],$
(89) $\displaystyle N^{4}_{nij}$
$\displaystyle=n(n+1)I_{c_{nij}}+2i(i+1)I_{c_{ijn}}-6I_{g_{nij}},$ (90)
$\displaystyle N^{5}_{nij}$
$\displaystyle=\\{6+i(i+1)+n(n+1)\\}I_{a_{nij}}-\left(\frac{n(n+1)}{2}+2\right)I_{c_{nij}}$
$\displaystyle\qquad\qquad-(i(i+1)+4)I_{c_{ijn}}+3I_{g_{nij}}.$ (91)

##  Appendix C Coefficients appearing in first-order differential equations

In the first order differential equations (57), the coefficients are given by:

$\displaystyle\beta_{n}$
$\displaystyle=\frac{(n+3)\Omega_{0,n}^{2}-(n+1)U+\left(n+\frac{3}{2}\right)\dot{\theta}_{n}^{2}-B(n+1)\left[n(n+1)-10\right]}{\mathcal{H}_{n}},$
(92) $\displaystyle\zeta_{n}$
$\displaystyle=3\Omega_{0,n}^{2}-(n-1)\Omega_{0}^{2}-3\dot{\theta}_{x}\dot{\theta}_{n}+2(n+1)\left[14-n(n+1)\right]B-(n+1)I_{\rm
p},$ (93) $\displaystyle{W_{nm}}$
$\displaystyle=\frac{\mathcal{H}_{n}}{4}\bigg{[}W_{0}\left(8I_{a_{nnm}}\right)+\Omega_{0,m}^{2}M_{c_{nnm}}+\Omega_{0,n}^{2}M_{c_{nmn}}-2\left(M_{d_{nnm}}+M_{d_{nmn}}\right)$
$\displaystyle\quad+8B\left(5I_{a_{nnm}}+I_{c_{nmn}}+\frac{I_{c_{nnm}}}{2}\right)+{4(N_{nnm}+N_{nmn})}\bigg{]},$
(94) $\displaystyle{O_{mn}}$
$\displaystyle=\frac{\mathcal{H}_{m}}{4}\bigg{[}W_{0}\left(4I_{a_{mnn}}\right)+\Omega_{0,n}^{2}M_{c_{mnn}}+M_{d_{mnn}}+4B\left(5I_{a_{mnn}}+I_{c_{nnm}}+\frac{I_{c_{mnn}}}{2}\right)+{4N_{mnn}}\bigg{]}.$
(95)

## References

  * Hoff (2001) L. Hoff,  _Acoustic characterization of contrast agents for medical ultrasound imaging_ (Springer Science & Business Media, 2001).
  * Stride and Saffari (2003) E. Stride and N. Saffari, “Microbubble ultrasound contrast agents: A review,” Proc. Inst. Mech. Eng., Part H: J. Eng. Med. 217, 429–447 (2003).
  * Postema  _et al._ (2004) M. Postema, A. Van Wamel, C. T. Lancée, and N. De Jong, “Ultrasound-induced encapsulated microbubble phenomena,” Ultrasound Med. Biol. 30, 827–840 (2004).
  * Lindner (2004) J. R. Lindner, “Microbubbles in medical imaging: current applications and future directions,” Nat. Rev. Drug Discovery 3, 527–533 (2004).
  * Klibanov (2006) A. L. Klibanov, “Microbubble contrast agents,” Invest. Radiol. 41, 354–362 (2006).
  * Errico  _et al._ (2015) C. Errico, J. Pierre, S. Pezet, Y. Desailly, Z. Lenkei, O. Couture, and M. Tanter, “Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging,” Nat. 527, 499–502 (2015).
  * Tachibana and Tachibana (1999) K. Tachibana and S. Tachibana, “Application of ultrasound energy as a new drug delivery system,” Jpn. J. Appl. Phys. 38, 3014 (1999).
  * Tsutsui, Xie, and Porter (2004) J. M. Tsutsui, F. Xie, and R. T. Porter, “The use of microbubbles to target drug delivery,” Cardiovasc. Ultrasound 2 (2004), 10.1186/1476-7120-2-23.
  * Hernot and Klibanov (2008) S. Hernot and A. L. Klibanov, “Microbubbles in ultrasound-triggered drug and gene delivery,” Adv. Drug Delivery Rev. 60, 1153 – 1166 (2008).
  * Hynynen (2008) K. Hynynen, “Ultrasound for drug and gene delivery to the brain,” Adv. Drug Delivery Rev. 60, 1209 – 1217 (2008).
  * Kooiman  _et al._ (2014) K. Kooiman, H. J. Vos, M. Versluis, and N. de Jong, “Acoustic behavior of microbubbles and implications for drug delivery,” Adv. Drug Delivery Rev. 72, 28 – 48 (2014).
  * Suslick (1990) K. S. Suslick, “Sonochemistry,” Sci. 247, 1439–1445 (1990).
  * Blake (1999) J. R. Blake, “Preface to acoustic cavitation and sonoluminescence,” Philos. Trans. R. Soc. London, Ser. A 357, 201–201 (1999).
  * Agarwal  _et al._ (2012) A. Agarwal, H. Xu, W. J. Ng, and Y. Liu, “Biofilm detachment by self-collapsing air microbubbles: a potential chemical-free cleaning technology for membrane biofouling,” J. Mater. Chem. 22, 2203–2207 (2012).
  * Seo  _et al._ (2018) Y. Seo, J. Leong, J. D. Park, Y. T. Hong, S. H. Chu, C. Park, D. H. Kim, Y. H. Deng, V. Dushnov, J. Soh, S. Rogers, Y. Y. Yang, and H. Kong, “Diatom microbubbler for active biofilm removal in confined spaces,” ACS Appl. Mater. Interfaces 10, 35685–35692 (2018).
  * Ohl  _et al._ (2006) C. D. Ohl, M. Arora, R. Dijkink, V. Janve, and D. Lohse, “Surface cleaning from laser-induced cavitation bubbles,” Appl. Phys. Lett. 89, 074102 (2006).
  * Reuter  _et al._ (2017) F. Reuter, S. Lauterborn, R. Mettin, and W. Lauterborn, “Membrane cleaning with ultrasonically driven bubbles,” Ultrason. Sonochem. 37, 542 – 560 (2017).
  * Lee  _et al._ (2015) M. Lee, E. Y. Lee, D. Lee, and B. J. Park, “Stabilization and fabrication of microbubbles: applications for medical purposes and functional materials,” Soft matter 11, 2067–2079 (2015).
  * de Jong  _et al._ (1992) N. de Jong, L. Hoff, T. Skotland, and N. Bom, “Absorption and scatter of encapsulated gas filled microspheres: Theoretical considerations and some measurements,” Ultrasonics 30, 95 – 103 (1992).
  * de Jong and Hoff (1993) N. de Jong and L. Hoff, “Ultrasound scattering properties of albunex microspheres,” Ultrasonics 31, 175 – 181 (1993).
  * Church (1995) C. C. Church, “The effects of an elastic solid surface layer on the radial pulsations of gas bubbles,” J. Acoust. Soc. Am. 97, 1510–1521 (1995).
  * Brenner, Lohse, and Dupont (1995) M. P. Brenner, D. Lohse, and T. F. Dupont, “Bubble shape oscillations and the onset of sonoluminescence,” Phys. Rev. Lett. 75, 954 (1995).
  * Hao and Prosperetti (1999) Y. Hao and A. Prosperetti, “The effect of viscosity on the spherical stability of oscillating gas bubbles,” Phys. Fluids 11, 1309–1317 (1999).
  * Versluis  _et al._ (2004) M. Versluis, S. M. van der Meer, D. Lohse, P. Palanchon, D. Goertz, C. T. Chin, and N. de Jong, “Microbubble surface modes [ultrasound contrast agents],” in  _IEEE Ultrasonics Symposium, 2004_, Vol. 1 (2004) pp. 207–209 Vol.1.
  * van der Meer  _et al._ (2006) S. M. van der Meer, B. Dollet, D. E. Goertz, N. de Jong, M. Versluis, and D. Lohse, “Surface modes of ultrasound contrast agent microbubbles,” in  _2006 IEEE Ultrasonics Symposium_ (2006) pp. 112–115.
  * Dollet  _et al._ (2008) B. Dollet, S. M. van der Meer, V. Garbin, N. de Jong, D. Lohse, and M. Versluis, “Nonspherical oscillations of ultrasound contrast agent microbubbles,” Ultrasound Med. Biol. 34, 1465–1473 (2008).
  * Vos  _et al._ (2011) H. J. Vos, B. Dollet, M. Versluis, and N. de Jong, “Nonspherical shape oscillations of coated microbubbles in contact with a wall,” Ultrasound Med. Biol. 37, 935–948 (2011).
  * Liu  _et al._ (2011) Y. Liu, K. Sugiyama, S. Takagi, and Y. Matsumoto, “Numerical study on the shape oscillation of an encapsulated microbubble in ultrasound field,” Phys. Fluids 23, 041904 (2011).
  * Tsiglifis and Pelekasis (2011) K. Tsiglifis and N. A. Pelekasis, “Parametric stability and dynamic buckling of an encapsulated microbubble subject to acoustic disturbances,” Phys. Fluids 23, 012102 (2011).
  * Tamadapu, Grishenkov, and Eriksson (2016) G. Tamadapu, D. Grishenkov, and A. Eriksson, “Modeling and parametric investigation of thick encapsulated microbubble’s nonspherical oscillations,” J. Acoust. Soc. Am. 140, 3884–3895 (2016).
  * Hilgenfeldt, Lohse, and Brenner (1996) S. Hilgenfeldt, D. Lohse, and M. P. Brenner, “Phase diagrams for sonoluminescing bubbles,” Phys. Fluids 8, 2808–2826 (1996).
  * Brenner, Hilgenfeldt, and Lohse (2002) M. P. Brenner, S. Hilgenfeldt, and D. Lohse, “Single-bubble sonoluminescence,” Rev. Mod. Phys. 74, 425–484 (2002).
  * Loughran, Eckersley, and Tang (2012) J. Loughran, R. J. Eckersley, and M.-X. Tang, “Modeling non-spherical oscillations and stability of acoustically driven shelled microbubbles,” J. Acoust. Soc. Am. 131, 4349–4357 (2012).
  * Doinikov (2004) A. A. Doinikov, “Translational motion of a bubble undergoing shape oscillations,” J. Fluid Mech. 501, 1–24 (2004).
  * Shaw (2006) S. J. Shaw, “Translation and oscillation of a bubble under axisymmetric deformation,” Phys. Fluids 18, 072104 (2006).
  * Shaw (2009) S. Shaw, “The stability of a bubble in a weakly viscous liquid subject to an acoustic traveling wave,” Phys. Fluids 21, 022104 (2009).
  * Shaw (2017) S. J. Shaw, “Nonspherical sub-millimeter gas bubble oscillations: Parametric forcing and nonlinear shape mode coupling,” Phys. Fluids 29, 122103 (2017).
  * Guédra and Inserra (2018) M. Guédra and C. Inserra, “Bubble shape oscillations of finite amplitude,” J. Fluid Mech. 857, 681–703 (2018).
  * Dash and Tamadapu (2022a) N. Dash and G. Tamadapu, “Radial dynamics of an encapsulated microbubble with interface energy,” J. Fluid Mech. 932, A26 (2022a).
  * Steigmann and Ogden (1999) D. J. Steigmann and R. W. Ogden, “Elastic surface-substrate interactions,” Proc. R. Soc. London, Ser. A 455, 437–474 (1999).
  * Gao  _et al._ (2014) X. Gao, Z. Huang, J. Qu, and D. Fang, “A curvature-dependent interfacial energy-based interface stress theory and its applications to nano-structured materials: (i) general theory,” J. Mech. Phys. Solids 66, 59 – 77 (2014).
  * Cattaneo and Supponen (2023) M. Cattaneo and O. Supponen, “Shell viscosity estimation of lipid-coated microbubbles,” Soft Matter 19, 5925–5941 (2023).
  * Dash and Tamadapu (2022b) N. Dash and G. Tamadapu, “Describing the dynamics of a nonlinear viscoelastic shelled microbubble with an interface energy model,” J. Appl. Phys. 132, 204702 (2022b).
  * Fung (1977) Y. C. Fung,  _A First Course in Continuum Mechanics_ (Prentice Hall, 1977).
  * Nayfeh (2008) A. H. Nayfeh,  _Perturbation methods_ (John Wiley & Sons, 2008).
  * Allen and Roy (2000a) J. S. Allen and R. A. Roy, “Dynamics of gas bubbles in viscoelastic fluids. ii. nonlinear viscoelasticity,” J. Acoust. Soc. Am. 108, 1640–1650 (2000a).
  * Allen and Roy (2000b) J. S. Allen and R. A. Roy, “Dynamics of gas bubbles in viscoelastic fluids. i. linear viscoelasticity,” J. Acoust. Soc. Am. 107, 3167–3178 (2000b).
  * Doinikov, Zhao, and Dayton (2009) A. A. Doinikov, S. Zhao, and P. A. Dayton, “Modeling of the acoustic response from contrast agent microbubbles near a rigid wall,” Ultrasonics 49, 195–201 (2009).
  * Suslov, Ooi, and Manasseh (2012) S. A. Suslov, A. Ooi, and R. Manasseh, “Nonlinear dynamic behavior of microscopic bubbles near a rigid wall,” Phys. Rev. E 85, 066309 (2012).
  * Stride  _et al._ (2009) E. Stride, C. Porter, A. G. Prieto, and Q. Pankhurst, “Enhancement of microbubble mediated gene delivery by simultaneous exposure to ultrasonic and magnetic fields,” Ultrasound Med. Biol. 35, 861–868 (2009).
  * Jamburidze  _et al._ (2019) A. Jamburidze, A. Huerre, D. Baresch, V. Poulichet, M. De Corato, and V. Garbin, “Nanoparticle-coated microbubbles for combined ultrasound imaging and drug delivery,” Langmuir 35, 10087–10096 (2019).
  * Malvar, Gontijo, and Cunha (2018) S. Malvar, R. Gontijo, and F. Cunha, “Nonlinear motion of an oscillating bubble immersed in a magnetic fluid,” J. Eng. Math. 108, 143–170 (2018).