1. Introduction
Bubbles widely exist in nature, such as the ocean, volcanoes (Lyons et al. Reference Lyons, Haney, Fee, Wech and Waythomas2019) and trees (Vincent et al. Reference Vincent, Marmottant, Quinto-Su and Ohl2012). Bubbles have many important applications in food science, medicine, chemical engineering, pharmaceutical engineering, etc. For instance, collapsing bubbles generate shock waves and high-speed microjets towards a solid boundary; these properties have been used to clean food, break ice (Cui et al. Reference Cui, Zhang, Wang and Khoo2018) and break the cell membrane to deliver drugs (Brennen Reference Brennen2015). As champagne is poured into a glass, bubbles rise to the air–champagne interface and burst. The bursting bubbles eject droplets, enhancing the champagne's fragrance (Liger-Belair et al. Reference Liger-Belair, Cilindre, Gougeon, Lucio, Gebefügi, Jeandet and Schmitt-Kopplin2009). Moreover, bubbles can be used to enhance mixing in microfluidic devices (Liu et al. Reference Liu, Yang, Pindera, Athavale and Grodzinski2002; Wang et al. Reference Wang, Jiao, Huang, Yang and Nguyen2009).
The dynamics of an uncoated spherical bubble in an incompressible Newtonian medium is well described by the Rayleigh–Plesset equation (Rayleigh Reference Rayleigh1917; Plesset Reference Plesset1949; Plesset & Prosperetti Reference Plesset and Prosperetti1977), which takes into account the surface tension and the liquid viscosity. Based on the Rayleigh–Plesset equation, various models have been developed to include liquid compressibility (Keller & Miksis Reference Keller and Miksis1980; Prosperetti & Lezzi Reference Prosperetti and Lezzi1986), thermal effects (Plesset & Zwick Reference Plesset and Zwick1954), non-equilibrium evaporation and condensation (Fujikawa & Akamatsu Reference Fujikawa and Akamatsu1980; Zhong et al. Reference Zhong, Eshraghi, Vlachos, Dabiri and Ardekani2020), viscoelastic media (Warnez & Johnsen Reference Warnez and Johnsen2015; Murakami et al. Reference Murakami, Yamakawa, Zhao, Johnsen and Ando2021) and surface rheology (Chatterjee & Sarkar Reference Chatterjee and Sarkar2003; Doinikov & Dayton Reference Doinikov and Dayton2007; Dollet, Marmottant & Garbin Reference Dollet, Marmottant and Garbin2019). Cavitation in protein solutions is also frequently encountered in the pharmaceutical (Veilleux et al. Reference Veilleux, Maeda, Colonius and Shepherd2018; Zhang et al. Reference Zhang, Dou, Veilleux, Shi, Collins, Vlachos, Dabiri and Ardekani2021) and food (Narsimhan Reference Narsimhan2016) industries. However, this problem is complex for the following reasons: (1) the bulk rheological properties of protein solutions are complicated (Sharma et al. Reference Sharma, Jaishankar, Wang and McKinley2011); (2) many proteins adsorb preferentially to the solution–air interface (Fainerman, Lucassen-Reynders & Miller Reference Fainerman, Lucassen-Reynders and Miller1998); (3) the surface tension varies with the protein concentration at the interface (Miller et al. Reference Miller, Fainerman, Makievski, Krägel, Grigoriev, Kazakov and Sinyachenko2000); (4) the surface of the protein-coated bubble is usually viscoelastic (Narsimhan Reference Narsimhan2016); and (5) protein conformation, aggregation and shedding may occur at the interface, which can lead to the buckling of the interface (Lin et al. Reference Lin, Pathak, Kim, Carlson, Riguero, Kim, Buff and Fuller2016).
Current theoretical studies on bubble dynamics have covered the varying surface tension (Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005) and the rheological behaviours of the bulk medium (Warnez & Johnsen Reference Warnez and Johnsen2015). Some works (Chatterjee & Sarkar Reference Chatterjee and Sarkar2003) further treat the bubble surface as a Newtonian interface and characterize it by the Boussinesq–Scriven constitutive equation (i.e. the surface dilatational viscosity has been included for spherical bubble dynamics). The surface dilatational elasticity is also an important parameter of the protein-coated bubble surface, but it has received much less attention. Thus, we aim to derive a theoretical description of a linear viscoelastic bubble surface and investigate the importance of the surface rheology. To simplify the problem, we will consider a spherical bubble shape and neglect the protein conformation, aggregation and shedding.
We will use two cases to illustrate the importance of surface rheology. The first case in which we are interested is bubble dissolution in a protein solution. For a bubble in water, the gas concentration in the bubble is higher than that in water due to the existence of surface tension; thus, the gas in the bubble continuously flows out. Epstein & Plesset (Reference Epstein and Plesset1950) adopted a quasi-static approximation and derived a simple solution for the rate of bubble growth/dissolution, which has been experimentally validated (Kapodistrias & Dahl Reference Kapodistrias and Dahl2012). An uncoated microbubble in fresh water is expected to dissolve in seconds, but a bubble coated by lipids (Kwan & Borden Reference Kwan and Borden2012), proteins (Khan & Dalvi Reference Khan and Dalvi2020), solid particles (Poulichet & Garbin Reference Poulichet and Garbin2015) or insoluble surfactants (Hanwright et al. Reference Hanwright, Zhou, Evans and Galvin2005) can stay in the liquid for a long time (Kloek, van Vliet & Meinders Reference Kloek, van Vliet and Meinders2001). We find that, for a bubble with a viscoelastic interface, the surface excess stress can balance the surface tension and thus stabilize the bubble against dissolution.
The second case in which we are interested is the response of a protein-coated bubble to an imposed fluctuating pressure field. Including both gas diffusion and the imposed pressure field – also known as rectified mass diffusion (Crum Reference Crum1984; Church Reference Church1988; Peñas-López et al. Reference Peñas-López, Parrales, Rodríguez-Rodríguez and van der Meer2016) – can lead to interesting results, such as the progressive growth or reduction of the oscillation. Still, the problem becomes complicated, which brings difficulties in understanding the role of surface rheology. Thus, we will neglect gas diffusion in this second case. We will focus on the resonance between the bubble and an imposed acoustic pressure. In this case, theoretical solutions for resonant frequency and the peak amplitude of the bubble radius will be given.
This paper is organized as follows. A model for bubble dynamics in a protein solution is developed in § 2. The dissolution process of a protein-coated bubble is investigated in § 3. The response of a protein-coated bubble to an imposed fluctuating pressure is studied in § 4. Conclusions are provided in § 5.
2. Model
We consider a spherical bubble of radius 
$R$ in an infinite domain of an incompressible protein solution. The origin of the spherical coordinate system is at the centre of the bubble; and 
$r$, 
$\theta$ and 
$\phi$ are the radial, azimuthal and polar coordinates, respectively. The system is spherically symmetric; thus, all physical quantities, such as the protein concentration in the bulk solution, are independent of 
$\theta$ and 
$\phi$.
We let an overdot ( 
$\dot {~}$ ) denote the derivative with respect to time 
$t$; and 
$v_r$ is the radial component of the liquid velocity. We consider a protein solution with a constant density 
$\rho$ that is independent of the protein concentration. Coupling the mass conservation equation 
$\partial (r^2v_r)/\partial r = 0$ and the boundary condition 
$v_r(r=R) = \dot {R}$, we get 
$v_r = \dot {R}R^2 / r^2$. We model the protein solution as a Newtonian fluid, with its viscosity 
$\mu (c)$ depending on the local protein concentration 
$c(t,r)$. We let 
$\boldsymbol {\tau }$ denote the viscous stress tensor in the protein solution, so that 
$\tau _{rr}=-4\mu R^2 \dot {R}/r^3$ and 
$\tau _{\theta \theta }= \tau _{\phi \phi } = -\tau _{rr}/2$. Substituting 
$v_r = \dot {R}R^2 / r^2$ into the radial component of the momentum equation 
$\rho (\partial v_r / \partial t + v_r \partial v_r / \partial r) = -\partial p / \partial r + (\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {\tau })_r$, we have
\begin{equation} \rho\left(\frac{R^2\ddot{R}+2R\dot{R}^2}{r^2} - \frac{2R^4\dot{R}^2}{r^5} \right) ={-}\frac{\partial p}{\partial r} + (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau})_r, \end{equation}
where 
$p$ is the pressure in the protein solution. We integrate the above equation from 
$r=R$ to 
$r \rightarrow \infty$, to obtain
\begin{equation} \rho\left(R \ddot{R} + \frac{3}{2}\dot{R}^2 \right) = p|_{r=R} - p_{\infty} + \int_R^\infty \frac{3\tau_{rr}}{r}\,\mathrm{d}r - \tau_{rr}|_{r=R}, \end{equation}
where 
$p_\infty$ is the pressure in the protein solution far from the bubble.
The balance of the normal stress across the bubble surface is
\begin{equation} p|_{r=R} - (\boldsymbol{\nabla}_s \boldsymbol{\cdot} {\boldsymbol{\mathsf{{P}}}}_s) \boldsymbol{\cdot} \boldsymbol{{n}} - \tau_{rr}|_{r=R} = p_b. \end{equation}
Here 
$\boldsymbol {{n}}$ is the unit normal vector to the bubble surface; 
$p_b = p_v + p_g$ is the pressure inside the bubble, where 
$p_v$ and 
$p_g$ are the vapour pressure and gas pressure, respectively; 
$\boldsymbol {\nabla }_s \equiv {\boldsymbol{\mathsf{{I}}}}_s \boldsymbol {\cdot } \boldsymbol {\nabla }$ is the surface gradient operator; 
${\boldsymbol{\mathsf{{I}}}}_s \equiv {\boldsymbol{\mathsf{{I}}}} - \boldsymbol {{n}}\boldsymbol {{n}}$ is the surface identity matrix, where 
${\boldsymbol{\mathsf{{I}}}}$ is the identity matrix of size three; and 
${\boldsymbol{\mathsf{{P}}}}_s = {\boldsymbol{\mathsf{{I}}}}_s\gamma + \boldsymbol {\tau }_s$ is the surface excess pressure tensor, where 
$\gamma$ and 
$\boldsymbol {\tau }_s$ are the surface tension and surface excess stress tensor, respectively.
The constitutive equation for a linear viscoelastic interface (Edwards, Brenner & Wasan Reference Edwards, Brenner and Wasan1991; Narsimhan Reference Narsimhan2016) is
\begin{equation} \boldsymbol{\tau}_s = \int_{-\infty}^t \{ G^d(t-t') - G^s(t-t') \} ({\boldsymbol{\mathsf{{I}}}}_s \boldsymbol{:} {\boldsymbol{\mathsf{{D}}}}_s(t') ) {\boldsymbol{\mathsf{{I}}}}_s \,\mathrm{d}t' + 2\int_{-\infty}^t G^s(t-t'){\boldsymbol{\mathsf{{D}}}}_s(t')\,\mathrm{d}t', \end{equation}
where 
$G^d$ and 
$G^s$ are the dilatational modulus and the shear modulus, respectively. The surface rate-of-deformation tensor 
${\boldsymbol{\mathsf{{D}}}}_s$ is
\begin{equation} {\boldsymbol{\mathsf{{D}}}}_s \equiv \frac{1}{2} [(\boldsymbol{\nabla}_s \boldsymbol{v}_s) \boldsymbol{\cdot} {\boldsymbol{\mathsf{{I}}}}_s +{\boldsymbol{\mathsf{{I}}}}_s \boldsymbol{\cdot} (\boldsymbol{\nabla}_s \boldsymbol{v}_s)^{\rm T} ] = \frac{\dot{R}}{r}{\boldsymbol{\mathsf{{I}}}}_s, \end{equation}
where 
$\boldsymbol {v}_s$ is the surface velocity. Substituting (2.5) into (2.4) and coupling equations (2.2) and (2.3), we have
\begin{align} \rho\left(R \ddot{R} + \frac{3}{2}\dot{R}^2 \right) &= p_v + p_g - p_{\infty} -12R^2 \dot{R} \int_R^\infty \frac{\mu(c) }{r^4}\,\mathrm{d}r - \frac{2\gamma}{R} \nonumber\\ & \quad - \frac{4}{R}\int_{-\infty}^t G^d (t-t') \frac{\dot{R}(t')}{R(t')}\,\mathrm{d}t'. \end{align}We consider a linear viscoelastic interface, which is characterized by the Maxwell model
\begin{equation} \frac{\partial}{\partial t}\left( \frac{\boldsymbol{\tau}_s}{E_s}\right) + \frac{\boldsymbol{\tau}_s}{\kappa_s} = {\boldsymbol{\mathsf{{D}}}}_s, \end{equation}
where 
$E_s$ and 
$\kappa _s$ are the surface dilatational elasticity and the surface dilatational viscosity, respectively. Moreover, we assume 
$E_s$ and 
$\kappa _s$ to be proportional to the protein concentration at the interface 
$\varGamma$ (Narsimhan Reference Narsimhan2016). Defining 
$E_s = E_{sl}\varGamma$ and 
$\kappa _s = \kappa _{sl}\varGamma$ and substituting them into (2.7), we have
\begin{equation} \boldsymbol{\tau}_s = \int_{-\infty}^t \{E_s(t)\exp[-(t-t')/\lambda_{s}]\}{\boldsymbol{\mathsf{{D}}}}_s(t')\,\mathrm{d}t', \end{equation}
where 
$\lambda _s = \kappa _s / E_s = \kappa _{sl}/E_{sl}$ is the (constant) relaxation time.
Thus, (2.6) for a linear viscoelastic bubble surface (characterized by the Maxwell model) can be rewritten as
\begin{align} \rho\left(R \ddot{R} + \frac{3}{2}\dot{R}^2 \right) &= p_v + p_g - p_{\infty} -12R^2 \dot{R} \int_R^\infty \frac{\mu(c) }{r^4}\,\mathrm{d}r - \frac{2\gamma}{R} \nonumber\\ & \quad - \frac{4E_s}{R}\int_{0}^t \exp[-(t-t')/\lambda_{s}] \frac{\dot{R}(t')}{R(t')} \,\mathrm{d}t', \end{align}or
\begin{equation} \left. \begin{aligned} & \rho\left(R \ddot{R} + \frac{3}{2}\dot{R}^2 \right) = p_v + p_g - p_{\infty} -12R^2 \dot{R} \int_R^\infty \frac{\mu(c) }{r^4}\,\mathrm{d}r - \frac{2\gamma}{R} - \frac{4E_s \mathcal{F}}{R},\\ & \dot{\mathcal{F}} + \frac{\mathcal{F}}{\lambda_s}=\frac{\dot{R}}{R}, \quad \mathcal{F}|_{t=0} = 0, \end{aligned} \right\} \end{equation}
where 
$\mathcal {F}=\int _0^t \exp [-(t-t')/\lambda _s]\dot {R}(t')/R(t')\,\mathrm {d}t'$.
If we consider a constant liquid viscosity 
$\mu$ and a Newtonian bubble surface, i.e. 
$\lambda _s \rightarrow 0$, then (2.9) reduces to the model for surfactant- or lipid-coated bubbles (Chatterjee & Sarkar Reference Chatterjee and Sarkar2003; Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005):
\begin{equation} \rho\left(R \ddot{R} + \frac{3}{2}\dot{R}^2 \right) = p_v + p_g - p_{\infty} -4\mu \frac{\dot{R}}{R} - \frac{2\gamma}{R} - 4\kappa_s \frac{\dot{R}}{R^2}. \end{equation} If we consider a constant liquid viscosity 
$\mu$ and an elastic bubble surface, i.e. 
$\lambda _s \rightarrow \infty$, then (2.9) reduces to
\begin{equation} \rho\left(R \ddot{R} + \frac{3}{2}\dot{R}^2 \right) = p_v + p_g - p_{\infty} -4\mu \frac{\dot{R}}{R} - \frac{2\gamma}{R} - \frac{4E_s}{R} \ln\left(\frac{R}{R_0} \right), \end{equation}
where 
$R_0$ is the bubble radius at 
$t=0$.
The protein concentration in the bulk of the liquid, 
$c$ with a unit 
$\mathrm {kg}~\mathrm {m}^{-3}$, is governed by the advection–diffusion equation
\begin{equation} \frac{\partial c}{\partial t} + v_r \frac{\partial c}{\partial r} = \frac{1}{r^2}\frac{\partial}{\partial r}\left(D_c r^2 \frac{\partial c}{\partial r} \right), \end{equation}which is subject to
\begin{equation} c|_{t=0} = c_\infty, \quad c|_{r\rightarrow \infty} = c_\infty, \quad \left. D_c\frac{\partial c}{\partial r}\right|_{r=R} = J, \end{equation}
where 
$D_c$ is the diffusion coefficient of protein in the bulk solution and 
$c_\infty$ is the protein concentration far from the bubble surface. There are several well-accepted models (Henry, Langmuir, Frumkin, van der Waals, etc.) to describe the adsorption isotherms of surfactants and proteins. The corresponding expressions for the surface tension and adsorption/desorption flux can be found in the literature (Edwards et al. Reference Edwards, Brenner and Wasan1991; Manikantan & Squires Reference Manikantan and Squires2020).
We model the net flux of protein to the interface 
$J$ as
\begin{equation} J = \left\{ \begin{array}{@{}ll} \theta^a c_s(\varGamma_m - \varGamma) - \theta^d \varGamma, & \text{when }\varGamma \leq \varGamma_m ,\\ - \theta^d \varGamma, & \text{when }\varGamma > \varGamma_m. \end{array} \right. \end{equation}
Here 
$\theta ^a$ and 
$\theta ^d$ are the adsorption and desorption constants, respectively; 
$c_s = c|_{r=R}$ represents the protein concentration in the bulk solution close to the bubble surface; and 
$\varGamma _m$ is the maximum packing density of protein at the bubble surface. If 
$\varGamma \leq \varGamma _m$, the first term in 
$J$, which is described by the Langmuir equation (Craster, Matar & Papageorgiou Reference Craster, Matar and Papageorgiou2009; Langevin Reference Langevin2014; Manikantan & Squires Reference Manikantan and Squires2020), denotes the influx that is proportional to the protein concentration near the interface and the available space on the interface (proportional to (
$\varGamma _m - \varGamma$)). The second term in 
$J$ denotes the desorption of protein from the interface. If 
$\varGamma > \varGamma _m$, we assume the influx to be zero such that 
$J = - \theta ^d \varGamma$.
The protein concentration at the bubble surface, 
$\varGamma$ with a unit 
$\mathrm {kg}~\mathrm {m}^{-2}$, is governed by (Stone Reference Stone1990)
\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} (\varGamma R^2) = JR^2. \end{equation}We use the modified Langmuir relation (Fainerman et al. Reference Fainerman, Lucassen-Reynders and Miller1998) to characterize the surface tension:
\begin{equation} \gamma = \left\{ \begin{array}{@{}ll} \gamma_f + R_g T \varGamma_{p}\ln \left( 1-\dfrac{\varGamma}{\varGamma_m}\right), & \text{when }\varGamma \leq \varGamma_c, \\ \gamma_c, & \text{when }\varGamma > \varGamma_c. \end{array} \right. \end{equation}
Here 
$R_g$, 
$T$, 
$\gamma _f$ and 
$\gamma _c$ are the gas constant, the temperature in the bulk solution near the interface, the surface tension of a clean solution–air interface, and the minimum surface tension, respectively; 
$\varGamma _{p}$ is a fitting parameter, with unit 
$\mathrm {mol}~\mathrm {m}^{-2}$; and 
$\varGamma _c = \varGamma _m - \varGamma _m \exp [(\gamma _c - \gamma _f) / (R_g T \varGamma _p)]$, beyond which the surface tension remains constant. Here 
$E_{\varGamma } \equiv - \varGamma \partial \gamma / \partial \varGamma$ is often called the Gibbs, Marangoni or simply dilatational modulus (Manikantan & Squires Reference Manikantan and Squires2020); and 
$E_\varGamma$ describes the rate of change of the surface tension with respect to the protein concentration on the bubble, while 
$E_s$ and 
$\kappa _s$ define the excess stress induced by the viscoelasticity of the interface. Thus, 
$E_\varGamma$ and 
$E_s$ (or 
$\kappa _s$) are two different parameters (Manikantan & Squires Reference Manikantan and Squires2020). We use Mooney's equation (Mooney Reference Mooney1951; Tomar et al. Reference Tomar, Kumar, Singh, Goswami and Li2016) to characterize the bulk viscosity,
\begin{equation} \mu = \mu_f \exp\left(\frac{2.5\phi}{1-k \phi}\right), \end{equation}
where 
$\mu _f$, 
$\phi =c/\rho$ and 
$k$ are the liquid viscosity at zero protein concentration, the volume fraction of protein in the bulk solution and the crowding factor (
$1.35 < k < 1.91$), respectively.
To close the model, we need an additional equation for the gas pressure inside the bubble. In the first case study in which the bubble dissolves in a protein solution due to gas diffusion, we assume the temperature in the bubble 
$T_g$ to remain constant and equal to the temperature in the bulk of the liquid, i.e. 
$T_g = T$. We let 
$\epsilon$, 
$\psi$, 
$V$, 
$n_g$, 
$M_g$, 
$m_g$ and 
$D_\psi$ denote the Henry's constant with a unit 
$\mathrm {kg}~\mathrm {m}^{-3}~\mathrm {Pa}^{-1}$, the gas concentration in the bulk of the liquid with a unit 
$\mathrm {kg}~\mathrm {m}^{-3}$, the volume of the bubble, the amount of gas in the bubble with a unit 
$\mathrm {mol}$, the molar weight of the gas, the mass of gas in the bubble, and the diffusion coefficient of the gas in the bulk of the liquid, respectively. The gas concentration in the liquid 
$\psi$ is governed by (Epstein & Plesset Reference Epstein and Plesset1950; Peñas-López et al. Reference Peñas-López, Parrales, Rodríguez-Rodríguez and van der Meer2016)
\begin{equation} \frac{\partial \psi}{\partial t} + v_r \frac{\partial \psi}{\partial r} = \frac{1}{r^2}\frac{\partial}{\partial r}\left(D_\psi r^2 \frac{\partial \psi}{\partial r} \right), \end{equation}which is subject to
\begin{equation} \psi|_{t=0} = \psi_\infty, \quad \psi|_{r\rightarrow \infty} = \psi_\infty, \quad \psi|_{t>0, r=R} = \epsilon p_g. \end{equation}
The far-field gas concentration is 
$\psi _\infty = \epsilon (p_a - p_v)$, where 
$p_a$ is the atmospheric pressure. The rate of change of the gas mass in the bubble is
\begin{equation} \frac{\mathrm{d}m_g}{\mathrm{d}t} = \left. 4{\rm \pi} R^2 D_\psi \frac{\partial \psi}{\partial r}\right|_{r=R}. \end{equation}
From the ideal gas law 
$p_g V = n_g R_g T_g$, we have
\begin{equation} \frac{\mathrm{d}m_g}{\mathrm{d}t} = \frac{\mathrm{d}(n_gM_g)}{\mathrm{d}t} = \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{4{\rm \pi} M_gR^3p_g}{3R_gT_g} \right). \end{equation}Coupling equations (2.21) and (2.22), we get
\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} (p_g R^3) = \left. 3 D_\psi R^2\frac{R_g T_g}{M_g}\frac{\partial \psi}{\partial r}\right|_{r=R}. \end{equation}In the second case study in which an imposed fluctuating pressure drives the bubble to oscillate, we assume the gas pressure to follow a polytropic process,
\begin{equation} p_g = p_{g0}\left(\frac{R_0}{R}\right)^{3\alpha}, \end{equation}
where 
$p_{g0}$ and 
$\alpha$ are the gas pressure at 
$t = 0$ and the polytropic index, respectively; 
$\alpha = 1$ represents an isothermal process, and 
$\alpha = 1.4$ denotes an adiabatic process.
In the following two case studies, unless otherwise stated, we will use the values of the physical parameters listed in table 1. Bovine serum albumin (BSA) solution is considered here. Values for the varying surface tension of the solution–air interface and dynamic adsorption/desorption of protein are obtained by fitting to the experimental data from the literature (Fainerman et al. Reference Fainerman, Lucassen-Reynders and Miller1998; Ybert & di Meglio Reference Ybert and di Meglio1998; Miller et al. Reference Miller, Fainerman, Makievski, Krägel, Grigoriev, Kazakov and Sinyachenko2000). At 
$t = 0$, we assume that: (1) the velocity of the bubble surface 
$\dot {R}(t=0)=0$; (2) unless otherwise stated, the adsorption/desorption of protein is in equilibrium, i.e. 
$J=0$, and thus the initial protein concentration at the bubble surface 
$\varGamma _0 = \theta ^a c_\infty \varGamma _m / (\theta ^a c_\infty + \theta ^d)$; and (3) the initial gas pressure 
$p_{g0} = p_a - p_v + 2\gamma _0/R_0$, where 
$\gamma _0$ is the initial surface tension. We introduce 
$\zeta = 1 - 2R/r$ such that 
$r\in [R, \infty )$ maps to 
$\zeta \in [-1, 1)$. We let 
$\zeta = 0.98$, corresponding to 
$r = 100R$, represent the far field, and then we divide 
$\zeta \in [-1, 0.98]$ into 300 uniform intervals. We use the backward Euler method to solve (2.10), (2.13), (2.16), and (2.23) or (2.24).
Table 1. Physical parameters and their values.
3. Bubble dissolution in a protein solution
3.1. Prediction of the dissolution rate
The dissolution of a bubble is caused by the imbalanced gas concentration across the bubble surface. Initially, the surface excess stress is negligible because 
$\mathcal {F}(t=0) = 0$. According to the Epstein–Plesset theory (Epstein & Plesset Reference Epstein and Plesset1950), the surface velocity is
\begin{equation} \frac{\mathrm{d}R}{\mathrm{d}t} ={-}\frac{ R_g T D_\psi\epsilon}{p_g M_g} \, \frac{2\gamma}{R} \left(\frac{1}{R} + \frac{1}{\sqrt{{\rm \pi} D_\psi t}}\right). \end{equation}
The above equation indicates that the dissolution rate increases with the gas diffusion coefficient but decreases with the bubble radius. Protein adsorption and desorption are in equilibrium at 
$t=0$. If the protein concentration in the bulk solution 
$c_\infty$ is low, increasing 
$c_\infty$ consequently yields a higher protein concentration at the bubble surface, and hence a smaller surface tension. Thus at small 
$t$, the dissolution rate decreases with the protein concentration in the bulk of the liquid, as shown in figure 1. If we define 
$R = R_0 \tilde {R}$ and 
$t = \lambda \tilde {t}$, where 
$\lambda = R_0^2 / D_\psi$ is the gas diffusion time scale, then (3.1) can be rewritten as
\begin{equation} \frac{\mathrm{d}\tilde{R}}{\mathrm{d}\tilde{t}} ={-}\frac{1}{{We}\,\tilde{R}}\left(\frac{1}{\tilde{R}} + \frac{1}{\sqrt{{\rm \pi} \tilde{t} }} \right), \quad {We} = \frac{\rho_g / \epsilon}{2\gamma/R_0}, \end{equation}
where 
$\rho _g$ is the instantaneous gas density. The Weber number, 
${We}$, compares the gas inertial force to the surface tension force. Thus, the surface tension drives the gas to leave the bubble (through diffusion), and the dissolution rate is initially determined by the ratio of the surface tension force to the gas inertial force. In other words, the dissolution of the bubble initially follows an inertio-capillary process. The dissolution rate decreases with the gas inertia, so, in medicine, where microbubbles are used to open cell membranes and deliver drugs, bubbles with a heavy gas core (e.g. 
$\mathrm {C}_3 \mathrm {F}_8$ with density 8.17 
$\mathrm {kg}~\mathrm {m}^{-3}$ or 
$\mathrm {C}_4 \mathrm {F}_{10}$ with density 11.2 
$\mathrm {kg}~\mathrm {m}^{-3}$) are generally used to extend the life of bubbles in the blood (Dauba et al. Reference Dauba, Delalande, Kamimura, Conti, Larrat, Tsapis and Novell2020).
Figure 1. Effect of protein concentration on the bubble dissolution process.
We define 
$Z_1$ as the ratio of surface excess stress to the surface tension and 
$Z_2$ as the ratio of bulk viscous stress to the surface tension, i.e.
$$\begin{gather} Z_1 =\left. \left|\frac{4E_s\mathcal{F}}{R} \right| \right/ \left(\frac{2\gamma}{R} \right) =\frac{2}{{Ec}}\int_{0}^t \exp[-(t-t')/\lambda_{s}] \frac{\dot{R}(t')}{R(t')} \,\mathrm{d}t', \end{gather}$$
$$\begin{gather}Z_2 =\left. \left|12R^2 \dot{R} \int_R^\infty \frac{\mu(c) }{r^4}\,\mathrm{d}r \right| \right/ \left(\frac{2\gamma}{R} \right)=\frac{3{Ca}}{4}\int_{{-}1}^{1}\exp\left(\frac{2.5\phi}{1-k \phi}\right)(1-\zeta)^2 \,\mathrm{d}\zeta . \end{gather}$$
Here 
${Ec} = \gamma /E_s$ is the elasto-capillary number, which compares the surface tension and surface dilatational elasticity; and 
${Ca}=\mu _f \dot {R}/\gamma$ is the capillary number, which measures the importance of the viscous stress in the bulk of the liquid compared to the surface tension.
As the bubble shrinks, the protein density at the bubble surface increases, which leads to the growth of the surface excess stress, as shown in figure 2(a). The rate of dissolution is determined by the gas inertial force, surface tension force and surface excess stress, during which the dissolution of the bubble follows an inertio-elasto-capillary process. As the bubble progressively shrinks, the surface excess stress eventually balances the surface tension force, and the dissolution of the bubble follows an elasto-capillary process.
Figure 2. (a) Time evolution of the total protein on the bubble surface 
$\varGamma R^2$, the ratio of surface excess stress to the surface tension 
$Z_1$, and the ratio of bulk viscous stress to the surface tension 
$Z_2$ in the case of 
$c_\infty = 10^{-4}~\mathrm {kg}~\mathrm {m}^{-3}$. (b) Comparison of the bubble behaviour in water and in a protein solution with 
$c_\infty = 10^{-4}~\mathrm {kg}~\mathrm {m}^{-3}$.
For a Newtonian interface (i.e. the surface relaxation time goes to zero), the balance of surface tension force and surface dilatational viscous stress yields 
$\dot {R} = -\gamma R / (2\kappa _s)$. Thus, the dissolution rate for a bubble with a Newtonian interface (e.g. a surfactant- or lipid-coated bubble) is inversely proportional to the surface dilatational viscosity.
For an elastic interface (i.e. the surface relaxation time goes to infinity), the balance of surface tension force and surface dilatational elastic stress yields 
$R/R_0 = \exp [-\gamma /(2E_s)]=\exp [-{Ec}/2]$. Thus the dimensionless bubble radius is solely determined by the elasto-capillary number. A value 
$Z_2 \ll 1$ also indicates that the viscous stress in the bulk of the liquid is negligible. Additionally, after the surface excess stress balances the surface tension force, the bubble dissolution rate is small such that the gas pressure inside the bubble 
$p_g$ remains almost constant (close to atmospheric pressure), as shown in figure 2(b).
Next, we aim to obtain the dissolution rate of a protein-coated bubble after the surface excess stress balances the surface tension, i.e. 
$Z_1 = 1$. Let us consider a protein solution with high concentration that satisfies 
$\theta ^a c_\infty \gg \theta ^d$; then 
$\varGamma _0 = \theta ^a c_\infty \varGamma _m / (\theta ^a c_\infty + \theta ^d) \approx \varGamma _m$, which is the maximum packing density of proteins at the bubble surface. Proteins typically have a low desorption ability. Thus, the protein concentration at the bubble surface increases as the bubble shrinks, i.e. 
$\varGamma (t) \geq \varGamma _m$, which yields 
$\gamma = \gamma _c$ (i.e. the surface tension remains constant as the bubble shrinks) and 
$J = -\theta ^d \varGamma$. Using (2.16), we obtain the total protein at the bubble surface 
$\varGamma R^2 = \varGamma _0 R_0^2 \exp (-\theta ^d t)$, which decreases exponentially with time, and the decay rate equals the protein desorption constant. Substituting the expression of 
$\varGamma$ into 
$Z_1 = 1$, i.e. the balance of surface tension and surface excess stress, we get
\begin{equation} \mathcal{F} ={-}\frac{\gamma_c}{2E_s} ={-}\frac{\gamma_c}{2E_{sl}\varGamma} ={-}\frac{\gamma_c R^2 \exp(\theta^d t)}{2E_{sl}\varGamma_0 R_0^2}. \end{equation}
Substituting (3.5) into 
$\dot {\mathcal {F}} + \mathcal {F}/\lambda _s=\dot {R}/R$ and after some manipulations, we obtain
\begin{equation} \frac{\mathrm{d}R}{\mathrm{d}t} ={-} \frac{R_0}{2\lambda_s} \, \frac{(\theta^d \lambda_s + 1)(R/R_0)^3} {(R/R_0)^2 + \exp(-\theta^d t)/{Ec}_0}, \end{equation}
where 
${Ec}_0 = \gamma _c / (E_{sl}\varGamma _0) = \gamma _c/E_{s0}$ denotes the ratio of the surface tension to the initial surface dilatational elasticity, 
$E_{s0}$.
The dissolution rate will not reduce to zero (i.e. the bubble will not remain static) because (1) the surface excess stress relaxes (decreases) when the strain on the surface remains constant and (2) protein desorbs from the bubble surface, which reduces the surface excess stress. If 
$\theta ^d \lambda _s \ll 1$, the variation of the surface excess stress with respect to time is mainly caused by the surface relaxation, and the dissolution rate is inversely proportional to the surface relaxation time 
$\lambda _s$ at large 
$t$, as indicated by (3.6). If 
$\theta ^d \lambda _s \gg 1$, the variation of the surface excess stress is induced by the protein desorption from the bubble surface, and hence the dissolution rate is independent of the surface relaxation time 
$\lambda _s$. Verification of the Epstein–Plesset theory (short-time asymptotic), i.e. (3.1), and the long-time asymptotic, i.e. (3.6), are shown in figure 3.
Figure 3. The time evolution of the dissolution rate.
During the elasto-capillary process, let us use the surface relaxation time 
$\lambda _s$ as the time scale and define 
$t = \lambda _s \check {t}$. Then (3.6) can be rewritten as
\begin{equation} \frac{\mathrm{d}\tilde{R}}{\mathrm{d}\check{t}} ={-}\frac{(\theta^d \lambda_s + 1 )\tilde{R}^3}{2(\exp(-\theta^d\lambda_s\check{t})/{Ec}_0 + \tilde{R}^2)}. \end{equation}
In the limiting case of 
$\theta ^d \lambda _s \check {t} \ll 1$ or 
$\theta ^d \lambda _s \check {t} \gg 1$, (3.7) reduces to
\begin{equation} \frac{\mathrm{d}\tilde{R}}{\mathrm{d}\check{t}} = \left\{ \begin{array}{@{}ll} -(\theta^d \lambda_s + 1)\tilde{R}^3/[2({Ec}_0^{{-}1} + \tilde{R}^2)], & \text{when}\ \theta^d \lambda_s \check{t} \ll 1,\\ -(\theta^d \lambda_s + 1)\tilde{R}/2, & \text{when}\ \theta^d \lambda_s \check{t} \gg 1. \end{array} \right. \end{equation}
The continuous decrease of 
$\tilde {R}$ yields a progressive reduction of the dissolution rate on both short and long time scales (
$\theta ^d \lambda _s \check {t}$ is much less or much greater than one). As the bubble shrinks, on the one hand, the reduction of the surface area increases the protein density on the bubble surface (i.e. decelerates dissolution); on the other hand, the protein desorbs from the bubble, which reduces the protein density on the bubble surface (i.e. accelerates dissolution). Thus, the bubble may accelerate or decelerate shrinking at intermediate scale (
$\theta ^d \lambda _s \check {t}$ of the order of one), depending on the magnitude of protein desorption ability and surface dilatational elasticity, as indicated by (3.7). Additionally, the surface relaxation time and surface dilatational elasticity help stabilize the bubble against dissolution, as shown in figure 4.
Figure 4. (a) Effects of surface relaxation time on the bubble dynamics; 
$E_{sl}=10^4~\mathrm {m}^2~\mathrm {s}^{-2}$ is used here. (b) Effects of surface dilatational elasticity on the bubble dynamics; 
$\lambda _s = 10~\mathrm {s}$ is used here.
Khan & Dalvi (Reference Khan and Dalvi2020) performed experiments to study the dissolution behaviour of microbubbles in aqueous BSA solution (
$3\,\%\,\mathrm {w/v}$). They observed that: (1) the protein-coated bubble shrinks at a slow rate (compared to the uncoated bubble) before the bubble buckles; and (2) once the bubble buckles, a large portion of the protein detaches from the bubble surface, and the microbubble quickly disappears. This is different from the behaviour of lipid-coated bubbles, which stay for a while after buckling because the surface tension also reduces to zero (Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005). Our model well captures the kinetics of the BSA-coated bubble before it buckles, as shown in figure 5. Additionally, a larger BSA-coated bubble has a longer relaxation time according to our calibration results. An uncoated microbubble can only stay in water for approximately 1 s, but the time for a protein-coated microbubble to disappear after buckling is of the order of 10 s (Khan & Dalvi Reference Khan and Dalvi2020). Two reasons may explain this: (1) the bubble is non-spherical after buckling; and/or (2) there is a small amount of protein left on the bubble surface; thus, the protein structure provides resistance.
Figure 5. Comparison between our model predictions and experimental data for the dissolution of microbubbles in aqueous BSA solution (
$3\,\%\,\mathrm {w/v}$) (Khan & Dalvi Reference Khan and Dalvi2020). All experiments were performed at room temperature. The BSA solution was preheated up to 
$70 \pm 1\,^\circ \mathrm {C}$ to produce a stable BSA microbubble formulation (Khan & Dalvi Reference Khan and Dalvi2020). Symbols: 
$\circ$, experimental data with initial bubble diameter 
$6.9~\mathrm {\mu }\mathrm {m}$; 
$\triangle$, experimental data with initial bubble diameter 
$10.8~\mathrm {\mu }\mathrm {m}$; 
$\Diamond$, experimental data with initial bubble diameter 
$14.7~\mathrm {\mu }\mathrm {m}$. Here 
$E_{sl}=4\times 10^5~\mathrm {m}^2~\mathrm {s}^{-2}$ is used in these three model predictions.
3.2. Protein adsorption to an initially clean bubble surface
The protein adsorption to a bubble surface consists of two steps: (1) proteins in the bulk solution advect/diffuse to the adjacent layer of the bubble surface; and (2) proteins in the adjacent layer adsorb to the bubble surface (Manikantan & Squires Reference Manikantan and Squires2020). For a static spherical bubble with a kinetic-limited flux (i.e. 
$c|_{r=R} = c_\infty$), the rate of change of the protein concentration on the bubble surface is governed by
\begin{equation} \frac{\mathrm{d}\varGamma}{\mathrm{d}t} = \theta^a c_\infty \varGamma_m - (\theta^a c_\infty + \theta^d) \varGamma, \end{equation}
which yields a characteristic time scale, 
$\tau _{k} = 1/(\theta ^a c_\infty + \theta ^d )$, for kinetic-limited adsorption (Jin, Balasubramaniam & Stebe Reference Jin, Balasubramaniam and Stebe2004). For a static spherical bubble with diffusion-limited adsorption, the protein adsorption/desorption on the bubble surface is rapid such that the protein concentration on the bubble surface satisfies the adsorption isotherm.
The mass of proteins on the bubble can be estimated by 
$4 {\rm \pi}R^2\varGamma _{{eq}}$, where 
$\varGamma _{{eq}} = \theta ^a c_\infty \varGamma _m / (\theta ^a c_\infty + \theta ^d)$. Assuming that there is a thin layer (around the bubble surface) in the bulk solution supplying proteins to the surface, i.e.
\begin{equation} 4 {\rm \pi}R^2\varGamma_{{eq}} = \frac{4}{3}{\rm \pi} c_\infty [(R+h_s)^3-R^3], \end{equation}we obtain the thickness of this thin layer (Alvarez, Walker & Anna Reference Alvarez, Walker and Anna2010b,Reference Alvarez, Walker and Annaa) as
\begin{equation} h_s = R\left[\left(\frac{3h_p}{R}+1 \right)^{1/3}-1 \right] . \end{equation}
Here 
$h_s$ is called the intrinsic spherical depletion depth; and 
$h_p = \varGamma _{{eq}} / c_\infty$ is the intrinsic planar depletion depth (Jin et al. Reference Jin, Balasubramaniam and Stebe2004; Alvarez et al. Reference Alvarez, Walker and Anna2010a).
The flux boundary condition
\begin{equation} \frac{\mathrm{d}\varGamma}{\mathrm{d}t} =\left. D_c\frac{\partial c}{\partial r}\right|_{r=R} \end{equation}
has a characteristic time scale 
$\tau _{{D}1} =h_s h_p/D_c$, derived from 
$\varGamma _{{eq}}/\tau _{{D}1} \sim D_c c_\infty / h_s$. The diffusion equation (2.13) for proteins in the bulk solution has another characteristic time scale 
$\tau _{{D}2}= h_s^2 / D_c$, derived from 
$c_\infty / \tau _{{D}2} \sim D_c c_\infty / h_s^2$. Alvarez et al. (Reference Alvarez, Walker and Anna2010a) used
\begin{equation} \tau_{{D}} = (\tau_{{D}1}\tau_{{D}2})^{1/2} = \frac{(h_s^3h_p)^{1/2}}{D_c} \end{equation}to represent the characteristic time scale for diffusion-limited adsorption and validated this time scale in experiments (Alvarez et al. Reference Alvarez, Walker and Anna2010b).
If the bubble radius is much smaller than the intrinsic planar depletion depth, i.e. 
$R \ll h_p$, then
\begin{equation} \tau_{D} = \frac{\sqrt{3}\,h_p R} { D_c} = \frac{\sqrt{3}\,\theta^a \varGamma_m}{\theta^a c_\infty + \theta^d} \, \frac{R}{D_c} = \tau_{k} \left(\frac{\sqrt{3}\,\theta^a R \varGamma_m}{D_c}\right). \end{equation}
For a small bubble (
$R \ll h_p$), the characteristic time scale of the diffusion-limited adsorption is proportional to the bubble radius and inversely proportional to the protein diffusion coefficient in the bulk solution. If the bubble radius is much greater than the intrinsic spherical depletion depth, i.e. 
$R \gg h_p$, the spherical depletion depth approaches the planar depletion depth, and the characteristic time scale for the diffusion-limited adsorption is
\begin{equation} \tau_{D} = \frac{h_p^2}{D_c} = \frac{1}{D_c}\left(\frac{\theta^a \varGamma_m}{\theta^a c_\infty + \theta^d} \right)^2 = \tau_{k} \left[ \frac{(\theta^a \varGamma_m)^2} {D_c(\theta^a c_\infty +\theta^d)} \right]. \end{equation}
Thus, the diffusion time scale is inversely proportional to the protein diffusion coefficient in the bulk solution but independent of the bubble radius when 
$R \gg h_p$.
If a clean bubble is formed in a protein solution, i.e. 
$\varGamma (t=0)=0$, proteins adsorb to the bubble's surface while the bubble shrinks simultaneously. Before the surface excess stress becomes significant, the dimensionless dissolution rate is governed by (3.2a,b). When 
$\tilde {t} \ll 1$, 
$\mathrm {d}\tilde {R}/\mathrm {d}\tilde {t} \sim 1 / ({We}\,\tilde {R}\sqrt {{\rm \pi} \tilde {t}}\,)$ yields a dimensionless time scale 
$\tilde {\tau }_{{dis1}} = {We}^2$, and hence the dimensional time scale is
\begin{equation} \tau_{{dis1}} = \lambda \tilde{\tau}_{{dis1}} \sim \frac{\rho_{g0}^2 R_0^4}{\epsilon^2 \gamma_f^2 D_\psi }, \end{equation}
where 
$\rho _{g0}$ is the initial gas density. When 
$\tilde {t} \gg 1$, 
$\mathrm {d}\tilde {R}/\mathrm {d}\tilde {t} \sim 1 / ({We}\,\tilde {R}^2)$ yields another dimensionless time scale 
$\tilde {\tau }_{{dis2}} = {We}$, and hence the dimensional time scale is
\begin{equation} \tau_{{dis2}} = \lambda \tilde{\tau}_{{dis2}} \sim \frac{\rho_{g0} R_0^3}{\epsilon \gamma_f D_\psi}. \end{equation}
Thus, the characteristic time scale 
$\tau _{{dis}}$ for the bubble dissolution is inversely proportional to the gas diffusion coefficient. Additionally, 
$\tau _{{dis}}$ follows 
$R_0^4$ when 
$t \ll R_0^2 / D_\psi$; as 
$t$ grows, 
$\tau _{{dis}}$ follows 
$R_0^3$ until the surface excess stress becomes significant.
The small bubble quickly dissolves in the protein solution (
$\tau _{{dis}}$ is much less than 
$\tau _{D}$ and 
$\tau _{k}$), and hence proteins have no time to adsorb to the bubble surface. For a large bubble, 
$\tau _{{dis}}$ is much greater than 
$\tau _{D}$ and 
$\tau _{k}$ such that the bubble can be treated as static. Figure 6 shows the protein adsorption to a large (static) bubble in three different cases: 
$\tau _{D} \ll \tau _{k}$, 
$\tau _{D} \sim \tau _{k}$ and 
$\tau _{D} \gg \tau _{k}$. In the case of kinetic-limited adsorption (
$\tau _{D} \ll \tau _{k}$), 
$c_s \approx c_\infty$ and the protein adsorption rate on the bubble is constrained by the protein adsorption/desorption ability (
$\theta ^a$ and 
$\theta ^d$), as shown in figure 6(a). In the case of diffusion-limited adsorption, i.e. 
$\tau _{D} \gg \tau _{k}$, the protein adsorption rate on the bubble is constrained by 
$c_s$ (i.e. the protein diffusion coefficient 
$D_c$), and 
$\varGamma \approx \varGamma _{{iso}} = \theta ^a c_s \varGamma _m / (\theta ^a c_s + \theta ^d)$, as shown in figure 6(c).
Figure 6. Protein adsorption to a large bubble for: (a) 
$D_c = 10^{-9}\ \mathrm {m}^2~\mathrm {s}^{-1}$, 
$\tau _{{D}}=35~\mathrm {s}$; (b) 
$D_c = 10^{-11}\ \mathrm {m}^2~\mathrm {s}^{-1}$, 
$\tau _{{D}}=3500~\mathrm {s}$; and (c) 
$D_c = 10^{-12}\ \mathrm {m}^2~\mathrm {s}^{-1}$, 
$\tau _{{D}}=35\,000~\mathrm {s}$. Here 
$R_0 = 1~\mathrm {mm}$ and 
$c_\infty = 0.01~\mathrm {kg}~\mathrm {m}^{-3}$ are used; and 
$\tau _{{dis}1} = 5\times 10^{12}~\mathrm {s}$ and 
$\tau _{{dis}2} = 5\times 10^{7}~\mathrm {s}$ are much larger than 
$\tau _{D}$ and 
$\tau _{k} = 1400~\mathrm {s}$ such that the bubble radius remains almost constant.
Figure 7 shows the transition from diffusion-limited to kinetic-limited behaviour as the bubble shrinks. If we use 
$R_0 = 0.2~\mathrm {mm}$, 
$D_c = 10^{-12}~\mathrm {m}^2~\mathrm {s}^{-1}$, 
$\theta ^a = 0.05\ \mathrm {m}^3~(\mathrm {s}~\mathrm {kg})^{-1}$ and 
$\theta ^d=0.0002~\mathrm {s}^{-1}$, then 
$\tau _{D} = 20\,000~\mathrm {s}$ and 
$\tau _{k} = 1400~\mathrm {s}$. The protein adsorption initially is diffusion-limited, and increasing the protein diffusion coefficient 
$D_c$ effectively improves the protein adsorption rate. The protein concentration on the bubble surface 
$\varGamma$ grows with the adsorption of proteins and the reduction of the surface area. When 
$\varGamma$ exceeds 
$\varGamma _m$, the net flux of proteins to the interface is constrained by the protein desorption ability (
$J = -\theta ^d \varGamma$); in other words, the protein diffusion coefficient 
$D_c$ no longer influences the bubble behaviour and the protein concentration on the bubble.
Figure 7. Transition from the diffusion-limited to the kinetic-limited behaviour. Here we use 
$R_0 = 0.2~\mathrm {mm}$, 
$c_\infty = 0.01~\mathrm {kg}~\mathrm {m}^{-3}$, 
$E_{sl} = 7 \times 10^{5}\ \mathrm {m}^2~\mathrm {s}^{-2}$ and 
$\lambda _s=50~\mathrm {s}$. The red dash-dotted line is used as the reference: 
$D_c = D_{c,{ref}} = 10^{-12}\ \mathrm {m}^2~\mathrm {s}^{-1}$, 
$\theta ^a = \theta ^{a}_{ref} = 0.05\ \mathrm {m}^3~(\mathrm {s}~\mathrm {kg})^{-1}$ and 
$\theta ^d=\theta ^d_{{ref}}=0.0002~\mathrm {s}^{-1}$. Blue solid line: 
$D_c = 5D_{c,{ref}}$, 
$\theta ^a = \theta ^{a}_{ref}$ and 
$\theta ^d=\theta ^d_{{ref}}$. Orange dashed line: 
$D_c = D_{c,{ref}}$, 
$\theta ^a =5 \theta ^{a}_{ref}$ and 
$\theta ^d=5\theta ^d_{{ref}}$.
The reduction of surface area increases the protein concentration on the bubble 
$\varGamma$, but proteins desorb from the bubble simultaneously. Therefore, after 
$\varGamma$ increases to 
$\varGamma _{eq}$, the value (trend) of 
$\varGamma$ depends on the magnitude of the bubble dissolution rate, the protein diffusion coefficient in the bulk solution and the protein adsorption/desorption ability:
(i) If the bubble dissolution rate is much less than the protein desorption rate and the protein diffusion rate in the bulk solution, the protein concentration on the bubble surface
$\varGamma$ remains at equilibrium as the bubble shrinks, i.e. $\varGamma = \varGamma _{{eq}} = \theta ^a c_\infty / (\theta ^a c_\infty + \theta ^d)$.(ii) If the protein diffusion rate in the bulk solution is much lower than the bubble dissolution rate, and also the dissolution rate is much less than the protein desorption rate, then in that case proteins in the bulk solution around the bubble surface will accumulate such that
$c_s \gg c_\infty$ and $\varGamma = \varGamma _m$.(iii) If the protein desorption rate is much lower than the bubble dissolution rate, then the mass of proteins on the bubble
$\varGamma R^2$ remains almost constant.
The bubble dissolution rate depends on the surface rheology, as discussed in § 3.1. The effects of the surface rheology on the protein adsorption/desorption rate are shown in figure 8.
Figure 8. Effects of the surface rheology on the protein adsorption/desorption rate for: (a) 
$E_{sl} = 10^{4}\ \mathrm {m}^2~\mathrm {s}^{-2}$, 
$\lambda _s=1~\mathrm {s}$; (b) 
$E_{sl} = 7 \times 10^{5}\ \mathrm {m}^2~\mathrm {s}^{-2}$, 
$\lambda _s=50~\mathrm {s}$; and (c) 
$E_{sl} = 10^{6}\ \mathrm {m}^2~\mathrm {s}^{-2}$, 
$\lambda _s=100~\mathrm {s}$. Here 
$R_0 = 0.1~\mathrm {mm}$, 
$c_\infty = 0.01\ \mathrm {kg}~\mathrm {m}^{-3}$ and 
$D_c = 10^{-11}\ \mathrm {m}^2~\mathrm {s}^{-1}$ are used; and 
$\tau _k = 1400~\mathrm {s}$, 
$\tau _{D} = 6600~\mathrm {s}$, 
$\tau _{{dis}1} = 5\times 10^{8}\ \mathrm {s}$ and 
$\tau _{{dis}2} = 5\times 10^{4}~\mathrm {s}$.
The interfacial parameters, even for BSA, depending on the experimental procedures, can vary by orders of magnitude. Thus, identifying interfacial parameters that have significant effects on the protein adsorption rate helps the design of experiments. Table 2 summarizes the bubble dynamics and the corresponding important interfacial parameters in different scenarios. Before the protein concentration on the bubble 
$\varGamma$ reaches 
$\varGamma _{{eq}}$, 
$\varGamma$ grows as the bubble shrinks due to the adsorption of proteins and the reduction of surface area. The rate of protein adsorption depends on the protein adsorption/desorption ability, protein diffusion coefficient and the bubble dissolution rate. Before the surface excess stress becomes significant and balances the surface tension force, the bubble dissolution rate depends on the gas density, gas diffusion coefficient and surface tension. However, after the surface excess stress balances the surface tension force, the surface dilatational elasticity 
$E_s$ and the surface relaxation time 
$\lambda _s$ become important, and the gas properties become insignificant. If the protein concentration on the bubble grows beyond 
$\varGamma _{{eq}}$ but smaller than 
$\varGamma _m$, the protein concentration on the bubble is still affected by the protein diffusion coefficient, protein adsorption/desorption ability and the bubble dissolution rate. Once the protein concentration on the bubble exceeds 
$\varGamma _m$, the net flux of proteins to the bubble only depends on the protein desorption ability 
$\theta ^d$; therefore, the protein adsorption ability 
$\theta ^a$ and the protein diffusion coefficient 
$D_c$ become unimportant.
Table 2. The bubble dynamics and the corresponding important interfacial parameters.
4. Response of a protein-coated bubble to an imposed fluctuating pressure
In this case, we consider a protein-coated bubble oscillating in a pressure field 
$p_\infty = p_a [1+\varepsilon \cos (\omega t)]$. Let us consider a small fluctuating pressure (
$\varepsilon \ll 1$), that the system is linear, and hence the bubble responds at the driving frequency at the steady state (Brennen Reference Brennen2013). To simplify this problem, we will neglect gas diffusion. We also ignore the adsorption/desorption of protein since its time scale (of the order of seconds) generally is much greater than the period of imposed pressure (of the order of microseconds or milliseconds). Additionally, we ignore the variation of the protein concentration in the bulk of the liquid such that the bulk viscosity 
$\mu$ is independent of radial position 
$r$.
Let subscript 
$E$ denote the equilibrium state. We write the complex form of the pressure at infinity as 
$p_\infty = p_a + \varepsilon p_a \mathrm {e}^{\mathrm {j}\omega t}$, where 
$\mathrm {j} = \sqrt {-1}$; then its physical form (projection onto the real axis) is 
$p_\infty = p_a + \mathrm {Re}\{\varepsilon p_a \mathrm {e}^{\mathrm {j}\omega t}\}$. Similarly, we express the complex form of the bubble radius 
$R$, protein concentration at the bubble surface 
$\varGamma$, gas pressure inside the bubble 
$p_g$ and surface tension 
$\gamma$ as
\begin{align} R =
R_E(1+\hat{R}\mathrm{e}^{\mathrm{j}\omega
t}), \varGamma =
\varGamma_E(1+\hat{\varGamma}\mathrm{e}^{\mathrm{j}\omega
t}), p_g = p_{gE}
(1+\hat{p}_{g}\mathrm{e}^{\mathrm{j}\omega
t}), \gamma = \gamma_E
(1+\hat{\gamma}\mathrm{e}^{\mathrm{j}\omega
t}).
\end{align}
Substituting them into (2.16), (2.17) and (2.24) and then linearizing, we get 
$\hat {\varGamma } = -2\hat {R}$, 
$\hat {p}_g=-3\alpha \hat {R}$ and 
$\hat {\gamma } = \chi \hat {R}$, where
\begin{equation} \chi = \left\{
\begin{array}{@{}ll} \left(\dfrac{R_g T
\varGamma_{p}}{\gamma_E}\right)
\left(\dfrac{2\varGamma_E}{\varGamma_m -
\varGamma_E}\right), & \text{when }\varGamma_E \leq
\varGamma_c,\\ 0, & \text{when }\varGamma_E > \varGamma_c.
\end{array} \right. \end{equation}
We assume 
$\mathcal {F} = \mathcal {F}_E + \hat {\mathcal {F}} \mathrm {e}^{\mathrm {j}\omega t}$ and substitute it into 
$\dot {\mathcal {F}} + \mathcal {F}/\lambda _s = \dot {R}/R$. After linearization we obtain 
$\mathcal {F}_E = 0$ and 
$\hat {\mathcal {F}} = (\mathrm {j}\omega \lambda _s + \omega ^2\lambda _s^2)\hat {R}/(1+\omega ^2 \lambda _s^2)$, and thus the surface excess stress is zero at the equilibrium state. Substituting the expressions of 
$p_\infty$, 
$R$, 
$\varGamma$, 
$\gamma$ and 
$\mathcal {F}$ into (2.9) and linearizing, we get the zeroth-order solution
\begin{equation} p_v + p_{gE} - p_a - 2\gamma_E/R_E = 0, \end{equation}and the first-order solution
\begin{align} \frac{\varepsilon p_a}{\rho R_E^2 \hat{R}} &= \omega^2 - \frac{4\mathrm{j}\omega }{\rho R_E^2}\left[\mu +\frac{\kappa_{sE}}{R_E} \, \frac{1}{1+(\omega \lambda_s)^2} \right] \nonumber\\ &\quad -\frac{1}{\rho R_E^2} \left[3\alpha (p_a - p_v)+ \frac{2\gamma_E}{R_E}(\chi + 3\alpha -1)+\frac{4E_{sE}}{R_E} \, \frac{(\omega \lambda_s)^2} {1+(\omega\lambda_s)^2}\right], \end{align}
where 
$\kappa _{sE}=\kappa _{sl}\varGamma _E$ and 
$E_{sE}=E_{sl}\varGamma _E$ are the surface dilatational viscosity and dilatational elasticity at the equilibrium state, respectively.
Equation (4.4) is a solution of a forced oscillation, in which a mass 
$m = \rho R_E^3$, attached to a spring and a dashpot in parallel arrangement, oscillates under an external force 
$F_e=\varepsilon p_a R_E^2 \cos (\omega t)$, as shown in figure 9. The stiffness of the spring, 
$s$, and the resistance from the dashpot, 
$R_m$, are
$$\begin{gather} s = m \omega_0^2 + 2\chi\gamma_E +4E_{sE} \frac{(\omega \lambda_s)^2} {1+(\omega\lambda_s)^2}, \end{gather}$$
$$\begin{gather}R_m = 2m \beta_0 + \frac{4\kappa_{sE}}{1+(\omega \lambda_s)^2}, \end{gather}$$
where 
$\omega _0 = \{[3\alpha (p_a - p_v)R_E + 2\gamma _E(3\alpha -1)]/(\rho R_E^3)\}^{1/2}$ and 
$\beta _0 = 2\mu /(\rho R_E^2)$ are the natural frequency and the damping coefficient of the uncoated bubble, respectively. The restoring force provided by the viscoelastic bubble surface includes an elastic force (due to the surface dilatational elasticity) and a viscous force (due to the surface dilatational viscosity). The surface dilatational viscosity contributes to the resistance term, 
$4\kappa _{sE}/(1+\omega ^2\lambda _s^2)$, which decreases with 
$\omega \lambda _s$. The surface dilatational elasticity contributes to the stiffness term, 
$4E_{sE}\omega ^2\lambda _s^2/(1+\omega ^2\lambda _s^2)$, which increases with 
$\omega \lambda _s$. Additionally, the variation of surface tension with protein concentration at the bubble surface, 
$\chi$, adds additional stiffness to the system, which is independent of the acoustic frequency.
Figure 9. Effects of surface dilatational viscosity and dilatational elasticity on the resistance and stiffness of the system. Here 
$\omega _0$ and 
$\beta _0$ are the natural frequency and damping coefficient of an uncoated bubble, respectively.
Let 
$\beta = R_m / (2m)$ and 
$\omega _n = \sqrt {s/m}$ denote the damping coefficient and the natural frequency, respectively. The amplitude of the dimensionless bubble radius is
\begin{equation} \varLambda = |\hat{R}| = \frac{\varepsilon p_a}{\rho R_E^2 \sqrt{4\beta^2 \omega^2 + (\omega_n^2 - \omega^2)^2}}. \end{equation}We find a critical frequency
\begin{equation} \omega_c = \sqrt{\omega_0^2 + \frac{2}{\rho R_E^2}\left(\frac{\chi\gamma_E}{R_E}+\frac{E_{sE}}{R_E}+\frac{2\mu}{\lambda_{s}} \right)}, \end{equation}
at which the viscoelasticity of the bubble surface does not affect the amplitude of oscillation: 
$\varLambda (\kappa _{sE}, E_{sE}) \geq \varLambda (\kappa _{sE}=0, E_{sE}=0)$ if 
$\omega \geq \omega _c$, otherwise 
$\varLambda (\kappa _{sE}, E_{sE})< \varLambda (\kappa _{sE}=0, E_{sE}=0)$. Thus, for a Newtonian interface (surface relaxation time 
$\lambda _s$ approaches zero and hence the critical frequency 
$\omega _c$ becomes infinitely large), including the surface dilatational viscosity reduces the oscillation amplitude. For a viscoelastic bubble surface, viscoelasticity enlarges the oscillation amplitude if and only if 
$\omega > \omega _c$ (i.e. the imposed acoustic pressure has a high frequency). We define 
$\mathcal {A}\equiv \varLambda / \varLambda _0$, where 
$\varLambda _0$ is the amplitude of the bubble radius at zero surface dilatational elasticity and dilatational viscosity. The effects of surface dilatational elasticity and dilatational viscosity on the dimensionless amplitude 
$\mathcal {A}$ are shown in figure 10.
Figure 10. Effects of surface dilatational viscosity and surface dilatational elasticity on the amplitude of the dimensionless bubble radius 
$\mathcal {A}$. Here 
$R_E = 10^{-5}~\mathrm {m}$, 
$\mu = 0.001\ \mathrm {Pa}~\mathrm {s}$, 
$\gamma _E = 0.052\ \mathrm {N}~\mathrm {m}^{-1}$ and 
$\chi = 0$ are used. The natural frequency of the uncoated bubble is 
$\omega _0 = 1.774\times 10^6 \mathrm {rad}~\mathrm {s}^{-1}$. (a) The frequency of the acoustic pressure is 
$\omega = 1.8\times 10^6\ \mathrm {rad}~\mathrm {s}^{-1}$, which is slightly larger than 
$\omega _0$; and (b) here 
$\omega =\omega _0$.
The resonant frequency, 
$\omega _r$, and the maximum amplitude of the dimensionless bubble radius, 
$\varLambda _{max} = |\hat {R}|_{max}$, satisfy
\begin{equation} \omega_r^2 = \omega_n^2 - 2\beta^2, \quad \varLambda_{max} = \left(\frac{\varepsilon p_a}{\rho R_E^2}\right) \left(\frac{1}{2\beta \sqrt{\omega_n^2 -\beta^2}} \right). \end{equation}
The natural frequency 
$\omega _n$ and the damping coefficient 
$\beta$ vary with the frequency of the imposed pressure 
$\omega$; thus we need to solve 
$\omega _r^2 = \omega _n^2 -2\beta ^2$ numerically to obtain the resonant frequency 
$\omega _r$. However, the explicit resonance frequency can be obtained in the following two limiting cases: (1) if 
$\omega \lambda _s \ll 1$, the bubble surface is Newtonian and hence it only contributes resistance, the natural frequency is 
$\omega _n = \{\omega _0^2+2\chi \gamma _E/m\}^{1/2}$ and the damping coefficient of the bubble is 
$\beta = \beta _0 + 2\kappa _{sE}/m$; and (2) if 
$\omega \lambda _s \gg 1$, the bubble surface is purely elastic, which only contributes stiffness, then the natural frequency is 
$\omega _n = \{\omega _0^2+(2\chi \gamma _E+4E_{sE})/m\}^{1/2}$ and the damping coefficient of the bubble is 
$\beta = \beta _0$. The effects of surface dilatational viscosity and surface dilatational elasticity on the resonant amplitude and frequency are shown in figure 11.
Figure 11. Effects of surface dilatational viscosity and dilatational elasticity on the resonant frequency and resonant amplitude. Here 
$\varepsilon = 0.001$, 
$\mu = 0.001\ \mathrm {Pa}~\mathrm {s}$, 
$\gamma _E = 0.052\ \mathrm {N}~\mathrm {m}^{-1}$ and 
$\chi = 0$ are used.
If we want to measure the dilatational elasticity and the relaxation time of the viscoelastic bubble surface, we can consider performing two experiments: (1) Impose a small fluctuating pressure with a frequency near the natural frequency of the uncoated bubble. In this case, 
$\omega \lambda _s$ generally is much greater than one, and hence the amplitude of oscillation is independent of the relaxation time; thus, we can use (4.7) to determine the surface dilatational elasticity. (2) After obtaining the surface dilatational elasticity, we can adjust the frequency to around 
$1/\lambda _s$ (generally in the range of 
$0.001 - 1~\mathrm {s}^{-1}$) and get the corresponding oscillation amplitude; next, we use (4.7) to determine the relaxation time. However, in the second experiment with a low frequency, we need to minimize the effect of gas diffusion on the bubble dynamics (e.g. we can do experiments in a saturated solution (Epstein & Plesset Reference Epstein and Plesset1950) or use gas that has a low diffusion coefficient in the liquid). Additionally, if the ratio of the surface dilatational elasticity to the equilibrium bubble radius, 
$E_{sE}/R_E$, is much less than atmospheric pressure or the surface tension 
$2\gamma _{E}/R_E$, the surface dilatational elasticity or the relaxation time has minor effects on oscillations, which will bring difficulties to experiments.
5. Conclusions
We have developed a model for bubble dynamics in a protein solution, including varying surface tension and dynamic adsorption/desorption of protein onto the bubble surface. We use the Maxwell model to describe the protein-coated bubble surface because the gas–liquid interface with an adsorbed protein layer is usually viscoelastic (Narsimhan Reference Narsimhan2016). If the surface relaxation time goes to zero, then our derived formula for a viscoelastic bubble surface reduces to the Boussinesq–Scriven constitutive equation for a Newtonian interface, which is suitable for a surfactant- or lipid-coated bubble (Chatterjee & Sarkar Reference Chatterjee and Sarkar2003; Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005).
We have studied two cases to understand the importance of surface rheology. In the first case study, in which the bubble dissolves in a protein solution, we find that the protein reduces the dissolution rate of the bubble. If we ignore the situation of buckling and protein shedding, the dissolution can be divided into three processes: (1) inertio-capillary process, in which the rate of dissolution is determined by the Weber number (the ratio of gas inertial force to the surface tension force); (2) inertio-elasto-capillary process, in which the surface excess stress counteracts the surface tension force such that the dissolution rate reduces; and (3) elasto-capillary process, in which the surface excess stress balances the surface tension force, and the dissolution rate is governed by both the protein desorption rate and the elasto-capillary number (the ratio of surface tension to the surface dilatational elasticity). If the protein-coated bubble surface has a high dilatational elasticity and a large relaxation time, the surface excess stress will quickly grow and balance the surface tension force such that the first two processes can be neglected; in that case, the dissolution rate can be predicted by (3.7). Suppose we have experimental data for the bubble radius with respect to time; we can use (3.7) to fit the experimental data and find the protein desorption rate, surface relaxation time and surface dilatational elasticity.
If a large clean bubble is formed in a protein solution, the dissolution rate is negligible. Hence the protein adsorption rate is determined by the protein diffusion coefficient in the bulk solution 
$D_c$ and the protein adsorption/desorption ability (
$\theta ^a$ and 
$\theta ^d$). If the bubble dissolution rate is not negligible compared with the protein adsorption rate, the reduction of the surface area increases the protein concentration on the bubble surface. After the protein concentration on the bubble surface 
$\varGamma$ reaches the equilibrium value, 
$\varGamma _{{eq}}=\theta ^a c_\infty / (\theta ^a c_\infty + \theta ^d)$. If the protein desorption rate is small compared with the bubble dissolution rate, the protein mass on the bubble will remain almost constant, i.e. the protein concentration on the bubble continues growing as the bubble shrinks. If the protein desorption rate is large, then the protein concentration on the bubble surface will satisfy the adsorption isotherm. In other words, 
$\varGamma = \varGamma _m$ if the protein diffusion coefficient 
$D_c$ is small; and 
$\varGamma = \varGamma _{{eq}}$ if the protein diffusion rate in the bulk solution is large compared with the dissolution rate.
In the second case study, in which a small fluctuating pressure drives a protein-coated bubble to oscillate, we have derived a formula for the oscillation at the steady state. We find that the surface dilatational viscosity contributes to bubble damping, while the surface dilatational elasticity contributes to bubble stiffness. Although we cannot obtain an explicit solution for the resonant frequency in the general case, we get the explicit solution for a bubble with a Newtonian interface (surface relaxation time goes to zero) or an elastic interface (surface relaxation time goes to infinity). We also find a critical frequency (greater than the natural frequency of an uncoated bubble), beyond which including the surface dilatational elasticity and the surface dilatational viscosity enlarges the amplitude of oscillation.
Our model is applicable for globular proteins. For the fibrillar structure of proteins, the molecules at the interface can have different conformations, and hence the protein adsorption isotherms (2.15) and (2.17) need to be modified or replaced with appropriate expressions. Additionally, our model is valid if the concentrations are below the critical micelle concentration of conventional non-ionic/ionic surfactants in the presence of salt, or below the critical aggregation concentration of proteins. If the concentrations are above the critical micelle concentration or the critical aggregation concentration, the solutions contain monomers and micelles/aggregates, and we need an additional equation to describe the kinetics of micelles/aggregates (Patist et al. Reference Patist, Kanicky, Shukla and Shah2002; Craster et al. Reference Craster, Matar and Papageorgiou2009; Morris, Watzky & Finke Reference Morris, Watzky and Finke2009). Moreover, conventional surfactants usually have a large desorption rate; thus, the bubble under shrinkage can remain spherical when the surfactant concentration on the bubble surface is close to the maximum possible concentration. However, if the protein adsorption is close to the maximum possible adsorption, the fast shrinkage of the bubble usually leads to a non-spherical bubble shape (the bubble volume decreases while the bubble surface area remains constant) since the protein desorption rate is low. Thus, our results are qualitative for relatively large concentrations of protein solutions.
Our findings in the first case study indicate that the Maxwell model can characterize the protein-coated bubble surface. Also, surface rheology can be important and needs to receive more attention. Additionally, our theoretical derivation for the dissolution rate and the bubble response to an acoustic pressure provide a way to predict the properties of a viscoelastic bubble surface. This work provides insight into protein-coated bubbles and helps guide the design of protein-coated bubbles, which are used as vehicles to deliver drugs, as active miniature tracers to probe the rheology of soft and biological materials (Dollet et al. Reference Dollet, Marmottant and Garbin2019), or as contrast agents to enhance the ultrasound backscatter in ultrasonic imaging (Dauba et al. Reference Dauba, Delalande, Kamimura, Conti, Larrat, Tsapis and Novell2020).
Funding
This research was supported by a grant from the Eli Lilly and Company.
Declaration of interests
The authors report no conflict of interest.
Open access
