# A unifying Rayleigh-Plesset-type equation for bubbles in viscoelastic media

Alexandros T. Oratis  Physics of Fluids Group, Max Planck Center for Complex
Fluid Dynamics, Faculty of Science and Technology, Technical Medical (TechMed)
Center, University of Twente, Enschede, The Netherlands  Kay Dijs  Physics of
Fluids Group, Max Planck Center for Complex Fluid Dynamics, Faculty of Science
and Technology, Technical Medical (TechMed) Center, University of Twente,
Enschede, The Netherlands  Guillaume Lajoinie  Physics of Fluids Group, Max
Planck Center for Complex Fluid Dynamics, Faculty of Science and Technology,
Technical Medical (TechMed) Center, University of Twente, Enschede, The
Netherlands  Michel Versluis  Physics of Fluids Group, Max Planck Center for
Complex Fluid Dynamics, Faculty of Science and Technology, Technical Medical
(TechMed) Center, University of Twente, Enschede, The Netherlands  Jacco H.
Snoeijer<EMAIL_ADDRESS>Physics of Fluids Group, Max Planck Center
for Complex Fluid Dynamics, Faculty of Science and Technology, Technical
Medical (TechMed) Center, University of Twente, Enschede, The Netherlands

###### Abstract

Understanding the ultrasound pressure-driven dynamics of microbubbles confined
in viscoelastic materials is relevant for multiple biomedical applications,
ranging from contrast-enhanced ultrasound imaging to ultrasound-assisted drug
delivery. The volumetric oscillations of spherical bubbles is analyzed using
the Rayleigh-Plesset equation, which describes the conservation of mass and
momentum in the surrounding medium. Several studies have considered an
extension of the Rayleigh-Plesset equation for bubbles embedded into
viscoelastic media, but these are restricted to a particular choice of
constitutive model and/or to small deformations. Here, we derive a unifying
equation applicable to bubbles in viscoelastic media with arbitrary complex
moduli and that can account for large bubble deformations. To derive this
equation, we borrow concepts from finite-strain theory. We validate our
approach by comparing the result of our model to previously published results
and extend it to show how microbubbles behave in arbitrary viscoelastic
materials. In particular, we use our viscoelastic Rayleigh-Plesset model to
compute the bubble dynamics in benchmarked viscoelastic liquids and solids.

##  I Introduction

The growth and collapse of bubbles in viscoelastic materials is a physical
process that features in a variety of environmental, food processing and
industrial settings. From the bubbles formed during volcanic eruptions [26] to
those used in the production of polymer foams [18, 4], understanding the
influence of viscoelasticity on the bubble dynamics is a vital step for
technological innovation. In particular, the crucial importance of
viscoelasticity has been a topic of focus for biomedical applications [61, 22,
15]. Indeed, micrometer-size bubbles are frequently used in ultrasound as
contrast agents to visualize organ perfusion and blood flow down to the
smallest capillary vessels [21]. They are also heavily investigated as
therapeutic agents for cancer treatment, blood-brain-barrier opening [10, 50,
3, 6, 45], sonoporation [32, 47, 14], lithotripsy, histotripsy [60], or even
sonothrombolysis [35, 13, 5]. These therapeutic applications often require
inducing controlled microbubble oscillations within tissue to impart treatment
while preventing collateral damage to the surrounding healthy tissue. However,
this embedding in a polymeric liquid or a soft solid has a dramatic effect on
bubble dynamics [15]. More specifically these materials have viscoelastic
properties and provide solid-like or liquid-like resistance to bubble motion,
depending on the deformation time scale. Therefore, the ability to control and
utilize bubbles in biomedical settings crucially relies on our understanding
of the micro and macro-rheology of the surrounding medium and its impact on
bubble oscillations.

Various biological materials exhibit viscoelastic properties that come with
intricate behavior of stress and strain. A simple method to understand such
responses consists of modelling the viscoelastic material as a combination of
Hookean springs and Newtonian dashpots [8, 7, 29]. The spring introduces
elastic characteristics to the material’s response, while the dashpot adds a
viscous contribution. Models built from a combination of springs and dashpots
typically lead to exponential relaxation of stress or strain. Even though the
spring-dashpot approach accurately describes various viscoelastic liquids and
solids, it fails to capture the rheology of multi-scale materials such as gels
or elastomeric networks [59, 33, 28], which often come along with power-law
stress relaxation rather than exponentials. An alternative and more general
approach consists of characterizing the material via its complex modulus
$\mu(\omega)=G^{\prime}(\omega)+iG^{\prime\prime}(\omega)$, which combines an
elastic storage modulus $G^{\prime}(\omega)$ and viscous loss modulus
$G^{\prime\prime}(\omega)$. These moduli can be experimentally obtained by
measuring the deformation of the material in response to a periodic
excitation, and are functions of the imposed angular frequency $\omega$. As an
example, the complex modulus of gels typically exhibits a non-integer scaling
behavior with respect to frequency, which cannot be represented by (a finite
number of) springs and dashpots. An equivalent formulation of the rheology in
the temporal domain can be achieved via the stress relaxation function
$\psi(t)$. This function describes how the stress relaxes after a material is
subjected to a step-strain, and can be related to the complex modulus as
$\mu(\omega)=i\omega\int_{0}^{\infty}\mathrm{d}t~{}\psi(t)e^{-i\omega t}$ [9];
hence, $\mu(\omega)$ and $\psi(t)$ contain the exact same rheological
information. Within this framework, $\psi\to 0$ at large times for
viscoelastic liquids, reflecting how the stress fully relaxes after stopping
the flow. Viscoelastic solids, by contrast, exhibit a finite $\psi$ at large
times, and thus a residual elastic stress, expressing a long-term memory of
the reference configuration at rest.

The dynamical behavior of spherical bubbles is typically described using the
Rayleigh-Plesset equation, initially developed for a spherical bubbles
oscillating in free-field. Possible ways of adding the effects of the
viscoelastic surrounding on the bubble dynamics in the context of the
Rayleigh-Plesset equation have been the focus of many studies, starting with
the pioneering work of Fogler and Goddard [19, 20]. Assuming a linear Maxwell
model, the authors investigated the effects of elasticity in the cavitation of
a spherical void. Subsequent studies extended this approach by considering
different constitutive models of the surrounding material [51, 46, 1, 61, 11,
24], including those applicable to large-amplitude bubble deformations [53,
41, 30, 2, 27, 39, 22, 57]. Yet, a unifying equation valid for large
deformations and applicable to materials with arbitrary complex rheology
(equivalently, arbitrary $\psi(t)$), is still lacking.

In this paper we derive a Rayleigh-Plesset-type equation for arbitrary
viscoelastic materials, i.e. materials described by an arbitrary stress
relaxation function $\psi(t)$, that is also valid for large bubble
deformations. For this, we resort to a modelling framework that combines
linear relaxation and finite strain. In the context of viscoelastic solids
this class of models is referred to as finite linear viscoelasticity [58],
while for viscoelastic liquids this corresponds to the K-BKZ model (Kaye-
Bernstein, Kearsley, Zapas) [52, 37]. Special cases of this modelling
framework include the neo-Hookean solid and the Oldroyd-B fluid. An overview
of commonly used constitutive models is provided in Tab. 1. As is known in
viscoelasticity, there is an ambiguity in the choice of upper-convected or
lower-convected tensor quantities, representing how tensors are transported by
the flow. We primarily focus on the more common upper-convected models, for
which, as we will demonstrate, the viscoelastic Rayleigh-Plesset equation for
a bubble with radius $R(t)$ takes the form:

$\begin{split}&\rho\left(R\ddot{R}+\frac{3}{2}\dot{R}^{2}\right)=\Delta
p-\frac{2\gamma}{R}\\\
&-2\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,\psi(t-t^{\prime})\,\frac{\dot{R}(t^{\prime})}{R(t)}\left\\{\left[\frac{R(t^{\prime})}{R(t)}\right]^{3}+1\right\\}.\end{split}$
(1)

Here, the dot denotes a differentiation with time, $\rho$ is the density of
the medium, $\gamma$ the interfacial tension, and $\Delta p=p_{g}-p_{\infty}$
the difference between the gas pressure $p_{g}$ and the far field pressure
$p_{\infty}$. Below, we will provide the detailed derivation of Eq. (1), as
well as some benchmarks with existing literature. Finally, as an example of
specific relevance, we will explore the resonance behavior of microbubbles
(with and without coating of phospholipid molecules) oscillating in
viscoelastic media.

Model | Schematic | Relaxation Function | Complex Modulus | Constitutive
---|---|---|---|---
|  | $\psi(t)$ | $\mu(\omega)=G^{\prime}(\omega)+iG^{\prime\prime}(\omega)$ | Differential Equation
Viscous Fluid |  | $\eta\delta(t)$ | $i\eta\omega$ | $\bm{\tau}=\eta\dot{\bm{\epsilon}}$
Maxwell Fluid |  | $Ge^{-t/\lambda}$ | $\frac{i\eta\omega}{1+i\lambda\omega}$ | $\bm{\tau}+\lambda\overset{\triangledown}{\bm{\tau}}=\eta\dot{\bm{\epsilon}}$
Oldroyd-B Fluid |  | $\eta\delta(t)+Ge^{-t/\lambda}$ | $\frac{-\eta\lambda\omega^{2}+i\omega(\eta+\eta_{p})}{1+i\lambda\omega}$ | $\bm{\tau}+\lambda\overset{\triangledown}{\bm{\tau}}=(\eta+\eta_{p})\dot{\bm{\epsilon}}+\eta\lambda\overset{\triangledown}{\dot{\bm{\epsilon}}}$
Critical Gel | - | $G(\lambda/t)^{1/2}$ | $G\left(i\pi\lambda\omega\right)^{1/2}$ | Not available
Neo-Hookean Solid |  | $G$ | $G$ | $\bm{\tau}=G\left(\mathrm{\textbf{B}}-\mathbf{I}\right)$
Kelvin-Voigt Solid |  | $G+\eta\delta(t)$ | $G+i\eta\omega$ | $\bm{\tau}=G\left(\mathrm{\textbf{B}}-\mathbf{I}\right)+\eta\dot{\bm{\epsilon}}$
Chasset-Thirion | - | $G[1+\Gamma(1-n)^{-1}\left(\frac{\lambda}{t}\right)^{n}]$ | $G\left[1+(i\lambda\omega)^{n}\right]$ | Not available
Table 1: Examples of viscoelastic fluid and solid models, defined by the
stress relaxation function $\psi(t)$ and complex modulus $\mu(\omega)$, that
are captured by Eq. (1). In some (but not all) cases the integral form of the
constitutive law (see Eq. 4 below) can be written as a differential equation
(constitutive DE), involving upper convected derivatives, the rate-of-strain
tensor $\dot{\mathbf{\epsilon}}$, and for solids involves the Finger tensor
$\mathbf{B}$. Note that the relaxation function of the Chasset-Thirion model
contains the gamma function
$\Gamma(t)=\int_{0}^{\infty}x^{t-1}e^{-x}\mathrm{d}x$.

##  II Derivation of the viscoelastic Rayleigh-Plesset equation

###  II.1 Constitutive equation

We start by considering a spherical bubble, whose radius $R(t)$ varies only in
time. We employ a spherical coordinate system
$(\hat{\mathbf{e}}_{r},\hat{\mathbf{e}}_{\theta},\hat{\mathbf{e}}_{\phi})$ at
the center of the bubble, with $r$, $\theta$, and $\phi$ denoting the radial,
azimuthal, and polar directions, respectively. Assuming a purely radial
incompressible flow, the flow velocity writes
$\mathbf{v}=(\dot{R}R^{2}/r^{2})\,\mathbf{\hat{e}}_{r}$. The momentum equation
combined with the appropriate boundary conditions can be used to obtain the
Rayleigh-Plesset equation [44]

$\begin{split}&\rho\left(R\ddot{R}+\frac{3}{2}\dot{R}^{2}\right)=\Delta
p-\frac{2\gamma}{R}+\\\
&+\int_{R}^{\infty}\mathrm{d}r\,\frac{1}{r}(2\tau_{rr}-\tau_{\theta\theta}-\tau_{\phi\phi}),\end{split}$
(2)

where $\tau_{rr}$, $\tau_{\theta\theta}$, and $\tau_{\phi\phi}$ represent the
radial, azimuthal and polar components of the deviatoric stress tensor
$\bm{\tau}$. In the context of small bubble deformations, the stress tensor is
usually assumed to be traceless, such that
$2\tau_{rr}-\tau_{\theta\theta}-\tau_{\phi\phi}=3\tau_{rr}$. This assumption,
often used to simplify the analysis, allows for an exact integration for
specific constitutive laws (e.g. Kelvin-Voigt model [61, 24]). This strategy,
however, is not suited to establish a generalized Rayleigh-Plesset equation
that encompasses arbitrary viscoelastic models for large deformations, which
is the present goal of this paper. Instead, we must maintain the integral as
in Eq. (2) and proceed with the so-called finite linear viscoelastic
formulation [58]. The essence of finite linear viscoelasticity is to combine a
linear relaxation in time, while admitting geometric nonlinearities associated
to large deformations. The corresponding constitutive relation involves an
integral over the history of deformation, and is of the form:

$\begin{split}\bm{\tau}=-\int_{-\infty}^{t}\mathrm{d}t^{\prime}&\,[\psi_{1}(t-t^{\prime})\frac{\partial\mathbf{B}(t,t^{\prime})}{\partial
t^{\prime}}\\\
&-\psi_{2}(t-t^{\prime})\frac{\partial\mathbf{B}^{-1}(t,t^{\prime})}{\partial
t^{\prime}}].\end{split}$ (3)

Here, we have introduced the Finger tensor $\mathbf{B}(t,t^{\prime})$, which
expresses the deformation between the states at time $t^{\prime}$ in the past
and the current time $t$. A precise definition of $\mathbf{B}(t,t^{\prime})$
is given below. Equation (3) resembles the K-BKZ model for viscoelastic
liquids [52, 37], which can be recovered via an integration by parts. The
appearance of two relaxations functions, respectively associated to
$\mathbf{B}$ and to its inverse $\mathbf{B}^{-1}$, reflects the freedom of
choosing upper-convected or lower-convected derivatives of tensors in
viscoelastic models (see Appendix C). We will focus on upper-convected
materials, which are based on $\mathbf{B}$ rather than its inverse [48].
Therefore, setting $\psi_{1}=\psi$ and $\psi_{2}=0$, we obtain:

$\bm{\tau}=-\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,\psi(t-t^{\prime})\frac{\partial\mathbf{B}(t,t^{\prime})}{\partial
t^{\prime}},$ (4)

which is the constitutive relation used in the remainder of this paper. For
completeness, the result obtained from using the lower convective derivative
is worked out in Appendix A.

Two remarks are in order here. First, in the limit of small deformations where
$t^{\prime}\to t$, the time derivative of the Finger tensor reduces to
$\frac{\partial\mathbf{B}(t,t^{\prime})}{\partial
t^{\prime}}=-\dot{\bm{\epsilon}}(t)$, where
$\dot{\bm{\epsilon}}(t)=\nabla\textbf{v}+(\nabla\textbf{v})^{\mathrm{T}}$ is
the rate of strain tensor (Appendix B, see also [17]). Inserting this
expression into Eq. (4), we obtain the conventional memory integral for the
stress at small deformations [7], as is frequently used in the context of
bubbles [19, 15]. In general, however,
$\frac{\partial\mathbf{B}(t,t^{\prime})}{\partial t^{\prime}}$ is not equal to
$\dot{\bm{\epsilon}}$, and this distinction is essential at large
deformations. Second, (4) reflects the deformation history of a certain
material point. By consequence, the integral must be carried out at constant
material point and calls for a Lagrangian description of the problem.

###  II.2 Lagrangian formulation

To define the Finger tensor $\textbf{B}(t,t^{\prime})$ we introduce a
Lagrangian description of the deformation, which features prominently in
finite-strain theory. Specifically, this description involves relating the
Eulerian position $\mathbf{x}=\xi(\mathbf{X},t)$ in the current configuration
at time $t$ to a Lagrangian material point $\mathbf{X}$. In solids, one
naturally defines $\mathbf{X}$ as the coordinates in the reference
configuration, but, in general, one can define $\mathbf{X}$ from the
configuration at some arbitrary time $t_{0}$. The mapping between the two
states $\xi$ can be used to evaluate the deformation gradient tensor
$\mathbf{F}(t,t_{0})=\frac{\partial\mathbf{x}}{\partial\mathbf{X}}$. This
tensor describes the material stretching through the transformation of a line
element $d\mathbf{X}$ (i.e. the distance between two material particles) in
the reference configuration at time $t_{0}$ to the same material line element
$d\mathbf{x}$ at time $t$. Using the deformation gradient, the Finger tensor
is defined as
$\mathbf{B}(t,t_{0})=\mathbf{F}(t,t_{0})\cdot\mathbf{F}^{\mathrm{T}}(t,t_{0})$.

During the purely spherical motion of the bubble, we only need to keep track
of the radial position of the materials points. In Fig. 1 we therefore denote
$\mathcal{R}$ as the reference radial position of any point in the medium. On
the interface of the bubble at rest, $\mathcal{R}=R_{0}=R(t_{0})$, which is
thus the bubble radius in the reference configuration. Since the deformation
is unidirectional, the tensor $\mathbf{F}$ is diagonal: the radial component
of the deformation gradient tensor reads $\frac{\partial
r}{\partial\mathcal{R}}$, while the azimuthal components are given by the
ratio $\frac{r}{\mathcal{R}}$ [25]. The latter represents the stretching of a
shell of constant material point. We thus find

Figure 1: Schematic illustrating the mapping between the reference and current
state for a bubble deforming with radius $R(t)$. The material coordinate $r$
in the current state can be related to the undeformed material coordinate
$\mathcal{R}$. To evaluate the stress in terms of the history of deformation,
it is instructive to introduce an intermediate past state at time
$t^{\prime}$. A new mapping can then be derived in terms of the material
coordinate in the intermediate state $r^{\prime}$. The mapping between each
state can be used to determine the deformation gradient tensor $\mathbf{F}$.

$\mathbf{F}(t,t_{0})=\begin{pmatrix}\dfrac{\partial
r}{\partial\mathcal{R}}&0&0\\\ 0&\dfrac{r}{\mathcal{R}}&0\\\
0&0&\dfrac{r}{\mathcal{R}}.\end{pmatrix}$ (5)

Incompressibility of the medium requires that $\mathrm{det}(\mathbf{F})=1$
[25], which implies $\frac{\partial
r}{\partial\mathcal{R}}=(\frac{\mathcal{R}}{r})^{2}$. Integrating this
relation, we identify the mapping between the two states:

$r^{3}(t,\mathcal{R})=\mathcal{R}^{3}+R(t)^{3}-R_{0}^{3},$ (6)

where the time-dependence is entirely encoded in the difference between the
bubble radius $R(t)$ compared to the reference radius $R_{0}$. This mapping
has indeed been used in many previous studies [19, 20, 53, 41, 1, 2, 22]. For
this spherical geometry, we simply recover the conservation of volume for
concentric spheres: $r^{3}-R^{3}=\mathcal{R}^{3}-R_{0}^{3}$. [2]

To further evaluate the memory integral in Eq. (4) we need to express the
deformation in terms of the entire history of the bubble motion, and not with
respect to the reference state. We thus introduce the position
$\mathbf{x}^{\prime}$ as the position of a material point at some past time
$t^{\prime}<t$. The deformation gradient between two arbitrary times $t$ and
$t^{\prime}$ then becomes
$\mathbf{F}(t,t^{\prime})=\frac{\partial\mathbf{x}}{\partial\mathbf{x^{\prime}}}=\mathbf{F}(t,t_{0})\cdot\mathbf{F}^{-1}(t^{\prime},t_{0})$.
As explained schematically in Fig. 1, the mapping between $t^{\prime}$ and $t$
can thus be obtained in two steps: first moving from the configuration at
$t^{\prime}$ to the reference configuration at $t_{0}$, and then going from
$t_{0}$ to $t$. Bearing in mind that
$\mathbf{B}(t,t^{\prime})=\mathbf{F}(t,t^{\prime})\cdot\mathbf{F}(t,t^{\prime})^{\mathrm{T}}$,
we thus obtain the components of the Finger tensor as

$\begin{split}B_{\phi\phi}(t,t^{\prime})=B_{\theta\theta}(t,t^{\prime})&=\left(\frac{r}{r^{\prime}}\right)^{2}\\\
&=\left(\frac{\mathcal{R}^{3}+R(t)^{3}-R_{0}^{3}}{\mathcal{R}^{3}+R(t^{\prime})^{3}-R_{0}^{3}}\right)^{2/3},\end{split}$
(7)
$\begin{split}B_{rr}(t,t^{\prime})=\frac{1}{B_{\theta\theta}(t,t^{\prime})^{2}}&=\left(\frac{r^{\prime}}{r}\right)^{4}\\\
&=\left(\frac{\mathcal{R}^{3}+R(t^{\prime})^{3}-R_{0}^{3}}{\mathcal{R}^{3}+R(t)^{3}-R_{0}^{3}}\right)^{4/3}.\end{split}$
(8)

In the final step, we have made use of the explicit radial mapping of Eq. (6).
This step is crucial, as it allows expressing $\mathbf{B}(t,t^{\prime})$ at a
constant material point $\mathcal{R}$, as is required for the evaluation of
the integral in Eq. (4). The components of the deviatoric stress tensor then
follow as

$\begin{split}\tau_{rr}=-4\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,&\\{\psi(t-t^{\prime})R^{2}(t^{\prime})\,\dot{R}(t^{\prime})\\\
&\frac{\left[\mathcal{R}^{3}+R^{3}(t^{\prime})-R_{0}^{3}\right]^{1/3}}{\left[\mathcal{R}^{3}+R^{3}(t)-R_{0}^{3}\right]^{4/3}}\\},\end{split}$
(9)
$\begin{split}\tau_{\theta\theta}=\tau_{\phi\phi}=2\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,&\\{\psi(t-t^{\prime})R^{2}(t^{\prime})\,\dot{R}(t^{\prime})\\\
&\frac{\left[\mathcal{R}^{3}+R^{3}(t)-R_{0}^{3}\right]^{2/3}}{\left[\mathcal{R}^{3}+R^{3}(t^{\prime})-R_{0}^{3}\right]^{5/3}}\\}.\end{split}$
(10)

The resulting components of $\bm{\tau}$ are now explicit functions of time,
encoded in $R(t)$, and the Lagrangian position $\mathcal{R}$.

The remaining step is to spatially integrate the stresses in the Rayleigh-
Plesset equation, as required in Eq. (2). These spatial integrals can be
carried out explicitly, so that we are only left with the temporal memory
integrals

$\begin{split}&\int_{R}^{\infty}\mathrm{d}r\,\frac{1}{r}(2\tau_{rr}-\tau_{\theta\theta}-\tau_{\phi\phi})=\\\
&-2\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,\psi(t-t^{\prime})\,\frac{\dot{R}(t^{\prime})}{R(t)}\left\\{\left[\frac{R(t^{\prime})}{R(t)}\right]^{3}+1\right\\}.\end{split}$
(11)

This concludes the derivation of the viscoelastic Rayleigh-Plesset equation
for arbitrary complex modulus as presented in Eq. (1).

##  III Applications

###  III.1 Special cases & benchmarking

We contextualize the derived Rayleigh-Plesset equation by considering a set of
special cases for $\psi(t)$, enabling a benchmark with existing literature.
The relaxation functions of various models encountered in the literature are
summarized in Table 1. Some of these can be represented schematically by a
spring-dashpot, in which case it is possible to articulate the constitutive
relation as a differential equation (see schematics in Table 1). Note that,
gel-like materials that exhibit a power-law relaxation function cannot be
represented by a differential constitutive law, and thus must be treated with
an integral formulation.

We first consider a Newtonian fluid of viscosity $\eta$, for which the
relaxation function takes the form $\psi(t)=\eta\,\delta(t)$, where
$\delta(t)$ is the Dirac delta function. Using the convolution of the delta
function, the integral in the Rayleigh-Plesset equation reduces to
$-4\eta\dot{R}/R$, which is the standard viscous contribution. Second, we turn
to the Neo-Hookean solid, which is a purely elastic medium, obtained when the
relaxation function is constant $\psi(t)=G$. In this case, we can carry out
explicitly the memory integral in Eq. (11). Applying the initial condition
$R(t^{\prime}\to-\infty)=R_{0}$, with $R_{0}$ the radius in the reference
configuration, yields an elastic contribution to the Rayleigh-Plesset equation
of the form $(G/2)[5-4(R_{0}/R)-(R_{0}/R)^{4}]$. This term exactly corresponds
to the one previously obtained for a bubble inside a Neo-Hookean solid [22,
57]. Third, we consider a bubble inside an Oldroyd-B fluid, for which the
relaxation function takes the form
$\psi(t)=\eta\,\delta(t)+\,G\mathrm{exp}(-t/\lambda)$, where $\eta$ is the
solvent viscosity, $\lambda$ the relaxation time, and $G=\eta_{p}/\lambda$ the
polymer’s elastic modulus (see also Appendix C). Here as well, the resulting
Rayleigh-Plesset equation is identical to that reported by Ting, who
considered the Oldroyd-B model in the context of bubble cavitation [53].

To further validate our model and its numerical implementation, we compare it
to the results obtained by Allen and Roy [2]. Specifically, Allen and Roy
considered large deformations of an acoustically driven micron-sized bubble.
The oscillatory pressure field is $p_{\infty}=p_{0}+p_{a}\sin(2\pi ft)$, where
$p_{0}$ is the ambient pressure, $p_{a}$ is the pressure amplitude, and $f$
the driving frequency. The gas pressure takes the form
$p_{g}=(p_{0}+2\gamma/R_{0})(R_{0}/R)^{3k}$, to account for thermal
dissipation through the polytropic constant $k$. The bubble is surrounded by
virtual Upper Convected Maxwell fluids (UCM) with relaxation times $\lambda=$
0, 0.5, and 1 $\mu$s. The ratio between the relaxation time and the
characteristic time scale of the flow is expressed by the Deborah number,
which takes the values $\mathrm{De}=2\pi f\lambda=$ 0, 1, and 2 for the three
relaxation times considered.

Our model is in excellent agreement with the results reported by Allen and
Roy, as shown in Fig. 2. Importantly, the numerical solution in [2] is not
based on the integral form of the constitutive equation, but on the
differential form. Consequently, the stress field outside the bubble must be
solved numerically, and the spatial integral of the stress needed in the
Rayleigh-Plesset equation must also be evaluated numerically. Hence, the
agreement in Fig. 2 offers a nontrivial validation of the proposed modelling
framework. We also wish to highlight that Eq. (1) has the advantage of
providing an autonomous equation for the bubble radius $R(t)$, which no longer
relies on a separate (numerical) evaluation of the stress outside the bubble.

Figure 2: Plot of the normalized radius $R/R_{0}$ against the acoustic cycles
$ft$ for a micron sized bubble in an Upper Convected Maxwell model. Our
results using the memory integral (solid lines) are overlayed on top of those
by Allen and Roy (circles) (see Fig. 12 in [2]). Here we take the pressure
amplitude as $p_{A}=0.4$ MPa and a driving frequency $f=3$ MHz. The properties
of the surrounding material are $\rho=1000$ $\mathrm{kg/m^{3}}$,
$\gamma=0.072$ N/m, and $\eta=0.03$ MPa$\cdot$s. The relaxation times are
varied such that the Deborah number $\mathrm{De}=2\pi f\lambda$ = 0, 1, 2 and
the corresponding shear modulus is computed as $G=\eta/\lambda$.

###  III.2 Resonance behavior of microbubbles

Figure 3: Response of the bubble radius $R$ when subjected to a pressure pulse
$p(t)$ with frequency $f=2$ MHz. (a) The pressure pulse $p(t)$ normalized by
its amplitude $p_{a}$ against time. (b) A bubble in a Newtonian liquid
oscillates periodically following the applied pressure pulse. The bubble
radius $R(t)$ achieves a maximum peak-to-peak amplitude $\Delta R\approx
0.4R_{0}$. If the bubble is surrounded by an Oldroyd-B liquid with relaxation
time $\lambda=1$ $\mu$s and shear modulus $G=500$ kPa, the resistance of the
surrounding fluid leads to smaller oscillation amplitudes. (c) Performing a
fast Fourier transform of the bubble oscillation amplitude $\Delta R$, we can
observe how it varies in the frequency domain, including higher and lower
harmonics.

Having established the validity of our model and its numerical implementation,
we now further explore the oscillatory motion of bubbles in viscoelastic
materials. We study the motion of a 2 $\mu$m sized bubble driven by a
Gaussian-tapered pressure pulse with 16 acoustic cycles (Fig. 3a), pressure
amplitude of 50 kPa and a frequency in the range $0.5\leq f\leq 8$ MHz. For
material properties, we select a liquid density $\rho=1,000$
$\mathrm{kg/m^{3}}$, and for most cases we fix the surface tension to a
constant value $\gamma=72$ mN/m. For viscoelastic liquids we consider the
bubble to be surrounded by an Oldroyd-B liquid or a Critical Gel, whose
relaxation functions are $\psi(t)=G\exp(-t/\lambda)$ and
$\psi(t)=G(\lambda/t)^{1/2}$ respectively (see Table 1). For the Oldroyd-B
liquid we set the solvent viscosity to $\eta=2$ $\mathrm{mPa\cdot s}$ and we
test the viscoelastic effects effects by varying the values of the shear
modulus $G$ and relaxation time $\lambda$. Specifically, we consider three
shear moduli $G=\text{10, 100, and 500}$ kPa and two relaxation times of
$\lambda=\text{0.01 and 1.00}$ $\mu$s. The same parameters are used for the
Critical Gel to enable a direct comparison. As an additional perspective, we
will also consider an example where the bubble is coated by a layer of lipids,
modelled via additional mechanical properties at the interface (see below). By
comparing the resonance behavior of coated and uncoated bubbles, one further
appreciates relative importance of the viscoelasticity of the surrounding
medium.

####  III.2.1 Oldroyd-B fluid

The results of our simulations for the evolution of the bubble radius $R(t)$
in a Newtonian and Oldroyd-B fluid are shown in Fig. 3. Considering first a
bubble in a Newtonian liquid, the oscillation amplitude of the bubble
increases initially as it follows the applied pressure pulse. The bubble
radius reaches stable oscillations with a peak-to-peak amplitude $\Delta
R=0.4R_{0}$, before decreasing to its initial radius, as the pressure pulse
decays (Fig. 3b). Switching to an Oldroyd-B fluid with relaxation time
$\lambda=1$ $\mu$s and shear modulus $G=500$ kPa, the amplitude decreases to
approximately 20% of $R_{0}$. We thus observe that the resistance of the
viscoelastic stresses can significantly affect the bubble oscillation. To
further examine the effects of the surrounding medium, it is more instructive
to evaluate the resonance behavior of the bubble. Indeed, bubbles are known to
act as damped harmonic oscillators when subjected to oscillatory pressure
driving [43, 36].

To this end, we obtain the response amplitude of the bubble through a fast
Fourier transform and plot it as a function of the driving frequency. The fast
Fourier transform of the oscillation amplitude is also indicative of the
presence of higher or lower harmonics (Fig. 3c). The resulting normalized
resonance curves of $\Delta R/R_{0}$ are shown in Fig. 4. For a 2 $\mu$m
radius bubble in a Newtonian liquid and a polytropic constant $k=1.4$, the
resonance frequency is very close to the Minnaert frequency $f_{0}\approx
1/(2\pi R_{0})(3kp_{0}/\rho)^{1/2}\approx 1.6$ MHz [36], where the bubble
displays a maximum relative amplitude of oscillation of 0.42 (Fig. 4a). The
viscoelastic effects become immediately apparent for the Oldroyd B liquid;
with a relaxation time $\lambda=0.01$ $\mu$s, the viscoelastic effects almost
exclusively translate into a reduction of the oscillation amplitude. Note that
$G=10$ kPa does not give rise to significant difference as compared to the
Newtonian case. By contrast, shear moduli of 100 and 500 kPa decrease the
response amplitude by 20 and 50%, respectively. In all three cases the
resonance frequency remains unaffected. The effects of viscoelasticity change
qualitatively when increasing the relaxation time to $\lambda=1$ $\mu$s (Fig.
4b). Even though $G=10$ kPa does not significantly change the bubble behavior,
setting $G=100$ and 500 kPa results in a drastic increase in resonance
frequency (20 and 100 %, respectively). In addition, there is a further
decrease in amplitude as compared to the case of $\lambda=0.01$ $\mu$s.

Figure 4: Effects of the shear modulus $G$ and relaxation time $\lambda$ on
the bubble amplitude $\Delta R$ of a microbubble in an Oldroyd-B liquid and a
Critical Gel. (a) For a relatively low relaxation time ($\lambda=0.01$ $\mu$s)
in an Oldroyd-B fluid, increasing the shear modulus leads to a decrease in the
bubble amplitude. Yet, the resonance frequency remains unaffected. (b)
Increasing the relaxation time to $\lambda=1$ $\mu$s now also changes the
resonance frequency of the bubble with the shear modulus. The bubble amplitude
also decreases but slightly compared to the smaller relaxation time. Inset:
The normalized storage $G^{\prime}/G$ and loss $G^{\prime\prime}/G$ moduli of
the Oldroyd-B fluid as a function of the normalized angular frequency
$\omega\lambda$, indicating a viscous behaviour at low time scales and an
elastic behavior at larger time scales. Note that the loss modulus has be
computed for a viscosity ratio $\eta/\eta_{p}=2\cdot 10^{-3}$. (c) For a
bubble in an Critical Gel with $\lambda=0.01$ $\mu$s, both and amplitude and
resonance frequency vary significantly from their Newtonian counterpart. (d)
An increase in the relaxation time has much more dramatic effects in a
Critical gel, where the amplitude sharply decreases for larger shear moduli.
Inset: Because the critical gel is a special case corresponding to the
gelation point, the storage and loss moduli are identical for all frequencies.

These results can be explained by calculating the Deborah number
$\mathrm{De}=\lambda f$, which compares the relaxation time $\lambda$ to the
characteristic time scale of the flow. For a relaxation time $\lambda=0.01$
$\mu$s, the maximum Deborah number is De = 0.1, allowing the elastic stresses
sufficient time to relax. As a result, the Oldroyd-B fluid at these low
extensional rates exhibits a viscous response rather than an elastic
contribution and $G^{\prime\prime}$ dominates (Fig. 4b inset). The effective
polymer viscosity $\eta_{p}=G\lambda$ works in conjunction with that of the
solvent $\eta$ to dampen the oscillation amplitude. For shear moduli of 10,
100, and 500 kPa, the polymer viscosity takes a value of $\eta_{p}=$ 0.1, 1,
and 5 mPa$\cdot$s, respectively. Only the latter two are high enough as
compared to the solvent viscosity $\eta$ = 2 mPa·s to affect the dynamics,
which is confirmed by our results. In contrast, for a relaxation time of
$\lambda=1$ $\mu$s the maximum Deborah number takes a value of De = 10. The
larger relaxation time does not provide sufficient time for the elastic
stresses to relax, and the surrounding medium behaves more like an elastic
solid than a viscous liquid. Indeed, for sufficiently small bubble
deformations, the natural frequency of a bubble in an elastic solid can be
computed as [22, 57, 15]

$\omega_{0}=\left\\{\frac{1}{\rho
R_{0}^{2}}\left[3kp_{0}+2(3k-1)\frac{\gamma}{R_{0}}+4G\right]\right\\}^{1/2}$
(12)

Thus, the elastic effects start affecting the resonance frequency when the
shear modulus reaches values close to the atmospheric pressure $p_{0}\approx
100$ kPa. Above this value, elasticity starts to dominate the response of the
gas. Even though Eq. (12) is only valid for small amplitudes, it qualitatively
shows how elasticity increases the resonance frequency. A fully quantitative
confirmation of this scenario is provided in Fig. 5. The figure shows that the
resonance curve for a bubble inside an Oldroyd-B fluid with $G=500$ kPa at
finite $\lambda$ switches from a Newtonian liquid behavior ($\lambda\to 0$,
blue dashed line) to that of a Kelvin-Voigt solid ($\lambda\to\infty$, orange
dashed line). Intriguingly, the transition is non-monotonic upon increasing
$\lambda$, as is clearly demonstrated in Fig. 5a.

Figure 5: Convergence of the Oldroyd-B fluid to a Kelvin-Voigt solid. (a) As
we gradually increase the relaxation time of an Oldroyd-B fluid ($G=500$ kPa)
its behaviour starts deviating from a Newtonian liquid and converges to that
of a Kelvin-Voigt solid. The amplitude at resonance follows a non-monotonic
trend with the relaxation time, highlighting the different effects of
viscosity and elasticity at different relaxation times. (b) Considering a
coated microbubble, the resonance curve qualitatively follows the same trends,
as it converges to a Kelvin-Voigt solid with an increase in the relaxation
time. Quantitatively, the additional viscous and elastic resistance of the
coating further decrease the oscillation amplitude and slightly increase the
resonance frequency.

As an additional perspective, we consider the resonance behavior of coated
microbubbles. Microbubbles can increase the contrast in ultrasound imaging or,
once embedded in tissue, used as controlled cavitation agents for therapy. By
contrast to the free gas bubble considered until now, these bubbles must be
coated with a shell to prevent gas diffusion driven by the excess Laplace
pressure. The coating has two additional contributions to the bubble dynamics
(see Appendix D). First, the dynamic change in surface area induces molecular
friction, which leads to a significantly increased damping contribution [54].
Second, the surface tension becomes a function of the packing density of the
phospholipid molecules and, therefore, for a certain initial packing fraction,
of the bubble size [34]. When the bubble size decreases, the lipids get
compressed leading to shell buckling like a spherical elastic shell. Using the
viscoelastic Rayleigh-Plesset equation for coated microbubbles (see Appendix
D), we test the resonance curve of a coated microbubble inside an Oldroyd-B
fluid with $\lambda=1$ $\mu$s and $G=500$ kPa. We take a surface elasticity
$\chi=1$ N/m and a shell dilatational viscosity of $\kappa_{s}$ = $10^{-8}$
kg/s. The resulting resonance behavior of the coated microbubble is shown in
Fig. 5b. Compared to the result for the uncoated bubbles (Fig. 5b), two
features stand out. First, the resonances shift to slightly higher
frequencies, which can be attributed to the additional stiffness induced by
surface elasticity of the coating. Second, the peaks are wider, and of smaller
amplitude, due to the additional damping induced by the shell’s viscosity.
Yet, the overall features of the resonance curves, specifically the trends
with $\lambda$, remain unaffected by the coating, showing that these are
dictated by the viscoelastic properties of the surrounding medium.

Figure 6: Resonance curve of the second harmonic response $\Delta
R_{\mathrm{2H}}$ normalized by the fundamental amplitude $\Delta R$. (a) For a
relaxation time of $\lambda=0.01$ $\mu$s, the trends of the second harmonics
follow those of the fundamental. (b) Similarly, increasing the relaxation time
to $\lambda=1$ $\mu$s we observe a shift of the resonance peak for higher
shear moduli. Yet, at approximately half the resonance frequency a new peak
appears. (c) The similar trends between fundamental and second harmonic
response can also be found for the Critical gel, where the amplitude of the
second harmonic reduces by a factor of approximately 5.

####  III.2.2 Critical Gel

Polymer networks typically exhibit a broad spectrum of relaxation times, and
therefore cannot be described by a simple exponential decay. For example,
silicone gels are well-described by the so-called Chasset-Thirion model (cf.
Table. 1), which exhibit an algebraic decay in the range $0<n<1$ [33, 28].
Here we focus on the special case of a material at the gelation point, for
which the equilibrium shear modulus vanishes. Following Winter & Chambon [59],
the Critical Gel corresponds to the case where $G^{\prime}=G^{\prime\prime}$
at all frequencies (Fig. 4d inset). This criterion is achieved only for the
special case of $\psi=G(\lambda/t)^{1/2}$, which gives
$G^{\prime}=G^{\prime\prime}\sim G\lambda^{1/2}\omega^{1/2}$. Such a material
at the gelation point has the peculiar property that, at low frequency, it has
a vanishing shear modulus ($\lim_{\omega\to 0}G^{\prime}=0$) and an infinite
steady viscosity ($\lim_{\omega\to 0}G^{\prime\prime}/\omega=\infty$).
Conversely, at large frequency, the storage modulus diverges
($\lim_{\omega\to\infty}G^{\prime}=\infty$) while the effective viscosity
vanishes ($\lim_{\omega\to\infty}G^{\prime\prime}/\omega=0$).

The resonance behavior of free gas microbubbles inside a Critical Gel is
indeed very different from that in an Oldroyd-B fluid. We find a stronger
reduction in the amplitude of oscillation. In addition, the Critical Gel also
induces a much stronger resonance frequency shift. This shift is not even
visible for the lowest relaxation time of $\lambda=0.01$ $\mu$s (Fig. 4c),
which did not induce a significant shift for the Oldroyd-B fluid. This is a
direct consequence of the absence of a characteristic timescale in the
Critical Gel. In fact, even though we separately defined $G$ and $\lambda$,
the gel is entirely characterized by the single combination $G\lambda^{1/2}$,
which is also known as the strength of the gel [59]. Therefore, the results
for a larger relaxation time of $\lambda=1$ $\mu$s (Fig. 4d), shows similar
trends as those for $\lambda=0.01$ $\mu$s. More quantitatively, we can see in
Fig. 4c and d that the combination of [$G=100$ kPa, $\lambda=0.01$ $\mu s$]
(green curve in Fig. 4c) and [$G=10$ kPa, $\lambda=1$ $\mu s$] (red curve in
Fig. 4c) result in identical curves, since the product of $G\lambda^{1/2}$ is
the equal. Thus, for the gel, in contrast to the Oldroyd-B fluid, increasing
$\lambda$ has a similar effect as increasing $G$.

The above results illustrate the importance of the rheology for the resonance
behavior of microbubbles; in particular when comparing exponentially decaying
relaxation functions to gel-like materials that exhibit a power-law spectrum.
The fact that the strength of a critical gel involves the combination
$G\lambda^{1/2}$ prevents the existence of an “elastic limit” like for the
Oldroyd-B fluid, since the equivalent solid would naturally have an infinite
stiffness. The absence of an elastic limit at $\lambda\to\infty$ is not
limited to the Critical Gel, but also applies to the Chasset-Thirion model
that is used to describe silicone gels.

####  III.2.3 Second harmonic

Microbubbles driven at sufficiently high pressure amplitudes can undergo
highly non-linear oscillations, generating higher order harmonic responses of
the bubble [42, 31], as well as subharmonic behavior [40]. In the context of
viscoelastic fluids, Allen and Roy showed that the presence of elasticity can
enhance the second harmonic response of an Upper Convected Maxwell fluid [2].
In this section we briefly comment on the second harmonic by examining the
fast Fourier transform of the bubble response (Fig. 3c). We record the
oscillation amplitude $\Delta R_{\mathrm{2H}}$ at twice the driving frequency
and directly compare the results with respect to the fundamental frequency in
the context of the Oldroyd-B and Critical Gel of Fig. 4. The resonance curve
of the second harmonic is shown in Fig. 6. For the Oldroyd-B fluid with a
relaxation time of $\lambda=0.01$ $\mu$s we observe the appearance of two
peaks for each curve (Fig. 6a). The main peak in the response occurs at
approximately half the fundamental frequency, and a smaller one at the
fundamental frequency. The higher peak at half the fundamental frequency
reflects the preference of the bubble to oscillate at its eigenfrequency. An
increase of the shear modulus dampens the curves and leads to smaller
amplitudes. The relaxation time is too short as compared to the oscillation
period, and thus excites a viscous response from the Oldroyd-B fluid. Upon
increasing the relaxation time to $\lambda=1$ $\mu$s, the resonance frequency
is shifted for the higher shear moduli (Fig. 6b). The two peaks still appear
at the fundamental and half the fundamental frequency but are shifted
depending on the elasticity. We thus find that the effects of the Oldroyd-B
fluid on the second harmonic are similar to those observed for the fundamental
frequency. The enhancement of the second harmonic response reported by Allen
and Roy could therefore be attributed to the non-monotonic behavior of the
Oldroyd-B model in the elastic limit. Finally, the second harmonic of the
Critical Gel also shares the same characteristics, as the trends follow the
behavior of the fundamental frequency (Fig. 6c). Two peaks also appear for
certain values of $G\lambda^{1/2}$, with the amplitude decreasing as the shear
modulus is increased.

##  IV Summary and outlook

In this article we have derived a unifying Rayleigh-Plesset-type equation for
bubbles in viscoelastic materials. Using concepts from finite linear
viscoelasticity, we express the viscoelastic stresses as a function of the
deformation history. This approach has the benefit of extending the Rayleigh-
Plesset equation to a broad class of viscoelastic materials with arbitrary
complex moduli through their relaxation function $\psi(t)$, while consistently
accounting for large deformations. Not only were we able to capture the bubble
dynamics in materials with well-defined constitutive relations, such as the
Upper Convected Maxwell fluid, but also to gel-like materials that can only be
characterized by their storage and loss moduli, and not via a constitutive
differential equation. With these ideas in mind, we tested the resonance
behavior of a microbubble in different types of viscoelastic materials when
driven by a tapered sinusoidal pressure pulse. The broad spectrum of
relaxation times of gel-like materials completely alter the resonance behavior
of a microbubble for sufficiently stiff shear moduli. On the other hand, the
Oldroyd-B fluid was found to behave both as a viscous liquid and elastic
solid, depending on the relaxation time. On the one hand, the Oldroyd-B fluid
acts as a viscous liquid when the relaxation time is much faster than the
typical oscillation period. Conversely, for relaxation times much slower than
the typical oscillation period, the resonance curve of an Oldroyd-B fluid
converges to that of a Kelvin-Voigt solid.

Our model is applicable to various viscoelastic materials; yet, it still has
its limitations. Even though the bubble deformation is nonlinear in terms of
the bubble radius $R$, our stress formulation of Eq. (3) can be classified as
finite linear viscoelastic, as the memory integral is still based on a
superposition principle, involving a linear operator in $\mathbf{B}$. A
possible extension to nonlinear response, while preserving the superposition
principle, is by taking the relaxation function $\psi(t)$ to be dependent on
the deformation [58]. For the derivation of Eq.(1), however, the assumption
that $\psi(t)$is independent of deformation is crucial in our approach as it
enabled us to separate the spatial and temporal integrals of the stress that
lead to Eq. (11). More generally, there are many constitute equations that are
not described by Eq. (3). Typical examples include models that exhibit finite
extensibility, which limits the extent to which the material can be stretched,
such as the Gent model for solids [23] or the FENE-P model for fluids [56].
Furthermore, viscoplastic materials also do not fall into the class of Eq. 3
[12]. The key advantage of Eq. (1) therefore lies in the applicability to
materials of arbitrary stress relaxation function, while preserving the
possibility of finite deformations.

Another limitation of our model is that we consider strictly radial bubble
motions, i.e. purely spherical volumetric oscillations. This assumption can be
consequential for sufficiently large oscillations, as the shape of a
microbubble is susceptible to parametric instabilities when driven at large
acoustic pressures [55]. The assumption of purely radial motions is also
crucial when extending our model to coated microbubbles, which are known to
exhibit non-spherical oscillations either as a result of a symmetry breaking
in the medium, or as surface modes develop [16, 15]. As a future perspective,
it would thus be worthwhile focusing on the non-spherical motion of coated and
uncoated microbubbles in viscoelastic materials. Indeed, our theoretical
approach lays a generalized framework for bubble dynamics in viscoelastic
media, which, as we have exemplified, could also be extended to model
viscoelastic dissipation at the interface of coated microbubbles.

###### Acknowledgements.

##  Appendix A Lower convected derivative

Our analysis focused on upper convected materials by setting $\psi_{2}=0$ in
Eq. (3). If we instead had chosen to work with the lower convected derivative,
it would require setting $\psi_{1}=0$ and using the inverse of
$\mathbf{B}(t,t^{\prime})$. The analysis is still straightforward, owing to
the purely radial flow that leads to a diagonal Finger tensor. The components
of the deviatoric stress tensor then adopt the form

$\begin{split}\tau_{rr}=4\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,&\\{\psi_{2}(t-t^{\prime})R^{2}(t^{\prime})\,\dot{R}(t^{\prime})\\\
&\frac{\left[\mathcal{R}^{3}+R^{3}(t)-R_{0}^{3}\right]^{4/3}}{\left[\mathcal{R}^{3}+R^{3}(t^{\prime})-R_{0}^{3}\right]^{7/3}}\\},\end{split}$
(13)
$\begin{split}\tau_{\theta\theta}=\tau_{\phi\phi}=-2\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,&\\{\psi_{2}(t-t^{\prime})R^{2}(t^{\prime})\,\dot{R}(t^{\prime})\\\
&\frac{\left[\mathcal{R}^{3}+R^{3}(t^{\prime})-R_{0}^{3}\right]^{-1/3}}{\left[\mathcal{R}^{3}+R^{3}(t)-R_{0}^{3}\right]^{2/3}}\\}.\end{split}$
(14)

The radial integration of the stresses in the Rayleigh-Plesset equation again
leaves only the temporal memory integrals such that

$\begin{split}&\int_{R}^{\infty}\mathrm{d}r\,\frac{1}{r}(2\tau_{rr}-\tau_{\theta\theta}-\tau_{\phi\phi})=\\\
&-2\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,\psi(t-t^{\prime})\,\frac{\dot{R}(t^{\prime})}{R(t)}\left[\frac{R(t)}{R(t^{\prime})}\right]^{2}\left\\{\left[\frac{R(t^{\prime})}{R(t)}\right]^{3}+1\right\\}.\end{split}$
(15)

Using the lower convected formulation leads to an additional factor
$[R(t)/R(t^{\prime})]^{2}$ in the memory integral.

##  Appendix B Kinematics

Here we provide an overview of the kinematic relations that can be used to
better interpret the stress Eq. (4) in certain limits. Recall that the
deformation gradient tensor is defined as
$\mathbf{F}(t,t^{\prime})=\frac{\partial\mathbf{x}}{\mathbf{\partial
x^{\prime}}}$. Taking the time derivative $\mathrm{d}/\mathrm{d}t$ at a
constant material point and using the chain rule, one can show that [38, 48,
17, 49]

$\displaystyle\frac{\mathrm{d}\mathbf{F}}{\mathrm{d}t}=(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}}\cdot\mathbf{F}$
(16a)
$\displaystyle\frac{\mathrm{d}\mathbf{F}^{-1}}{\mathrm{d}t}=-\mathbf{F}^{-1}\cdot(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}}$
(16b)

where $\mathbf{v}(\mathbf{x})$ is the Eulerian velocity field.

Considering now the deformations with respect to past time $t^{\prime}$,
passing via the reference state at $t_{0}$, the deformation tensor can be
expressed as
$\mathbf{F}(t,t^{\prime})=\mathbf{F}(t)\cdot\mathbf{F}^{-1}(t^{\prime})$,
where for notational convenience we use a single argument when the reference
state is involved, i.e. $\mathbf{F}(t)=\mathbf{F}(t,t_{0})$. As a result, the
Finger tensor becomes

$\begin{split}\mathbf{B}(t,t^{\prime})&=\mathbf{F}(t,t^{\prime})\cdot\mathbf{F}^{\mathrm{T}}(t,t^{\prime})\\\
&=\mathbf{F}(t)\cdot\mathbf{F}^{-1}(t^{\prime})\cdot\mathbf{F}^{-\mathrm{T}}(t^{\prime})\cdot\mathbf{F}^{\mathrm{T}}(t),\end{split}$
(17)

and its inverse

$\begin{split}\mathbf{B}^{-1}(t,t^{\prime})&=\mathbf{F}^{-\mathrm{T}}(t,t^{\prime})\cdot\mathbf{F}^{-1}(t,t^{\prime})\\\
&=\mathbf{F}^{-\mathrm{T}}(t)\cdot\mathbf{F}^{\mathrm{T}}(t^{\prime})\cdot\mathbf{F}(t^{\prime})\cdot\mathbf{F}^{-1}(t).\end{split}$
(18)

We now wish to take the time derivative of $\mathbf{B}(t,t^{\prime})$ with
respect to $t^{\prime}$ as expressed in the stress in Eq. (4). Making use of
the kinematic relations Eq. (16), we obtain [17]

$\begin{split}&\frac{\partial\mathbf{B}(t,t^{\prime})}{\partial
t^{\prime}}=\\\
&\mathbf{F}(t)\cdot\left[\frac{\mathrm{d}\mathbf{F}^{-1}(t^{\prime})}{\mathrm{d}t^{\prime}}\cdot\mathbf{F}^{-\mathrm{T}}(t^{\prime})+\mathbf{F}^{-1}(t^{\prime})\cdot\frac{\mathrm{d}\mathbf{F}^{-\mathrm{T}}(t^{\prime})}{\mathrm{d}t^{\prime}}\right]\cdot\mathbf{F}^{\mathrm{T}}(t)\\\
&=-\mathbf{F}(t)\cdot\mathbf{F}^{-1}(t^{\prime})\cdot\left[(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}}(t^{\prime})+\mathbf{\nabla}\mathbf{v}(t^{\prime})\right]\cdot\mathbf{F}^{\mathrm{-T}}(t^{\prime})\cdot\mathbf{F}^{\mathrm{T}}(t)\\\
&=-\mathbf{F}(t,t^{\prime})\cdot\dot{\bm{\epsilon}}(t^{\prime})\cdot\mathbf{F}^{\mathrm{T}}(t),\end{split}$
(19)

where we recognize the rate of strain tensor
$\dot{\bm{\epsilon}}=(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}}+\mathbf{\nabla}\mathbf{v}$.
This proves that, in general, one cannot replace
$\frac{\partial\mathbf{B}(t,t^{\prime})}{\partial t^{\prime}}$ by
$\dot{\bm{\epsilon}}$ in the constitutive relation. Similarly for the inverse
of the Finger tensor we get [17]

$\frac{\partial\mathbf{B}^{-1}(t,t^{\prime})}{\partial
t^{\prime}}=\mathbf{F}^{-\mathrm{T}}(t,t^{\prime})\cdot\dot{\bm{\epsilon}}(t^{\prime})\cdot\mathbf{F}^{-1}(t,t^{\prime}).$
(20)

Only when taking the limit $t^{\prime}\to t$, or when deformations are small,
we recover by definition that $\mathbf{F}(t,t^{\prime})\simeq\mathbf{I}$, and
we observe that the time derivatives of $\mathbf{B}(t,t^{\prime})$ and its
inverse reduce to $\mp\dot{\bm{\epsilon}}(t)$.

##  Appendix C Convected Maxwell Models

We show how, for special choices of $\psi(t)$, the integral formulation for
the stress can be expressed as a constitutive differential equations (Tab. 1).
We first introduce the upper convected derivative

$\overset{\triangledown}{\mathbf{A}}=\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}t}-(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}}\cdot\mathbf{A}-\mathbf{A}\cdot\mathbf{\nabla}\mathbf{v},$
(21)

which ensures that the dynamics of a second rank tensor $\mathbf{A}$ remain
frame invariant. The first term is a time derivative evaluated at a constant
material point, while the last two terms ensure that $\mathbf{A}$ transforms
appropriately as it gets deformed by the flow. Applying the upper convected
derivative to Eq. (4), we get

$\begin{split}\overset{\triangledown}{\bm{\tau}}&=\frac{\mathrm{d}\bm{\tau}}{\mathrm{d}t}-(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}}\cdot\bm{\tau}-\bm{\tau}\cdot\mathbf{\nabla}\mathbf{v}\\\
&=-\left.\psi(t-t^{\prime})\frac{\partial\mathbf{B}(t,t^{\prime})}{\partial
t^{\prime}}\right|_{t^{\prime}=t}-(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}}\cdot\bm{\tau}-\bm{\tau}\cdot\mathbf{\nabla}\mathbf{v}\\\
&-\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,\left[\frac{\mathrm{d}\psi(t-t^{\prime})}{\mathrm{d}t}\frac{\partial\mathbf{B}(t,t^{\prime})}{\partial
t^{\prime}}+\psi(t-t^{\prime})\frac{\partial^{2}\mathbf{B}(t,t^{\prime})}{\partial
t\partial t^{\prime}}\right]\\\
&=\psi(0)\dot{\bm{\epsilon}}(t)-\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,\frac{\mathrm{d}\psi(t-t^{\prime})}{\mathrm{d}t}\frac{\partial\mathbf{B}(t,t^{\prime})}{\partial
t^{\prime}}.\end{split}$ (22)

Note that we have utilized the relation
$\frac{\partial^{2}\mathbf{B}(t,t^{\prime})}{\partial t\partial
t^{\prime}}=(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}}(t)\cdot\frac{\partial\mathbf{B}(t,t^{\prime})}{\partial
t^{\prime}}+\frac{\partial\mathbf{B}(t,t^{\prime})}{\partial
t^{\prime}}\cdot\mathbf{\nabla}\mathbf{v}(t)$, which exactly cancels the last
two terms from the upper convected derivative.

Upon taking the upper convected derivative, we thus transformed the integral
equation to a integro-differential equation. So, not much progress is made,
except for special choices for the stress relaxation function. Specifically,
we recover the Upper Convected Maxwell model for $\psi(t)=G\exp(-t/\lambda)$.
Using this relaxation function, we obtain the differential equation

$\bm{\tau}+\lambda\overset{\triangledown}{\bm{\tau}}=\eta\dot{\bm{\epsilon}},$
(23)

where we have introduced the Maxwell viscosity $\eta=G\lambda$. Indeed, this
equation is the conventional form of the Upper Convected Maxwell model for
large deformations [7], which was also used by Allen and Roy [1] to study
large amplitude bubble oscillations in viscoelastic liquids.

For completeness, we also apply the same analysis for the Lower Convected
Maxwell Model. We first introduce the lower convected derivative

$\overset{\triangle}{\mathbf{A}}=\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}t}+(\mathbf{\nabla}\mathbf{v})\cdot\mathbf{A}+\mathbf{A}\cdot(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}},$
(24)

which has the same property of preserving frame invariance. We set
$\psi_{1}(t)=0$ in Eq. (3) and retain the $\psi_{2}(t)$ relaxation function.
Applying the lower convected derivative to the stress we get

$\begin{split}\overset{\triangle}{\bm{\tau}}&=\frac{\mathrm{d}\bm{\tau}}{\mathrm{d}t}+(\mathbf{\nabla}\mathbf{v})\cdot\bm{\tau}+\bm{\tau}\cdot(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}}\\\
&=\left.\psi_{2}(t-t^{\prime})\frac{\partial\mathbf{B}^{-1}(t,t^{\prime})}{\partial
t^{\prime}}\right|_{t^{\prime}=t}+(\mathbf{\nabla}\mathbf{v})\cdot\bm{\tau}+\bm{\tau}\cdot(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}}\\\
&+\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,\left[\frac{\mathrm{d}\psi_{2}(t-t^{\prime})}{\mathrm{d}t}\frac{\partial\mathbf{B}^{-1}(t,t^{\prime})}{\partial
t^{\prime}}+\psi_{2}(t-t^{\prime})\frac{\partial^{2}\mathbf{B}^{-1}(t,t^{\prime})}{\partial
t\partial t^{\prime}}\right]\\\
&=\psi_{2}(0)\dot{\bm{\epsilon}}(t)+\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,\frac{\mathrm{d}\psi_{2}(t-t^{\prime})}{\mathrm{d}t}\frac{\partial\mathbf{B}^{-1}(t,t^{\prime})}{\partial
t^{\prime}}.\end{split}$ (25)

Similar to the upper convected case, the relation
$\frac{\partial^{2}\mathbf{B}^{-1}(t,t^{\prime})}{\partial t\partial
t^{\prime}}=-\mathbf{\nabla}\mathbf{v}(t)\cdot\frac{\partial\mathbf{B}^{-1}(t,t^{\prime})}{\partial
t^{\prime}}-\frac{\partial\mathbf{B}^{-1}(t,t^{\prime})}{\partial
t^{\prime}}\cdot(\mathbf{\nabla}\mathbf{v})^{\mathrm{T}}(t)$ yields exact
cancellations with the last two terms of the lower convected derivative. When
we apply the relaxation function of the lower convected derivative
$\psi_{2}(t)=G\exp{(-t/\lambda)}$, we get the standard form of the Lower
Convected Maxwell model

$\bm{\tau}+\lambda\overset{\triangle}{\bm{\tau}}=\eta\dot{\bm{\epsilon}}.$
(26)

##  Appendix D Coated microbubbles

We show how the Rayleigh-Plesset equation gets modified when considering the
dynamics of coated microbubbles. Introducing a coating around the bubble
introduces an additional dilatational viscosity $\kappa_{s}$ and a surface
elasticity $\chi$. The effect of the surface elasticity is well-described by
the Marmottant model, which consists of a piece-wise function for the surface
tension $\gamma(R)$ [34]:

$\gamma(R)=\begin{cases}0\quad&\text{if}\,\,R\leq R_{\mathrm{buck}}\\\
\gamma(R_{0})+\chi\left(\frac{R^{2}}{R_{\mathrm{buck}}^{2}}-1\right)\quad&\text{if}\,\,R_{\mathrm{buck}}\leq
R\leq R_{\mathrm{rup}}\\\ \gamma_{l}\quad&\text{if}\,R\geq R_{\mathrm{rup}}\\\
\end{cases}$

Here, $R_{\mathrm{buck}}$ denotes the radius at which the coating buckles and
below which the surface tension is zero, $R_{\mathrm{rup}}$ the radius at
which the coating ruptures. Above this radius, the free gas interface is
exposed to the surrounding medium and the surface tension is equal to that of
the medium $\gamma_{l}$. The parameter $\chi$ denotes the shell elasticity of
the coating, i.e. the rate of change in surface tension with respect to the
bubble surface area $A$: $\chi=A\frac{\mathrm{d}\gamma}{\mathrm{d}A}$.
Molecular dissipation is introduced through a dilatational viscosity
$\kappa_{s}$, associated with the phospholipid monolayer. As a result, the
Rayleigh-Plesset equation now becomes

$\begin{split}&\rho\left(R\ddot{R}+\frac{3}{2}\dot{R}^{2}\right)=\left(p_{0}+\frac{\gamma(R_{0})}{R_{0}}\right)\left(\frac{R_{0}}{R}\right)^{3k}-\frac{2\gamma(R)}{R}-p_{\infty}\\\
&-\frac{4\kappa_{s}\dot{R}}{R^{2}}-2\int_{-\infty}^{t}\mathrm{d}t^{\prime}\,\psi(t-t^{\prime})\,\frac{\dot{R}(t^{\prime})}{R(t)}\left\\{\left[\frac{R(t^{\prime})}{R(t)}\right]^{3}+1\right\\}.\end{split}$
(27)

## References

  * Allen and Roy [2000a] Allen, J. S., and Roy, R. A. (2000a). “Dynamics of gas bubbles in viscoelastic fluids. i. linear viscoelasticity,” The Journal of the Acoustical Society of America 107(6), 3167–3178.
  * Allen and Roy [2000b] Allen, J. S., and Roy, R. A. (2000b). “Dynamics of gas bubbles in viscoelastic fluids. ii. nonlinear viscoelasticity,” The Journal of the Acoustical Society of America 108(4), 1640–1650.
  * Alonso [2015] Alonso, A. (2015). “Ultrasound-induced blood-brain barrier opening for drug delivery,” Translational Neurosonology 36, 106–115.
  * Andrieux _et al._ [2018] Andrieux, S., Quell, A., Stubenrauch, C., and Drenckhan, W. (2018). “Liquid foam templating–a route to tailor-made polymer foams,” Advances in Colloid and Interface Science 256, 276–290.
  * Bader _et al._ [2016] Bader, K. B., Bouchoux, G., and Holland, C. K. (2016). “Sonothrombolysis,” Ther. Ultrasound 339–362.
  * Beccaria _et al._ [2013] Beccaria, K., Canney, M., Goldwirt, L., Fernandez, C., Adam, C., Piquet, J., Autret, G., Clément, O., Lafon, C., Chapelon, J.-Y. _et al._ (2013). “Opening of the blood-brain barrier with an unfocused ultrasound device in rabbits,” J. Neurosurg. 119(4), 887–898.
  * Bird _et al._ [1987] Bird, R. B., Armstrong, R. C., and Hassager, O. (1987). _Dynamics of polymeric liquids. Vol. 1: Fluid mechanics_ (John Wiley and Sons Inc., New York, NY).
  * Bland [2016] Bland, D. R. (2016). _The theory of linear viscoelasticity_ (Courier Dover Publications).
  * Carcione _et al._ [2019] Carcione, J. M., Cavallini, F., Ba, J., Cheng, W., and Qadrouh, A. N. (2019). “On the kramers-kronig relations,” Rheologica Acta 58, 21–28.
  * Choi _et al._ [2011] Choi, J. J., Selert, K., Vlachos, F., Wong, A., and Konofagou, E. E. (2011). “Noninvasive and localized neuronal delivery using short ultrasonic pulses and microbubbles,” Proc. Natl. Acad. Sci. 108(40), 16539–16544.
  * Cunha and Albernaz [2013] Cunha, F., and Albernaz, D. (2013). “Oscillatory motion of a spherical bubble in a non-newtonian fluid,” Journal of Non-Newtonian fluid mechanics 191, 35–44.
  * De Corato _et al._ [2019] De Corato, M., Saint-Michel, B., Makrigiorgos, G., Dimakopoulos, Y., Tsamopoulos, J., and Garbin, V. (2019). “Oscillations of small bubbles and medium yielding in elastoviscoplastic fluids,” Physical Review Fluids 4(7), 073301\.
  * de Saint Victor _et al._ [2014] de Saint Victor, M., Crake, C., Coussios, C.-C., and Stride, E. (2014). “Properties, characteristics and applications of microbubbles for sonothrombolysis,” Expert Opin. Drug Deliv. 11(2), 187–209.
  * Deprez _et al._ [2021] Deprez, J., Lajoinie, G., Engelen, Y., De Smedt, S., and Lentacker, I. (2021). “Opening doors with ultrasound and microbubbles: Beating biological barriers to promote drug delivery,” Adv. Drug Deliv. Rev. 172, 9–36.
  * Dollet _et al._ [2019] Dollet, B., Marmottant, P., and Garbin, V. (2019). “Bubble dynamics in soft and biological matter,” Annu. Rev. Fluid Mech. 51, 331–355.
  * Dollet _et al._ [2008] Dollet, B., van Der Meer, S. M., Garbin, V., de Jong, N., Lohse, D., and Versluis, M. (2008). “Nonspherical oscillations of ultrasound contrast agent microbubbles,” Ultrasound in medicine & biology 34(9), 1465–1473.
  * Essink [2022] Essink, M. H. (2022). “Soft contact: from wetting to adhesion” .
  * Everitt _et al._ [2003] Everitt, S., Harlen, O., Wilson, H., and Read, D. (2003). “Bubble dynamics in viscoelastic fluids with application to reacting and non-reacting polymer foams,” Journal of Non-Newtonian Fluid Mechanics 114(2-3), 83–107.
  * Fogler and Goddard [1970] Fogler, H. S., and Goddard, J. D. (1970). “Collapse of spherical cavities in viscoelastic fluids,” Phys. Fluids 13(5), 1135–1141.
  * Fogler and Goddard [1971] Fogler, H. S., and Goddard, J. D. (1971). “Oscillations of a gas bubble in viscoelastic liquids subject to acoustic and impulsive pressure variations,” J. Appl. Phys. 42(1), 259–263.
  * Frinking _et al._ [2020] Frinking, P., Segers, T., Luan, Y., and Tranquart, F. (2020). “Three decades of ultrasound contrast agents: a review of the past, present and future improvements,” Ultrasound Med. Biol. 46(4), 892–908.
  * Gaudron _et al._ [2015] Gaudron, R., Warnez, M., and Johnsen, E. (2015). “Bubble dynamics in a viscoelastic medium with nonlinear elasticity,” Journal of Fluid Mechanics 766, 54–75.
  * Gent [1996] Gent, A. N. (1996). “A new constitutive relation for rubber,” Rubber chemistry and technology 69(1), 59–61.
  * Hamaguchi and Ando [2015] Hamaguchi, F., and Ando, K. (2015). “Linear oscillation of gas bubbles in a viscoelastic material under ultrasound irradiation,” Physics of Fluids 27(11), 113103\.
  * Holzapfel [2002] Holzapfel, G. A. (2002). “Nonlinear solid mechanics: a continuum approach for engineering science” .
  * Ichihara [2008] Ichihara, M. (2008). “Dynamics of a spherical viscoelastic shell: Implications to a criterion for fragmentation/expansion of bubbly magma,” Earth and Planetary Science Letters 265(1-2), 18–32.
  * Jiménez-Fernández and Crespo [2005] Jiménez-Fernández, J., and Crespo, A. (2005). “Bubble oscillation and inertial cavitation in viscoelastic fluids,” Ultrasonics 43(8), 643–651.
  * Karpitschka _et al._ [2015] Karpitschka, S., Das, S., van Gorcum, M., Perrin, H., Andreotti, B., and Snoeijer, J. H. (2015). “Droplets move over viscoelastic substrates by surfing a ridge,” Nat. Commun. 6(1), 7891.
  * Kelly [2013] Kelly, P. (2013). “Solid mechanics part i: An introduction to solid mechanics,” A Creative Commons Attributions, Mountain View, CA 94042.
  * Kim [1994] Kim, C. (1994). “Collapse of spherical bubbles in maxwell fluids,” Journal of non-newtonian fluid mechanics 55(1), 37–58.
  * Lauterborn [1976] Lauterborn, W. (1976). “Numerical investigation of nonlinear oscillations of gas bubbles in liquids,” J. Acoust. Soc. Am. 59(2), 283–293.
  * Lentacker _et al._ [2014] Lentacker, I., De Cock, I., Deckers, R., De Smedt, S., and Moonen, C. (2014). “Understanding ultrasound induced sonoporation: definitions and underlying mechanisms,” Adv. Drug Deliv. Rev. 72, 49–64.
  * Long _et al._ [1996] Long, D., Ajdari, A., and Leibler, L. (1996). “Static and dynamic wetting properties of thin rubber films,” Langmuir 12(21), 5221–5230.
  * Marmottant _et al._ [2005] Marmottant, P., Van Der Meer, S., Emmer, M., Versluis, M., De Jong, N., Hilgenfeldt, S., and Lohse, D. (2005). “A model for large amplitude oscillations of coated bubbles accounting for buckling and rupture,” The Journal of the Acoustical Society of America 118(6), 3499–3505.
  * Miller _et al._ [2018] Miller, D. L., Lu, X., Dou, C., Zhu, Y. I., Fuller, R., Fields, K., Fabiilli, M. L., Owens, G. E., Gordon, D., and Kripfgans, O. D. (2018). “Ultrasonic cavitation-enabled treatment for therapy of hypertrophic cardiomyopathy: Proof of principle,” Ultrasound Med. Biol. 44(7), 1439–1450.
  * Minnaert [1933] Minnaert, M. (1933). “On musical air-bubbles and the sounds of running water, london edinburgh dublin philos. mag,” Philosophical Magazine 16, 235–248.
  * Mitsoulis [2013] Mitsoulis, E. (2013). “50 years of the k-bkz constitutive relation for polymers,” International Scholarly Research Notices 2013.
  * Morozov and Spagnolie [2015] Morozov, A., and Spagnolie, S. E. (2015). “Introduction to complex fluids,” Complex Fluids in Biological Systems: Experiment, Theory, and Computation 3–52.
  * Naude and Mendez [2008] Naude, J., and Mendez, F. (2008). “Periodic and chaotic acoustic oscillations of a bubble gas immersed in an upper convective maxwell fluid,” Journal of non-newtonian fluid mechanics 155(1-2), 30–38.
  * Neppiras [1969] Neppiras, E. (1969). “Subharmonic and other low-frequency emission from bubbles in sound-irradiated liquids,” J. Acoust. Soc. Am. 46(3B), 587–601.
  * Papanastasiou _et al._ [1984] Papanastasiou, A., Scriven, L., and Macosko, C. (1984). “Bubble growth and collapse in viscoelastic liquids analyzed,” Journal of Non-Newtonian Fluid Mechanics 16(1-2), 53–75.
  * Prosperetti [1974] Prosperetti, A. (1974). “Nonlinear oscillations of gas bubbles in liquids: steady-state solutions,” J. Acoust. Soc. Am. 56(3), 878–885.
  * Prosperetti [1977] Prosperetti, A. (1977). “Thermal effects and damping mechanisms in the forced radial oscillations of gas bubbles in liquids,” The Journal of the Acoustical Society of America 61(1), 17–27.
  * Prosperetti [1982] Prosperetti, A. (1982). “A generalization of the rayleigh–plesset equation of bubble dynamics,” Phys. Fluids 25(3), 409–410.
  * Sheikov _et al._ [2008] Sheikov, N., McDannold, N., Sharma, S., and Hynynen, K. (2008). “Effect of focused ultrasound applied with an ultrasound contrast agent on the tight junctional integrity of the brain microvascular endothelium,” Ultrasound Med. Biol. 34(7), 1093–1104.
  * Shima _et al._ [1986] Shima, A., Tsujino, T., and Nanjo, H. (1986). “Nonlinear oscillations of gas bubbles in viscoelastic fluids,” Ultrasonics 24(3), 142–147.
  * Sirsi and Borden [2012] Sirsi, S. R., and Borden, M. A. (2012). “Advances in ultrasound mediated gene therapy using microbubble contrast agents,” Theranostics 2(12), 1208.
  * Snoeijer _et al._ [2020] Snoeijer, J., Pandey, A., Herrada, M., and Eggers, J. (2020). “The relationship between viscoelasticity and elasticity,” Proc. Roy. Soc. A 476(2243), 20200419.
  * Stone _et al._ [2023] Stone, H. A., Shelley, M. J., and Boyko, E. (2023). “A note about convected time derivatives for flows of complex fluids,” arXiv preprint arXiv:2304.06449 .
  * Sulheim _et al._ [2019] Sulheim, E., Mørch, Y., Snipstad, S., Borgos, S. E., Miletic, H., Bjerkvig, R., de Lange Davies, C., and Åslund, A. K. (2019). “Therapeutic effect of cabazitaxel and blood-brain barrier opening in a patient-derived glioblastoma model,” Nanotheranostics 3(1), 103\.
  * Tanasawa and Yang [1970] Tanasawa, I., and Yang, W.-J. (1970). “Dynamic behavior of a gas bubble in viscoelastic liquids,” Journal of Applied Physics 41(11), 4526–4531.
  * Tanner [1988] Tanner, R. (1988). “From a to (bk) z in constitutive relations,” Journal of Rheology 32(7), 673–702.
  * Ting [1975] Ting, R. Y. (1975). “Viscoelastic effect of polymers on single bubble dynamics,” AIChE Journal 21(4), 810–813.
  * Van der Meer _et al._ [2007] Van der Meer, S. M., Dollet, B., Voormolen, M. M., Chin, C. T., Bouakaz, A., de Jong, N., Versluis, M., and Lohse, D. (2007). “Microbubble spectroscopy of ultrasound contrast agents,” The Journal of the Acoustical Society of America 121(1), 648–656.
  * Versluis _et al._ [2010] Versluis, M., Goertz, D. E., Palanchon, P., Heitman, I. L., van Der Meer, S. M., Dollet, B., de Jong, N., and Lohse, D. (2010). “Microbubble shape oscillations excited through ultrasonic parametric driving,” Phys. Rev. E 82(2), 026321.
  * Warner Jr [1972] Warner Jr, H. R. (1972). “Kinetic theory and rheology of dilute suspensions of finitely extendible dumbbells,” Industrial & Engineering Chemistry Fundamentals 11(3), 379–387.
  * Warnez and Johnsen [2015] Warnez, M., and Johnsen, E. (2015). “Numerical modeling of bubble dynamics in viscoelastic media with relaxation,” Physics of Fluids 27(6), 063103.
  * Wineman [2009] Wineman, A. (2009). “Nonlinear viscoelastic solids—a review,” Mathematics and mechanics of solids 14(3), 300–366.
  * Winter and Chambon [1986] Winter, H. H., and Chambon, F. (1986). “Analysis of linear viscoelasticity of a crosslinking polymer at the gel point,” J. Rheol. 30(2), 367–382.
  * Xu _et al._ [2021] Xu, Z., Hall, T. L., Vlaisavljevich, E., and Lee Jr, F. T. (2021). “Histotripsy: the first noninvasive, non-ionizing, non-thermal ablation technique based on ultrasound,” nt. J. Hyperthermia 38(1), 561–575.
  * Yang and Church [2005] Yang, X., and Church, C. C. (2005). “A model for the dynamics of gas bubbles in soft tissue,” The Journal of the Acoustical Society of America 118(6), 3595–3606.