Hostname: page-component-7b7567c8c-fvjsw Total loading time: 0 Render date: 2026-11-14T09:16:51.576Z Has data issue: false hasContentIssue false
  • English
  • Français

Mathematical modelling for acoustic microstreaming produced by a gas bubble undergoing asymmetric oscillations

Published online by Cambridge University Press:  08 May 2025

Claude InserraOpen the ORCID record for Claude Inserra[Opens in a new window],
Cyril Mauger,
Philippe Blanc-Benon and
Alexander A. Doinikov
Claude Inserra*
Affiliation:
University of Lyon, Université Claude Bernard Lyon 1, Centre Léon Bérard, INSERM, UMR 1032, LabTAU, Lyon F-69003, France
Cyril Mauger
Affiliation:
INSA Lyon, CNRS, École Centrale de Lyon, Université Claude Bernard Lyon 1, UMR 5509, LMFA, Villeurbanne 69621, France
Philippe Blanc-Benon
Affiliation:
INSA Lyon, CNRS, École Centrale de Lyon, Université Claude Bernard Lyon 1, UMR 5509, LMFA, Villeurbanne 69621, France
Alexander A. Doinikov
Affiliation:
INSA Lyon, CNRS, École Centrale de Lyon, Université Claude Bernard Lyon 1, UMR 5509, LMFA, Villeurbanne 69621, France
*
Corresponding author: Claude Inserra, claude.inserra@inserm.fr
Rights & Permissions [Opens in a new window]

Abstract

An exact solution is developed for bubble-induced acoustic microstreaming in the case of a gas bubble undergoing asymmetric oscillations. The modelling is based on the decomposition of the solenoidal, first- and second-order, vorticity fields into poloidal and toroidal components. The result is valid for small-amplitude bubble oscillations without restriction on the size of the viscous boundary layer $(2\nu /\omega )^{1/2}$ in comparison to the bubble radius. The non-spherical distortions of the bubble interface are decomposed over the set of orthonormal spherical harmonics $Y_{n}^{m}(\theta , \phi )$ of degree $n$ and order $m$. The present theory describes the steady flow produced by the non-spherical oscillations $(n,\pm m)$ that occur at a frequency different from that of the spherical oscillation, as in the case of a parametrically excited surface oscillation. The three-dimensional aspect of the streaming pattern is revealed as well as the particular flow signatures associated with different asymmetric oscillations.

JFM classification

Drops and Bubbles: Drops and Bubbles

Information

Type
JFM Papers
Information
Journal of Fluid Mechanics , Volume 1010 , 10 May 2025, A48
Check for updates
Copyright
© The Author(s), 2025. Published by Cambridge University Press

1. Introduction

The high-frequency interfacial motion of ultrasound-driven gas bubbles can generate a steady flow in the surrounding fluid, called acoustic microstreaming. Interest in the bubble-induced streaming and the resulting stresses exerted on surrounding bodies is motivated by a variety of technological, chemical and biomedical applications. Oscillating bubbles close to or in contact with a wall can actuate the transport of particles within the viscous boundary layer at solid surfaces, with applications in the removal of contamination particles and surface cleaning (Kim et al. Reference Kim, Kim, Choi and Kim2009; Reuter et al. Reference Reuter, Lauterborn, Mettin and Lauterborn2017), the amplification of chemical treatments (Mason Reference Mason1999), the handwashing and cleaning of sensitive surfaces (Birkin, Offin & Leighton Reference Birkin, Offin and Leighton2016), and the micromixing of fluids (Wang, Rallabandi & Hilgenfeldt Reference Wang, Rallabandi and Hilgenfeldt2013). Cavitation microstreaming is also a propelling source for small water-floating objects and bubble-driven microrobots (Dijkink et al. Reference Dijkink, van der Dennen, Ohl and Prosperetti2006; Ahmed et al. Reference Ahmed, Lu, Nourhani, Lammert, Stratton, Muddana, Crespi and Huang2015). These artificial, acoustically driven microswimmers act as propulsion devices that can carry a payload, with a great potential for microfluidic applications and targeted drug delivery. The commonly used geometry for these microrobots includes a gas bubble trapped in the robot body cavity, whose oscillations create an axisymmetric flow that orients and propels the microrobot (Zhou, Dai & Jiao Reference Zhou, Dai and Jiao2022). The orientation of the propeller can be controlled by inducing one or several holes on armoured bubbles (Bertin et al. Reference Bertin, Spelman, Stephan, Gredy, Bouriau, Lauga and Marmottant2015) by combining acoustic powering and magnetic steering (Aghakhani et al. Reference Aghakhani, Yasa, Wrede and Sitti2020), or by designing gear-like bubble-based propellers in order to actuate rotational motion along closed trajectories (Mohanty et al. Reference Mohanty, Zhang, McNeill, Kuenen, Linde, Rouwkema and Misra2021). Biomedical applications based on bubble-induced flows include cell detachment (Ohl & Wolfrum Reference Ohl and Wolfrum2003), the sorting and manipulation of biological materials (Volk et al. Reference Volk, Rossi, Rallabandi, Kähler, Hilgenfeldt and Marin2020), the lysis of vesicles (Marmottant, Biben & Hilgenfeldt Reference Marmottant, Biben and Hilgenfeldt2008), and ultrasound-mediated targeted drug delivery (Lajoinie et al. Reference Lajoinie, De Cock, Coussios, Lentacker, Le Gac, Stride and Versluis2016; Pereno et al. Reference Pereno2018). The ultrasound-mediated delivery of a drug is based on the action of oscillating microbubbles nearby biological barriers that increase their permeability and allow drug and genes to penetrate into individual cells without serious consequence for the cell viability (Fan et al. Reference Fan, Liu, Ting, Lee, Huang, Ma, Wei, Yen and Yeh2014). The temporary permeabilisation of biological barriers is caused by shear stresses exerted on cell tissues by the bubble-induced flows (Wu & Nyborg Reference Wu and Nyborg2008), which are responsible for a ‘massage’ effect on cells and the creation of transient pores on the cell membrane. Even if no consensus exists on the exact mechanism responsible for cell poration and the required bubble activity (collapsing/inertial regime or stably oscillating regime), determining the bubble-induced flows and resulting shear stress is mandatory.

The first theoretical investigation of acoustic microstreaming was performed by Nyborg (Reference Nyborg1958). By investigating the near-boundary streaming induced by a compressible body (a gas bubble) resting on a surface, he showed how resonant bubbles produced a pronounced microstreaming in the surface vicinity. Davidson & Riley (Reference Davidson and Riley1971) have considered the case of a spherical bubble oscillating laterally in an unbounded fluid. Their investigation covers a wide range of situations with respect to the orders of magnitude of two dimensionless parameters. The first, $\epsilon =U/R_{0}\omega$ , is the ratio of the vibration amplitude of the velocity $U$ occuring at the angular frequency $\omega$ to the bubble radius $R_{0}$ . The second parameter quantifies the ratio of the thickness of the Stokes layer to the bubble radius, $\gamma =(2\nu /\omega )^{1/2}/R_{0}=\delta _{v}/R_{0}$ , where $\nu$ is the kinematic viscosity, and $\delta _{v}$ is the thickness of the oscillatory shear layer. A great theoretical achievement of their work is the introduction of a matching approximation between the solution within the inner boundary layer and the one in the outer boundary, performed in the case of large bubbles ( $\gamma \ll 1$ ) driven at relatively low frequencies and with small amplitudes of lateral oscillations ( $\epsilon \ll 1$ ). This approximation has been used in later decades by several authors, adding the contribution of small radial oscillation (Longuet-Higgins Reference Longuet-Higgins1998), of small-amplitude axisymmetric shape oscillations in the case of the $n$ th distortion mode with $n\gg 1$ (Maksimov Reference Maksimov2007), and of any arbitrary combination of axisymmetric shape oscillations (Spelman & Lauga Reference Spelman and Lauga2017). All the above-mentioned theoretical works are based on the matching of the inner/outer solution that assumes a small viscous penetration depth in comparison to the bubble radius. This assumption limits the findings to the case of large bubbles in low-viscosity fluids. Recently, Doinikov et al. (Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a ) have overcome these limitations by calculating exactly the second-order mean flow induced by all possible interactions between axisymmetric shape oscillations (including the spherical and translational ones).

Yet the similarity between all these theoretical derivations is to consider initially (at rest) spherical bubbles far from any boundary. While this scenario has been recovered experimentally by using levitating, acoustically trapped bubbles (Cleve et al. Reference Cleve, Guedra, Mauger, Inserra and Blanc-Benon2019), the majority of experimental work on acoustic microstreaming is performed on substrate-attached microbubbles. The positional stability of the bubble is therefore ensured, hence facilitating the capture of the interface dynamics as well as the surrounding fluid motion. Marmottant et al. (Reference Marmottant, Versluis, de Jong, Hilgenfeldt and Lohse2006) resolved the acoustic streaming surrounding a wall-attached bubble experiencing spherical oscillations and a translational one occurring perpendicularly to the wall. Tho, Manasseh & Ooi (Reference Tho, Manasseh and Ooi2007) performed an extensive study of streaming patterns surrounding a substrate-attached bubble from the top view. In addition to the cases of varying translational and/or oscillating motion of the bubble interface, the authors have also extended the analysis to the case of shape oscillations, without mentioning the triggered shape instability. Marin et al. (Reference Marin, Rossi, Rallabandi, Wang, Hilgenfeldt and Kähler2015) revealed the three-dimensional nature of the acoustic streaming flow surrounding a wall-attached cylindrical bubble by using an astigmatism particle tracking velocimetry (APTV) technique. Using the same APTV technique, Bolanos-Jimenez et al. (Reference Bolanos-Jimenez, Rossi, Fernandez Rivas, Kähler and Marin2017) captured the three-dimensional axisymmetric, ‘fountain-like’ flow pattern surrounding a hemispherical bubble. Interestingly, the strength of the flow was used as a way of finding the bubble’s lowest resonant frequency. When summarizing the literature, it becomes obvious that experimental and theoretical works on bubbles experiencing asymmetric (i.e. non-axisymmetric) oscillations are scarce, while the ease of triggering asymmetric deformations is facilitated by the contact line dynamics and the breaking of the spherical symmetry. Usually, the complexity of the asymmetric oscillations is disregarded experimentally, where undetermined shape modes are sometimes reported (Saint-Michel & Garbin Reference Saint-Michel and Garbin2020). The landscape of emergence of specific spherical harmonics for high-amplitude-driven microbubbles was proposed by Fauconnier, Bera & Inserra (Reference Fauconnier, Bera and Inserra2020). The ease of triggering specific asymmetric modes near the minimum of the instability threshold, as well as the partitioning of the triggered modes within the instability region of existence for shape oscillations, was demonstrated. The only theoretical work discussing the emergence and triggering of asymmetric oscillations was performed by Maksimov (Reference Maksimov2020) in the case of a spherical bubble located in front of a wall. The splitting (partitioning) of the shape modes was recovered theoretically. Concerning the induced flows, a systematic study of the acoustic microstreaming surrounding a wall-attached bubble has been performed by Fauconnier et al. (Reference Fauconnier, Mauger, Bera and Inserra2022). High-amplitude acoustic driving allowed the triggering of non-spherical shape oscillations, including asymmetric ones. The flow signatures were correlated to the bubble interface dynamics decomposed over the set of orthonormal spherical harmonics $Y_{n}^{m}(\theta , \phi )$ of degree $n$ and order $m$ . The self-interaction of a given asymmetric shape oscillation results in various flower-like patterns whose numbers of lobes are associated with the degree $n$ and order $m$ of the triggered harmonics. Only top view observations were performed, and the three-dimensional nature of these asymmetric patterns has still to be determined.

The present paper provides the mathematical modelling of the second-order mean flow surrounding a bubble experiencing asymmetric oscillations at a frequency different from that of the spherical oscillations, so that the microstreaming comes from the interaction of the shape oscillation with itself (self-interaction). This case is relevant to parametrically excited shape oscillations. Section 2 describes the derivation of the second-order mean flow, with a particular focus on the decomposition of the solenoidal, first- and second-order, vorticity fields into poloidal and toroidal components. The derivation is valid for small-amplitude bubble oscillations with no restriction on the size of the viscous boundary layer $\delta _{v}=(2\nu /\omega )^{1/2}$ in comparison to the bubble radius. This means that our theory is valid for any value of the liquid viscosity $\nu$ , which in turn means that our theory is valid for any value of the ratio $\delta _{v}/R_{0}$ , where $R_{0}$ is the equilibrium bubble radius. We emphasise this fact because in many previous studies on acoustic streaming, the ratio $\delta _{v}/R_{0}$ is assumed to be small, which means that the liquid viscosity is assumed to be very low and/or $R_{0}$ to be big. Section 3 presents various numerical examples obtained by the proposed model, highlighting the three-dimensional nature of the streaming patterns, the reversal of the flow atop the bubble by means of the generation of an anti-fountain behaviour, and the possibility of creating flows with a strong azimuthal component when a travelling surface wave propagates at the interface of the bubble.

2. Theory

We consider a gas bubble surrounded by an infinite viscous incompressible liquid. We assume that the bubble, which is spherical at rest, undergoes asymmetric oscillations in response to an external acoustic excitation. Figure 1 shows a spherical coordinate system, originated at the equilibrium centre of the bubble, that is used in our calculations. Our derivation follows the conventional procedure for calculating acoustic streaming. We assume that the amplitudes of the bubble oscillation modes are small compared to the equilibrium bubble radius. This assumption allows us to linearise the equations of liquid motion (Navier–Stokes equations) and to find their solutions, assuming the amplitudes of the bubble oscillation modes to be given quantities. These solutions give us a first-order velocity field generated by the bubble in the liquid. Then the equations of liquid motion are written with accuracy up to terms of second order of smallness with respect to the first-order solutions and averaged over time. As a result, we obtain equations that describe the velocity field of acoustic streaming produced by the bubble.

Figure 1.

Coordinate system used in calculations.

2.1. First-order solutions

The surface of a bubble undergoing the oscillation modes $(n,m)$ and $(n,-m)$ can be represented by

(2.1) \begin{equation}r_{s}(\theta ,\phi ,t)=R_{0}+e^{-i\omega t}\,s^{(n,m)}\,Y_{n}^{m}(\theta ,\phi )+e^{-i\omega t}\,s^{(n,-m)}\,Y_{n}^{-m}(\theta ,\phi ),\end{equation}

where $R_{0}$ is the bubble radius at rest, $\omega$ is the angular frequency of the bubble oscillations, $s^{(n,m)}$ and $s^{(n,-m)}$ are the complex amplitudes of the modes $(n,m)$ and $(n,-m)$ , respectively, and $Y_{n}^{m}(\theta ,\phi )$ are spherical harmonics of degree $n$ and order $-n\leq m\leq n$ , which are defined by (D1) in Appendix D. Here, the azimuthal dependency as a complex exponential has been preferred to the cosine form (known as the real spherical harmonics) for the following reason. The cosine representation of the azimuthal part of the angular function implies an established, standing wave at the interface of the bubble. In order to allow for the most general case of two oppositely moving travelling waves at the interface of the bubble, the two components $(n,\pm m)$ are considered. When the amplitudes of the two components are equal, $s^{(n,m)}=s^{(n,-m)}$ , a stationary azimuthal wave is established. When they differ, a quasi-stationary wave exists that is composed of a partially stationary and a partially travelling component. In the general case, any three-dimensional bubble deformation can be decomposed over the set of orthonormal spherical harmonics, which means that several modes can coexist. The interaction of several modes oscillating at the same frequency (mixed-mode streaming) has been described as a dominant contribution to acoustic microstreaming in many situations; see, for example, Longuet-Higgins (Reference Longuet-Higgins1998) and Rallabandi, Wang & Hilgenfeldt (Reference Rallabandi, Wang and Hilgenfeldt2014). However, the consideration of two interacting modes $(n_{1},m_{1})$ and $(n_{2}\neq n_{1},m_{2}\neq m_{1})$ will lead to serious mathematical complexity in calculating the resulting flow. According to the work of Fauconnier et al. (Reference Fauconnier, Mauger, Bera and Inserra2022), the possibility of triggering a single, well-identified, asymmetric mode was demonstrated. The associated flow, resulting from the interaction of the triggered asymmetric mode with itself, was captured. It is therefore worth considering a single asymmetric mode interacting with itself for the description of the induced acoustic microstreaming, as in (2.1).

The values of $s^{(n,\pm m)}$ and $\omega$ are considered as known quantities. They are assumed to be measured experimentally and serve as input data in our study. We also assume that $|s^{(n,\pm m)}|/R_{0}\ll 1$ , which allows us to linearise the equations of liquid motion. The linearised equations of a viscous incompressible liquid are given by (Landau & Lifshitz Reference Landau and Lifshitz1987)

(2.2) \begin{align} \boldsymbol{\nabla }\cdot \boldsymbol{v}_{1}=0 ,\qquad\quad\end{align}
(2.3) \begin{align} \frac{\partial \boldsymbol{v}_{1}}{\partial t}=-\frac{1}{\rho }\,{\boldsymbol \nabla }p_{1}+\nu\, \Delta \boldsymbol{v}_{1} ,\end{align}

where $\boldsymbol{v}_{1}$ and $p_{1}$ are the first-order liquid velocity and pressure, respectively, $\rho$ is the constant liquid density, $\nu =\eta /\rho$ is the kinematic liquid viscosity, and $\eta$ is the dynamic liquid viscosity.

The first-order velocity field generated by modes $(n,m)$ and $(n,-m)$ can be written as

(2.4) \begin{equation} \boldsymbol{v}_{1}=\boldsymbol{v}_{1}^{\left(n,m\right)}+\boldsymbol{v}_{1}^{\left(n,-m\right)} ,\end{equation}

where $\boldsymbol{v}_{1}^{(n,m)}$ and $\boldsymbol{v}_{1}^{(n,-m)}$ are the first-order liquid velocities generated by the modes $(n,m)$ and $(n,-m)$ , respectively. Both $\boldsymbol{v}_{1}^{(n,m)}$ and $\boldsymbol{v}_{1}^{(n,-m)}$ obey (2.2) and (2.3).

The commonly used procedure to solve (2.2) and (2.3) consists in decomposing the first-order liquid velocity $\boldsymbol{v}_{1}$ into the scalar $\varphi$ and the vector $\boldsymbol{\psi }$ velocity potentials using the Helmholtz decomposition $\boldsymbol{v}_{1}=\nabla \varphi +\nabla \times \boldsymbol{\psi }$ . In spherical coordinates associated with the basis set of vectors ( $\boldsymbol{e}_{\boldsymbol{r}},\boldsymbol{e}_{\boldsymbol{\theta } },\boldsymbol{e}_{\boldsymbol{\phi }}$ ), the general forms of these potentials are $\varphi =\varphi (r,\theta ,\phi ,t)$ and $\boldsymbol{\psi }=\psi _{r}(r,\theta ,\phi ,t)\,\boldsymbol{e}_{\boldsymbol{r}}+\psi _{\theta }(r,\theta ,\phi ,t)\,\boldsymbol{e}_{\theta }+\psi _{\phi }(r,\theta ,\phi ,t)\,\boldsymbol{e}_{\phi }$ . The Helmholtz decomposition is suitable for solving axisymmetric cases because the azimuthal invariance reduces the decomposition of the velocity field to only two unknowns, $\varphi (r,\theta ,\phi ,t)$ and $\psi _{\phi }(r,\theta ,\phi ,t)$ (Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a ). In the present problem, the general form of the scalar and vector potentials should be kept, leading to four unknown functions. Using the divergence-free property of the velocity field in (2.2), an alternative decomposition can be proposed to reduce the number of unknown functions, such as the poloidal–toroidal decomposition sometimes called the Mie representation of a vector field. Indeed, Mie (Reference Mie1908) was the first to introduce this decomposition in the investigation of Maxwell’s equations in a spherical geometry, after Lamb (Reference Lamb1881), who introduced the toroidal part only in hydrodynamics. The application of the poloidal–toroidal decomposition of a solenoidal field concerns mainly terrestrial magnetism (Elsasser Reference Elsasser1946; Bullard & Gellman Reference Bullard and Gellman1954; Backus Reference Backus1986). Chandrasekhar (Reference Chandrasekhar1961) adapted this representation to treat a convective-flow stability problem in a spherical cavity. The validity of the poloidal–toroidal decomposition for any solenoidal fields has been demonstrated by Chadwick & Trowbridge (Reference Chadwick and Trowbridge1967) in a bounded annular region, by Backus (Reference Backus1958) in a closed ball, and by Padmavati & Amaranath (Reference Padmavati and Amaranath2002) for unbounded regions. In the context of bubble physics, Prosperetti (Reference Prosperetti1977) proposed this decomposition for the vorticity field induced by bubble non-spherical oscillations in order to determine the equations of motion of the bubble interface accounting for non-spherical perturbations in viscous fluids. Therefore, we follow the approach of Prosperetti (Reference Prosperetti1977), and calculate the curl of both sides of (2.3), which results in

(2.5) \begin{equation} \frac{\partial \boldsymbol{\omega }_{1}^{\left(n,m\right)}}{\partial t}=\nu\, \Delta \boldsymbol{\omega }_{1}^{\left(n,m\right)} ,\end{equation}

where $\boldsymbol{\omega }_{1}^{(n,m)}={\boldsymbol \nabla }\times \boldsymbol{v}_{1}^{(n,m)}$ is called the vorticity of the velocity field $\boldsymbol{v}_{1}^{(n,m)}$ . We then apply the poloidal–toroidal decomposition of the vorticity (Backus Reference Backus1986)

(2.6) \begin{align} \boldsymbol{\omega }_{1}^{\left(n,m\right)}=\boldsymbol{P}_{1}^{\left(n,m\right)}+\boldsymbol{T}_{1}^{\left(n,m\right)} ,\end{align}

where, in view of the equation of the bubble surface (2.1), the poloidal $\boldsymbol{P}_{1}^{(n,m)}$ and toroidal $\boldsymbol{T}_{1}^{(n,m)}$ fields are written by

(2.7) \begin{align} \boldsymbol{P}_{1}^{\left(n,m\right)}=e^{-i\omega t}\,{\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\,P_{nm}(r)\,Y_{n}^{m}(\theta ,\phi )\right] ,\end{align}
(2.8) \begin{align} \boldsymbol{T}_{1}^{\left(n,m\right)}=e^{-i\omega t}\,{\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\,T_{nm}(r)\,Y_{n}^{m}(\theta ,\phi )\right] ,\quad \end{align}

with $\boldsymbol{e}_{r}=\boldsymbol{r}/r$ being the unit radial vector.

Substitution of (2.6) into (2.5) yields

(2.9) \begin{equation} \frac{\partial \boldsymbol{P}_{1}^{\left(n,m\right)}}{\partial t}+\frac{\partial \boldsymbol{T}_{1}^{\left(n,m\right)}}{\partial t}=\nu\, \Delta \boldsymbol{P}_{1}^{\left(n,m\right)}+\nu\, \Delta \boldsymbol{T}_{1}^{\left(n,m\right)} .\end{equation}

In view of the orthogonality of the poloidal and toroidal fields, (2.9) can be split into two equations:

(2.10) \begin{align} \frac{\partial \boldsymbol{P}_{1}^{\left(n,m\right)}}{\partial t}&=\nu\, \Delta \boldsymbol{P}_{1}^{\left(n,m\right)} ,\end{align}
(2.11) \begin{align} \frac{\partial \boldsymbol{T}_{1}^{\left(n,m\right)}}{\partial t}&=\nu\, \Delta \boldsymbol{T}_{1}^{\left(n,m\right)} .\end{align}

Substituting (2.7) and (2.8) into (2.10) and (2.11), and using the identity

(2.12) \begin{align} \Delta \left[\boldsymbol{e}_{r}\,F(r)\,Y_{n}^{m}(\theta ,\phi )\right]&=\boldsymbol{e}_{r}\,Y_{n}^{m}(\theta ,\phi )\left[\frac{{\rm d}^{2}F(r)}{{\rm d}r^{2}}-\frac{n\left(n+1\right)}{r^{2}}F(r)\right]\nonumber\\ &\quad +2\,{\boldsymbol \nabla }\left[\frac{F(r)}{r}Y_{n}^{m}(\theta ,\phi )\right] ,\end{align}

where $F(r)$ is an arbitrary function, one obtains

(2.13) \begin{align} &-i\omega\, {\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\,P_{nm}(r)\,Y_{n}^{m}(\theta ,\phi )\right]=\nu\, {\boldsymbol \nabla }\times {\boldsymbol \nabla }\nonumber\\ &\quad \times \left\{\boldsymbol{e}_{r}\,Y_{n}^{m}(\theta ,\phi )\left[\frac{{\rm d}^{2}P_{nm}(r)}{{\rm d}r^{2}}-\frac{n\left(n+1\right)}{r^{2}}P_{nm}(r)\right]\right\} ,\end{align}
(2.14) \begin{align} -i\omega\, {\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\,T_{nm}(r)\,Y_{n}^{m}(\theta ,\phi )\right]=\nu\, {\boldsymbol \nabla }\times \left\{\boldsymbol{e}_{r}\,Y_{n}^{m}(\theta ,\phi )\!\left[\frac{{\rm d}^{2}T_{nm}(r)}{{\rm d}r^{2}}-\frac{n\left(n+1\right)}{r^{2}}T_{nm}(r)\right]\right\} . \nonumber\\[3pt] \end{align}

It follows from (2.13) and (2.14) that $P_{nm}(r)$ and $T_{nm}(r)$ obey the equations

(2.15) \begin{align} \left[\frac{{\rm d}^{2}}{{\rm d}r^{2}}-\frac{n\left(n+1\right)}{r^{2}}+k_{v}^{2}\right]P_{nm}(r)=0 ,\end{align}
(2.16) \begin{align} \left[\frac{{\rm d}^{2}}{{\rm d}r^{2}}-\frac{n\left(n+1\right)}{r^{2}}+k_{v}^{2}\right]T_{nm}(r)=0 ,\end{align}

where $k_{v}=(1+i)/\delta _{v}$ is the viscous wavenumber, and $\delta _{v}=\sqrt{2\nu /\omega }$ is the viscous penetration depth.

It is easy to check that both (2.15) and (2.16) are transformed to the Riccati–Bessel equation (Abramowitz & Stegun Reference Abramowitz and Stegun1972) by multiplying by $r^{2}$ . Solutions to the above equation are given by $k_{v}r\,z_{n}(k_{v}r)$ , where $z_{n}$ is the spherical Bessel function of the first or second kind, or the spherical Hankel function of the first or second kind. Since we are looking for solutions in the form of an outgoing wave, we get

(2.17) \begin{align} P_{nm}(r)=a_{nm}k_{v}r\,h_{n}^{\left(1\right)}(k_{v}r) ,\end{align}
(2.18) \begin{align} T_{nm}(r)=b_{nm}k_{v}r\,h_{n}^{\left(1\right)}(k_{v}r) ,\end{align}

where $a_{nm}$ and $b_{nm}$ are constants, called the linear scattering coefficients, that are determined by boundary conditions at the bubble surface, and $h_{n}^{(1)}$ is the spherical Hankel function of the first kind.

Substituting (2.7) and (2.8) into (2.6), and using (2.17) and (2.18), one obtains

(2.19) \begin{align} \boldsymbol{\omega }_{1}^{\left(n,m\right)}&=e^{-i\omega t}\left\{{\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\,a_{nm}k_{v}r\,h_{n}^{\left(1\right)}(k_{v}r)\,Y_{n}^{m}(\theta ,\phi )\right]\right. \nonumber\\ &\quad \left. +{\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\,b_{nm}k_{v}r\,h_{n}^{\left(1\right)}(k_{v}r)\,Y_{n}^{m}(\theta ,\phi )\right]\right\}. \end{align}

It follows from (2.19) that $\boldsymbol{v}_{1}^{(n,m)}$ can be written as

(2.20) \begin{align} \boldsymbol{v}_{1}^{\left(n,m\right)}&=e^{-i\omega t}\left\{{\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\,a_{nm}k_{v}r\,h_{n}^{\left(1\right)}(k_{v}r)\,Y_{n}^{m}(\theta ,\phi )\right]\right. \nonumber\\ &\quad \left.+\boldsymbol{e}_{r}\,b_{nm}k_{v}r\,h_{n}^{\left(1\right)}(k_{v}r)\,Y_{n}^{m}(\theta ,\phi )-{\boldsymbol \nabla }\varphi \right\} ,\end{align}

where a scalar function $\varphi$ is introduced in order to satisfy (2.2).

Substitution of (2.20) into (2.2) yields

(2.21) \begin{equation} \Delta \varphi ={\boldsymbol \nabla }\cdot \left[\boldsymbol{e}_{r}\,b_{nm}k_{v}r\,h_{n}^{\left(1\right)}(k_{v}r)\,Y_{n}^{m}(\theta ,\phi )\right] .\end{equation}

Equation (2.21) is the inhomogeneous Laplace equation, which is also called Poisson’s equation. Since we are looking for a solution for $r\gt R_{0}$ , the solution to (2.21) can be represented by

(2.22) \begin{equation} \varphi \left(r,\theta ,\phi \right)=c_{nm}r^{-\left(n+1\right)}\,Y_{n}^{m}(\theta ,\phi )+\varphi _{nm}(r)\,Y_{n}^{m}(\theta ,\phi ) ,\end{equation}

where the first term on the right-hand side is the solution to the Laplace equation ( $\Delta \varphi =0$ ), $c_{nm}$ being a constant coefficient, and the second term is a particular solution to Poisson’s equation (2.21), $\varphi _{nm}(r)$ being a function to be found.

To find $\varphi _{nm}(r)$ , we substitute (2.22) into (2.21), and calculate the right-hand side of (2.21). As a result, we obtain the equation

(2.23) \begin{align} \left[\frac{{\rm d}^{2}}{{\rm d}x^{2}}+\frac{2}{x}\frac{\rm d}{{\rm d}x}-\frac{n\left(n+1\right)}{x^{2}}\right]\varphi _{nm}(x)&=\frac{b_{nm}}{k_{v}}\left[3\,h_{n}^{\left(1\right)}(x)+x\,h_{n}^{\left(1\right)\prime}(x)\right]\nonumber\\&=\frac{b_{nm}}{k_{v}}\left[\left(n+3\right)\,h_{n}^{\left(1\right)}(x)-x\,h_{n+1}^{\left(1\right)}(x)\right] ,\end{align}

where $x=k_{v}r$ , and the prime denotes the derivative with respect to an argument in brackets.

To solve (2.23), we use the identities (Abramowitz & Stegun Reference Abramowitz and Stegun1972)

(2.24) \begin{align} \frac{n}{x}h_{n}^{\left(1\right)}(x)-h_{n}^{\left(1\right)\prime}(x)=h_{n+1}^{\left(1\right)}(x) ,\quad\end{align}
(2.25) \begin{align} \frac{n+1}{x}h_{n}^{\left(1\right)}(x)+h_{n}^{\left(1\right)\prime}(x)=h_{n-1}^{\left(1\right)}(x) .\end{align}

With the help of (2.24) and (2.25), we obtain

(2.26) \begin{equation} \left[\frac{{\rm d}^{2}}{{\rm d}x^{2}}+\frac{2}{x}\frac{\rm d}{{\rm d}x}-\frac{n\left(n+1\right)}{x^{2}}\right]\!\left[\alpha\, h_{n}^{\left(1\right)}(x)+\beta x\,h_{n+1}^{\left(1\right)}(x)\right]=(2\beta -\alpha )\,h_{n}^{\left(1\right)}(x)-\beta x\,h_{n+1}^{\left(1\right)}(x) .\end{equation}

Comparison of (2.23) with (2.26) reveals that the solution to (2.23) is obtained by setting $\alpha =-(n+1)b_{nm}/k_{v}$ and $\beta =b_{nm}/k_{v}$ :

(2.27) \begin{equation} \varphi _{nm}\left(x\right)=\frac{b_{nm}}{k_{v}}\left[x\,h_{n+1}^{\left(1\right)}(x)-(n+1)\,h_{n}^{\left(1\right)}(x)\right] .\end{equation}

It follows that the general solution to (2.21) is given by

(2.28) \begin{equation} \varphi \left(r,\theta ,\phi \right)=\left\{c_{nm}r^{-\left(n+1\right)}+\frac{b_{nm}}{k_{v}}\left[k_{v}r\,h_{n+1}^{\left(1\right)}\left(k_{v}r\right)-\left(n+1\right)h_{n}^{\left(1\right)}\left(k_{v}r\right)\right]\right\}Y_{n}^{m}\left(\theta ,\phi \right) .\end{equation}

The coefficients $a_{nm} , b_{nm}, c_{nm}$ are calculated in Appendix A. In the process of this calculation, the following boundary conditions at the bubble surface have been used: the condition that the normal component of the first-order liquid velocity is equal to the normal component of the bubble surface velocity, and the condition that the tangential stress generated by the liquid motion vanishes at the bubble surface because the gas viscosity is much lower than the liquid viscosity. By using the results obtained in Appendix A, $\boldsymbol{v}_{1}^{(n,m)}$ is represented by

(2.29) \begin{align} \boldsymbol{v}_{1}^{\left(n,m\right)}=v_{1r}^{\left(n,m\right)}\boldsymbol{e}_{r}+v_{1\theta }^{\left(n,m\right)}\boldsymbol{e}_{\theta }+v_{1\phi }^{\left(n,m\right)}\boldsymbol{e}_{\phi } ,\qquad\quad\end{align}
(2.30) \begin{align} v_{1r}^{\left(n,m\right)}=e^{-i\omega t}\,s^{\left(n,m\right)}\,V_{n}\left(r\right)\,Y_{n}^{m}\left(\theta ,\phi \right) ,\qquad\quad\end{align}
(2.31) \begin{align} v_{1\theta }^{\left(n,m\right)}=e^{-i\omega t}\,s^{\left(n,m\right)}\,W_{n}\left(r\right)\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta },\quad n\geq 1 ,\quad\end{align}
(2.32) \begin{align} v_{1\phi }^{\left(n,m\right)}=e^{-i\omega t}\,s^{\left(n,m\right)}\,\frac{W_{n}\left(r\right)}{\sin \theta }\,\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \phi },\quad n,m\geq 1, \end{align}

where the functions $V_{n}(r)$ and $W_{n}(r)$ are calculated by

(2.33) \begin{align} V_{n}\left(r\right)=\frac{\left(n+1\right)\alpha _{n}}{r^{n+2}}+\frac{n\left(n+1\right)\beta _{n}h_{n}^{\left(1\right)}\left(k_{v}r\right)}{k_{v}r} ,\qquad\quad\end{align}
(2.34) \begin{align} W_{n}\left(r\right)=-\frac{\alpha _{n}}{r^{n+2}}+\frac{\beta _{n}}{k_{v}r}\left[\left(n+1\right)h_{n}^{\left(1\right)}\left(k_{v}r\right)-k_{v}rh_{n+1}^{\left(1\right)}\left(k_{v}r\right)\right] ,\end{align}

and the coefficients $\alpha _{n}$ and $\beta _{n}$ are given by

(2.35) \begin{align} \alpha _{n}=\frac{i\omega R_{0}^{n+2}\left[\left(2-n-n^{2}\right)h_{n}^{\left(1\right)}\left(\overline{x}\right)-\overline{x}^{2}h_{n}^{\left(1\right)\prime\prime}\left(\overline{x}\right)\right]}{\left(n+1\right)\left[\overline{x}^{2}h_{n}^{\left(1\right)\prime\prime}\left(\overline{x}\right)-\left(n^{2}+3n+2\right)h_{n}^{\left(1\right)}\left(\overline{x}\right)\right]},\quad n\geq 0 ,\end{align}
(2.36) \begin{align} \beta _{n}=\frac{2i\left(n+2\right)\overline{x}\omega }{\left(n+1\right)\left[\overline{x}^{2}h_{n}^{\left(1\right)\prime\prime}\left(\overline{x}\right)-\left(n^{2}+3n+2\right)h_{n}^{\left(1\right)}\left(\overline{x}\right)\right]},\quad n\geq 1 ,\end{align}

where $\overline{x}=k_{v}R_{0}$ . Note that (2.33)–(2.36) follow from (A8), (A9), (A22) and (A23) by setting $c_{nm}=s^{(n,m)}\alpha _{n}$ and $b_{nm}=s^{(n,m)}\beta _{n}$ .

The first-order velocity generated by the mode $(n,-m)$ is obtained by replacing $m$ with $-m$ in (2.29)–(2.32):

(2.37) \begin{align} \boldsymbol{v}_{1}^{\left(n,-m\right)}=v_{1r}^{\left(n,-m\right)}\boldsymbol{e}_{r}+v_{1\theta }^{\left(n,-m\right)}\boldsymbol{e}_{\theta }+v_{1\phi }^{\left(n,-m\right)}\boldsymbol{e}_{\phi } ,\qquad\end{align}
(2.38) \begin{align} v_{1r}^{\left(n,-m\right)}=e^{-i\omega t}s^{\left(n,-m\right)}V_{n}\left(r\right)Y_{n}^{-m}\left(\theta ,\phi \right) ,\qquad\quad\end{align}
(2.39) \begin{align} v_{1\theta }^{\left(n,-m\right)}=e^{-i\omega t}s^{\left(n,-m\right)}W_{n}\left(r\right)\frac{\partial Y_{n}^{-m}\left(\theta ,\phi \right)}{\partial \theta },\quad n\geq 1 ,\,\,\end{align}
(2.40) \begin{align} v_{1\phi }^{\left(n,-m\right)}=e^{-i\omega t}s^{\left(n,-m\right)}\frac{W_{n}\left(r\right)}{\sin \theta }\frac{\partial Y_{n}^{-m}\left(\theta ,\phi \right)}{\partial \phi },\quad n,m\geq 1 .\end{align}

With the help of (D1) and (D3) from Appendix D, (2.38)–(2.40) are represented by

(2.41) \begin{align} v_{1r}^{\left(n,-m\right)}=e^{-i\omega t}s^{\left(n,-m\right)}V_{n}\left(r\right)\left(-1\right)^{m}Y_{n}^{m\ast }\left(\theta ,\phi \right) ,\qquad\quad\end{align}
(2.42) \begin{align} v_{1\theta }^{\left(n,-m\right)}=e^{-i\omega t}s^{\left(n,-m\right)}W_{n}\left(r\right)\left(-1\right)^{m}\frac{\partial Y_{n}^{m\ast }\left(\theta ,\phi \right)}{\partial \theta },\quad n\geq 1 ,\quad\end{align}
(2.43) \begin{align} v_{1\phi }^{\left(n,-m\right)}=e^{-i\omega t}s^{\left(n,-m\right)}\frac{W_{n}\left(r\right)}{\sin \theta }\left(-1\right)^{m}\frac{\partial Y_{n}^{m\ast }\left(\theta ,\phi \right)}{\partial \phi },\quad n,m\geq 1 ,\end{align}

where the asterisk denotes the complex conjugate.

2.2. Solutions of the equations of acoustic streaming

The equations of acoustic streaming are given by (Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a )

(2.44) \begin{align} {\boldsymbol \nabla }\cdot \boldsymbol{v}_{E}=0 ,\qquad\qquad\end{align}
(2.45) \begin{align} \Delta {\boldsymbol \nabla }\times \boldsymbol{v}_{E}=\frac{1}{\nu }{\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}\right\rangle ,\end{align}

where $\boldsymbol{v}_{E}$ is the Eulerian streaming velocity, i.e. the time-averaged second-order liquid velocity, and $\langle\ \rangle$ means the time average.

Substitution of (2.4) into (2.45) yields

(2.46) \begin{align} &\Delta {\boldsymbol \nabla }\times \boldsymbol{v}_{E}=\frac{1}{\nu }{\boldsymbol \nabla }\nonumber\\&\quad \times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}+\boldsymbol{v}_{1}^{\left(n,-m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}+\boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}+\boldsymbol{v}_{1}^{\left(n,-m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle .\end{align}

This equation can be divided into three equations:

(2.47) \begin{align} \Delta {\boldsymbol \nabla }\times \boldsymbol{v}_{E}^{\left(n,m\right)}=\frac{1}{\nu }{\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle ,\qquad\qquad\end{align}
(2.48) \begin{align} \Delta {\boldsymbol \nabla }\times \boldsymbol{v}_{E}^{\left(n,-m\right)}=\frac{1}{\nu }{\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,-m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}\right\rangle ,\qquad\quad\end{align}
(2.49) \begin{align} \Delta {\boldsymbol \nabla }\times \boldsymbol{v}_{E}^{\left(\times \right)}=\frac{1}{\nu }{\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}+\boldsymbol{v}_{1}^{\left(n,-m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle , \end{align}

where (2.47) describes the acoustic streaming generated by the mode $(n,m)$ alone, (2.48) describes the acoustic streaming generated by the mode $(n,-m)$ alone, and (2.49) describes the acoustic streaming due to the interaction of the above modes (streaming due to cross-terms).

Let us first consider (2.47). Since ${\boldsymbol \nabla }\times \boldsymbol{v}_{E}^{(n,m)}$ is a divergence-free vector field, it can be represented by a poloidal–toroidal decomposition (Backus Reference Backus1986)

(2.50) \begin{align} {\boldsymbol \nabla }\times \boldsymbol{v}_{E}^{\left(n,m\right)}&={\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}P_{kl}^{\left(n,m\right)}\left(r\right)Y_{k}^{l}\left(\theta ,\phi \right)\right]\nonumber\\&\quad +{\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}T_{kl}^{\left(n,m\right)}\left(r\right)Y_{k}^{l}\left(\theta ,\phi \right)\right]\nonumber\\ &=\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\left\{\boldsymbol{e}_{r}\frac{k\left(k+1\right)P_{kl}^{\left(n,m\right)}\left(r\right)}{r^{2}}Y_{k}^{l}\left(\theta ,\phi \right)\right.\nonumber\\&\quad +\frac{P_{kl}^{\left(n,m\right)\prime}\left(r\right)}{r}\left[\boldsymbol{e}_{\theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }+\frac{\boldsymbol{e}_{\phi }}{\sin \theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \phi }\right]\nonumber\\ &\left.\quad{}+\frac{T_{kl}^{\left(n,m\right)}\left(r\right)}{r}\left[\frac{\boldsymbol{e}_{\theta }}{\sin \theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \phi }-\boldsymbol{e}_{\phi }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }\right]\right\}. \end{align}

Recall that the prime in the superscript denotes the derivative with respect to an argument in brackets. Substitution of (2.50) into (2.47) and the calculation of $\Delta {\boldsymbol \nabla }\times \boldsymbol{v}_{E}^{(n,m)}$ yield

(2.51) \begin{align} &\Delta {\boldsymbol \nabla }\times \boldsymbol{v}_{E}^{\left(n,m\right)}=\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\left\{\boldsymbol{e}_{r}\frac{k\left(k+1\right)}{r^{2}}\left[P_{kl}^{\left(n,m\right)\prime\prime}\left(r\right)-\frac{k\left(k+1\right)P_{kl}^{\left(n,m\right)}\left(r\right)}{r^{2}}\right]\right.Y_{k}^{l}\left(\theta ,\phi \right)\nonumber\\ &\quad{}+\frac{1}{r}\left[P_{kl}^{\left(n,m\right)\prime\prime\prime}\left(r\right)-\frac{k\left(k+1\right)P_{kl}^{\left(n,m\right)/}\left(r\right)}{r^{2}}+\frac{2k\left(k+1\right)P_{kl}^{\left(n,m\right)}\left(r\right)}{r^{3}}\right]\nonumber\\ &\qquad{}\times \left[\boldsymbol{e}_{\theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }+\frac{\boldsymbol{e}_{\phi }}{\sin \theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \phi }\right]\nonumber\\ &\quad{}+\left.\frac{1}{r}\left[T_{kl}^{\left(n,m\right)\prime\prime}\left(r\right)-\frac{k\left(k+1\right)T_{kl}^{\left(n,m\right)}\left(r\right)}{r^{2}}\right]\left[\frac{\boldsymbol{e}_{\theta }}{\sin \theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \phi }-\boldsymbol{e}_{\phi }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }\right]\right\}\nonumber\\ &=\frac{1}{\nu }{\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle. \end{align}

Equating the $r$ -components of both sides of (2.51), one obtains

(2.52) \begin{align} &\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\frac{k\left(k+1\right)}{r^{2}}\left[P_{kl}^{\left(n,m\right)\prime\prime}\left(r\right)-\frac{k\left(k+1\right)P_{kl}^{\left(n,m\right)}\left(r\right)}{r^{2}}\right]Y_{k}^{l}\left(\theta ,\phi \right)\nonumber\\&\quad=\frac{1}{\nu }\boldsymbol{e}_{r}\cdot \left[{\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle \right] .\end{align}

Calculating the curl of both sides of (2.51) and taking the $r$ -component, one finds

(2.53) \begin{align} &\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\frac{k\left(k+1\right)}{r^{2}}\left[T_{kl}^{\left(n,m\right)\prime\prime}\left(r\right)-\frac{k\left(k+1\right)T_{kl}^{\left(n,m\right)}\left(r\right)}{r^{2}}\right]Y_{k}^{l}\left(\theta ,\phi \right)\nonumber\\&\quad=\frac{1}{\nu }\boldsymbol{e}_{r}\cdot \left[{\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle \right] .\end{align}

Equations (2.52) and (2.53) make it possible to calculate $P_{kl}^{(n,m)}(r)$ and $T_{kl}^{(n,m)}(r)$ , and hence $\boldsymbol{v}_{E}^{(n,m)}$ . In view of the cumbersome nature of these calculations, they are performed in Appendix B. As a result, we obtain the components of the Eulerian streaming velocity in the following form:

(2.54) \begin{align} v_{Er}^{\left(n,m\right)}\left(r,\theta \right)=\frac{1}{2\sqrt{\pi}}\sum _{k=1}^{\infty }\sqrt{2k+1}\left[T_{k0}^{\left(n,m\right)}\left(r\right)+\Phi _{k0}^{\left(n,m\right)/}\left(r\right)\right]P_{k}\left(\cos \theta \right) , \end{align}
(2.55) \begin{align} v_{E\theta }^{\left(n,m\right)}\left(r,\theta \right)=\frac{1}{2\sqrt{\pi}\,r}\sum _{k=1}^{\infty }\sqrt{2k+1}\,\Phi _{k0}^{\left(n,m\right)}\left(r\right)P_{k}^{1}\left(\cos \theta \right) ,\qquad\end{align}
(2.56) \begin{align} v_{E\phi }^{\left(n,m\right)}\left(r,\theta \right)=-\frac{1}{2\sqrt{\pi}\,r}\sum _{k=1}^{\infty }\sqrt{2k+1}\,P_{k0}^{\left(n,m\right)}\left(r\right)P_{k}^{1}\left(\cos \theta \right) ,\qquad\end{align}

where $P_{k}$ is the Legendre polynomial of degree $k$ , $P_{k}^{1}$ is the associated Legendre polynomial of first order and degree $k$ , and the functions $P_{k0}^{(n,m)}(r) , T_{k0}^{(n,m)}(r) , \Phi _{k0}^{(n,m)}(r), \Phi _{k0}^{(n,m)\prime}(r)$ are calculated by (B18), (B24), (B33) and (B45).

In the process of the calculation of $\boldsymbol{v}_{E}^{(n,m)}$ in Appendix B, we also calculate the Stokes drift velocity $\boldsymbol{v}_{S}^{(n,m)}$ (Longuet-Higgins Reference Longuet-Higgins1998), whose components are given by

(2.57) \begin{align} v_{Sr}^{\left(n,m\right)}\left(r,\theta \right)&=\frac{|s^{\left(n,m\right)}|^{2}}{4\sqrt{\pi}\,\omega }\sum _{k=1}^{\infty }\sqrt{2k+1}\,A_{k0}^{nmnm}\,\nonumber\\&\quad\times \mathrm{Re}\left\{i\,V_{n}^{\ast }\left(r\right)\left[\frac{2n\left(n+1\right)-k\left(k+1\right)}{2r}W_{n}\left(r\right)-V_{n}^{\prime}\left(r\right)\right]\right\}P_{k}\left(\cos \theta \right) ,\end{align}
(2.58) \begin{align} v_{S\theta }^{\left(n,m\right)}\left(r,\theta \right)=\frac{|s^{\left(n,m\right)}|^{2}}{8\sqrt{\pi}\,\omega }\,\mathrm{Re}\!\left\{i\,V_{n}^{\ast }\left(r\right)\!\left[\frac{W_{n}\left(r\right)}{r}-W_{n}^{\prime}\left(r\right)\right]\right\}\!\sum _{k=1}^{\infty }\sqrt{2k+1}\,A_{k0}^{nmnm}\,P_{k}^{1}\!\left(\cos \theta \right) ,\end{align}
(2.59) \begin{align} v_{S\phi }^{\left(n,m\right)}\left(r,\theta \right)=\sum _{k=1}^{\infty }E_{k}^{\left(n,m\right)}\left(r\right)P_{k}^{1}\left(\cos \theta \right) ,\qquad\qquad\qquad\qquad\qquad\quad\end{align}

where Re means ‘the real part of’, the coefficients $A_{k0}^{nmnm}$ are calculated by (D20), and the function $E_{k}^{(n,m)}$ is defined by (B63). Note that (2.57) and (2.58) follow from (B54) and (B56) using $f_{nm}(r)=s^{(n,m)}V_{n}(r)$ and $g_{nm}(r)=s^{(n,m)}W_{n}(r)$ .

Knowing $\boldsymbol{v}_{E}^{(n,m)}$ and $\boldsymbol{v}_{S}^{(n,m)}$ , we can calculate the Lagrangian streaming velocity produced by the mode $(n,m)$ , which is defined by $\boldsymbol{v}_{L}^{(n,m)}=\boldsymbol{v}_{E}^{(n,m)}+\boldsymbol{v}_{S}^{(n,m)}$ .

Due to symmetry, solutions to (2.48) are expressed in terms of the solutions to (2.47) as follows:

(2.60) \begin{align} v_{Lr}^{\left(n,-m\right)}\left(r,\theta \right)=\frac{|s^{\left(n,-m\right)}|^{2}}{|s^{\left(n,m\right)}|^{2}}v_{Lr}^{\left(n,m\right)}\left(r,{\pi} -\theta \right) , \end{align}
(2.61) \begin{align} v_{L\theta }^{\left(n,-m\right)}\left(r,\theta \right)=-\frac{|s^{\left(n,-m\right)}|^{2}}{|s^{\left(n,m\right)}|^{2}}v_{L\theta }^{\left(n,m\right)}\left(r,{\pi} -\theta \right) ,\end{align}
(2.62) \begin{align} v_{L\phi }^{\left(n,-m\right)}\left(r,\theta \right)=-\frac{|s^{\left(n,-m\right)}|^{2}}{|s^{\left(n,m\right)}|^{2}}v_{L\phi }^{\left(n,m\right)}\left(r,{\pi} -\theta \right) .\end{align}

To solve (2.49), by analogy with the solution of (2.47), ${\boldsymbol \nabla }\times \boldsymbol{v}_{E}^{(\times )}$ is represented by a poloidal–toroidal decomposition,

(2.63) \begin{align} {\boldsymbol \nabla }\times \boldsymbol{v}_{E}^{\left(\times \right)}&=\mathrm{Re}\left\{{\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}P_{kl}^{\left(\times \right)}\left(r\right)Y_{k}^{l}\left(\theta ,\phi \right)\right]\right.\nonumber\\&\quad+\left. {\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}T_{kl}^{\left(\times \right)}\left(r\right)Y_{k}^{l}\left(\theta ,\phi \right)\right]\right\}\nonumber\\&=\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\left\{\boldsymbol{e}_{r}\frac{k\left(k+1\right)P_{kl}^{\left(\times \right)}\left(r\right)}{r^{2}}Y_{k}^{l}\left(\theta ,\phi \right)\right.\nonumber\\&\quad+\frac{P_{kl}^{\left(\times \right)\prime}\left(r\right)}{r}\left[\boldsymbol{e}_{\theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta } +\frac{\boldsymbol{e}_{\phi }}{\sin \theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \phi }\right]\nonumber\\&\quad\left.+\frac{T_{kl}^{\left(\times \right)}\left(r\right)}{r}\left[\frac{\boldsymbol{e}_{\theta }}{\sin \theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \phi }-\boldsymbol{e}_{\phi }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }\right]\right\}, \end{align}

where the functions $P_{kl}^{(\times )}(r)$ and $T_{kl}^{(\times )}(r)$ , by analogy with (2.52) and (2.53), obey the equations

(2.64) \begin{align} &\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\frac{k\left(k+1\right)}{r^{2}}\left[P_{kl}^{\left(\times \right)\prime\prime}\left(r\right)-\frac{k\left(k+1\right)P_{kl}^{\left(\times \right)}\left(r\right)}{r^{2}}\right]Y_{k}^{l}\left(\theta ,\phi \right)\nonumber\\&\quad=\frac{1}{\nu }\boldsymbol{e}_{r}\cdot \left[{\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}+\boldsymbol{v}_{1}^{\left(n,-m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle \right], \end{align}
(2.65) \begin{align} &\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\frac{k\left(k+1\right)}{r^{2}}\left[T_{kl}^{\left(\times \right)\prime\prime}\left(r\right)-\frac{k\left(k+1\right)T_{kl}^{\left(\times \right)}\left(r\right)}{r^{2}}\right]Y_{k}^{l}\left(\theta ,\phi \right)\nonumber\\&\quad=\frac{1}{\nu }\boldsymbol{e}_{r}\cdot \left[{\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}+\boldsymbol{v}_{1}^{\left(n,-m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle \right]. \end{align}

Equations (2.64) and (2.65) are solved in Appendix C. As a result, we obtain the components of $\boldsymbol{v}_{E}^{(\times )}$ in the form

(2.66) \begin{align} &\qquad\qquad v_{Er}^{\left(\times \right)}\left(r,\theta ,\phi \right)=\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\left[T_{kl}^{\left(\times \right)}\left(r\right)+\Phi _{kl}^{\left(\times \right)\prime}\left(r\right)\right]Y_{k}^{l}\left(\theta ,\phi \right) ,\end{align}
(2.67) \begin{align} &v_{E\theta }^{\left(\times \right)}\left(r,\theta ,\phi \right)=\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\frac{\Phi _{kl}^{\left(\times \right)}\left(r\right)}{2r} \nonumber\\&\quad\times \left[\sqrt{k\left(k+1\right)-l\left(l+1\right)}Y_{k}^{l+1}\left(\theta ,\phi \right)e^{-i\phi }-\sqrt{k\left(k+1\right)-l\left(l-1\right)}Y_{k}^{l-1}\left(\theta ,\phi \right)e^{i\phi }\right], \end{align}
(2.68) \begin{align} & v_{E\phi }^{\left(\times \right)}\left(r,\theta ,\phi \right)=-\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\frac{il\Phi _{kl}^{\left(\times \right)}\left(r\right)e^{i\phi }}{r}\sqrt{\frac{\left(2k+1\right)\left(k-l\right)!}{\left(k+l\right)!}}\nonumber\\& \quad\qquad\qquad\qquad \times \sum _{s=1}^{\left[\left(k-l+2\right)/2\right]}\sqrt{\frac{\left(2k-4s+3\right)\left(k+l-2s\right)!}{\left(k-l-2s+2\right)!}}Y_{k-2s+1}^{l-1}\left(\theta ,\phi \right), \end{align}

where the functions $T_{kl}^{(\times )}(r) , \Phi _{kl}^{(\times )}(r)$ and $\Phi _{kl}^{(\times )/}(r)$ are calculated by (C13), (C22) and (C33).

The components of the Stokes drift velocity produced by the interaction of the modes $(n,m)$ and $(n,-m)$ are given by

(2.69) \begin{equation} v_{Sr}^{\left(\times \right)}\left(r,\theta ,\phi \right)=\mathrm{Re}\sum _{k=1}^{\infty }S_{k}^{\left(\times \right)}\left(r\right)\sum _{l=-k}^{k}D_{kl}^{nmnm}Y_{k}^{l}\left(\theta ,\phi \right) ,\end{equation}
(2.70) \begin{align} v_{S\theta }^{\left(\times \right)}\left(r,\theta ,\phi \right)&=\frac{1}{2}\mathrm{Re}\left\{U_{nm}\left(r\right)\times \sum _{k=1}^{\infty }\sum _{l=-k}^{k}D_{kl}^{nmnm}\!\left[\!\sqrt{k\left(k+1\right)-l\left(l+1\right)}Y_{k}^{l+1}\!\left(\theta ,\phi \right)\!e^{-i\phi }\right.\right.\nonumber\\&\quad\left.\left.\vphantom{\sum _{k=1}^{\infty }\sum _{l=-k}^{k}D_{kl}^{nmnm}} -\sqrt{k\left(k+1\right)-l\left(l-1\right)}\,Y_{k}^{l-1}\left(\theta ,\phi \right)e^{i\phi }\right]\right\}, \end{align}
(2.71) \begin{align} v_{S\phi }^{\left(\times \right)}\left(r,\theta ,\phi \right)&=-\mathrm{Re}\left\{iU_{nm}\left(r\right)e^{i\phi }\sum _{k=1}^{\infty }\sum _{l=-k}^{k}l\sqrt{\frac{\left(2k+1\right)\left(k-l\right)!}{\left(k+l\right)!}}D_{kl}^{nmnm}\right.\nonumber\\&\quad\left.\times \sum _{s=1}^{\left[\left(k-l+2\right)/2\right]}\sqrt{\frac{\left(2k-4s+3\right)\left(k+l-2s\right)!}{\left(k-l-2s+2\right)!}}Y_{k-2s+1}^{l-1}\left(\theta ,\phi \right)\right\}, \end{align}

where the coefficients $D_{kl}^{nmnm}$ are calculated by (D27), and the functions $S_{k}^{(\times )}(r)$ and $U_{nm}(r)$ are defined by (C43) and (C44).

Knowing $\boldsymbol{v}_{E}^{(\times )}$ and $\boldsymbol{v}_{S}^{(\times )}$ , we can calculate the corresponding Lagrangian streaming velocity: $\boldsymbol{v}_{L}^{(\times )}=\boldsymbol{v}_{E}^{(\times )}+\boldsymbol{v}_{S}^{(\times )}$ .

Finally, summing the three components of the acoustic streaming, we get the total flow:

(2.72) \begin{equation} \boldsymbol{v}_{L}^{\left(\textit{total}\right)}\left(r,\theta ,\phi \right)=\boldsymbol{v}_{L}^{\left(n,m\right)}\left(r,\theta \right)+\boldsymbol{v}_{L}^{\left(n,-m\right)}\left(r,\theta \right)+\boldsymbol{v}_{L}^{\left(\times \right)}\left(r,\theta ,\phi \right) .\end{equation}

In the process of the calculations described in this section, we have applied the boundary conditions for the acoustic streaming at the bubble surface that assume that the normal velocity and the tangential stress produced by the Lagrangian streaming vanish at the equilibrium bubble surface.

It is worth noting that we do not provide expressions for the first- and second-order pressures because they are not used in the calculation of acoustic streaming. The above pressures are needed, for example, for the calculation of the acoustic radiation force on the bubble. However, this problem is beyond the scope of our paper. It requires individual consideration somewhere else because, in particular, the calculation of the second-order pressure is not a trivial mathematical problem.

3. Results and discussion

3.1. Three classes of spherical harmonics

Let us first consider how the asymmetric shape deformations are described mathematically. In the general case, the shape of a deformed three-dimensional body can always be decomposed over the set of orthonormal spherical harmonics $Y_{n}^{m}(\theta ,\phi )$ of degree $n$ and order $m$ . These two indexes are related to the spatial evolution of the spherical harmonics along the spherical angular coordinates, the colatitude $\theta \in [0,{\pi} ]$ and the longitude $\phi \in [0,2{\pi} ]$ . For a given degree $n$ , three classes of spherical harmonics are usually considered: the zonal harmonics when $m=0\lt n$ , the tesseral harmonics when $0\lt m\lt n$ , and the sectoral harmonics when $m=n$ . The three classes of spherical harmonics are represented in figure 2 for degree $n=5$ .

Figure 2.

Numerical shape deformations of a bubble experiencing (a) the zonal harmonic $Y_{5}^{0}$ , (b) the tesseral harmonic $Y_{5}^{3}$ , and (c) the sectoral harmonic $Y_{5}^{5}$ . Every isometric view of the spherical harmonics is displayed in top (upper) and side (lower) views.

In figure 2(a), the bubble undergoes the zonal harmonic $Y_{5}^{0}$ . The bubble interface is axisymmetric, the bubble contour looks spherical from the top view, and the shape from the side view is a Legendre polynomial ( $m=0$ ). Figure 2(b) illustrates the case of the tesseral harmonic $Y_{5}^{3}$ . The bubble interface is complex to describe at first, as angular deviations from the sphere appear along both the elevation and azimuthal directions. A closer look reveals that the bubble interface possesses $n-m=2$ nodal lines along the elevation (see figure 2(b), side view), and $m=3$ meridian nodal lines (see figure 2(b), top view). Along the azimuthal direction, and equivalently along the elevation, two successive extrema are out of phase from both side of the nodal lines. Figure 2(c) displays a bubble experiencing the sectoral harmonic $Y_{5}^{5}$ . The top view contour of a sectoral harmonic of degree $n$ corresponds to an $n$ -lobe deformation and is easily recognisable. The side view of a sectoral harmonic is close to the spherical shape. In fact, the side view contour resembles a bell-shape function, with a higher amplitude at the equator $\theta ={\pi} /2$ , and a decreasing oscillation amplitude when reaching the poles. It is worth mentioning that the so-called top and side views are ruled experimentally by the existence of a preferential direction for the triggering of a given shape oscillation at the bubble interface. For instance, in the case of a substrate-attached deformable body, the normal to the substrate is always the preferential direction for the triggering of asymmetric oscillations, as observed in Fauconnier et al. (Reference Fauconnier, Mauger, Bera and Inserra2022) for wall-attached bubbles, and in Chang et al. (Reference Chang, Bostwick, Daniel and Steen2015) for sessile drops.

In the present derivation of the bubble-induced microstreaming, only self-interacting asymmetric oscillations are considered. We recall here that the main mechanism for generating shape oscillations at the bubble interface is the Faraday instability occurring at the subharmonic of the driving frequency. Therefore, interactions of shape oscillations with the radial (spherical) oscillations of the bubble cannot lead to fluid mean flows, as the radial mode oscillates at the driving frequency. In addition, near the instability threshold of a given shape oscillation, we also assume that a single instability is triggered without considering the energy transfer to secondary shape oscillations, as theoretically discussed by Shaw (Reference Shaw2006) and experimentally observed by Guédra et al. (Reference Guédra, Inserra, Mauger and Gilles2016). An analysis of the intermodal interactions and their impact on the microstreaming pattern is performed in Regnault et al. (Reference Regnault, Mauger, Blanc-Benon, Doinikov and Inserra2021). The analysis of self-interacting axisymmetric modes is well-known theoretically (Inserra et al. Reference Inserra, Regnault, Cleve, Mauger and Doinikov2020a ) and has been observed from a side view in a plane containing the bubble symmetry $z$ axis using acoustically trapped bubbles, being far from any boundary (Cleve et al. Reference Cleve, Guedra, Mauger, Inserra and Blanc-Benon2019). The case of self-interacting asymmetric oscillations was considered by Fauconnier et al. (Reference Fauconnier, Mauger, Bera and Inserra2022). The authors investigated experimentally the top view microstreaming flow induced by a wall-attached bubble experiencing various asymmetric oscillations. The predominant shape oscillations were selected using a spectral analysis of the displacement of the bubble interface, and the modal content of the bubble could safely be associated with the observed fluid pattern. For some experimental cases, the secondary excited shape oscillations were so weak that the microstreaming pattern could be confidently related to the self-interaction of the main triggered surface oscillation. Such patterns resulting from self-interacting asymmetric oscillations are reproduced in figure 3 (adapted from Fauconnier et al. Reference Fauconnier, Mauger, Bera and Inserra2022). Figure 3(a) shows the top view microstreaming pattern induced by the self-interacting zonal harmonic $Y_{4}^{0}$ . The pattern is exclusively radial, since the oscillating bubble interface exhibits no azimuthal dependence (figure 2 a, top view). Figure 3(b) represents the top view microstreaming pattern induced by the self-interacting tesseral harmonic $Y_{4}^{2}$ . The pattern is characterised by $4m=8$ lobes, arranged by pairs around the $z$ axis orthogonal to the image, each pair being located between two displacement nodes of the bubble interface (the location of the four meridian lines is clearly visible through the light shades inside the bubble). Figure 3(c) shows the top view microstreaming pattern induced by the self-interacting sectoral harmonic $Y_{3}^{3}$ . The microstreaming pattern displays $4m=12$ lobes that are also arranged in pairs. Because sectoral modes are devoid of nodal parallels and exhibit an azimuthal shape that corresponds to the $\cos (m\phi )$ function along the equator, $2m=6$ extrema of displacements are observed on the bubble interface. Therefore, two vortices are expected between two successive nodal meridians, leading to the observed 12-lobe pattern.

Figure 3.

Microstreaming patterns associated with the three classes of spherical harmonics, i.e. (a) zonal, (b) tesseral and (c) sectoral, experienced by a wall-attached bubble and observed with microscope in a top view configuration. Images are reproduced from Fauconnier et al. (Reference Fauconnier, Mauger, Bera and Inserra2022).

For all these experimental microstreaming patterns, the top view configuration avoids the analysis of the full three-dimensional signature of fluid flows induced by asymmetric oscillations, as the seeded particles tracking the fluid motion can escape from the focal plane of observation. The dependence in elevation for the above-mentioned pattern is therefore not ensured and can only be guessed using physical arguments. The three-dimensional nature of the fluid flow was investigated by Marin et al. (Reference Marin, Rossi, Rallabandi, Wang, Hilgenfeldt and Kähler2015) and Bolanos-Jimenez et al. (Reference Bolanos-Jimenez, Rossi, Fernandez Rivas, Kähler and Marin2017) using an APTV technique. However, the analysis was performed on a wall-attached bubble exhibiting an axisymmetric flow around the perpendicular to the wall. In the following, we illustrate the three-dimensional nature of the flow surrounding asymmetric bubbles using the present modelling.

3.2. Signature for axisymmetric, zonal ( $m=0$ ) harmonics

Axisymmetric shape oscillations are the widely investigated cases because of the ease of triggering them in various experimental configurations. The lack of the dependence of the zonal harmonics in the azimuthal direction significantly facilitates their mathematical analysis, as already discussed in the Introduction. Therefore, we first analyse the case of a zonal (axisymmetric) harmonic, namely the case $n=3,\ m=0$ . The bubble equilibrium radius is set to 73.8 µm, as in the experimental case (figure 3 c) displaying a shape oscillation of degree $n=3$ . This value is close to the resonant radius of the $n$ th shape oscillation, as derived by Lamb (Reference Lamb1916):

(2.73) \begin{equation}R_{res}^{n}=\sqrt[3]{\frac{\left(n-1\right)\left(n+1\right)\left(n+2\right)\sigma }{\rho \omega _{d}^{2}/4}},\end{equation}

Figure 4.

Presentation of the microstreaming pattern observed for a bubble (equilibrium radius 73.8 µm) exhibiting zonal oscillations on the mode $(n=3, m=0$ ): (a) numerical contours of the bubble interface in the top, side and three-dimensional views; (b) top, (c) side and (d) isometric views of the resulting flow.

where $\sigma$ is the surface tension, $\rho$ is the liquid density, and $\omega _{d}=2{\pi} f_{d}$ is the angular frequency of acoustic driving. The value of $\omega _{d}$ differs from the frequency of the shape oscillation when the surface deformation is induced by a parametric instability. It should hence be noted that in all the following numerical simulations of microstreaming induced by a shape oscillation, the angular frequency $\omega$ used in the theoretical derivation is assumed equal to half the driving frequency: here, $\omega =\omega _{d}/2$ .

For the axisymmetric shape oscillations of degree $n=3$ , the resonant radius is $R_{res}^{n}=65.9$ µm for an air bubble in water with $\rho =1000$ kg m−3, $\sigma =0.0727$ N m $^{-1}$ at $f_{d}=30.5$ kHz. The value of the frequency corresponds to the experimental driving frequency in Fauconnier et al. (Reference Fauconnier, Mauger, Bera and Inserra2022). We recall that the resonant radius (2.73) is theoretically independent of order $m$ of the spherical harmonics. The amplitudes of the two modes $s^{(n,\pm m)}$ are set equal to 15 µm ( $|s^{(n,\pm m)}/R_{0}|\sim 0.2$ ), so that the case of a stationary surface oscillation is investigated. It is worth mentioning that we have chosen the value of the oscillation amplitude that corresponds approximately to measured amplitudes of asymmetric oscillations in Fauconnier et al. (Reference Fauconnier, Mauger, Bera and Inserra2022). This value leads to a not very small value of the dimensionless parameter $|s^{(n,\pm m)}/R_{0}|$ , which, however, in practice is small enough to ensure the convergence of the first-order solutions (see § 2.1). What is more, from the point of view of the qualitative behaviour of acoustic streaming, the value of the oscillation amplitude is of no importance in our simulations because mathematically the streaming velocity is directly proportional to the oscillation amplitude squared. Hence a change in the oscillation amplitude does not change, qualitatively, the obtained results. Numerical simulations are performed for an asymmetric deformation oscillating at half the driving frequency, i.e. for $f=15.25$ kHz, as the shape oscillations are supposed to be triggered by the Faraday parametric instability. The dimensionless thickness of the viscous boundary layer is therefore $\gamma =\delta _{v}/R_{0}\sim 0.06\ll 1.$

The bubble interface and the resulting microstreaming are presented in figure 4. This display will be used in all investigated numerical cases of standing-wave patterns. Figure 4(a) shows the deformation of the bubble interface on the axisymmetric mode $Y_{3}^{0}$ , with a display similar to that in figure 2 (top, side and three-dimensional views). Keeping the same orientation of the Cartesian frame, figures 4(b)–4(d) present the microstreaming pattern derived from the Lagrangian velocity field. The streamlines have been obtained using the library pyvista from Python programming. The top and side views of the microstreaming are obtained by calculating the streamlines in a narrow area (in depth) along the whole equator (for the top view, figure 4 b) or along a whole meridian (side view, figure 4 c). The procedure is as follows. The considered equator/meridian is segmented in 48 particle clouds, termed as source, homogeneously distributed along the perimeter. A source particle cloud is located at the radial distance $1.01\,R_{0}$ from the centre of the bubble, for a given angular position $(\theta , \phi)$ . Each source contains 10 particles. The streamlines are generated by using each particle as a starting location. For the sake of readability, the three-dimensional streamlines are computed in figure 4(d) only inside a constrained range of the longitude $\phi$ of the three-dimensional space surrounding the bubble, and by calculating 400 streamlines in the investigated region. Note that because of these different numerical calculations for the streamlines, the velocity range displayed in the colour bars may differ from one view to another. The derived streamlines are in full agreement with the theoretical predictions (Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a ) and experimental observations (Cleve et al. Reference Cleve, Guedra, Mauger, Inserra and Blanc-Benon2019; Fauconnier et al. Reference Fauconnier, Mauger, Bera and Inserra2022) of the microstreaming induced by an axisymmetric oscillation. The side view exhibits a cross-like shape with small recirculation loops in the vicinity of the bubble interface. It is worth mentioning that while theoretically all streamlines are closed vortical trajectories, some of the vortices are so extended spatially that they appear as straight lines in figure 4(c). In fact, these streamlines will close at a large distance from the centre of the bubble. Increasing the number of calculated streamlines would fill the region between the two orthogonal branches of the cross-like shape in figure 4(c), for instance. Due to the absence of azimuthal dependence, the streamlines from the top view are purely radial. Note that the small recirculation loops along the meridian are not observed because of the cutting of the microstreaming field along the equator. This situation is similar to that encountered when looking at a streaming field in the focal plane of a microscope (figure 2 a). The rotational direction of the flow can be perceived from a close look at the streamlines in the side and three-dimensional views: it is known that the fluid particles are propelled away from the anti-nodes (where the streamlines start) and attracted back to the nodes of displacement of the bubble interface (where the streamlines finish). This fact will be very useful in the following to discern the three-dimensional pattern of the flow, in the case of asymmetric modes.

3.3. Signature for various asymmetric modes with the same degree $n$

For the same bubble equilibrium radius and mode amplitudes, figures 5–7 show the resulting microstreaming patterns for the class of tesseral ( $Y_{3}^{1}$ , figure 5, and $Y_{3}^{2}$ , figure 6) and sectoral ( $Y_{3}^{3}$ , figure 7) harmonics of the same degree $n=3$ .

Figure 5.

Presentation of the microstreaming pattern observed for a bubble (equilibrium radius 73.8 µm) exhibiting tesseral oscillations on the mode $(n=3, m=1$ ): (a) numerical contours of the bubble interface in the top, side and three-dimensional views; (b) top, (c) side and (d) isometric views of the resulting flow.

Figure 6.

Presentation of the microstreaming pattern observed for a bubble (equilibrium radius 73.8 µm) exhibiting tesseral oscillations on the mode $(n=3, m=2$ ): (a) numerical contours of the bubble interface in the top, side and three-dimensional views; (b) top, (c) side and (d) isometric views of the resulting flow.

Figure 7.

Presentation of the microstreaming pattern observed for a bubble (equilibrium radius 73.8 µm) exhibiting sectoral oscillations on the mode $(n=3, m=3$ ): (a) numerical contours of the bubble interface in the top, side and three-dimensional views; (b) top, (c) side and (d) isometric views of the resulting flow.

The tesseral mode $Y_{3}^{1}$ possesses $m=1$ meridian nodal lines (figure 5 a, top view) and $n-m=2$ nodal parallels (figure 5 a, side view). From the top view configuration, the bubble contour exhibits two nodes and two anti-nodes of displacement. The resulting top view pattern is characterised by a cross shape (figure 5 b), similar to the quadrupole pattern generated by a solid-body translation oscillation without shape deformation (Longuet-Higgins Reference Longuet-Higgins1998; Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019b ). The main difference with this axisymmetric quadrupole pattern lies in the out-of-plane component of the velocity field, as streamlines go away from the equatorial plane and from a closed loop that reaches the nearest displacement node on the bubble interface. The closeness of this top view asymmetric pattern to the translation-induced quadrupole shape confirms the importance of capturing both the acoustic (high-frequency) dynamics of the bubble interface and the fluid flow, in order to avoid misinterpretation of the physical origin of the streaming pattern. In particular, from the top view, the tesseral $m=1$ harmonic, whatever the degree $n$ , always exhibits an up-and-down oscillation around the $x$ axis that resembles the translational oscillation of the bubble along the $y$ axis. From the side view, the bubble interface displays six extrema of displacement. It results in a recognisable 12-lobe flow pattern (figure 5 c), where the lobes are arranged in pairs and are oriented from the extremum of the bubble interface towards the nearest displacement node. Combining these two views, and the fact that the fluid particles move from the anti-node towards the nearest displacement node on the bubble interface, the isometric view of the pattern can be guessed in figure 5(d).

The tesseral mode $Y_{3}^{2}$ possesses $m=2$ meridian nodal lines (figure 6 a, top view) and $n-m=1$ nodal parallels (figure 6 a, side view). From the top view configuration, the bubble contour exhibits four nodes and four anti-nodes of displacement. Note that because the nodal parallel is located on the equator of the bubble interface, the so-called top view nodes and anti-nodes are located out of the equatorial plane. The resulting top view pattern (figure 6 b) is characterised by a lobe-type pattern with $4m=8$ recirculation loops arranged in pairs. The overall structure of the top view streaming is in good agreement with the experimental top view pattern displayed in figure 3(b) for another tesseral oscillation whose order also equals $m=2.$ The microstreaming pattern from the side view (figure 6 c) is harder to interpret as the bubble interface displays four extrema of displacement whose location is dependent on the orientation of the bubble around the $z$ axis. Here, eight lobes are arranged in pairs and are superimposed on a cross-type pattern. A closer look at the pairs of lobes located at the poles ( $\theta =0$ and $\theta ={\pi}$ ) reveals that the recirculation loops are in fact split into two loops. This artefact is due to the slicing of the streamlines along the plane containing meridians. In fact, the isometric view (figure 6 d) indicates how the streamlines start from the displacement anti-node located in the upper right corner in figure 6(d), and reach the nearest displacement nodes that form a triangular shape surrounding this anti-node. It results in a seemingly axisymmetric vortex around the normal at the bubble interface passing through the displacement anti-node.

The results for an asymmetric bubble experiencing sectoral oscillations on the mode $Y_{3}^{3}$ are shown in figure 7. Being devoid of nodal parallels, the maximal longitudinal displacement occurs at the equatorial plane, where the bubble interface exhibits as many azimuthal lobes as the modal degree $n$ to which it belongs. The amplitude of the oscillation displacement at a given acoustic phase decreases along the line going from the equator to the poles of the bubble. For $n=3$ , it results in $2n=6$ extrema for the displacement of the interface along the equator (figure 7 a), an azimuthal shape that is easily recognisable from the top view observation. From the side view, the bubble contour is almost spherical, with an outgrowth located at the equator. Due to the (relative) simplicity of the bubble deformation along the elevation and the azimuthal direction, the resulting microstreaming pattern is easy to analyse from the top, side, and even isometric view. The top view is characterised by a $4n=12$ -lobe flower-like pattern, where the lobes are arranged in pairs around all extrema of the bubble interface displacement. This pattern is in good qualitative agreement with the experimental microstreaming pattern shown in figure 3(c), both for the number of recirculation loops and for the spatial extension of the vortices. The side view resembles a cross-like shape (figure 7 c), typical of a quadrupole pattern. This pattern is the analogue of the top view streaming resulting from the tesseral mode $m=1$ (figure 5 b) when switching from the top to the side view observation. Indeed, sectoral oscillations from the side view exhibit an out-of-phase motion along a meridian from one side of the bubble to the opposite one around the $z$ axis, so that the oscillation seems like that of a solid-body translation oscillation. Because the sectoral harmonic $Y_{3}^{3}$ possesses only meridian nodal lines, the vortices observed from the top view configuration spread with a decreasing amplitude towards the poles of the bubble (figure 7 d).

The numerical streaming patterns obtained for the self-interacting tesseral (figure 6) or sectoral (figure 7) oscillations are found to be in good agreement with the patterns (figure 3) obtained by Fauconnier et al. (Reference Fauconnier, Mauger, Bera and Inserra2022). Unlike the experimental observations of the latter, the present theory considers neither the presence of a wall in the bubble vicinity, nor the tethering of the bubble on a surface. It was shown in Fauconnier et al. (Reference Fauconnier, Bera and Inserra2020) that the contact angle between the bubble interface and the substrate lies in the range 40°–60°, with no dependence on the bubble equilibrium radius. It is also indicated that the contact angle remains the same after the activation of the ultrasound driving sequence. For this range of contact angle values, the bubble tethering does not influence strongly the dynamics of the oscillation, and the majority of the existing asymmetric oscillations for a given degree $n$ are triggered. Concerning the resulting acoustic microstreaming, it is known that the mathematical modelling of the flow surrounding a free bubble, being far from any boundary, can be in fairly good agreement with the experimental flows in the form of a large-scale streaming induced by a substrate-attached bubble (Marmottant & Hilgenfledt Reference Marmottant and Hilgenfeldt2003). Here, in the case of asymmetric oscillations, the vortices with the highest velocity magnitude are confined in the close vicinity of the bubble interface (figure 7 b) and extend much less than the large-scale quadrupole-like streaming. It can then be inferred that the influence of the bubble tethering on the resulting flow can be less pronounced in comparison to the case of an axisymmetric flow. In addition, it should be mentioned that our theory does not consider the effect of buoyancy, which may cause the deviation of experimental streamlines from theoretical Lagrangian streamlines. This fact should be taken into account when comparing experimental results with our theory.

3.4. Extension to other spherical harmonics of any degree $n$

The previously investigated case of the microstreaming induced by the three classes of spherical harmonics for the same degree $n=3$ can help us in predicting the pattern for any asymmetric modes, whatever the degree $n$ and order $m$ . The simplest case concerns bubbles exhibiting zonal oscillations $Y_{n}^{0}$ , for which theoretical predictions exist (Inserra et al. Reference Inserra, Regnault, Cleve, Mauger and Doinikov2020a ). Due to the absence of azimuthal dependence for the bubble displacement, the top view pattern of a zonal-induced microstreaming is unique and consists in a purely radial flow. Due to the axisymmetry of the zonal harmonics, slicing the three-dimensional microstreaming pattern in a plane containing the $z$ axis always provides the same pattern, which consists of a large-scale cross pattern with small recirculation zones in the vicinity of the bubble interface. This general signature is identical whatever the degree $n$ of the investigated zonal harmonics. Only the number of small vortices in the vicinity of the bubble interface increases as the degree $n$ increases.

In general, a bubble exhibiting an asymmetric oscillation on a tesseral harmonic of order $m$ , whatever the degree $n$ , will exhibit the same top view patterns as shown in figure 5 when $m=1$ , or figure 6 when $m=2$ . Indeed, these patterns are uniquely linked to the number of meridian nodal lines, which equals $m$ . For tesseral harmonics with higher order $m$ , top view streaming will be characterised by a lobe-type pattern with $4m$ lobes arranged in pairs. This top view signature of the pattern resulting from a self-interacting tesseral mode is similar to the shape obtained from another self-interacting spherical harmonics (obtained from a side or top view). This means that without knowing the bubble equilibrium radius and hence the degree $n$ ruled by (2.73), the identification of the triggered mode at the bubble interface cannot be inferred when viewing the resulting streaming pattern. The analysis from the side view only is even more complicated, as the orientation of the focal plane around the $z$ axis significantly affects the conformation of the captured flow.

The prediction of the microstreaming pattern induced by the class of sectoral harmonics is far easier to infer in comparison to the tesseral harmonics. All the sectoral harmonics $m=n$ possess the same spatial conformation for the bubble displacement, with a maximal longitudinal displacement occurring at the equatorial plane. Along the equator, as many azimuthal lobes as the modal degree $n$ to which they belong exist. The bubble displacement at the poles always tends to zero, so that the side view microstreaming pattern always resembles a cross-like shape, whatever the degree $n$ of the sectoral oscillations. This is evidenced in figure 8, where the top and side view patterns resulting from self-interacting sectoral oscillations from $n=2$ to $n=5$ are shown. Here, the bubble equilibrium radii are chosen so that they correspond to the resonant radius of a given degree $n$ : $R_{res}^{2}=46\ \unicode{x03BC}\rm m$ ( $\gamma =\delta _{v}/R_{0}\sim 0.1)$ , $R_{res}^{3}=68\ \unicode{x03BC}\rm m$ ( $\gamma =\delta _{v}/R_{0}\sim 0.06)$ , $R_{res}^{4}=90\ \unicode{x03BC}\rm m$ ( $\gamma =\delta _{v}/R_{0}\sim 0.05)$ and $R_{res}^{5}=110\ \unicode{x03BC}\rm m$ ( $\gamma =\delta _{v}/R_{0}\sim 0.04$ ). The amplitudes of the two modes $s^{(n,\pm n)}$ are set equal to 15 µm, resulting in a stationary sectoral oscillation. The top view microstreaming is characterised by a lobe-type pattern where the lobes are arranged in pairs, and their number equals twice the number of extrema of the bubble displacement along the equator. Therefore, a $4n$ -lobe shape is expected from the top view, as shown in figure 8. This signature is in full agreement with experimental observations of top view microstreaming patterns induced by sectoral oscillation for various degrees $n$ in Fauconnier et al. (Reference Fauconnier, Mauger, Bera and Inserra2022).

Figure 8.

Presentation of the microstreaming pattern observed for various sectoral oscillations, with the degree $n$ ranging from 2 to 5. Every line corresponds to a given degree $n$ . (a,d,g,j) The projected side view microstreaming pattern. (b,e,h,k) The deformation of the bubble interface exhibiting the sectoral oscillation. (c,f,i,l) The top view microstreaming pattern.

3.5. Influence of the thickness of the viscous boundary layer on the streaming pattern

The acoustic microstreaming is driven by the streaming inside the oscillatory boundary layer around the bubble, while nonlinear second-order effects of the hydrodynamic equations are responsible for extending the streaming field much further than the viscous boundary layer. This explains why Elder (Reference Elder1959) reported different types of streaming patterns whose appearance depended mostly on the bubble surface velocity and the fluid viscosity. For a given bubble equilibrium radius and driving frequency, the influence of the thickness of the viscous boundary layer was investigated for varying values of the liquid viscosity in Doinikov et al. (Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a ) and Inserra et al. (Reference Inserra, Regnault, Cleve, Mauger and Doinikov2020b ) in the case of axisymmetric shape oscillations. For the translational oscillation of a solid body, Li et al. (Reference Li, Collis, Brumley, Schneiders and Sader2023) reveal changes in the structure of the streaming pattern for varying driving frequencies while keeping the body size and the liquid viscosity constant. Here, we investigate the influence of the value of the dimensionless thickness of the viscous boundary layer on the microstreaming pattern. We have previously considered the case (figure 7) of a sectoral oscillation $n=m=3$ at frequency $f=15.25$ kHz for a bubble with equilibrium radius 73.8 µm in water (dynamic viscosity $\eta =10^{-6}$ Pa s). It corresponds to the dimensionless thickness of the viscous boundary layer $\gamma =\delta _{v}/R_{0}\sim 0.06.$ In order to modify the thickness of the viscous boundary layer $\delta _{v}$ , the liquid viscosity is varied while keeping the bubble equilibrium radius and the oscillation frequency constant, so that the triggered shape oscillation given by (2.73) remains unchanged. For the same parameters as in figure 7, figure 9 presents the microstreaming patterns when the liquid viscosity is varied from $\eta =10^{-6}$ Pa s ( $\gamma =0.06$ , figure 9 a) to $\eta =4\times10^{-6}$ Pa s ( $\gamma =0.12$ , figure 9 b) and $\eta =16\times10^{-6}$ Pa s ( $\gamma =0.24$ , figure 9 c).

Figure 9.

Evolution of the microstreaming pattern for the sectoral oscillation $n=m=3$ for varying values of the dimensionless thickness of the viscous boundary layer. The liquid dynamic viscosities and resulting dimensionless thickness of the viscous boundary layer are (a) $\eta =10^{-6}$ Pa s, $\gamma =0.06$ ; (b) $\eta =4\times10^{-6}$ Pa s, $\gamma =0.12$ , and (c) $\eta =16\times10^{-6}$ Pa s, $\gamma =0.24$ . Only the top and side views are displayed.

When the size of the viscous boundary layer is varied by four times, the overall structure of the streaming pattern remains unchanged: the top view is characterised by a 12-lobe flower-like pattern, and the side view looks like a cross (quadrupolar) shape. When the ratio $\gamma =\delta _{v}/R_{0}$ is increased, the magnitude of the maximum streaming velocity decreases, while the size of the small vortices visible in the top view increases slightly. More pronounced modifications of the microstreaming pattern can be expected for large values of the dimensionless parameter $\gamma$ , but such values would correspond to unrealistic experimental conditions (too large kinematic viscosity or a bubble size or driving frequency no longer matching (2.73)).

3.6. Sign reversal of the flow at the bubble zenith

Figure 10.

Investigation of the orientation of the flow at the zenith of the bubble ( $\theta =0$ ). (a) Three-dimensional view of the bubble interface oscillating on asymmetric modes with $n\in [2,5]$ . (b) Evolution of the value of the colatitude $\theta _{max}$ at which the displacement of the bubble interface is maximum, as a function of the order $m$ for various degrees $n$ . (c) Evolution of the amplitude of the radial velocity at the zenith ( $\theta =0$ ) of the bubble, at the distance $r=3R_{0}$ as a function of the order $m$ for various degrees $n$ . A positive value corresponds to the fountain-like streaming, while a negative value is typical of the anti-fountain behaviour.

Substrate-attached microbubbles are commonly used as microactuators for the trapping and manipulation of small particles or biological bodies in their vicinity. The hydrodynamic flow field generated by the bubble oscillation can rotate single microparticles within microchannels (Ahmed et al. Reference Ahmed, Ozcelik, Bojanala, Nama, Upadhyay, Chen, Hanna-Rose and Huang2016), precisely manipulate millimetric fish eggs (Lee et al. Reference Lee, Lee, Won, Rhee and Chung2012), or manipulate and rupture vesicles (Marmottant et al. Reference Marmottant, Biben and Hilgenfeldt2008). In most microfluidic applications, a substrate-attached microbubble experiences a combination of spherical and translational oscillations, at the origin of a large-scale microstreaming pattern constituted of two symmetric vortices (in the plane of observation, from a side view) atop the bubble. When particles are injected near the oscillating bubble, they are expelled from the top bubble interface and move along a closed loop back to the substrate. This motion is often called fountain-like streaming, as the fluid particles are expelled from the bubble interface to the surrounding medium at the zenith ( $\theta =0$ ) of the bubble. A similar fountain behaviour is obtained when axisymmetric modes are triggered at the bubble interface, because the zenith of the bubble interface is always a displacement anti-node, whatever the degree $n$ of the zonal harmonics (see the first column of figure 10 a). In the case of wall-bounded semi-cylindrical bubbles, Wang et al. (Reference Wang, Rallabandi and Hilgenfeldt2013) and Rallabandi et al. (Reference Rallabandi, Wang and Hilgenfeldt2014) have observed an inversion of the flow direction atop the bubble, identifying an anti-fountain behaviour of the vortices near the pole. This flow inversion was shown to result from the interaction between the radial oscillation and an axisymmetric shape oscillation (mixed-mode streaming). When asymmetric modes are triggered, the location of the maximum displacement on the bubble interface evolves as the order $m$ increases, tending to a maximum at the equator of the bubble in the case of a sectoral mode. Figure 10(a) illustrates the three-dimensional bubble shapes for various degrees $n$ ranging from 2 to 5. In addition, the value of the colatitude $\theta _{max}$  at which the location of the maximum displacement on the bubble interface occurs, whatever the value of the azimuthal angle $\phi$ , is shown in figure 10(b). When the order $m$ is smaller than $(n+1)/2$ , $\theta _{max}$ ranges from $0$ to ${\pi} /4$ , and the maximum displacement at the bubble interface is located in the upper part of the bubble, between the zenith and the bisector $\theta _{max}={\pi} /4$ . For higher values of $m$ , $\theta _{max}$ gets closer to ${\pi} /2$ , and the maximum displacement tends to the equatorial plane.

For the investigated asymmetric modes, the direction of the resulting microstreaming has been evaluated atop the bubble, at the radial distance $r=3R_{0}$ and for $\theta =0$ . Figure 10(c) shows the amplitude of the radial component of the Lagrangian velocity at this location. A sign reversal of the velocity amplitude occurs when the order $m$ reaches $(n+1)/2$ , in agreement with the evolution of the colatitude $\theta _{max}$ . In the case of zonal oscillations, the fluid particles are expelled (positive radial component of the velocity) from the zenith of the bubble, resulting in a fountain-like behaviour. In the case of sectoral oscillations, the analysis of the deformation of the bubble interface reveals that the maximum displacement occurs at the equator, while the zenith ( $\theta =0$ ) of the bubble becomes a displacement node. The fluid particles are therefore expelled from the equator and are attracted from the zenith of the bubble to the bubble interface. This results in an anti-fountain behaviour, which has not been observed experimentally so far in the case of self-interacting mode streaming. It is worth noting that the observation of the flow direction atop the bubble indirectly provides information on the triggered shape oscillation, as the direction of the flow is correlated to the order $m$ of the asymmetric oscillation. If the bubble equilibrium radius is known, then the degree $n$ of the shape oscillations is obtained from the Lamb spectrum (2.73), and the range of possible triggered order $m$ can be deduced from the observation of a fountain or anti-fountain flow.

3.7. The case of travelling surface waves

Up to now, only the stationary shape oscillations were considered, meaning that the amplitudes of the two modes $s^{(n,\pm m)}$ were equal. A slight change in the amplitude of the surface perturbation propagating in one direction along the azimuthal angle $\phi$ in comparison to the one propagating in the opposite direction will result in a quasi-stationary surface wave. Such a bubble will experience a rotating asymmetric surface wave, whose angular velocity will depend on the difference in the two amplitudes of modes $(n,\pm m)$ . This has already been reported once in the literature by Mekki-Berrada, Thibault & Marmottant (Reference Mekki-Berrada, Thibault and Marmottant2016), who measured a constant angular velocity of approximately 0.5 revolutions per second for every modal degree $n$ , in the case of bubbles flattened between two elastic walls. The underlying cause for the triggering of a rotating shape oscillation so far remains unknown.

The case of travelling surface waves is investigated in figure 11, for different ratios of the amplitudes of the two propagating modes $(n,\pm m)$ . The case of the sectoral oscillation $n=m=3$ is investigated, for which the microstreaming pattern obtained in the case of a stationary wave is presented in figure 7. Only the top and isometric views of the streaming pattern are displayed in figure 11, and the amplitude of the $(n,m)$ wave is kept equal to $15$ µm. When the amplitude of the $(n,-m)$ wave is modified to $14.9$ µm, i.e. when $s^{(n,-m)}/s^{(n,m)}=0.99$ , streamlines are curved along the azimuthal direction (figure 11(a), top view) and the maximum amplitude of the Lagrangian velocity reaches 8.3 cm s−1. From figure 7, it is worth noting that the maximum amplitude of the Lagrangian velocity reaches 7.65 cm s−1 in the case of a stationary wave, i.e. when $s^{(n,-m)}/s^{(n,m)}=1$ . The increase in difference between the amplitudes of the two counter-propagating surface waves results in a stronger azimuthal component of the velocity field. The amplitude of the streaming flow can increase by a factor of 2 (figure 11 b) when $s^{(n,-m)}/s^{(n,m)}=0.93$ , or a factor of 4 (figure 11 c) when $s^{(n,-m)}/s^{(n,m)}=0.66$ , in comparison to the stationary case.

Figure 11.

Evolution of the microstreaming pattern in the case of a quasi-stationary wave, for the sectoral oscillation $n=m=3$ . Three ratios of the amplitudes of the two modes $s^{(n,-m)}/s^{(n,m)}$ are investigated for (a) $s^{(n,-m)}/s^{(n,m)}=0.99$ , (b) $s^{(n,-m)}/s^{(n,m)}=0.93$ and (c) $s^{(n,-m)}/s^{(n,m)}=0.66$ . Only the top and isometric views are displayed.

The morphology of the swirling flow can be explained as follows. When the amplitudes $s^{(n,\pm m)}$ of a given shape oscillation are set equal, a standing shape oscillation is formed on the bubble interface. This means that displacement anti-nodes (starting points of the streamlines) and nodes (end points of the streamlines) remain at the same location with time. The resulting flow is constituted of vortical lines that are spatially stationary. When the amplitudes $s^{(n,\pm m)}$ of the two components of a shape oscillation differ, a quasi-propagating wave is triggered in a preferential direction on the bubble interface. ‘Quasi-propagating’ refers to the coexistence of a purely standing wave at the bubble interface and a purely propagating wave. Hence the locations of the displacement nodes and anti-nodes constantly move along the bubble interface with time. This prevents the spatial stationarity of the flow, so a convective-like flow is generated around the bubble. This flow follows the direction of propagation of the propagative surface wave, and is characterised by a strong azimuthal component of the flow velocity and a large-scale extension at the equator. However, why this flow is stronger in amplitude in comparison to the one induced by the standing wave alone is difficult to explain. Nevertheless, these results highlight the possibility to drastically enhance the streaming efficiency around bubbles if the bubble rotation can be controlled.

4. Conclusion

We present for the first time a theoretical modelling of acoustic microstreaming induced by a bubble experiencing asymmetric oscillations. The modelling is based on the decomposition of the first- and second-order vorticity field into the poloidal and toroidal fields, allowing the exact analytical derivation of the Lagrangian streaming induced by the bubble. This derivation is valid whatever the bubble size and the liquid viscosity, without any limitation on the thickness of the viscous boundary layer at the bubble interface. By representing the shape of the bubble interface in terms of the set of the orthonormal spherical harmonics, the signatures of the microstreaming patterns are investigated as a function of different classes of triggered spherical harmonics (zonal, tesseral or sectoral oscillations). Keeping in mind that the streamlines start from the displacement anti-node at the bubble interface, and close at the displacement nodes, the three-dimensional asymmetric patterns are shown to be easily predicted once the asymmetric shape deformation of the bubble is known. For a given class of spherical harmonics of degree $n$ , it is shown that the sign reversal in the flow orientation atop the bubble is obtained when the order $m$ evolves between $0$ and $n$ : zonal harmonics result in the fountain-like streaming, while the anti-fountain flow occurs in the case of sectoral oscillations. In between, the flow reverses for a typical value of the order $m$ that helps in the identification of the triggered asymmetric mode by means of the indirect observation of the flow orientation. The case of travelling surface waves along the bubble interface reveals the possibility of a drastic enhancement of the streaming strength. This observation highlights the interest in the controlling of oscillations of rotating bubbles for practical use.

Funding.

This work was funded by l’Agence Nationale de la Recherche under the project ANR-22-CE 92-0062, and was supported by the LabEx CeLyA of the University of Lyon (grant nos ANR-10-LABX-0060 and ANR-11-IDEX-0007).

Declaration of interests.

The authors report no conflict of interest.

Appendix A. Calculation of $\boldsymbol{a}_{\boldsymbol{nm}}$ , $\boldsymbol{b}_{\boldsymbol{nm}}$ and $\boldsymbol{c}_{\boldsymbol{nm}}$

To find $a_{nm}$ , $b_{nm}$ and $c_{nm}$ , we apply boundary conditions at the bubble surface.

The first boundary condition requires that the normal component of $\boldsymbol{v}_{1}^{(n,m)}$ at $r=R_{0}$ be equal to the respective normal component of the bubble surface velocity. This condition is written as

(A1) \begin{equation} \left.v_{1r}^{\left(n,m\right)}\right| _{r={R_{0}}}=\frac{{\rm d}r_{s}^{\left(n,m\right)}}{{\rm d}t}=-i\omega s^{\left(n,m\right)}\,e^{-i\omega t}\,Y_{n}^{m}\left(\theta ,\phi \right) ,\end{equation}

where $v_{1r}^{(n,m)}$ is the normal component of $\boldsymbol{v}_{1}^{(n,m)}$ , and $r_{s}^{(n,m)}$ is given by the second term of (2.1).

The second boundary condition assumes slippage on the bubble surface, which means that the tangential stress generated by the liquid motion should be zero on the bubble surface. This condition is written as

(A2) \begin{align} \sigma _{r\theta }=\eta \left(\frac{1}{r}\frac{\partial v_{1r}^{\left(n,m\right)}}{\partial \theta }+\frac{\partial v_{1\theta }^{\left(n,m\right)}}{\partial r}-\frac{v_{1\theta }^{\left(n,m\right)}}{r}\right)=0\quad \text{at }r=R_{0} ,\quad \end{align}
(A3) \begin{align} \sigma _{r\phi }=\eta \left(\frac{\partial v_{1\phi }^{\left(n,m\right)}}{\partial r}+\frac{1}{r\sin \theta }\frac{\partial v_{1r}^{\left(n,m\right)}}{\partial \phi }-\frac{v_{1\phi }^{\left(n,m\right)}}{r}\right)=0\quad \text{at }r=R_{0} , \end{align}

where $\sigma _{r\theta }$ and $\sigma _{r\phi }$ are the tangential components of the liquid stress in spherical coordinates (Landau & Lifshitz Reference Landau and Lifshitz1987). It should be mentioned that in the case when there are impurities and surfactants on the bubble surface, and in the case of a contrast agent bubble having a shell, the boundary condition of slippage should be replaced by the no-slip boundary condition, i.e. by the condition of adhesion of liquid particles to the bubble surface. In the case of a bubble with a shell, the behaviour of the shell material should also be taken into consideration. In view of its complexity, this case requires a separate study. The general algorithm of calculations remains the same, but all constants related to the boundary conditions at the bubble surface should be recalculated.

Substitution of (2.28) into (2.20) yields

(A4) \begin{align} \boldsymbol{v}_{1}^{\left(n,m\right)}=v_{1r}^{\left(n,m\right)}\boldsymbol{e}_{r}+v_{1\theta }^{\left(n,m\right)}\boldsymbol{e}_{\theta }+v_{1\phi }^{\left(n,m\right)}\boldsymbol{e}_{\phi } ,\qquad\qquad\quad\end{align}
(A5) \begin{align} v_{1r}^{\left(n,m\right)}=e^{-i\omega t}\,f_{nm}\left(r\right)\,Y_{n}^{m}\left(\theta ,\phi \right) ,\qquad\qquad\qquad\end{align}
(A6) \begin{align} v_{1\theta }^{\left(n,m\right)}=e^{-i\omega t}\left[\frac{a_{nm}k_{v}h_{n}^{\left(1\right)}\left(k_{v}r\right)}{\sin \theta }\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \phi }+g_{nm}\left(r\right)\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta }\right] ,\end{align}
(A7) \begin{align} v_{1\phi }^{\left(n,m\right)}=e^{-i\omega t}\left[-a_{nm}k_{v}h_{n}^{\left(1\right)}\left(k_{v}r\right)\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta }+\frac{g_{nm}\left(r\right)}{\sin \theta }\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \phi }\right] ,\end{align}

where $f_{nm}(r)$ and $g_{nm}(r)$ are given by

(A8) \begin{align} f_{nm}\left(r\right)=\frac{\left(n+1\right)c_{nm}}{r^{n+2}}+\frac{n\left(n+1\right)b_{nm}h_{n}^{\left(1\right)}\left(k_{v}r\right)}{k_{v}r} ,\qquad\quad \end{align}
(A9) \begin{align} g_{nm}\left(r\right)=-\frac{c_{nm}}{r^{n+2}}+\frac{b_{nm}}{k_{v}r}\left[\left(n+1\right)h_{n}^{\left(1\right)}\left(k_{v}r\right)-k_{v}rh_{n+1}^{\left(1\right)}\left(k_{v}r\right)\right] . \end{align}

Substituting (A5) and (A8) into (A1), one obtains

(A10) \begin{equation} c_{nm}+\frac{nR_{0}^{n+2}h_{n}^{\left(1\right)}\left(\overline{x}\right)}{\overline{x}}b_{nm}=-\frac{i\omega R_{0}^{n+2}s^{\left(n,m\right)}}{n+1} ,\end{equation}

where $\overline{x}=k_{v}R_{0}$ .

Substitution of (A5)–(A7) into (A2) and (A3) yields

(A11) \begin{align} &a_{nm}k_{v}\left[\overline{x}h_{n}^{\left(1\right)\prime}\left(\overline{x}\right)-h_{n}^{\left(1\right)}\left(\overline{x}\right)\right]\frac{1}{\sin \theta }\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \phi }\nonumber\\&\quad+\left[f_{nm}\left(R_{0}\right)-g_{nm}\left(R_{0}\right)+R_{0}g_{nm}^{\prime}\left(R_{0}\right)\right]\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta }=0, \end{align}
(A12) \begin{align} &a_{nm}k_{v}\left[h_{n}^{\left(1\right)}\left(\overline{x}\right)-\overline{x}h_{n}^{\left(1\right)\prime}\left(\overline{x}\right)\right]\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta }\nonumber\\&\quad+\left[f_{nm}\left(R_{0}\right)-g_{nm}\left(R_{0}\right)+R_{0}g_{nm}^{\prime}\left(R_{0}\right)\right]\frac{1}{\sin \theta }\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \phi }=0, \\[6pt] \nonumber \end{align}

where

(A13) \begin{equation} g_{nm}^{\prime}\!\left(r\right)=\frac{\left(n+2\right)c_{nm}}{r^{n+3}}+\frac{b_{nm}}{k_{v}r^{2}}\!\left\{\left(n+1\right)\!\left[k_{v}rh_{n}^{\left(1\right)/}\!\left(k_{v}r\right)-h_{n}^{\left(1\right)}\!\left(k_{v}r\right)\right]-\left(k_{v}r\right)^{2}h_{n+1}^{\left(1\right)\prime}\left(k_{v}r\right)\right\}\! .\end{equation}

Multiplying (A11) by $(\sin ^{-1}\theta)\, \partial Y_{n}^{m}(\theta ,\phi )/\partial \phi$ and subtracting (A12) multiplied by $\partial Y_{n}^{m}(\theta ,\phi )/\partial \theta$ , one obtains

(A14) \begin{equation} a_{nm}\left[\overline{x}h_{n}^{\left(1\right)/}\left(\overline{x}\right)-h_{n}^{\left(1\right)}\left(\overline{x}\right)\right]\left\{\left[\frac{1}{\sin \theta }\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \phi }\right]^{2}+\left[\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta }\right]^{2}\right\}=0 .\end{equation}

It follows from (A14) that

(A15) \begin{equation} a_{nm}=0 .\end{equation}

Accordingly, (A6) and (A7) take the form

(A16) \begin{align} v_{1\theta }^{\left(n,m\right)}=e^{-i\omega t}g_{nm}\left(r\right)\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta } , \end{align}
(A17) \begin{align} v_{1\phi }^{\left(n,m\right)}=e^{-i\omega t}\frac{g_{nm}\left(r\right)}{\sin \theta }\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \phi } .\end{align}

Multiplying (A11) by $\partial Y_{n}^{m}(\theta ,\phi )/\partial \theta$ and adding to (A12) multiplied by $(\sin ^{-1}\theta)\, \partial Y_{n}^{m}(\theta ,\phi )/\partial \phi$ , one obtains

(A18) \begin{equation} \left[f_{nm}\left(R_{0}\right)-g_{nm}\left(R_{0}\right)+R_{0}g_{nm}^{\prime}\left(R_{0}\right)\right]\left\{\left[\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta }\right]^{2}+\left[\frac{1}{\sin \theta }\frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \phi }\right]^{2}\right\}=0 .\end{equation}

It follows from (A18) that

(A19) \begin{equation} f_{nm}\left(R_{0}\right)-g_{nm}\left(R_{0}\right)+R_{0}g_{nm}^{\prime}\left(R_{0}\right)=0 .\end{equation}

Substituting (A8) and (A9) into (A19) and using the identity

(A20) \begin{equation} h_{n+1}^{\left(1\right)}\left(x\right)=\frac{n}{x}h_{n}^{\left(1\right)}\left(x\right)-h_{n}^{\left(1\right)\prime}\left(x\right) ,\end{equation}

one obtains

(A21) \begin{equation} \frac{2\left(n+2\right)}{R_{0}^{n+2}}c_{nm}+\frac{b_{nm}}{\overline{x}}\left[(n^{2}+n-2)\,h_{n}^{\left(1\right)}\left(\overline{x}\right)+\overline{x}^{2}h_{n}^{\left(1\right)\prime\prime}\left(\overline{x}\right)\right]=0 .\end{equation}

Combining (A10) and (A21), one finds

(A22) \begin{equation} c_{nm}=s^{\left(n,m\right)}\frac{i\omega R_{0}^{n+2}\left[\left(2-n-n^{2}\right)h_{n}^{\left(1\right)}\left(\overline{x}\right)-\overline{x}^{2}h_{n}^{\left(1\right)\prime\prime}\left(\overline{x}\right)\right]}{\left(n+1\right)\left[\overline{x}^{2}h_{n}^{\left(1\right)\prime\prime}\left(\overline{x}\right)-\left(n^{2}+3n+2\right)h_{n}^{\left(1\right)}\left(\overline{x}\right)\right]} ,\end{equation}
(A23) \begin{equation} b_{nm}=s^{\left(n,m\right)}\frac{2i\left(n+2\right)\overline{x}\omega }{\left(n+1\right)\left[\overline{x}^{2}h_{n}^{\left(1\right)\prime\prime}\left(\overline{x}\right)-\left(n^{2}+3n+2\right)h_{n}^{\left(1\right)}\left(\overline{x}\right)\right]},\quad n\geq 1 .\end{equation}

Appendix B. Solving the equations of acoustic streaming for the mode $\boldsymbol{(n,m)}$

To solve (2.52) and (2.53), we first need to calculate their right-hand sides.

The calculation of $\boldsymbol{v}_{1}^{(n,m)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{(n,m)}$ yields

(B1) \begin{align} &\boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}=\left(\!v_{1r}^{\left(n,m\right)}\frac{\partial }{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial }{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial }{\partial \phi }\!\right)\!\left(\!v_{1r}^{\left(n,m\right)}\boldsymbol{e}_{r}+v_{1\theta }^{\left(n,m\right)}\boldsymbol{e}_{\theta }+v_{1\phi }^{\left(n,m\right)}\boldsymbol{e}_{\phi }\!\right)\nonumber\\&=\boldsymbol{e}_{r}\left(v_{1r}^{\left(n,m\right)}\frac{\partial v_{1r}^{\left(n,m\right)}}{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial v_{1r}^{\left(n,m\right)}}{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial v_{1r}^{\left(n,m\right)}}{\partial \phi }-\frac{\left(v_{1\theta }^{\left(n,m\right)}\right)^{2}+\left(v_{1\phi }^{\left(n,m\right)}\right)^{2}}{r}\right)\nonumber\\&+\boldsymbol{e}_{\theta }\!\left(\!v_{1r}^{\left(n,m\!\right)}\frac{\partial v_{1\theta }^{\left(n,m\!\right)}}{\partial r}+\frac{v_{1\theta }^{\left(n,m\!\right)}}{r}\frac{\partial v_{1\theta }^{\left(n,m\!\right)}}{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\!\right)}}{r\sin \theta }\frac{\partial v_{1\theta }^{\left(n,m\!\right)}}{\partial \phi }+\frac{v_{1r}^{\left(n,m\!\right)}v_{1\theta }^{\left(n,m\!\right)}}{r}-\frac{\cos \theta }{r\sin \theta }\!\left(\!v_{1\phi }^{\left(n,m\!\right)}\!\right)^{2}\!\right)\nonumber\\&+\boldsymbol{e}_{\phi }\!\left(\!v_{1r}^{\left(n,m\!\right)}\frac{\partial v_{1\phi }^{\left(n,m\!\right)}}{\partial r}\!+\!\frac{v_{1\theta }^{\left(n,m\!\right)}}{r}\frac{\partial v_{1\phi }^{\left(n,m\!\right)}}{\partial \theta }\!+\!\frac{v_{1\phi }^{\left(n,m\!\right)}}{r\sin \theta }\frac{\partial v_{1\phi }^{\left(n,m\!\right)}}{\partial \phi }\!+\!\frac{v_{1r}^{\left(n,m\!\right)}v_{1\phi }^{\left(n,m\!\right)}}{r}\!+\!\frac{\cos \theta }{r\sin \theta }v_{1\theta }^{\left(n,m\!\right)}v_{1\phi }^{\left(n,m\!\right)}\!\right)\!. \end{align}

Substituting (A5), (A16) and (A17) into (B 1), and averaging it over time, one obtains

(B2) \begin{align} &\left\langle \boldsymbol{v}_{1}^{\left(n,m\!\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\!\right)}\right\rangle \nonumber\\&=\frac{1}{2}\,\mathrm{Re}\left\{\boldsymbol{e}_{r}\left[f_{nm}^{\prime}f_{nm}^{\ast }Y_{n}^{m}Y_{n}^{m\ast }+\frac{\left(f_{nm}-g_{nm}\right)g_{nm}^{\ast }}{r}\left(\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial Y_{n}^{m\ast }}{\partial \theta }+\frac{1}{\sin ^{2}\theta }\frac{\partial Y_{n}^{m}}{\partial \phi }\frac{\partial Y_{n}^{m\ast }}{\partial \phi }\right)\right]\right. \nonumber\\&+\boldsymbol{e}_{\theta }\left[f_{nm}^{\ast }\left(g_{nm}^{\prime}+\frac{g_{nm}}{r}\right)Y_{n}^{m\ast }\frac{\partial Y_{n}^{m}}{\partial \theta }+\frac{|g_{nm}|^{2}}{r}\left(\frac{\partial Y_{n}^{m\ast }}{\partial \theta }\frac{\partial ^{2}Y_{n}^{m}}{\partial \theta ^{2}}+\frac{1}{\sin ^{2}\theta }\frac{\partial Y_{n}^{m\ast }}{\partial \phi }\frac{\partial ^{2}Y_{n}^{m}}{\partial \theta \partial \phi }\right.\right. \nonumber\\&-\left.\left.\frac{\cos \theta }{\sin ^{3}\theta }\frac{\partial Y_{n}^{m\ast }}{\partial \phi }\frac{\partial Y_{n}^{m}}{\partial \phi }\right)\right]+\frac{\boldsymbol{e}_{\phi }}{\sin \theta }\left[f_{nm}^{\ast }\left(g_{nm}^{\prime}+\frac{g_{nm}}{r}\right)Y_{n}^{m\ast }\frac{\partial Y_{n}^{m}}{\partial \phi }\right. \nonumber\\&\left.\left.+\frac{|g_{nm}|^{2}}{r}\left(\frac{\partial Y_{n}^{m\ast }}{\partial \theta }\frac{\partial ^{2}Y_{n}^{m}}{\partial \theta \partial \phi }+\frac{1}{\sin ^{2}\theta }\frac{\partial Y_{n}^{m\ast }}{\partial \phi }\frac{\partial ^{2}Y_{n}^{m}}{\partial \phi ^{2}}\right)\right]\right\}. \end{align}

For simplicity, we drop the arguments of the functions in (B2).

The calculation of ${\boldsymbol \nabla }\times \langle \boldsymbol{v}_{1}^{(n,m)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{(n,m)}\rangle$ , after some rearrangements, results in

(B3) \begin{align} &{\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\!\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\!\right)}\right\rangle =\frac{1}{2}\,\mathrm{Re}\left\{\boldsymbol{e}_{r}\frac{f_{nm}^{\ast }\left(g_{nm}+rg_{nm}^{\prime}\right)}{r^{2}\sin \theta }\left(\frac{\partial Y_{n}^{m}}{\partial \phi }\frac{\partial Y_{n}^{m\ast }}{\partial \theta }-\frac{\partial Y_{n}^{m\ast }}{\partial \phi }\frac{\partial Y_{n}^{m}}{\partial \theta }\right)\right.\nonumber\\&-\frac{\left(f_{nm}^{\ast }g_{nm}+rf_{nm}^{\ast }g_{nm}^{\prime}\right)^{\prime}}{r}\left(\frac{\boldsymbol{e}_{\theta }}{\sin \theta }Y_{n}^{m\ast }\frac{\partial Y_{n}^{m}}{\partial \phi }-\boldsymbol{e}_{\phi }Y_{n}^{m\ast }\frac{\partial Y_{n}^{m}}{\partial \theta }\right)\nonumber\\&+\frac{f_{nm}f_{nm}^{\prime\ast }+f_{nm}^{\ast }f_{nm}^{\prime}}{2r}\left[\frac{\boldsymbol{e}_{\theta }}{\sin \theta }\frac{\partial }{\partial \phi }\left(Y_{n}^{m}Y_{n}^{m\ast }\right)-\boldsymbol{e}_{\phi }\frac{\partial }{\partial \theta }\left(Y_{n}^{m}Y_{n}^{m\ast }\right)\right]\nonumber\\&+\frac{1}{2r}\left[\frac{f_{nm}g_{nm}^{\ast }+f_{nm}^{\ast }g_{nm}-2g_{nm}g_{nm}^{\ast }}{r}-\left(g_{nm}g_{nm}^{\ast }\right)^{\prime}\right]\nonumber\\&\times \left[\frac{\boldsymbol{e}_{\theta }}{\sin \theta }\frac{\partial }{\partial \phi }\left(\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial Y_{n}^{m\ast }}{\partial \theta }+\frac{1}{\sin ^{2}\theta }\frac{\partial Y_{n}^{m}}{\partial \phi }\frac{\partial Y_{n}^{m\ast }}{\partial \phi }\right)\right.\nonumber\\&\left.\left.-\boldsymbol{e}_{\phi }\frac{\partial }{\partial \theta }\left(\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial Y_{n}^{m\ast }}{\partial \theta }+\frac{1}{\sin ^{2}\theta }\frac{\partial Y_{n}^{m}}{\partial \phi }\frac{\partial Y_{n}^{m\ast }}{\partial \phi }\right)\right]\right\}. \end{align}

Comparison of (2.52) and (B3) leads to

(B4) \begin{align} &\sum _{k=1}^{\infty }\sum _{l=-k}^{k}k\left(k+1\right)\left[P_{kl}^{\left(n,m\right)\prime\prime}\left(r\right)-\frac{k\left(k+1\right)P_{kl}^{\left(n,m\right)}\left(r\right)}{r^{2}}\right]Y_{k}^{l}\left(\theta ,\phi \right)\nonumber\\&\quad=\frac{1}{2\nu }\,\mathrm{Re}\left\{\frac{f_{nm}^{\ast }\left(g_{nm}+rg_{nm}^{\prime}\right)}{\sin \theta }\left(\frac{\partial Y_{n}^{m}}{\partial \phi }\frac{\partial Y_{n}^{m\ast }}{\partial \theta }-\frac{\partial Y_{n}^{m\ast }}{\partial \phi }\frac{\partial Y_{n}^{m}}{\partial \theta }\right)\right\}. \end{align}

With the help of (D8) and (D21), the right-hand side of (B4) is transformed so that (B4) takes the form

(B5) \begin{align} &\sum _{k=1}^{\infty }\sum _{l=-k}^{k}k\left(k+1\right)\left[P_{kl}^{\left(n,m\right)\prime\prime}\left(r\right)-\frac{k\left(k+1\right)P_{kl}^{\left(n,m\right)}\left(r\right)}{r^{2}}\right]Y_{k}^{l}\left(\theta ,\phi \right)\nonumber\\&\quad=\frac{1}{2\nu }\,\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}imf_{nm}^{\ast }\left(g_{nm}+rg_{nm}^{\prime}\right)\times \left[nC_{\left(n+1\right)m}\left(B_{kl}^{\left(n+1\right)mnm}+B_{kl}^{nm\left(n+1\right)m}\right)\right.\nonumber\\&\qquad\left.-\left(n+1\right)C_{nm}\left(B_{kl}^{\left(n-1\right)mnm}+B_{kl}^{nm\left(n-1\right)m}\right)\right]Y_{k}^{l}\left(\theta ,\phi \right), \end{align}

where $C_{nm}$ and $B_{kl}^{n_{1}m_{1}n_{2}m_{2}}$ are constants that are calculated by (D11) and (D22), respectively. In fact, (B5) involves only terms with $l=0$ , as $m_{1}=m_{2} B_{kl}^{n_{1}m_{1}n_{2}m_{2}}=0$ for $l\neq 0$ . This fact is in agreement with the right-hand side of (B4), which is independent of $\phi$ .

Keeping only non-zero terms, one obtains from (B5),

(B6) \begin{equation} P_{k0}^{\left(n,m\right)\prime\prime}\left(r\right)-\frac{k\left(k+1\right)}{r^{2}}P_{k0}^{\left(n,m\right)}\left(r\right)=F_{k}^{\left(n,m\right)}\left(r\right) ,\end{equation}
(B7) \begin{align} &F_{k}^{\left(n,m\right)}\left(r\right)=\frac{m}{2k\left(k+1\right)\nu }\,\mathrm{Re}\left\{if_{nm}^{\ast }\left(r\right)\left[g_{nm}\left(r\right)+rg_{nm}^{\prime}\left(r\right)\right]\right\}\nonumber\\&\quad\times \left[nC_{\left(n+1\right)m}\left(B_{k0}^{\left(n+1\right)mnm}+B_{k0}^{nm\left(n+1\right)m}\right)-\left(n+1\right)C_{nm}\left(B_{k0}^{\left(n-1\right)mnm}+B_{k0}^{nm\left(n-1\right)m}\right)\right]. \end{align}

The calculation of the $r$ -component of ${\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \langle \boldsymbol{v}_{1}^{(n,m)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{(n,m)}\rangle$ gives

(B8) \begin{align} &\boldsymbol{e}_{r}\cdot \left[{\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle \right]=\mathrm{Re}\left\{\frac{f_{nm}f_{nm}^{\prime\ast }+f_{nm}^{\ast }f_{nm}^{\prime}}{4r^{2}}L^{2}\left[Y_{n}^{m}Y_{n}^{m\ast }\right]\right.\nonumber\\&+\frac{f_{nm}g_{nm}^{\ast }+f_{nm}^{\ast }g_{nm}-2g_{nm}g_{nm}^{\ast }-r\left(g_{nm}g_{nm}^{\ast }\right)^{\prime}}{4r^{3}}L^{2}\left[\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial Y_{n}^{m\ast }}{\partial \theta }+\frac{m^{2}Y_{n}^{m}Y_{n}^{m\ast }}{\sin ^{2}\theta }\right]\nonumber\\&\left.+\frac{\left(f_{nm}^{\ast }g_{nm}+rf_{nm}^{\ast }g_{nm}^{\prime}\right)^{{\prime}}}{2r^{2}}\left[\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial Y_{n}^{m\ast }}{\partial \theta }+\frac{m^{2}Y_{n}^{m}Y_{n}^{m\ast }}{\sin ^{2}\theta }-n\left(n+1\right)Y_{n}^{m}Y_{n}^{m\ast }\right]\right\}, \end{align}

where $L^{2}$ stands for the square of the orbital angular momentum operator (Varshalovich, Moskalev & Khersonskii Reference Varshalovich, Moskalev and Khersonskii1988),

(B9) \begin{equation} L^{2}=-\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\sin \theta\, \frac{\partial }{\partial \theta }-\frac{1}{\sin ^{2}\theta }\frac{\partial ^{2}}{\partial \phi ^{2}} .\end{equation}

By using (D5), it can be shown that the following identity holds:

(B10) \begin{equation} -\frac{1}{2}L^{2}\left[Y_{n}^{m}Y_{n}^{m\ast }\right]=\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial Y_{n}^{m\ast }}{\partial \theta }+\frac{m^{2}Y_{n}^{m}Y_{n}^{m\ast }}{\sin ^{2}\theta }-n\left(n+1\right)Y_{n}^{m}Y_{n}^{m\ast } .\end{equation}

Substitution of (B10) into (B8) yields

(B11) \begin{align} &\boldsymbol{e}_{r}\cdot \left[{\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle \right]\nonumber\\&\quad=\frac{1}{8r^{3}}\,\mathrm{Re}\left\{L^{4}\left[Y_{n}^{m}Y_{n}^{m\ast }\right]\left[2g_{nm}g_{nm}^{\ast }-f_{nm}g_{nm}^{\ast }-f_{nm}^{\ast }g_{nm}+r\left(g_{nm}g_{nm}^{\ast }\right)^{\prime}\right]\right.\nonumber\\&\qquad+2L^{2}\left[Y_{n}^{m}Y_{n}^{m\ast }\right]\left\{r\left[f_{nm}f_{nm}^{\prime\ast }+f_{nm}^{\ast }f_{nm}^{\prime}-\left(f_{nm}^{\ast }g_{nm}+rf_{nm}^{\ast }g_{nm}^{\prime}\right)^{{\prime}}\right]\right.\nonumber\\&\qquad\left.\left. +n\left(n+1\right)\left[f_{nm}g_{nm}^{\ast }+f_{nm}^{\ast }g_{nm}-2g_{nm}g_{nm}^{\ast }-r\left(g_{nm}g_{nm}^{\ast }\right)^{\prime}\right]\right\}\right\}. \end{align}

By using (D17) and the identity $L^{2}Y_{k}^{l}=k(k+1)Y_{k}^{l}$ , one obtains

(B12) \begin{align} &\boldsymbol{e}_{r}\cdot \left[{\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle \right]=\frac{1}{8r^{3}}\,\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}k\left(k+1\right)A_{kl}^{nmnm}Y_{k}^{l}\left(\theta ,\phi \right)\nonumber\\&\quad\times \left\{2r\left[f_{nm}f_{nm}^{\prime\ast }+f_{nm}^{\ast }f_{nm}^{\prime}-\left(f_{nm}^{\ast }g_{nm}+rf_{nm}^{\ast }g_{nm}^{\prime}\right)^{\prime}\right]\right.\nonumber\\&\quad\left.+\left[k\left(k+1\right)-2n\left(n+1\right)\right]\left[2g_{nm}g_{nm}^{\ast }-f_{nm}g_{nm}^{\ast }-f_{nm}^{\ast }g_{nm}+r\left(g_{nm}g_{nm}^{\ast }\right)^{\prime}\right]\right\}, \end{align}

where the constant coefficients $A_{kl}^{nmnm}$ are calculated by (D20). In fact, (B12) involves only terms with $l=0$ as $A_{kl}^{nmnm}=0$ for $l\neq 0$ .

Comparing (2.53) with (B12), and keeping only non-zero terms, one obtains

(B13) \begin{equation} T_{k0}^{\left(n,m\right)\prime\prime}\left(r\right)-\frac{k\left(k+1\right)}{r^{2}}T_{k0}^{\left(n,m\right)}\left(r\right)=G_{k}^{\left(n,m\right)}\left(r\right) ,\end{equation}
(B14) \begin{align} G_{k}^{\left(n,m\right)}\left(r\right)&=\frac{A_{k0}^{nmnm}}{4\nu r}\,\mathrm{Re}\left\{rf_{nm}^{\ast }\left(r\right)\left[2f_{nm}^{\prime}\left(r\right)-2g_{nm}^{\prime}\left(r\right)-rg_{nm}^{\prime\prime}\left(r\right)\right]\right.\nonumber\\&\quad-rf_{nm}^{\prime\ast }\left(r\right)\left[g_{nm}\left(r\right)+rg_{nm}^{\prime}\left(r\right)\right]\nonumber\\&\quad\left.+\left[k\left(k+1\right)-2n\left(n+1\right)\right]g_{nm}^{\ast }\left(r\right)\left[g_{nm}\left(r\right)+rg_{nm}^{\prime}\left(r\right)-f_{nm}\left(r\right)\right]\right\}. \end{align}

The derivatives that appear in (B14) are calculated by

(B15) \begin{equation} f_{nm}^{\prime}\left(r\right)=-\frac{\left(n+1\right)\left(n+2\right)c_{nm}}{r^{n+3}}+\frac{n\left(n+1\right)b_{nm}}{k_{v}r^{2}}\left[\left(n-1\right)h_{n}^{\left(1\right)}\left(k_{v}r\right)-k_{v}rh_{n+1}^{\left(1\right)}\left(k_{v}r\right)\right] ,\end{equation}
(B16) \begin{equation} g_{nm}^{\prime}\left(r\right)=\frac{\left(n+2\right)c_{nm}}{r^{n+3}}+\frac{b_{nm}}{k_{v}r^{2}}\left\{\left[n^{2}-1-\left(k_{v}r\right)^{2}\right]h_{n}^{\left(1\right)}\left(k_{v}r\right)+k_{v}rh_{n+1}^{\left(1\right)}\left(k_{v}r\right)\right\} ,\end{equation}
(B17) \begin{align} g_{nm}^{\prime\prime}\left(r\right)&=-\frac{\left(n+2\right)\left(n+3\right)c_{nm}}{r^{n+4}}+\frac{b_{nm}}{k_{v}r^{3}}\left\{2\left(1-n^{2}\right)h_{n}^{\left(1\right)}\left(k_{v}r\right)\right.\nonumber\\&\quad\left.+k_{v}r\left[n^{2}-1-\left(k_{v}r\right)^{2}\right]h_{n}^{\left(1\right)\prime}\left(k_{v}r\right)-k_{v}rh_{n+1}^{\left(1\right)}\left(k_{v}r\right)+\left(k_{v}r\right)^{2}h_{n+1}^{\left(1\right)\prime}\left(k_{v}r\right)\right\}. \end{align}

Equations (B6) and (B13) are solved by the method of variation of parameters (Boyce & DiPrima Reference Boyce and DiPrima2001). According to this method, a solution to (B6) is given by

(B18) \begin{equation} P_{k0}^{\left(n,m\right)}\left(r\right)=r^{k+1}C_{1k}^{\left(n,m\right)}\left(r\right)+r^{-k}C_{2k}^{\left(n,m\right)}\left(r\right) ,\end{equation}

where $C_{1k}^{(n,m)}(r)$ and $C_{2k}^{(n,m)}(r)$ obey the equations

(B19) \begin{align} r^{k+1}C_{1k}^{\left(n,m\right)\prime}\left(r\right)+r^{-k}C_{2k}^{\left(n,m\right)\prime}\left(r\right)=0 ,\qquad\qquad\end{align}
(B20) \begin{align} \left(k+1\right)r^{k}C_{1k}^{\left(n,m\right)\prime}\left(r\right)-kr^{-k-1}C_{2k}^{\left(n,m\right)\prime}\left(r\right)=F_{k}^{\left(n,m\right)}\left(r\right) . \end{align}

Solving (B19) and (B20) gives

(B21) \begin{align} C_{1k}^{\left(n,m\right)}\left(r\right)=\overline{C}_{1k}^{\left(n,m\right)}+\frac{1}{2k+1}\int _{R_{0}}^{r}s^{-k}F_{k}^{\left(n,m\right)}\left(s\right){\rm d}s ,\end{align}
(B22) \begin{align} C_{2k}^{\left(n,m\right)}\left(r\right)=\overline{C}_{2k}^{\left(n,m\right)}-\frac{1}{2k+1}\int _{R_{0}}^{r}s^{k+1}F_{k}^{\left(n,m\right)}\left(s\right){\rm d}s ,\end{align}

where $\overline{C}_{1k}^{(n,m)}$ and $\overline{C}_{2k}^{(n,m)}$ are constants. From the condition that the acoustic streaming vanishes at infinity, one finds

(B23) \begin{equation} \overline{C}_{1k}^{\left(n,m\right)}=-\frac{1}{2k+1}\int _{R_{0}}^{\infty }s^{-k}F_{k}^{\left(n,m\right)}\left(s\right){\rm d}s ,\end{equation}

while $\overline{C}_{2k}^{(n,m)}$ is calculated by boundary conditions at the bubble surface; see below.

Equation (B13) is similar to (B6), so its solution can be written by analogy,

(B24) \begin{align} T_{k0}^{\left(n,m\right)}\left(r\right)=r^{k+1}C_{3k}^{\left(n,m\right)}\left(r\right)+r^{-k}C_{4k}^{\left(n,m\right)}\left(r\right) , \quad\end{align}
(B25) \begin{align} C_{3k}^{\left(n,m\right)}\left(r\right)=\overline{C}_{3k}^{\left(n,m\right)}+\frac{1}{2k+1}\int _{R_{0}}^{r}s^{-k}G_{k}^{\left(n,m\right)}\left(s\right){\rm d}s , \end{align}
(B26) \begin{align} C_{4k}^{\left(n,m\right)}\left(r\right)=\overline{C}_{4k}^{\left(n,m\right)}-\frac{1}{2k+1}\int _{R_{0}}^{r}s^{k+1}G_{k}^{\left(n,m\right)}\left(s\right){\rm d}s , \end{align}
(B27) \begin{align} \overline{C}_{3k}^{\left(n,m\right)}=-\frac{1}{2k+1}\int _{R_{0}}^{\infty }s^{-k}G_{k}^{\left(n,m\right)}\left(s\right){\rm d}s , \qquad \end{align}

where $\overline{C}_{4k}^{(n,m)}$ is a constant to be calculated by boundary conditions at the bubble surface; see below.

Continuing the calculation of $\boldsymbol{v}_{E}^{(n,m)}$ , one obtains from (2.50),

(B28) \begin{align} \boldsymbol{v}_{E}^{\left(n,m\right)}&={\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}P_{kl}^{\left(n,m\right)}\left(r\right)Y_{k}^{l}\left(\theta ,\phi \right)\right]\nonumber\\&\quad+\boldsymbol{e}_{r}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}T_{kl}^{\left(n,m\right)}\left(r\right)Y_{k}^{l}\left(\theta ,\phi \right)+{\boldsymbol \nabla }\Phi ^{\left(n,m\right)}\left(r,\theta ,\phi \right) .\end{align}

It follows from (2.44) that

(B29) \begin{align} \Delta \Phi ^{(n,m)}(r,\theta ,\phi )&=-\displaystyle\sum _{k=1}^{\infty }\sum _{l=-k}^{k}{\boldsymbol \nabla }\cdot [T_{kl}^{(n,m)}(r)Y_{k}^{l}(\theta ,\phi )\boldsymbol{e}_{r}]\nonumber\\&=-\sum _{k=1}^{\infty }\sum _{l=-k}^{k}[T_{kl}^{(n,m)/}(r)+2r^{-1}T_{kl}^{(n,m)}(r)]Y_{k}^{l}(\theta ,\phi ). \end{align}

In fact, (B29) involves only terms with $l=0$ as $T_{kl}^{(n,m)}(r)=0$ for $l\neq 0$ .

We assume that

(B30) \begin{equation} \Phi ^{\left(n,m\right)}\left(r,\theta ,\phi \right)=\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\Phi _{kl}^{\left(n,m\right)}\left(r\right)Y_{k}^{l}\left(\theta ,\phi \right) .\end{equation}

Substituting (B30) into (B29), using (B24) and keeping only non-zero terms, one has

(B31) \begin{equation} \Phi _{k0}^{\left(n,m\right)\prime\prime}\left(r\right)+\frac{2}{r}\Phi _{k0}^{\left(n,m\right)\prime}\left(r\right)-\frac{k\left(k+1\right)}{r^{2}}\Phi _{k0}^{\left(n,m\right)}\left(r\right)=H_{k}^{\left(n,m\right)}\left(r\right) ,\end{equation}

where

(B32) \begin{align} H_{k}^{\left(n,m\right)}\left(r\right)&=-T_{k0}^{\left(n,m\right)\prime}\left(r\right)-2r^{-1}T_{k0}^{\left(n,m\right)}\left(r\right)\nonumber\\&=-\left(k+3\right)r^{k}C_{3k}^{\left(n,m\right)}\left(r\right)+\left(k-2\right)r^{-k-1}C_{4k}^{\left(n,m\right)}\left(r\right) .\end{align}

Equation (B31) is solved by the method of variation of parameters, which results in

(B33) \begin{equation} \Phi _{k0}^{\left(n,m\right)}\left(r\right)=r^{k}C_{5k}^{\left(n,m\right)}\left(r\right)+r^{-k-1}C_{6k}^{\left(n,m\right)}\left(r\right) ,\end{equation}

where $C_{5k}^{(n,m)}(r)$ and $C_{6k}^{(n,m)}(r)$ obey the equations

(B34) \begin{align} r^{k}C_{5k}^{\left(n,m\right)^{\prime}}\left(r\right)+r^{-k-1}C_{6k}^{\left(n,m\right)\prime}\left(r\right)=0 ,\qquad\qquad\end{align}
(B35) \begin{align} kr^{k-1}C_{5k}^{\left(n,m\right)\prime}\left(r\right)-\left(k+1\right)r^{-k-2}C_{6k}^{\left(n,m\right)\prime}\left(r\right)=H_{k}^{\left(n,m\right)}\left(r\right) . \end{align}

Solving (B34) and (B35) gives

(B36) \begin{align} C_{5k}^{\left(n,m\right)}\left(r\right)=\overline{C}_{5k}^{\left(n,m\right)}+\frac{1}{2k+1}\int _{R_{0}}^{r}s^{1-k}H_{k}^{\left(n,m\right)}\left(s\right){\rm d}s ,\end{align}
(B37) \begin{align} C_{6k}^{\left(n,m\right)}\left(r\right)=\overline{C}_{6k}^{\left(n,m\right)}-\frac{1}{2k+1}\int _{R_{0}}^{r}s^{k+2}H_{k}^{\left(n,m\right)}\left(s\right){\rm d}s .\end{align}

The constant $\overline{C}_{6k}^{(n,m)}$ is calculated by boundary conditions at the bubble surface; see below.

From the condition that the acoustic streaming vanishes at infinity, one finds

(B38) \begin{equation} \overline{C}_{5k}^{\left(n,m\right)}=-\frac{1}{2k+1}\int _{R_{0}}^{\infty }s^{1-k}H_{k}^{\left(n,m\right)}\left(s\right){\rm d}s .\end{equation}

Substituting (B32) into (B38) and using (B26), one obtains

(B39) \begin{align} \overline{C}_{5k}^{\left(n,m\right)}=-\frac{\left(k-2\right)\overline{C}_{4k}^{\left(n,m\right)}}{\left(4k^{2}-1\right)R_{0}^{2k-1}}+I_{k}^{\left(n,m\right)} ,\end{align}
(B40) \begin{align} I_{k}^{\left(n,m\right)}=\frac{k+3}{2k+1}\int _{R_{0}}^{\infty }rC_{3k}^{\left(n,m\right)}\left(r\right){\rm d}r+\frac{k-2}{\left(2k+1\right)^{2}}\int _{R_{0}}^{\infty }\left[\int _{R_{0}}^{r}s^{k+1}G_{k}^{\left(n,m\right)}\left(s\right){\rm d}s\right]\frac{{\rm d}r}{r^{2k}} .\end{align}

Substitution of (B30) into (B28) yields

(B41) \begin{align} &\boldsymbol{v}_{E}^{\left(n,m\right)}=\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\left\{\boldsymbol{e}_{r}\left[T_{kl}^{\left(n,m\right)}\left(r\right)+\Phi _{kl}^{\left(n,m\right)\prime}\left(r\right)\right]Y_{k}^{l}\left(\theta ,\phi \right)\right.\nonumber\\&+\frac{\Phi _{kl}^{\left(n,m\right)}\left(r\right)}{r}\left[\boldsymbol{e}_{\theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }+\frac{\boldsymbol{e}_{\phi }}{\sin \theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \phi }\right]\nonumber\\&\left.+\frac{P_{kl}^{\left(n,m\right)}\left(r\right)}{r}\left[\frac{\boldsymbol{e}_{\theta }}{\sin \theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \phi }-\boldsymbol{e}_{\phi }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }\right]\right\}. \end{align}

Keeping only non-zero terms (those with $l=0$ ) and using (D1) and (D9), one obtains

(B42) \begin{align} v_{Er}^{\left(n,m\right)}=\frac{1}{2\sqrt{\pi}}\sum _{k=1}^{\infty }\sqrt{2k+1}\left[T_{k0}^{\left(n,m\right)}\left(r\right)+\Phi _{k0}^{\left(n,m\right)/}\left(r\right)\right]P_{k}\left(\cos \theta \right) , \end{align}
(B43) \begin{align} v_{E\theta }^{\left(n,m\right)}=\frac{1}{2\sqrt{\pi}\,r}\sum _{k=1}^{\infty }\sqrt{2k+1}\,\Phi _{k0}^{\left(n,m\right)}\left(r\right)P_{k}^{1}\left(\cos \theta \right) ,\qquad\quad\end{align}
(B44) \begin{align} v_{E\phi }^{\left(n,m\right)}=-\frac{1}{2\sqrt{\pi}\,r}\sum _{k=1}^{\infty }\sqrt{2k+1}\,P_{k0}^{\left(n,m\right)}\left(r\right)P_{k}^{1}\left(\cos \theta \right) .\qquad \end{align}

Equation (B42) involves $\Phi _{k0}^{(n,m)/}(r)$ , which is calculated by (B33) and (B34) to be

(B45) \begin{equation} \Phi _{k0}^{\left(n,m\right)/}\left(r\right)=kr^{k-1}C_{5k}^{\left(n,m\right)}\left(r\right)-\left(k+1\right)r^{-k-2}C_{6k}^{\left(n,m\right)}\left(r\right). \end{equation}

In order to go on with the calculation, we need to apply boundary conditions for acoustic streaming at the bubble surface. To do this, we need to know the Stokes drift velocity (Longuet-Higgins Reference Longuet-Higgins1998), which is calculated by (Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a )

(B46) \begin{equation} \boldsymbol{v}_{S}^{\left(n,m\right)}=\frac{1}{2\omega }\,\mathrm{Re}\left\{i\boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)\ast }\right\}. \end{equation}

Equation (B46) gives

(B47) \begin{equation} v_{Sr}^{\left(n,m\right)}=-\frac{1}{2\omega }\,\mathrm{Im}\left\{v_{1r}^{\left(n,m\right)}\frac{\partial v_{1r}^{\left(n,m\right)\ast }}{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial v_{1r}^{\left(n,m\right)\ast }}{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial v_{1r}^{\left(n,m\right)\ast }}{\partial \phi }\right\}, \end{equation}
(B48) \begin{equation} v_{S\theta }^{\left(n,m\right)}=-\frac{1}{2\omega }\,\mathrm{Im}\!\left\{v_{1r}^{\left(n,m\right)}\frac{\partial v_{1\theta }^{\left(n,m\right)\ast }}{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial v_{1\theta }^{\left(n,m\right)\ast }}{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial v_{1\theta }^{\left(n,m\right)\ast }}{\partial \phi }+\frac{v_{1r}^{\left(n,m\right)\ast }v_{1\theta }^{\left(n,m\right)}}{r}\!\right\}\!, \end{equation}
(B49) \begin{align} v_{S\phi }^{\left(n,m\right)}&=-\frac{1}{2\omega }\,\mathrm{Im}\left\{v_{1r}^{\left(n,m\right)}\frac{\partial v_{1\phi }^{\left(n,m\right)\ast }}{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial v_{1\phi }^{\left(n,m\right)\ast }}{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial v_{1\phi }^{\left(n,m\right)\ast }}{\partial \phi }+\frac{v_{1\phi }^{\left(n,m\right)}v_{1r}^{\left(n,m\right)\ast }}{r}\right.\nonumber\\&\quad\left.+\frac{\cos \theta\, v_{1\phi }^{\left(n,m\right)}v_{1\theta }^{\left(n,m\right)\ast }}{r\sin \theta }\right\}. \end{align}

Substituting (A5), (A16) and (A17) into (B47)–(B49) and using (D6), one obtains

(B50) \begin{equation} v_{Sr}^{\left(n,m\right)}=-\frac{1}{2\omega }\,\mathrm{Re}\left\{if_{nm}^{\ast }f_{nm}^{\prime}Y_{n}^{m}Y_{n}^{m\ast }+\frac{if_{nm}g_{nm}^{\ast }}{r}\left(\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial Y_{n}^{m\ast }}{\partial \theta }+\frac{m^{2}Y_{n}^{m}Y_{n}^{m\ast }}{\sin ^{2}\theta }\right)\right\} ,\end{equation}
(B51) \begin{equation} v_{S\theta }^{\left(n,m\right)}=-\frac{1}{4\omega }\,\mathrm{Re}\left\{i\left[f_{nm}^{\ast }g_{nm}^{\prime}-\frac{\left(f_{nm}^{\ast }-n\left(n+1\right)g_{nm}^{\ast }\right)g_{nm}}{r}\right]\frac{\partial \left(Y_{n}^{m}Y_{n}^{m\ast }\right)}{\partial \theta }\right\} ,\end{equation}
(B52) \begin{align} v_{S\phi }^{\left(n,m\right)}&=\frac{m}{2\omega }\,\mathrm{Re}\left\{f_{nm}^{\ast }\left(g_{nm}^{\prime}-\frac{g_{nm}}{r}\right)\frac{Y_{n}^{m}Y_{n}^{m\ast }}{\sin \theta }\right.\nonumber\\&\quad\left.+\frac{g_{nm}g_{nm}^{\ast }}{r\sin \theta }\left(\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial Y_{n}^{m\ast }}{\partial \theta }+\frac{m^{2}Y_{n}^{m}Y_{n}^{m\ast }}{\sin ^{2}\theta }-\frac{2\cos \theta }{\sin \theta }Y_{n}^{m}\frac{\partial Y_{n}^{m\ast }}{\partial \theta }\right)\right\}. \end{align}

With the help of (B10) and (D17), (B50) is represented by

(B53) \begin{align} v_{Sr}^{\left(n,m\right)}&=-\frac{1}{2\omega }\,\mathrm{Re}\left\{i\left[f_{nm}^{\ast }f_{nm}^{\prime}+\frac{n\left(n+1\right)f_{nm}g_{nm}^{\ast }}{r}\right]Y_{n}^{m}Y_{n}^{m\ast }-\frac{if_{nm}g_{nm}^{\ast }}{2r}L^{2}\left[Y_{n}^{m}Y_{n}^{m\ast }\right]\right\}\nonumber\\&=\frac{1}{2\omega }\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}iA_{kl}^{nmnm}\left[\frac{k\left(k+1\right)-2n\left(n+1\right)}{2r}f_{nm}g_{nm}^{\ast }-f_{nm}^{\ast }f_{nm}^{\prime}\right]Y_{k}^{l}\left(\theta ,\phi \right). \end{align}

Keeping only non-zero terms and using (D1), one obtains

(B54) \begin{align} v_{Sr}^{\left(n,m\right)}&=\frac{1}{4\sqrt{\pi}\,\omega }\sum _{k=1}^{\infty }\sqrt{2k+1}\,A_{k0}^{nmnm}\nonumber\\&\quad\times\mathrm{Re}\left\{if_{nm}^{\ast }\left(r\right)\left[\frac{2n\left(n+1\right)-k\left(k+1\right)}{2r}g_{nm}\left(r\right)-f_{nm}^{\prime}\left(r\right)\right]\right\}P_{k}\left(\cos \theta \right). \end{align}

Substitution of (D17) into (B51) yields

(B55) \begin{equation} v_{S\theta }^{\left(n,m\right)}=\frac{1}{4\omega }\,\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}iA_{kl}^{nmnm}\left\{\frac{g_{nm}\left[f_{nm}^{\ast }-n\left(n+1\right)g_{nm}^{\ast }\right]}{r}-f_{nm}^{\ast }g_{nm}^{\prime}\right\}\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta } .\end{equation}

Keeping only non-zero terms and using (D1) and (D9), one obtains

(B56) \begin{equation} v_{S\theta }^{\left(n,m\right)}=\frac{1}{8\sqrt{\pi}\omega }\,\mathrm{Re}\left\{if_{nm}^{\ast }\left(r\right)\left[\frac{g_{nm}\left(r\right)}{r}-g_{nm}^{\prime}\left(r\right)\right]\right\}\sum _{k=1}^{\infty }\sqrt{2k+1}\,A_{k0}^{nmnm}P_{k}^{1}\left(\cos \theta \right) .\end{equation}

With the help of (B10), (B52) is represented by

(B57) \begin{align} v_{S\phi }^{\left(n,m\right)}&=\frac{m}{2\omega }\,\mathrm{Re}\left\{\left[f_{nm}^{\ast }\left(g_{nm}^{\prime}-\frac{g_{nm}}{r}\right)+\frac{n\left(n+1\right)|g_{nm}|^{2}}{r}\right]\frac{Y_{n}^{m}Y_{n}^{m\ast }}{\sin \theta }\right.\nonumber\\&\quad\left.+\frac{|g_{nm}|^{2}}{r}\frac{\partial }{\partial \theta }\left(\frac{Y_{n}^{m\ast }}{\sin \theta }\frac{\partial Y_{n}^{m}}{\partial \theta }\right)\right\} .\end{align}

Note that $v_{S\phi }^{(n,m)}=0$ for $m=0$ . For $n=m=1$ , with the help of (D1), (B57) is represented by

(B58) \begin{equation} v_{S\phi }^{\left(1,1\right)}=E^{\left(1,1\right)}\left(r\right)P_{1}^{1}\left(\cos \theta \right) ,\end{equation}
(B59) \begin{equation} E^{\left(1,1\right)}\left(r\right)=\frac{3}{16{\pi} \omega r}\,\mathrm{Re}\left\{f_{11}^{\ast }\left(r\right)\left[g_{11}\left(r\right)-rg_{11}^{\prime}\left(r\right)\right]-|g_{11}\left(r\right)|^{2}\right\} .\end{equation}

For $n\gt 1$ , by using (D23), (D24), (D1), (D3) and (D15), keeping only non-zero terms, (B57) is represented by

(B60) \begin{equation} v_{S\phi }^{\left(n,m\right)}=\sum _{k=1}^{\infty }E_{k}^{\left(n\gt 1,m\right)}\left(r\right)P_{k}^{1}\left(\cos \theta \right) ,\end{equation}

where

(B61) \begin{align}& E_{k}^{(n \gt 1,m)} (r) = \frac{m\sqrt{2k+1}}{4\sqrt{\pi\omega r}} \left\{\left[nC_{(n+1)m}B_{k0}^{(n+1)mnm} - (n+1)C_{nm}B_{k0}^{(n+1)mnm}\right]|g_{nm}(r)|^{2}\vphantom{\sum _{s=1}^{\left[\frac{(n-m+2)}{2}\right]}}\right.\nonumber\\&-\sqrt{\frac{(2n+1)(n-m)!}{(n+m)!}}\,\mathrm{Re} \left\{f_{nm}^{*}(r)\left[g_{nm}(r)-rg_{nm}^{\prime}(r)\right]-n(n+1)|g_{nm}(r)|^{2}\right\}\nonumber\\&\left.\times \frac{1}{\sqrt{k(k+1)}}\sum _{s=1}^{\left[\frac{(n-m+2)}{2}\right]}\sqrt{\frac{(2n-4s+3)(n+m-2s)!}{(n-m-2s+2)!}}A_{k(-1)}^{(n-2s+1)(m-1)nm}\right\}.\end{align}

Equations (B58) and (B60) can be combined as follows:

(B62) \begin{align} v_{S\phi }^{\left(n,m\right)}=\sum _{k=1}^{\infty }E_{k}^{\left(n,m\right)}\left(r\right)P_{k}^{1}\left(\cos \theta \right) ,\end{align}

where

(B63) \begin{align} E_{k}^{\left(n,m\right)}\left(r\right)=\left\{\begin{array}{ll} E^{\left(1,1\right)}\left(r\right)\delta _{1k}, & n=1,\\ E_{k}^{\left(n\gt 1,m\right)}\left(r\right), & n\gt 1. \end{array}\right. \end{align}

We can now apply the boundary conditions for the acoustic streaming at the bubble surface. They are given by

(B64) \begin{align} v_{Lr}^{\left(n,m\right)}=0\quad\text{at }r=R_{0} ,\qquad\qquad\qquad\qquad\end{align}
(B65) \begin{align} \sigma _{Lr\theta }^{\left(n,m\right)}=\eta \left(\frac{1}{r}\frac{\partial v_{Lr}^{\left(n,m\right)}}{\partial \theta }+\frac{\partial v_{L\theta }^{\left(n,m\right)}}{\partial r}-\frac{v_{L\theta }^{\left(n,m\right)}}{r}\right)=0\quad\text{at }r=R_{0} ,\quad \end{align}
(B66) \begin{align} \sigma _{Lr\phi }^{\left(n,m\right)}=\eta \left(\frac{\partial v_{L\phi }^{\left(n,m\right)}}{\partial r}+\frac{1}{r\sin \theta }\frac{\partial v_{Lr}^{\left(n,m\right)}}{\partial \phi }-\frac{v_{L\phi }^{\left(n,m\right)}}{r}\right)=0\quad\text{at }r=R_{0} ,\end{align}

where $\boldsymbol{v}_{L}^{(n,m)}=\boldsymbol{v}_{E}^{(n,m)}+\boldsymbol{v}_{S}^{(n,m)}$ is the Lagrangian streaming velocity, and $\sigma _{Lr\theta }^{(n,m)}$ and $\sigma _{Lr\phi }^{(n,m)}$ are the tangential components of the stress produced by $\boldsymbol{v}_{L}^{(n,m)}$ . We use (B64)–(B66) in order to calculate the constants $\overline{C}_{2k}^{(n,m)}$ , $\overline{C}_{4k}^{(n,m)}$ and $\overline{C}_{6k}^{(n,m)}$ .

Substituting (B42) and (B54) into (B64), one obtains

(B67) \begin{align} T_{k0}^{\left(n,m\right)}\left(R_{0}\right)+\Phi _{k0}^{\left(n,m\right)\prime}\left(R_{0}\right)=Q_{k}^{\left(n,m\right)} ,\end{align}

where

(B68) \begin{align} Q_{k}^{\left(n,m\right)}=\frac{A_{k0}^{nmnm}}{2\omega }\,\mathrm{Re}\left\{if_{nm}^{\ast }\left(R_{0}\right)\left[f_{nm}^{\prime}\left(R_{0}\right)+\frac{k\left(k+1\right)-2n\left(n+1\right)}{2R_{0}}g_{nm}\left(R_{0}\right)\right]\right\} . \end{align}

To calculate $\sigma _{Lr\theta }^{(n,m)}$ , we use (B42), (B43), (B54) and (B56):

(B69) \begin{align} &\sigma _{Lr\theta }^{\left(n,m\right)}=\frac{\eta }{2\sqrt{\pi}}\sum _{k=1}^{\infty }\sqrt{2k+1}\,P_{k}^{1}\left(\cos \theta \right)\left\{\frac{T_{k0}^{\left(n,m\right)}\left(r\right)}{r}+\frac{2\Phi _{k0}^{\left(n,m\right)\prime}\left(r\right)}{r}-\frac{2\Phi _{k0}^{\left(n,m\right)}\left(r\right)}{r^{2}}\right.\nonumber\\&+\frac{A_{k0}^{nmnm}}{4\omega r}\,\mathrm{Re}\left\{if_{nm}^{\prime\ast }\left(r\right)\left[g_{nm}\left(r\right)-rg_{nm}^{\prime}\left(r\right)\right]\right.\nonumber\\&\left.\left.+if_{nm}^{\ast }\left(r\right)\left[\frac{2n\left(n+1\right)-k\left(k+1\right)-2}{r}g_{nm}\left(r\right)+2g_{nm}^{\prime}\left(r\right)-rg_{nm}^{\prime\prime}\left(r\right)-2f_{nm}^{\prime}\left(r\right)\right]\right\}\right\}. \end{align}

Substitution of (B69) into (B65) yields

(B70) \begin{equation} T_{k0}^{\left(n,m\right)}\left(R_{0}\right)+2\Phi _{k0}^{\left(n,m\right)\prime}\left(R_{0}\right)-\frac{2}{R_{0}}\Phi _{k0}^{\left(n,m\right)}\left(R_{0}\right)=S_{k}^{\left(n,m\right)} ,\end{equation}

where

(B71) \begin{align} &S_{k}^{\left(n,m\right)}=\frac{A_{k0}^{nmnm}}{4\omega }\,\mathrm{Re}\left\{if_{nm}^{/\ast }\left(R_{0}\right)\left[R_{0}g_{nm}^{\prime}\left(R_{0}\right)-g_{nm}\left(R_{0}\right)\right]+if_{nm}^{\ast }\left(R_{0}\right)\right.\nonumber\\&\left.\times\left[\frac{2+k\left(k+1\right)-2n\left(n+1\right)}{R_{0}}g_{nm}\left(R_{0}\right)-2g_{nm}^{\prime}\left(R_{0}\right)+R_{0}g_{nm}^{//}\left(R_{0}\right)+2f_{nm}^{\prime}\left(R_{0}\right)\right]\right\}. \end{align}

Combining (B67) and (B70), and using (B24), (B33), (B39) and (B45), one finds

(B72) \begin{align} \overline{C}_{4k}^{\left(n,m\right)}&=\frac{\left(2k-1\right)R_{0}^{k}}{\left(k+1\right)\left(2k+1\right)}\left[2\left(k+2\right)Q_{k}^{\left(n,m\right)}-\left(k+1\right)S_{k}^{\left(n,m\right)}-\left(k+3\right)R_{0}^{k+1}\overline{C}_{3k}^{\left(n,m\right)}\right.\nonumber\\&\quad\left.-2\left(2k+1\right)R_{0}^{k-1}I_{k}^{\left(n,m\right)}\right], \end{align}
(B73) \begin{equation} \overline{C}_{6k}^{\left(n,m\right)}=\frac{R_{0}^{2}}{k+1}\left[\frac{\left(k+1\right)\left(3k-1\right)}{4k^{2}-1}\overline{C}_{4k}^{\left(n,m\right)}+R_{0}^{2k+1}\overline{C}_{3k}^{\left(n,m\right)}+kR_{0}^{2k-1}I_{k}^{\left(n,m\right)}-R_{0}^{k}Q_{k}^{\left(n,m\right)}\right] .\end{equation}

To calculate $\sigma _{Lr\phi }^{(n,m)}$ , we use (B44) and (B62):

(B74) \begin{align} \sigma _{Lr\phi }^{\left(n,m\right)}&=\frac{\eta }{r^{2}}\sum _{k=1}^{\infty }P_{k}^{1}\left(\cos \theta \right)\left\{\frac{\sqrt{2k+1}}{2\sqrt{\pi}}\left[2P_{k0}^{\left(n,m\right)}\left(r\right)-rP_{k0}^{\left(n,m\right)\prime}\left(r\right)\right]\right.\nonumber\\&\quad\left. -rE_{k}^{\left(n,m\right)}\left(r\right)+r^{2}E_{k}^{\left(n,m\right)\prime}\left(r\right)\right\} ,\end{align}

where $P_{k0}^{(n,m)\prime}(r)$ and $E_{k}^{(n,m)\prime}(r)$ are calculated by

(B75) \begin{equation} P_{k0}^{\left(n,m\right)\prime}\left(r\right)=\left(k+1\right)r^{k}C_{1k}^{\left(n,m\right)}\left(r\right)-kr^{-k-1}C_{2k}^{\left(n,m\right)}\left(r\right) ,\end{equation}
(B76) \begin{equation} E_{k}^{\left(n,m\right)\prime}\left(r\right)=\left\{\begin{array}{ll} E^{\left(1,1\right)\prime}\left(r\right)\delta _{1k}, & n=1,\\ E_{k}^{\left(n\gt 1,m\right)\prime}\left(r\right), & n\gt 1, \end{array}\right. \end{equation}
(B77) \begin{align} E^{\left(1,1\right)\prime}\left(r\right)&=\frac{3}{16{\pi} \omega r^{2}}\,\mathrm{Re}\left\{f_{11}^{\ast }\left(r\right)\left[rg_{11}^{\prime}\left(r\right)-r^{2}g_{11}^{\prime\prime}\left(r\right)-g_{11}\left(r\right)\right]\right.\nonumber\\&\quad\left.+rf_{11}^{\prime\ast }\left(r\right)\left[g_{11}\left(r\right)-rg_{11}^{\prime}\left(r\right)\right]+g_{11}^{\ast }\left(r\right)\left[g_{11}\left(r\right)-2rg_{11}^{\prime}\left(r\right)\right]\right\}, \end{align}
(B78) \begin{align} & E_{k}^{(n \gt 1,m)}/ (r) = -\frac{E_{k}^{(n \gt 1,m)}(r)}{r}\nonumber\\& \quad +\frac{m\sqrt{2k+1}}{4\sqrt{\pi\omega r}} \left\{2\ \textrm{Re}\{g_{nm}^{\prime} (r)g_{nm}^{*}(r)\}\left[nC_{(n+1)m}B_{k0}^{(n+1)mnm} - (n+1)C_{nm}B_{k0}^{(n+1)mnm}\right]\vphantom{\sum _{s=1}^{\left[\frac{(n-m+2)}{2}\right]}}\right.\nonumber\\& \quad -\sqrt{\frac{(2n+1)(n-m)!}{(n+m)!}}\,\mathrm{Re} \left\{f_{nm}^{\prime*}(r)\left[g_{nm}(r)-rg_{nm}^{\prime}(r)\right] -rf_{nm}^{*}(r)g_{nm}^{\prime\prime}(r)\right.\nonumber\\& \qquad \left. -2n(n+1)g_{nm}^{/}(r)g_{nm}^{*}(r)\right\}\nonumber\\&\left.\times \frac{1}{\sqrt{k(k+1)}}\sum _{s=1}^{\left[\frac{(n-m+2)}{2}\right]}\sqrt{\frac{(2n-4s+3)(n+m-2s)!}{(n-m-2s+2)!}}A_{k(-1)}^{(n-2s+1)(m-1)nm}\right\}.\end{align}

Substituting (B74) into (B66), and using (B18) and (B75), one finds

(B79) \begin{equation} \overline{C}_{2k}^{\left(n,m\right)}=\frac{R_{0}^{k+1}}{k+2}\left\{\left(k-1\right)R_{0}^{k}\overline{C}_{1k}^{\left(n,m\right)}+\frac{2\sqrt{\pi}}{\sqrt{2k+1}}\left[E_{k}^{\left(n,m\right)}\left(R_{0}\right)-R_{0}E_{k}^{\left(n,m\right)/}\left(R_{0}\right)\right]\right\} .\end{equation}

Appendix C. Solving the equations of acoustic streaming for the cross terms

To solve (2.64) and (2.65), we first need to calculate their right-hand sides.

The calculation of $\boldsymbol{v}_{1}^{(n,m)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{(n,-m)}$ yields

(C1) \begin{align} &\boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}\nonumber\\&=\left(v_{1r}^{\left(n,m\right)}\frac{\partial }{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial }{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial }{\partial \phi }\right)\left(v_{1r}^{\left(n,-m\right)}\boldsymbol{e}_{r}+v_{1\theta }^{\left(n,-m\right)}\boldsymbol{e}_{\theta }+v_{1\phi }^{\left(n,-m\right)}\boldsymbol{e}_{\phi }\right)\nonumber\\&=\boldsymbol{e}_{r}\left(v_{1r}^{\left(n,m\right)}\frac{\partial v_{1r}^{\left(n,-m\right)}}{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial v_{1r}^{\left(n,-m\right)}}{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial v_{1r}^{\left(n,-m\right)}}{\partial \phi }\right.\nonumber\\&\quad\left.-\frac{v_{1\theta }^{\left(n,m\right)}v_{1\theta }^{\left(n,-m\right)}+v_{1\phi }^{\left(n,m\right)}v_{1\phi }^{\left(n,-m\right)}}{r}\right)\nonumber\\&+\boldsymbol{e}_{\theta }\left(v_{1r}^{\left(n,m\right)}\frac{\partial v_{1\theta }^{\left(n,-m\right)}}{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial v_{1\theta }^{\left(n,-m\right)}}{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial v_{1\theta }^{\left(n,-m\right)}}{\partial \phi }+\frac{v_{1\theta }^{\left(n,m\right)}v_{1r}^{\left(n,-m\right)}}{r}\right.\nonumber\\&\quad\left.-\frac{\cos \theta\, v_{1\phi }^{\left(n,m\right)}v_{1\phi }^{\left(n,-m\right)}}{r\sin \theta }\right)\nonumber\\&+\boldsymbol{e}_{\phi }\left(v_{1r}^{\left(n,m\right)}\frac{\partial v_{1\phi }^{\left(n,-m\right)}}{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial v_{1\phi }^{\left(n,-m\right)}}{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial v_{1\phi }^{\left(n,-m\right)}}{\partial \phi }+\frac{v_{1\phi }^{\left(n,m\right)}v_{1r}^{\left(n,-m\right)}}{r}\right.\nonumber\\&\quad\left.+\frac{\cos \theta\, v_{1\phi }^{\left(n,m\right)}v_{1\theta }^{\left(n,-m\right)}}{r\sin \theta }\right). \end{align}

Substituting (2.30)–(2.32) and (2.41)–(2.43) into (C1), and averaging over time, one obtains

(C2) \begin{align} & \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}\right\rangle \nonumber\\& =\displaystyle\frac{\left(-1\right)^{m}}{2}\,\mathrm{Re}\left\{\boldsymbol{e}_{r}\varepsilon \left\{V_{n}V_{n}^{\prime\ast }\left(Y_{n}^{m}\right)^{2}+\displaystyle\frac{W_{n}\left(V_{n}^{\ast }-W_{n}^{\ast }\right)}{r}\left[\left(\frac{\partial Y_{n}^{m}}{\partial \theta }\right)^{2}+\frac{1}{\sin ^{2}\theta }\left(\frac{\partial Y_{n}^{m}}{\partial \phi }\right)^{2}\right]\right\}\right.\nonumber\\&+\boldsymbol{e}_{\theta }\varepsilon \left\{\left(V_{n}W_{n}^{\prime\ast }+\frac{W_{n}V_{n}^{\ast }}{r}\right)Y_{n}^{m}\frac{\partial Y_{n}^{m}}{\partial \theta }\right.\nonumber\\&\left.+\frac{|W_{n}|^{2}}{r}\left[\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial ^{2}Y_{n}^{m}}{\partial \theta ^{2}}+\frac{1}{\sin ^{2}\theta }\frac{\partial Y_{n}^{m}}{\partial \phi }\frac{\partial ^{2}Y_{n}^{m}}{\partial \theta\, \partial \phi }-\frac{\cos \theta }{\sin ^{3}\theta }\left(\frac{\partial Y_{n}^{m}}{\partial \phi }\right)^{2}\right]\right\}\nonumber\\&\left.+\frac{\boldsymbol{e}_{\phi }\varepsilon }{\sin \theta }\left\{\left(V_{n}W_{n}^{\prime\ast }+\frac{W_{n}V_{n}^{\ast }}{r}\right)Y_{n}^{m}\frac{\partial Y_{n}^{m}}{\partial \phi }+\frac{|W_{n}|^{2}}{r}\left(\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial ^{2}Y_{n}^{m}}{\partial \theta\, \partial \phi }+\frac{1}{\sin ^{2}\theta }\frac{\partial Y_{n}^{m}}{\partial \phi }\frac{\partial ^{2}Y_{n}^{m}}{\partial \phi ^{2}}\right)\right\}\right\}, \end{align}

where $\varepsilon =s^{(n,m)}s^{(n,-m)\ast }$ . For simplicity, we drop the arguments of the functions in (C2).

The calculation of the curl of (C2) results in

(C3) \begin{align} &{\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}\right\rangle =\frac{\left(-1\right)^{m}}{2}\,\mathrm{Re}\left\{\frac{\boldsymbol{e}_{\theta }\varepsilon }{r\sin \theta }\left\{\left[2V_{n}V_{n}^{\prime\ast }-\left(rV_{n}W_{n}^{\prime\ast }+W_{n}V_{n}^{\ast }\right)^{\prime}\right]Y_{n}^{m}\frac{\partial Y_{n}^{m}}{\partial \phi }\right.\right.\nonumber\\&\left.+\left[\frac{2W_{n}\left(V_{n}^{\ast }-W_{n}^{\ast }\right)}{r}-\left(W_{n}W_{n}^{\ast }\right)^{\prime}\right]\left(\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial ^{2}Y_{n}^{m}}{\partial \theta\, \partial \phi }+\frac{1}{\sin ^{2}\theta }\frac{\partial Y_{n}^{m}}{\partial \phi }\frac{\partial ^{2}Y_{n}^{m}}{\partial \phi ^{2}}\right)\right\}\nonumber\\&+\frac{\boldsymbol{e}_{\phi }\varepsilon }{r}\left\{\left[\left(rV_{n}W_{n}^{/\ast }+W_{n}V_{n}^{\ast }\right)^{\prime}-2V_{n}V_{n}^{\prime\ast }\right]Y_{n}^{m}\frac{\partial Y_{n}^{m}}{\partial \theta }+\left[\left(W_{n}W_{n}^{\ast }\right)^{\prime}-\frac{2W_{n}\left(V_{n}^{\ast }-W_{n}^{\ast }\right)}{r}\right]\right.\nonumber\\&\left.\left.\times\left[\frac{\partial Y_{n}^{m}}{\partial \theta }\frac{\partial ^{2}Y_{n}^{m}}{\partial \theta ^{2}}+\frac{1}{\sin ^{2}\theta }\frac{\partial Y_{n}^{m}}{\partial \phi }\frac{\partial ^{2}Y_{n}^{m}}{\partial \theta\, \partial \phi }-\frac{\cos \theta }{\sin ^{3}\theta }\left(\frac{\partial Y_{n}^{m}}{\partial \phi }\right)^{2}\right]\right\}\right\}. \end{align}

Since (C3) does not contain a radial component, (2.64) leads to

(C4) \begin{equation} P_{kl}^{\left(\times \right)\prime\prime}\left(r\right)-\frac{k\left(k+1\right)P_{kl}^{\left(\times \right)}\left(r\right)}{r^{2}}=0 .\end{equation}

A solution to (C4) is given by

(C5) \begin{equation} P_{kl}^{\left(\times \right)}\left(r\right)=\overline{C}_{1k}^{\left(\times \right)}r^{k+1}+\overline{C}_{2k}^{\left(\times \right)}r^{-k} ,\end{equation}

where $\overline{C}_{1k}^{(\times )}$ and $\overline{C}_{2k}^{(\times )}$ are constants. From the condition that the acoustic streaming vanishes at infinity, it follows that

(C6) \begin{equation} \overline{C}_{1k}^{\left(\times \right)}=0 ,\end{equation}

while $\overline{C}_{2k}^{(\times )}$ is calculated by boundary conditions at the bubble surface; see below.

The calculation of the $r$ -component of ${\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \langle \boldsymbol{v}_{1}^{(n,m)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{(n,-m)}\rangle$ results in

(C7) \begin{align} &\boldsymbol{e}_{r}\cdot \left[{\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}\right\rangle \right]\nonumber\\&=\frac{\left(-1\right)^{m}}{4r^{2}}\,\mathrm{Re}\left\{\varepsilon \left[2V_{n}V_{n}^{\prime\ast }-\left(rV_{n}W_{n}^{\prime\ast }+W_{n}V_{n}^{\ast }\right)^{\prime}\right]L^{2}\left[\left(Y_{n}^{m}\right)^{2}\right]\right.\nonumber\\&\left.+\varepsilon \left[\frac{2W_{n}\left(V_{n}^{\ast }-W_{n}^{\ast }\right)}{r}-\left(W_{n}W_{n}^{\ast }\right)^{\prime}\right]L^{2}\left[\left(\frac{\partial Y_{n}^{m}}{\partial \theta }\right)^{2}+\frac{1}{\sin ^{2}\theta }\left(\frac{\partial Y_{n}^{m}}{\partial \phi }\right)^{2}\right]\right\}, \end{align}

where the operator $L^{2}$ is given by (B9).

The $r$ -component of ${\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \langle \boldsymbol{v}_{1}^{(n,-m)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{(n,m)}\rangle$ is calculated from (C7) by swapping $m$ with $-m$ . Doing so and using (D3), one obtains

(C8) \begin{align} &\boldsymbol{e}_{r}\cdot \left[{\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}+\boldsymbol{v}_{1}^{\left(n,-m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle \right]=\frac{\left(-1\right)^{m}}{4r^{2}}\nonumber\\&\times\mathrm{Re}\left\{\varepsilon \!\left[2V_{n}V_{n}^{/\ast }+2V_{n}^{\ast }V_{n}^{\prime}\!-\left(rV_{n}W_{n}^{/\ast }+W_{n}V_{n}^{\ast }\right)^{\prime}\!-\left(rV_{n}^{\ast }W_{n}^{\prime}+W_{n}^{\ast }V_{n}\right)^{\prime}\right]L^{2}\left[\left(Y_{n}^{m}\right)^{2}\right]\right.\nonumber\\&\left.+2\varepsilon \left[\frac{V_{n}W_{n}^{\ast }+V_{n}^{\ast }W_{n}-2W_{n}W_{n}^{\ast }}{r}-\left(W_{n}W_{n}^{\ast }\right)^{\prime}\right]L^{2}\left[\left(\frac{\partial Y_{n}^{m}}{\partial \theta }\right)^{2}-\frac{m^{2}\left(Y_{n}^{m}\right)^{2}}{\sin ^{2}\theta }\right]\right\}. \end{align}

By using (D5), it can be shown that the following identity holds:

(C9) \begin{equation} \left(\frac{\partial Y_{n}^{m}}{\partial \theta }\right)^{2}-\frac{m^{2}\left(Y_{n}^{m}\right)^{2}}{\sin ^{2}\theta }=n\left(n+1\right)\left(Y_{n}^{m}\right)^{2}-\frac{1}{2}L^{2}\left[\left(Y_{n}^{m}\right)^{2}\right] .\end{equation}

Substitution of (C9) into (C8) yields

(C10) \begin{align} &\boldsymbol{e}_{r}\cdot \left[{\boldsymbol \nabla }\times {\boldsymbol \nabla }\times \left\langle \boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)}+\boldsymbol{v}_{1}^{\left(n,-m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)}\right\rangle \right]=\frac{\left(-1\right)^{m}}{4r^{2}}\nonumber\\&\times\mathrm{Re}\left\{\varepsilon L^{2}\left[\left(Y_{n}^{m}\right)^{2}\right]\right.\left\{2V_{n}V_{n}^{/\ast }+2V_{n}^{\ast }V_{n}^{\prime}-\left(rV_{n}W_{n}^{/\ast }+V_{n}^{\ast }W_{n}+rV_{n}^{\ast }W_{n}^{\prime}+V_{n}W_{n}^{\ast }\right)^{\prime}\right.\nonumber\\&\left.+2n\left(n+1\right)\left[\frac{V_{n}W_{n}^{\ast }+V_{n}^{\ast }W_{n}-2W_{n}W_{n}^{\ast }}{r}-\left(W_{n}W_{n}^{\ast }\right)^{\prime}\right]\right\}\nonumber\\&\left.+\varepsilon L^{4}\left[\left(Y_{n}^{m}\right)^{2}\right]\left[\left(W_{n}W_{n}^{\ast }\right)^{\prime}-\frac{V_{n}W_{n}^{\ast }+V_{n}^{\ast }W_{n}-2W_{n}W_{n}^{\ast }}{r}\right]\right\}. \end{align}

Substituting (C10) into (2.65), expressing $(Y_{n}^{m})^{2}$ by (D25), and using the identity $L^{2}Y_{k}^{l}=k(k+1)Y_{k}^{l}$ , one obtains

(C11) \begin{equation} T_{kl}^{\left(\times \right)//}\left(r\right)-\frac{k\left(k+1\right)}{r^{2}}T_{kl}^{\left(\times \right)}\left(r\right)=D_{kl}^{nmnm}G_{k}^{\left(\times \right)}\left(r\right) ,\end{equation}

where the constant coefficients $D_{kl}^{nmnm}$ are calculated by (D27), and $G_{k}^{(\times )}(r)$ is given by

(C12) \begin{align} &G_{k}^{\left(\times \right)}\left(r\right)\nonumber\\&=\frac{\left(-1\right)^{m}\varepsilon }{2\nu }\mathrm{Re}\left\{V_{n}^{\ast }\left(r\right)\left[2V_{n}^{\prime}\left(r\right)-2W_{n}^{\prime}\left(r\right)-rW_{n}^{//}\left(r\right)\right]-V_{n}^{/\ast }\left(r\right)\left[W_{n}\left(r\right)+rW_{n}^{\prime}\left(r\right)\right]\right.\nonumber\\&\left.+\frac{2n\left(n+1\right)-k\left(k+1\right)}{r}W_{n}^{\ast }\left(r\right)\left[V_{n}\left(r\right)-W_{n}\left(r\right)-rW_{n}^{\prime}\left(r\right)\right]\right\}. \end{align}

A solution to (C11) is given by

(C13) \begin{align} T_{kl}^{\left(\times \right)}\left(r\right)=D_{kl}^{nmnm}\left[r^{k+1}C_{3k}^{\left(\times \right)}\left(r\right)+r^{-k}C_{4k}^{\left(\times \right)}\left(r\right)\right] ,\end{align}
(C14) \begin{align} C_{3k}^{\left(\times \right)}\left(r\right)=\overline{C}_{3k}^{\left(\times \right)}+\frac{1}{2k+1}\int _{R_{0}}^{r}s^{-k}G_{k}^{\left(\times \right)}\left(s\right){\rm d}s ,\quad\end{align}
(C15) \begin{align} C_{4k}^{\left(\times \right)}\left(r\right)=\overline{C}_{4k}^{\left(\times \right)}-\frac{1}{2k+1}\int _{R_{0}}^{r}s^{k+1}G_{k}^{\left(\times \right)}\left(s\right){\rm d}s. \quad\end{align}

The constant $\overline{C}_{3k}^{(\times )}$ is calculated from the condition that the acoustic streaming vanishes at infinity,

(C16) \begin{equation} \overline{C}_{3k}^{\left(\times \right)}=-\frac{1}{2k+1}\int _{R_{0}}^{\infty }s^{-k}G_{k}^{\left(\times \right)}\left(s\right){\rm d}s ,\end{equation}

and the constant $\overline{C}_{4k}^{(\times )}$ is calculated by boundary conditions at the bubble surface; see below.

From (2.63), it follows that

(C17) \begin{align} \boldsymbol{v}_{E}^{\left(\times \right)}&=\mathrm{Re}\left\{{\boldsymbol \nabla }\times \left[\boldsymbol{e}_{r}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}P_{kl}^{\left(\times \right)}\left(r\right)Y_{k}^{l}\left(\theta ,\phi \right)\right]\right.\nonumber\\&\quad\left.+\boldsymbol{e}_{r}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}T_{kl}^{\left(\times \right)}\left(r\right)Y_{k}^{l}\left(\theta ,\phi \right)+{\boldsymbol \nabla }\Phi ^{\left(\times \right)}\left(r,\theta ,\phi \right)\right\} .\end{align}

Substitution of (C17) into (2.44) yields

(C18) \begin{align} \Delta \Phi ^{\left(\times \right)}\left(r,\theta ,\phi \right)&=-\sum _{k=1}^{\infty }\sum _{l=-k}^{k}{\boldsymbol \nabla }\cdot \left[T_{kl}^{\left(\times \right)}\left(r\right)Y_{k}^{l}\left(\theta ,\phi \right)\boldsymbol{e}_{r}\right]\nonumber\\&=-\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\left[T_{kl}^{\left(\times \right)/}\left(r\right)+2r^{-1}T_{kl}^{\left(\times \right)}\left(r\right)\right]Y_{k}^{l}\left(\theta ,\phi \right) .\end{align}
(C19) \begin{equation} \Phi ^{\left(\times \right)}\left(r,\theta ,\phi \right)=\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\Phi _{kl}^{\left(\times \right)}\left(r\right)Y_{k}^{l}\left(\theta ,\phi \right) .\end{equation}

Substituting (C19) into (C18), and using (C13), one obtains

(C20) \begin{align} \Phi _{kl}^{\left(\times \right)//}\left(r\right)+\frac{2}{r}\Phi _{kl}^{\left(\times \right)/}\left(r\right)-\frac{k\left(k+1\right)}{r^{2}}\Phi _{kl}^{\left(\times \right)}\left(r\right)=D_{kl}^{nmnm}H_{k}^{\left(\times \right)}\left(r\right) , \end{align}
(C21) \begin{align} H_{k}^{\left(\times \right)}\left(r\right)=-\left(k+3\right)r^{k}C_{3k}^{\left(\times \right)}\left(r\right)+\left(k-2\right)r^{-k-1}C_{4k}^{\left(\times \right)}\left(r\right) . \quad\end{align}

Equation (C20) is solved by the method of variation of parameters, which results in

(C22) \begin{equation} \Phi _{kl}^{\left(\times \right)}\left(r\right)=D_{kl}^{nmnm}\left[r^{k}C_{5k}^{\left(\times \right)}\left(r\right)+r^{-k-1}C_{6k}^{\left(\times \right)}\left(r\right)\right] ,\end{equation}

where $C_{5k}^{(\times )}(r)$ and $C_{6k}^{(\times )}(r)$ obey the equations

(C23) \begin{align} r^{k}C_{5k}^{\left(\times \right)/}\left(r\right)+r^{-k-1}C_{6k}^{\left(\times \right)/}\left(r\right)=0 ,\qquad\qquad\end{align}
(C24) \begin{align} kr^{k-1}C_{5k}^{\left(\times \right)/}\left(r\right)-\left(k+1\right)r^{-k-2}C_{6k}^{\left(\times \right)/}\left(r\right)=H_{k}^{\left(\times \right)}\left(r\right) . \end{align}

Solutions to (C23) and (C24) are given by

(C25) \begin{align} C_{5k}^{\left(\times \right)}\left(r\right)=\overline{C}_{5k}^{\left(\times \right)}+\frac{1}{2k+1}\int _{R_{0}}^{r}s^{1-k}H_{k}^{\left(\times \right)}\left(s\right){\rm d}s ,\end{align}
(C26) \begin{align} C_{6k}^{\left(\times \right)}\left(r\right)=\overline{C}_{6k}^{\left(\times \right)}-\frac{1}{2k+1}\int _{R_{0}}^{r}s^{k+2}H_{k}^{\left(\times \right)}\left(s\right){\rm d}s ,\end{align}

where $\overline{C}_{5k}^{(\times )}$ and $\overline{C}_{6k}^{(\times )}$ are constants. The constant $\overline{C}_{6k}^{(\times )}$ is calculated by boundary conditions at the bubble surface; see below. The constant $\overline{C}_{5k}^{(\times )}$ is calculated from the condition that the acoustic streaming vanishes at infinity, which gives

(C27) \begin{equation} \overline{C}_{5k}^{\left(\times \right)}=-\frac{1}{2k+1}\int _{R_{0}}^{\infty }s^{1-k}H_{k}^{\left(\times \right)}\left(s\right){\rm d}s .\end{equation}

Substituting (C21) into (C27), and using (C15), one obtains

(C28) \begin{align} \overline{C}_{5k}^{\left(\times \right)}=-\frac{\left(k-2\right)\overline{C}_{4k}^{\left(\times \right)}}{\left(4k^{2}-1\right)R_{0}^{2k-1}}+I_{k}^{\left(\times \right)} ,\qquad\qquad\qquad\qquad\qquad\end{align}
(C29) \begin{align} I_{k}^{\left(\times \right)}=\frac{k+3}{2k+1}\int _{R_{0}}^{\infty }rC_{3k}^{\left(\times \right)}\left(r\right){\rm d}r+\frac{k-2}{\left(2k+1\right)^{2}}\int _{R_{0}}^{\infty }r^{-2k}\left[\int _{R_{0}}^{r}s^{k+1}G_{k}^{\left(\times \right)}\left(s\right){\rm d}s\right]{\rm d}r . \end{align}

Substitution of (C19) into (C17) yields

(C30) \begin{align} v_{Er}^{\left(\times \right)}=\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\left[T_{kl}^{\left(\times \right)}\left(r\right)+\Phi _{kl}^{\left(\times \right)/}\left(r\right)\right]Y_{k}^{l}\left(\theta ,\phi \right) ,\qquad\quad\end{align}
(C31) \begin{align} v_{E\theta }^{\left(\times \right)}=\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\left[\frac{\Phi _{kl}^{\left(\times \right)}\left(r\right)}{r}\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }+\frac{P_{kl}^{\left(\times \right)}\left(r\right)}{r\sin \theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \phi }\right] , \end{align}
(C32) \begin{align} v_{E\phi }^{\left(\times \right)}=\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\left[\frac{\Phi _{kl}^{\left(\times \right)}\left(r\right)}{r\sin \theta }\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \phi }-\frac{P_{kl}^{\left(\times \right)}\left(r\right)}{r}\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }\right] . \end{align}

The function $\Phi _{kl}^{(\times )/}(r)$ , which appears in (C30), is calculated by (C22) and (C23) to be

(C33) \begin{equation} \Phi _{kl}^{\left(\times \right)/}\left(r\right)=D_{kl}^{nmnm}\left[kr^{k-1}C_{5k}^{\left(\times \right)}\left(r\right)-\left(k+1\right)r^{-k-2}C_{6k}^{\left(\times \right)}\left(r\right)\right] .\end{equation}

In order to go on with the calculation, we need to apply the boundary conditions for the acoustic streaming at the bubble surface. To do this, we need to know the Stokes drift velocity, which is calculated by

(C34) \begin{equation} \boldsymbol{v}_{S}^{\left(\times \right)}=\frac{1}{2\omega }\mathrm{Re}\left\{i\boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)\ast }+i\boldsymbol{v}_{1}^{\left(n,-m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,m\right)\ast }\right\} .\end{equation}

The expression $\mathrm{Re}\{i\boldsymbol{v}_{1}^{(n,m)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{(n,-m)\ast }\}$ is calculated by

(C35) \begin{align} & \mathrm{Re}\left\{i\boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)\ast }\right\}\nonumber\\ &=\mathrm{Re}\left\{\boldsymbol{e}_{r}i\left(v_{1r}^{\left(n,m\right)}\frac{\partial v_{1r}^{\left(n,-m\right)\ast }}{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial v_{1r}^{\left(n,-m\right)\ast }}{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial v_{1r}^{\left(n,-m\right)\ast }}{\partial \phi }-\frac{v_{1\theta }^{\left(n,m\right)}v_{1\theta }^{\left(n,-m\right)\ast }}{r}\right.\right.\nonumber\\ &\quad\left.-\frac{v_{1\phi }^{\left(n,m\right)}v_{1\phi }^{\left(n,-m\right)\ast }}{r}\right)\nonumber\\ &+\boldsymbol{e}_{\theta }i\left(v_{1r}^{\left(n,m\right)}\frac{\partial v_{1\theta }^{\left(n,-m\right)\ast }}{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial v_{1\theta }^{\left(n,-m\right)\ast }}{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial v_{1\theta }^{\left(n,-m\right)\ast }}{\partial \phi }+\frac{v_{1\theta }^{\left(n,m\right)}v_{1r}^{\left(n,-m\right)\ast }}{r}\right.\nonumber\\ &\quad\left.-\frac{\cos \theta\, v_{1\phi }^{\left(n,m\right)}v_{1\phi }^{\left(n,-m\right)\ast }}{r\sin \theta }\right)\nonumber\\ &+\boldsymbol{e}_{\phi }i\left(v_{1r}^{\left(n,m\right)}\frac{\partial v_{1\phi }^{\left(n,-m\right)\ast }}{\partial r}+\frac{v_{1\theta }^{\left(n,m\right)}}{r}\frac{\partial v_{1\phi }^{\left(n,-m\right)\ast }}{\partial \theta }+\frac{v_{1\phi }^{\left(n,m\right)}}{r\sin \theta }\frac{\partial v_{1\phi }^{\left(n,-m\right)\ast }}{\partial \phi }+\frac{v_{1\phi }^{\left(n,m\right)}v_{1r}^{\left(n,-m\right)\ast }}{r}\right.\nonumber\\ &\quad\left.\left.+\frac{\cos \theta\, v_{1\phi }^{\left(n,m\right)}v_{1\theta }^{\left(n,-m\right)\ast }}{r\sin \theta }\right)\right\}. \end{align}

Substituting (2.30)–(2.32) and (2.41)–(2.43) into (C35), and using (C9), one obtains

(C36) \begin{align} & \mathrm{Re}\left\{i\boldsymbol{v}_{1}^{\left(n,m\right)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{\left(n,-m\right)\ast }\right\}=\left(-1\right)^{m}\notag\\ &\times\mathrm{Re}\!\left\{\boldsymbol{e}_{r}i\varepsilon \!\left\{\left[\!V_{n}V_{n}^{/\ast }\!+\frac{n\left(n+1\right)W_{n}\left(V_{n}^{\ast }-W_{n}^{\ast }\right)}{r}\right]\!\left(Y_{n}^{m}\right)^{2}-\frac{W_{n}\left(V_{n}^{\ast }-W_{n}^{\ast }\right)}{2r}L^{2}\left[\left(Y_{n}^{m}\right)^{2}\right]\right\}\right.\nonumber\\&+\boldsymbol{e}_{\theta }\frac{i\varepsilon }{2}\frac{\partial }{\partial \theta }\left\{\left[V_{n}W_{n}^{/\ast }+\frac{W_{n}\left(V_{n}^{\ast }+n\left(n+1\right)W_{n}^{\ast }\right)}{r}\right]\left(Y_{n}^{m}\right)^{2}-\frac{W_{n}W_{n}^{\ast }}{2r}L^{2}\left[\left(Y_{n}^{m}\right)^{2}\right]\right\}\nonumber\\&\left.+\boldsymbol{e}_{\phi }\frac{i\varepsilon }{2\sin \theta }\frac{\partial }{\partial \phi }\left\{\left[V_{n}W_{n}^{/\ast }+\frac{W_{n}\left(V_{n}^{\ast }+n\left(n+1\right)W_{n}^{\ast }\right)}{r}\right]\!\left(Y_{n}^{m}\right)^{2}-\frac{W_{n}W_{n}^{\ast }}{2r}L^{2}\left[\left(Y_{n}^{m}\right)^{2}\right]\right\}\!\right\}\!. \end{align}

The expression $\mathrm{Re}\{i\boldsymbol{v}_{1}^{(n,-m)}\cdot {\boldsymbol \nabla }\boldsymbol{v}_{1}^{(n,m)\ast }\}$ is calculated from (C36) by swapping $m$ with $-m$ . Doing so, and using (D3), one obtains

(C37) \begin{align} v_{Sr}^{\left(\times \right)}&=\frac{\left(-1\right)^{m}}{2\omega }\mathrm{Re}\left\{i\varepsilon \left\{\left[V_{n}V_{n}^{/\ast }-V_{n}^{\ast }V_{n}^{\prime}-\frac{n\left(n+1\right)\left(V_{n}W_{n}^{\ast }-V_{n}^{\ast }W_{n}\right)}{r}\right]\left(Y_{n}^{m}\right)^{2}\right.\right.\nonumber\\&\quad\left.\left.+\frac{V_{n}W_{n}^{\ast }-V_{n}^{\ast }W_{n}}{2r}L^{2}\left[\left(Y_{n}^{m}\right)^{2}\right]\right\}\right\} \end{align}
(C38) \begin{equation} v_{S\theta }^{\left(\times \right)}=\frac{\left(-1\right)^{m}}{4\omega }\mathrm{Re}\left\{i\varepsilon \left[V_{n}W_{n}^{/\ast }-V_{n}^{\ast }W_{n}^{\prime}-\frac{V_{n}W_{n}^{\ast }-V_{n}^{\ast }W_{n}}{r}\right]\frac{\partial \left(Y_{n}^{m}\right)^{2}}{\partial \theta }\right\} ,\end{equation}
(C39) \begin{equation} v_{S\phi }^{\left(\times \right)}=\frac{\left(-1\right)^{m}}{4\omega }\mathrm{Re}\left\{i\varepsilon \left[V_{n}W_{n}^{/\ast }-V_{n}^{\ast }W_{n}^{\prime}-\frac{V_{n}W_{n}^{\ast }-V_{n}^{\ast }W_{n}}{r}\right]\frac{1}{\sin \theta }\frac{\partial \left(Y_{n}^{m}\right)^{2}}{\partial \phi }\right\} .\end{equation}

With the help of (D25) and the identity $L^{2}Y_{k}^{l}=k(k+1)Y_{k}^{l}$ , (C37)–(C39) are transformed to

(C40) \begin{align} v_{Sr}^{\left(\times \right)}=\mathrm{Re}\sum _{k=1}^{\infty }S_{k}^{\left(\times \right)}\left(r\right)\sum _{l=-k}^{k}D_{kl}^{nmnm}Y_{k}^{l}\left(\theta ,\phi \right) ,\qquad\qquad\qquad\qquad\qquad\end{align}
(C41) \begin{align} v_{S\theta }^{\left(\times \right)}=\mathrm{Re}\left\{U_{nm}\left(r\right)\sum _{k=1}^{\infty }\sum _{l=-k}^{k}D_{kl}^{nmnm}\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }\right\} , \qquad\qquad\qquad\qquad\end{align}

(C42) \begin{align} v_{S\phi }^{\left(\times \right)}=\mathrm{Re}\left\{iU_{nm}\left(r\right)\sum _{k=1}^{\infty }\sum _{l=-k}^{k}lD_{kl}^{nmnm}\frac{Y_{k}^{l}\left(\theta ,\phi \right)}{\sin \theta }\right\} , \end{align}
(C43) \begin{align} S_{k}^{\left(\times \right)}\left(r\right)=\frac{\left(-1\right)^{m}\varepsilon }{\omega }\mathrm{Im}\left\{V_{n}^{\ast }\left(r\right)\left[V_{n}^{\prime}\left(r\right)-\frac{n\left(n+1\right)}{r}W_{n}\left(r\right)\right]+\frac{k\left(k+1\right)}{2r}V_{n}^{\ast }\left(r\right)W_{n}\left(r\right)\right\} , \end{align}
(C44) \begin{align} U_{nm}\left(r\right)=\frac{\left(-1\right)^{m}\varepsilon }{2\omega }\mathrm{Im}\left\{V_{n}^{\ast }\left(r\right)\left[W_{n}^{\prime}\left(r\right)-\frac{W_{n}\left(r\right)}{r}\right]\right\} .\end{align}

Equations (C41) and (C42) are convenient to use in the boundary conditions at the bubble surface. However, for the numerical calculation of $v_{S\theta }^{(\times )}$ and $v_{S\phi }^{(\times )}$ , it is convenient to transform (C41) and (C42).

With the help of (D9), (C41) is transformed to

(C45) \begin{align} v_{S\theta }^{\left(\times \right)}&=\frac{1}{2}\mathrm{Re}\left\{U_{nm}\left(r\right)\times \sum _{k=1}^{\infty }\sum _{l=-k}^{k}D_{kl}^{nmnm}\left[\sqrt{k\left(k+1\right)-l\left(l+1\right)}\,Y_{k}^{l+1}\left(\theta ,\phi \right)\,e^{-i\phi }\right.\right.\nonumber\\&\quad\left.\left.\vphantom{\sum _{k=1}^{\infty }}-\sqrt{k\left(k+1\right)-l\left(l-1\right)}\,Y_{k}^{l-1}\left(\theta ,\phi \right)\,e^{i\phi }\right]\right\}. \end{align}

With the help of (D16), (C42) is transformed to

(C46) \begin{align} v_{S\phi }^{\left(\times \right)}&=-\mathrm{Re}\left\{iU_{nm}\left(r\right)e^{i\phi }\sum _{k=1}^{\infty }\sum _{l=-k}^{k}l\sqrt{\frac{\left(2k+1\right)\left(k-l\right)!}{\left(k+l\right)!}}D_{kl}^{nmnm}\right.\nonumber\\&\quad\left.\times \sum _{s=1}^{\left[\left(k-l+2\right)/2\right]}\sqrt{\frac{\left(2k-4s+3\right)\left(k+l-2s\right)!}{\left(k-l-2s+2\right)!}}Y_{k-2s+1}^{l-1}\left(\theta ,\phi \right)\right\}. \end{align}

We can now apply the boundary conditions for the acoustic streaming at the bubble surface.

They are given by

(C47) \begin{equation} v_{Lr}^{\left(\times \right)}=0\quad\text{at }r=R_{0} ,\end{equation}
(C48) \begin{equation} \sigma _{Lr\theta }^{\left(\times \right)}=\eta \left(\frac{1}{r}\frac{\partial v_{Lr}^{\left(\times \right)}}{\partial \theta }+\frac{\partial v_{L\theta }^{\left(\times \right)}}{\partial r}-\frac{v_{L\theta }^{\left(\times \right)}}{r}\right)=0\quad\text{at }r=R_{0} ,\end{equation}
(C49) \begin{equation} \sigma _{Lr\phi }^{\left(\times \right)}=\eta \left(\frac{\partial v_{L\phi }^{\left(\times \right)}}{\partial r}+\frac{1}{r\sin \theta }\frac{\partial v_{Lr}^{\left(\times \right)}}{\partial \phi }-\frac{v_{L\phi }^{\left(\times \right)}}{r}\right)=0\quad\text{at }r=R_{0} ,\end{equation}

where $\boldsymbol{v}_{L}^{(\times )}=\boldsymbol{v}_{E}^{(\times )}+\boldsymbol{v}_{S}^{(\times )}$ is the Lagrangian streaming velocity, and $\sigma _{Lr\theta }^{(\times )}$ and $\sigma _{Lr\phi }^{(\times )}$ are the tangential components of the stress produced by $\boldsymbol{v}_{L}^{(\times )}$ . We use (C47)–(C49) in order to calculate the constants $\overline{C}_{2k}^{(\times )}$ , $\overline{C}_{4k}^{(\times )}$ and $\overline{C}_{6k}^{(\times )}$ .

Substituting (C30) and (C40) into (C47), and using (C13), (C28) and (C33), one obtains

(C50) \begin{equation} \frac{\left(k+1\right)\left(3k-1\right)}{4k^{2}-1}\overline{C}_{4k}^{\left(\times \right)}-\left(k+1\right)R_{0}^{-2}\overline{C}_{6k}^{\left(\times \right)}=-R_{0}^{2k+1}\overline{C}_{3k}^{\left(\times \right)}-kR_{0}^{2k-1}I_{k}^{\left(\times \right)}-R_{0}^{k}S_{k}^{\left(\times \right)}\left(R_{0}\right) .\end{equation}

The calculation of $\sigma _{Lr\theta }^{(\times )}$ and $\sigma _{Lr\phi }^{(\times )}$ results in

(C51) \begin{align} \sigma _{Lr\theta }^{\left(\times \right)}&=\eta\, \mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\left\{D_{kl}^{nmnm}\left[r^{k}C_{3k}^{\left(\times \right)}\left(r\right)+r^{-k-1}C_{4k}^{\left(\times \right)}\left(r\right)+2\left(k-1\right)r^{k-2}C_{5k}^{\left(\times \right)}\left(r\right)\vphantom{\frac{\partial Y_{k}^{l}}{\partial \theta }}\right.\right.\nonumber\\[5pt]&\quad\left.-2\left(k+2\right)r^{-k-3}C_{6k}^{\left(\times \right)}\left(r\right)+\frac{S_{k}^{\left(\times \right)}\left(r\right)-U_{nm}\left(r\right)+rU_{nm}^{\prime}\left(r\right)}{r}\right]\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }\nonumber\\[5pt]&\quad\left.-il\left(k+2\right)\overline{C}_{2k}^{\left(\times \right)}r^{-k-2}\frac{Y_{k}^{l}\left(\theta ,\phi \right)}{\sin \theta }\right\}, \end{align}
(C52) \begin{align} \sigma _{Lr\phi }^{\left(\times \right)}&=\eta\, \mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\left\{D_{kl}^{nmnm}\left[r^{k}C_{3k}^{\left(\times \right)}\left(r\right)+r^{-k-1}C_{4k}^{\left(\times \right)}\left(r\right)+2\left(k-1\right)r^{k-2}C_{5k}^{\left(\times \right)}\left(r\right)\vphantom{\frac{\partial Y_{k}^{l}}{\partial \theta }}\right.\right.\nonumber\\&\quad\left.-2\left(k+2\right)r^{-k-3}C_{6k}^{\left(\times \right)}\left(r\right)+\frac{S_{k}^{\left(\times \right)}\left(r\right)-U_{nm}\left(r\right)+rU_{nm}^{\prime}\left(r\right)}{r}\right]\frac{ilY_{k}^{l}}{\sin \theta }\nonumber\\&\quad\left.+\left(k+2\right)\overline{C}_{2k}^{\left(\times \right)}r^{-k-2}\frac{\partial Y_{k}^{l}}{\partial \theta }\right\}. \end{align}

Equations (C51) and (C52) show that (C48) and (C49) are satisfied if $\overline{C}_{2k}^{(\times )}=0$ . In this case, (C48) and (C49) give

(C53) \begin{align} &\frac{2k^{2}+6k-5}{4k^{2}-1}\overline{C}_{4k}^{\left(\times \right)}-2\left(k+2\right)R_{0}^{-2}\overline{C}_{6k}^{\left(\times \right)}=-R_{0}^{2k+1}\overline{C}_{3k}^{\left(\times \right)}-2\left(k-1\right)R_{0}^{2k-1}I_{k}^{\left(\times \right)}\nonumber\\[5pt]&\quad -R_{0}^{k}\left[S_{k}^{\left(\times \right)}\left(R_{0}\right)-U_{nm}\left(R_{0}\right)+R_{0}U_{nm}^{\prime}\left(R_{0}\right)\right]. \end{align}

Combining (C50) and (C53), one obtains

(C54) \begin{align} \overline{C}_{4k}^{\left(\times \right)}&= -\frac{\left(k+3\right)\left(2k-1\right)}{\left(k+1\right)\left(2k+1\right)}R_{0}^{2k+1}\overline{C}_{3k}^{\left(\times \right)}-\frac{2\left(2k-1\right)}{k+1}R_{0}^{2k-1}I_{k}^{\left(\times \right)} \nonumber\\[5pt] &\quad-\frac{2k-1}{2k+1}R_{0}^{k}\left[\frac{k+3}{k+1}S_{k}^{\left(\times \right)}\left(R_{0}\right)+U_{nm}\left(R_{0}\right)-R_{0}U_{nm}^{\prime}\left(R_{0}\right)\right] ,\end{align}
(C55) \begin{align} \overline{C}_{6k}^{\left(\times \right)}=\frac{3k-1}{4k^{2}-1}R_{0}^{2}\overline{C}_{4k}^{\left(\times \right)}+\frac{R_{0}^{k+2}}{k+1}\left[R_{0}^{k+1}\overline{C}_{3k}^{\left(\times \right)}+kR_{0}^{k-1}I_{k}^{\left(\times \right)}+S_{k}^{\left(\times \right)}\left(R_{0}\right)\right] ,\end{align}

where

(C56) \begin{align} U_{nm}^{\prime}\left(r\right)&=\frac{\left(-1\right)^{m}\varepsilon }{2\omega }\mathrm{Im}\left\{V_{n}^{/\ast }\left(r\right)\left[W_{n}^{\prime}\left(r\right)-\frac{W_{n}\left(r\right)}{r}\right]\right.\nonumber\\[5pt] &\quad\left.+V_{n}^{\ast }\left(r\right)\left[W_{n}^{//}\left(r\right)-\frac{W_{n}^{\prime}\left(r\right)}{r}+\frac{W_{n}\left(r\right)}{r^{2}}\right]\right\} .\end{align}

Since $\overline{C}_{1k}^{(\times )}=\overline{C}_{2k}^{(\times )}=0$ and hence $P_{kl}^{(\times )}(r)\equiv 0$ , $v_{E\theta }^{(\times )}$ and $v_{E\phi }^{(\times )}$ , given by (C31) and (C32), are represented by using (D9) and (D16) as follows:

(C57) \begin{align} v_{E\theta }^{\left(\times \right)}&=\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\frac{\Phi _{kl}^{\left(\times \right)}\left(r\right)}{r}\frac{\partial Y_{k}^{l}\left(\theta ,\phi \right)}{\partial \theta }=\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\frac{\Phi _{kl}^{\left(\times \right)}\left(r\right)}{2r}\nonumber\\&\times\left[\sqrt{k\left(k+1\right)-l\left(l+1\right)}Y_{k}^{l+1}\left(\theta ,\phi \right)e^{-i\phi }-\sqrt{k\left(k+1\right)-l\left(l-1\right)}Y_{k}^{l-1}\left(\theta ,\phi \right)e^{i\phi }\right], \end{align}
(C58) \begin{align} v_{E\phi }^{\left(\times \right)}&=\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\frac{il\Phi _{kl}^{\left(\times \right)}\left(r\right)}{r}\frac{Y_{k}^{l}\left(\theta ,\phi \right)}{\sin \theta }\nonumber\\&=-\mathrm{Re}\sum _{k=1}^{\infty }\sum _{l=-k}^{k}\frac{il\Phi _{kl}^{\left(\times \right)}\left(r\right)e^{i\phi }}{r}\sqrt{\frac{\left(2k+1\right)\left(k-l\right)!}{\left(k+l\right)!}}\nonumber\\&\quad\times\sum _{s=1}^{\left[\left(k-l+2\right)/2\right]}\sqrt{\frac{\left(2k-4s+3\right)\left(k+l-2s\right)!}{\left(k-l-2s+2\right)!}}Y_{k-2s+1}^{l-1}\left(\theta ,\phi \right). \end{align}

Appendix D. Mathematical identities used in calculations

In our derivation, the function $Y_{n}^{m}(\theta ,\phi )$ is defined by the following equations (Varshalovich et al. Reference Varshalovich, Moskalev and Khersonskii1988):

(D1) \begin{align} Y_{n}^{m}\left(\theta ,\phi \right)=\sqrt{\frac{\left(2n+1\right)\left(n-m\right)!}{4{\pi} \left(n+m\right)!}}\,e^{im\phi }\,P_{n}^{m}\left(\cos \theta \right),\quad 0\leq m\leq n ,\end{align}
(D2) \begin{align} P_{n}^{m}\left(\mu \right)=\left(-1\right)^{m}\left(1-\mu ^{2}\right)^{{m}/{2}}\frac{{\rm d}^{m}P_{n}\left(\mu \right)}{{\rm d}\mu ^{m}},\quad m\geq 0 ,\qquad\end{align}
(D3) \begin{align} Y_{n}^{-m}\left(\theta ,\phi \right)=\left(-1\right)^{m}Y_{n}^{m}\left(\theta ,-\phi \right)=\left(-1\right)^{m}Y_{n}^{m\ast }\left(\theta ,\phi \right) ,\qquad\end{align}
(D4) \begin{align} \int _{0}^{2{\pi} }{\rm d}\phi \int _{0}^{\pi}{\rm d}\theta \sin \theta\, Y_{n_{1}}^{m_{1}}\left(\theta ,\phi \right)\,Y_{n_{2}}^{m_{2}\ast }\left(\theta ,\phi \right)=\delta _{{n_{1}}{n_{2}}}\delta _{{m_{1}}{m_{2}}} ,\quad\end{align}

where $P_{n}(\mu )$ is the Legendre polynomial of degree $n$ , $P_{n}^{m}(\mu )$ is the associated Legendre polynomial of order $m$ and degree $n$ , $\mu =\cos \theta$ , $\delta _{nm}$ is the Kronecker delta, and the asterisk denotes the complex conjugate.

In the process of calculations, the following identities are used (Varshalovich et al. Reference Varshalovich, Moskalev and Khersonskii1988):

(D5) \begin{equation} \frac{\partial ^{2}Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta ^{2}}+\cot \theta\, \frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta }=\frac{m^{2}Y_{n}^{m}\left(\theta ,\phi \right)}{\sin ^{2}\theta }-n\left(n+1\right)Y_{n}^{m}\left(\theta ,\phi \right) ,\end{equation}
(D6) \begin{equation} \frac{{\rm d}P_{n}\left(\mu \right)}{{\rm d}\mu }=\sum _{k=1}^{\left[\left(n+1\right)/2\right]}\left(2n-4k+3\right)P_{n-2k+1}\left(\mu \right),\quad n\geq 1 ,\end{equation}
(D7) \begin{align} & \int _{0}^{2{\pi} }{\rm d}\phi \int _{0}^{\pi}{\rm d}\theta \sin \theta\, Y_{n_{1}}^{m_{1}}\left(\theta ,\phi \right)\,Y_{n_{2}}^{m_{2}}\left(\theta ,\phi \right)\,Y_{n_{3}}^{m_{3}\ast }\left(\theta ,\phi \right)\nonumber\\&\quad=\sqrt{\frac{\left(2n_{1}+1\right)\left(2n_{2}+1\right)}{4{\pi} \left(2n_{3}+1\right)}}\,C_{n_{1}0n_{2}0}^{n_{3}0}\,C_{n_{1}m_{1}n_{2}m_{2}}^{n_{3}m_{3}} ,\end{align}

(D8) \begin{equation} \sin \theta\, \frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta }=nC_{\left(n+1\right)m}Y_{n+1}^{m}\left(\theta ,\phi \right)-\left(n+1\right)C_{nm}Y_{n-1}^{m}\left(\theta ,\phi \right),\quad n\geq 1 ,\end{equation}
(D9) \begin{align} \frac{\partial Y_{n}^{m}\left(\theta ,\phi \right)}{\partial \theta }&=\frac{1}{2}\sqrt{n\left(n+1\right)-m\left(m+1\right)}\,Y_{n}^{m+1}\left(\theta ,\phi \right)\,e^{-i\phi }\nonumber\\&\quad-\frac{1}{2}\sqrt{n\left(n+1\right)-m\left(m-1\right)}\,Y_{n}^{m-1}\left(\theta ,\phi \right)\,e^{i\phi } ,\end{align}
(D10) \begin{equation} \cos \theta\, Y_{n}^{m}\left(\theta ,\phi \right)=C_{\left(n+1\right)m}Y_{n+1}^{m}\left(\theta ,\phi \right)+C_{nm}Y_{n-1}^{m}\left(\theta ,\phi \right) ,\end{equation}

where $C_{nm}$ is defined by

(D11) \begin{equation} C_{nm}=\sqrt{\frac{n^{2}-m^{2}}{\left(2n-1\right)\left(2n+1\right)}}. \end{equation}

Here, [ ] means the integer part of an expression in brackets, and $C_{n_{1}m_{1}n_{2}m_{2}}^{n_{3}m_{3}}=\langle n_{1}m_{1}n_{2}m_{2}\,|\,n_{3}m_{3}\rangle$ are the Clebsch–Gordan coefficients. The Clebsch–Gordan coefficients are zero unless the following conditions are satisfied: $m_{3}=m_{1}+m_{2} , n_{1}+n_{2}-n_{3}\geq 0 , n_{1}-n_{2}+n_{3}\geq 0 , -n_{1}+n_{2}+n_{3}\geq 0 , |m_{1}|\leq n_{1} , |m_{2}|\leq n_{2} , |m_{3}|\leq n_{3}$ .

In our derivation, we use the fact that an arbitrary function $f(\theta ,\phi )$ can be expanded in spherical harmonics by

(D12) \begin{equation} f\left(\theta ,\phi \right)=\sum _{n=0}^{\infty }\sum _{m=-n}^{n}a_{nm}Y_{n}^{m}\left(\theta ,\phi \right) ,\end{equation}

where the expansion coefficients are calculated by

(D13) \begin{equation} a_{nm}=\int _{0}^{2{\pi} }{\rm d}\phi \int _{0}^{\pi}{\rm d}\theta \sin \theta\, Y_{n}^{m\ast }\left(\theta ,\phi \right)f\left(\theta ,\phi \right) .\end{equation}

From (D13), by using (D3), one obtains

(D14) \begin{equation} a_{n\left(-m\right)}=\left(-1\right)^{m}\left[\int _{0}^{2{\pi} }{\rm d}\phi \int _{0}^{\pi}{\rm d}\theta \sin \theta\, Y_{n}^{m\ast }\left(\theta ,\phi \right)f^{\ast }\left(\theta ,\phi \right)\right]^{\ast } .\end{equation}

If $f(\theta ,\phi )$ is a real function, then it follows from (D14) that

(D15) \begin{equation} a_{n\left(-m\right)}=\left(-1\right)^{m}a_{nm}^{\ast } .\end{equation}

With the help of (D1), (D2) and (D6), one obtains

(D16) \begin{align} &\frac{Y_{n}^{m}\left(\theta ,\phi \right)}{\sin \theta }=\frac{1}{\sqrt{1-\mu ^{2}}}\sqrt{\frac{\left(2n+1\right)\left(n-m\right)!}{4{\pi} \left(n+m\right)!}}e^{im\phi }\left(-1\right)^{m}\left(1-\mu ^{2}\right)^{\frac{m}{2}}\frac{d^{m}P_{n}\left(\mu \right)}{{\rm d}\mu ^{m}}\nonumber\\[8pt]&=-e^{i\phi }\sqrt{\frac{\left(2n+1\right)\left(n-m\right)!}{4{\pi} \left(n+m\right)!}}e^{i\left(m-1\right)\phi }\left(-1\right)^{m-1}\left(1-\mu ^{2}\right)^{\frac{m-1}{2}}\frac{d^{m-1}}{{\rm d}\mu ^{m-1}}\left(\frac{{\rm d}P_{n}\left(\mu \right)}{{\rm d}\mu }\right)\nonumber\\&=-e^{i\phi }\sqrt{\frac{\left(2n+1\right)\left(n-m\right)!}{4{\pi} \left(n+m\right)!}}e^{i\left(m-1\right)\phi }\left(-1\right)^{m-1}\left(1-\mu ^{2}\right)^{\frac{m-1}{2}}\nonumber\\&\quad\times\sum _{k=1}^{\left[\left(n+1\right)/2\right]}\left(2n-4k+3\right)\frac{d^{m-1}P_{n-2k+1}\left(\mu \right)}{{\rm d}\mu ^{m-1}}\nonumber\\&=-e^{i\phi }\sqrt{\frac{\left(2n+1\right)\left(n-m\right)!}{\left(n+m\right)!}}\sum _{k=1}^{\left[\left(n-m+2\right)/2\right]}\sqrt{\frac{\left(2n-4k+3\right)\left(n+m-2k\right)!}{\left(n-m-2k+2\right)!}}\nonumber\\&\quad\times Y_{n-2k+1}^{m-1}\left(\theta ,\phi \right), n,m\geq 1. \end{align}

Let us expand $Y_{n_{1}}^{m_{1}}(\theta ,\phi )Y_{n_{2}}^{m_{2}\ast }(\theta ,\phi )$ in spherical harmonics. According to (D12) and (D13), we obtain

(D17) \begin{align} Y_{n_{1}}^{m_{1}}\left(\theta ,\phi \right)Y_{n_{2}}^{m_{2}\ast }\left(\theta ,\phi \right)=\sum _{n=0}^{\infty }\sum _{m=-n}^{n}A_{nm}^{n_{1}m_{1}n_{2}m_{2}}Y_{n}^{m}\left(\theta ,\phi \right) ,\end{align}

where

(D18) \begin{align} A_{nm}^{n_{1}m_{1}n_{2}m_{2}}=\int _{0}^{2{\pi} }{\rm d}\phi \int _{0}^{\pi}{\rm d}\theta \sin \theta\, Y_{n}^{m\ast }\left(\theta ,\phi \right)Y_{n_{1}}^{m_{1}}\left(\theta ,\phi \right)Y_{n_{2}}^{m_{2}\ast }\left(\theta ,\phi \right) . \end{align}

It follows from (D7) that

(D19) \begin{align} &\int _{0}^{2{\pi} }{\rm d}\phi \int _{0}^{\pi}{\rm d}\theta \sin \theta\, Y_{n}^{m\ast }\left(\theta ,\phi \right)Y_{n_{1}}^{m_{1}}\left(\theta ,\phi \right)Y_{n_{2}}^{m_{2}\ast }\left(\theta ,\phi \right)\nonumber\\[4pt]&\quad=\sqrt{\frac{\left(2n+1\right)\left(2n_{2}+1\right)}{4{\pi} \left(2n_{1}+1\right)}}C_{n0n_{2}0}^{n_{1}0}C_{nmn_{2}m_{2}}^{n_{1}m_{1}} .\end{align}

This equation is non-zero only if $m+m_{2}=m_{1}$ and $| n_{1}-n_{2}| \leq n\leq n_{1}+n_{2}$ . Therefore, $A_{nm}^{n_{1}m_{1}n_{2}m_{2}}$ is calculated by

(D20) \begin{equation} A_{nm}^{n_{1}m_{1}n_{2}m_{2}}=\left\{\begin{array}{ll} \sqrt{\dfrac{\left(2n+1\right)\left(2n_{2}+1\right)}{4{\pi} \left(2n_{1}+1\right)}}C_{n0n_{2}0}^{n_{1}0}C_{nmn_{2}m_{2}}^{n_{1}m_{1}}, & m=m_{1}-m_{2}\text{ and }\left| n_{1}-n_{2}\right| \\[4pt] \leq n\leq n_{1}+n_{2},\\[4pt] 0, & m\neq m_{1}-m_{2}\\[4pt] \text{ or }n\lt \left| n_{1}-n_{2}\right| \text{ or }n\gt n_{1}+n_{2}. \end{array}\right. \end{equation}

Note that $A_{nm}^{n_{1}m_{1}n_{2}m_{2}}$ is real. It also follows from the properties of the Clebsch–Gordan coefficients that $A_{nm}^{n_{1}m_{1}n_{2}m_{2}}=0$ unless $|m_{1}|\leq n_{1}$ , $|m_{2}|\leq n_{2}$ and $|m|\leq n$ . This fact should be kept in mind when implementing a numerical code.

The use of (D16) and (D17) gives the expansion

(D21) \begin{equation} \frac{Y_{n_{1}}^{m_{1}}\left(\theta ,\phi \right)Y_{n_{2}}^{m_{2}\ast }\left(\theta ,\phi \right)}{\sin ^{2}\theta }=\sum _{n=0}^{\infty }\sum _{m=-n}^{n}B_{nm}^{n_{1}m_{1}n_{2}m_{2}}Y_{n}^{m}\left(\theta ,\phi \right),\quad n_{1,2},m_{1,2}\geq 1 ,\end{equation}

where

(D22) \begin{align} &B_{nm}^{n_{1}m_{1}n_{2}m_{2}}=\sqrt{\frac{\left(2n_{1}+1\right)\left(2n_{2}+1\right)\left(n_{1}-m_{1}\right)!\left(n_{2}-m_{2}\right)!}{\left(n_{1}+m_{1}\right)!\left(n_{2}+m_{2}\right)!}}\times \sum _{k_{1}=1}^{\left[\left(n_{1}-m_{1}+2\right)/2\right]}\sum _{k_{2}=1}^{\left[\left(n_{2}-m_{2}+2\right)/2\right]}\nonumber\\&\quad\times\sqrt{\frac{\left(2n_{1}-4k_{1}+3\right)\left(2n_{2}-4k_{2}+3\right)\left(n_{1}+m_{1}-2k_{1}\right)!\left(n_{2}+m_{2}-2k_{2}\right)!}{\left(n_{1}-m_{1}-2k_{1}+2\right)!\left(n_{2}-m_{2}-2k_{2}+2\right)!}}\nonumber\\&\quad\times A_{nm}^{\left(n_{1}-2k_{1}+1\right)\left(m_{1}-1\right)\left(n_{2}-2k_{2}+1\right)\left(m_{2}-1\right)}. \end{align}

Note that $B_{nm}^{n_{1}m_{1}n_{2}m_{2}}$ is real. It should be emphasised that (D21) and (D22) are not valid if $m_{1}=0$ or $m_{2}=0$ . In our derivation, such terms vanish. However, in order to avoid problems in numerical simulations, one should specify $B_{nm}^{n_{1}m_{1}n_{2}m_{2}}=0$ if $m_{1}=0$ or $m_{2}=0$ . It also follows from (D21) that $B_{nm}^{n_{1}m_{1}n_{2}m_{2}}=0$ if $|m_{1}|\gt n_{1}$ or $|m_{2}|\gt n_{2}$ or $|m|\gt n$ .

In order to transform the angle-dependent functions that appear in (B64), we apply the following equations.

With the help of (D16) and (D17), we obtain the identity

(D23) \begin{align} \frac{Y_{n_{1}}^{m_{1}}\left(\theta ,\phi \right)Y_{n_{2}}^{m_{2}\ast }\left(\theta ,\phi \right)}{\sin \theta }&=-\sqrt{\frac{\left(2n_{1}+1\right)\left(n_{1}-m_{1}\right)!}{\left(n_{1}+m_{1}\right)!}}\sum _{n=0}^{\infty }\sum _{m=-n}^{n}Y_{n}^{m}\left(\theta ,\phi \right)e^{i\phi }\nonumber\\&\quad\times \sum _{k=1}^{\left[\left(n_{1}-m_{1}+2\right)/2\right]}\sqrt{\frac{\left(2n_{1}-4k+3\right)\left(n_{1}+m_{1}-2k\right)!}{\left(n_{1}-m_{1}-2k+2\right)!}}\nonumber\\&\quad\times A_{nm}^{\left(n_{1}-2k+1\right)\left(m_{1}-1\right)n_{2}m_{2}},\quad n_{1},m_{1}\geq 1. \end{align}

The use of (D8) and (D21) gives

(D24) \begin{align} &\frac{1}{\sin \theta }\frac{\partial Y_{n_{1}}^{m_{1}}\left(\theta ,\phi \right)}{\partial \theta }Y_{n_{2}}^{m_{2}\ast }\left(\theta ,\phi \right)\nonumber\\&=\sum _{n=0}^{\infty }\sum _{m=-n}^{n}\left[n_{1}C_{\left({n_{1}}+1\right){m_{1}}}B_{nm}^{\left(n_{1}+1\right)m_{1}n_{2}m_{2}}-\left(n_{1}+1\right)C_{{n_{1}}{m_{1}}}B_{nm}^{\left(n_{1}-1\right)m_{1}n_{2}m_{2}}\right]Y_{n}^{m}\left(\theta ,\phi \right). \end{align}

Let us expand $Y_{n_{1}}^{m_{1}}(\theta ,\phi )Y_{n_{2}}^{m_{2}}(\theta ,\phi )$ in spherical harmonics. According to (D12) and (D13), we obtain

(D25) \begin{equation} Y_{n_{1}}^{m_{1}}\left(\theta ,\phi \right)Y_{n_{2}}^{m_{2}}\left(\theta ,\phi \right)=\sum _{k=0}^{\infty }\sum _{l=-k}^{k}D_{kl}^{n_{1}m_{1}n_{2}m_{2}}Y_{k}^{l}\left(\theta ,\phi \right) ,\end{equation}
(D26) \begin{equation} D_{kl}^{n_{1}m_{1}n_{2}m_{2}}=\int _{0}^{2{\pi} }{\rm d}\phi \int _{0}^{\pi}{\rm d}\theta \sin \theta\, Y_{k}^{l\ast }\left(\theta ,\phi \right)Y_{n_{1}}^{m_{1}}\left(\theta ,\phi \right)Y_{n_{2}}^{m_{2}}\left(\theta ,\phi \right) .\end{equation}
(D27) \begin{equation} D_{kl}^{n_{1}m_{1}n_{2}m_{2}}=\sqrt{\frac{\left(2n_{1}+1\right)\left(2n_{2}+1\right)}{4{\pi} \left(2k+1\right)}}C_{n_{1}0n_{2}0}^{k0}C_{n_{1}m_{1}n_{2}m_{2}}^{kl} .\end{equation}

The coefficients $D_{kl}^{n_{1}m_{1}n_{2}m_{2}}$ are non-zero only if the following conditions are satisfied: $l=m_{1}+m_{2} , | n_{1}-n_{2}| \leq k\leq n_{1}+n_{2} , |m_{1}|\leq n_{1} , |m_{2}|\leq n_{2} , |l|\leq k$ .

References

Abramowitz, M. & Stegun, I.A. 1972 Handbook of Mathematical Functions.Dover.Google Scholar
Aghakhani, A., Yasa, O., Wrede, P. & Sitti, M. 2020 Acoustically powered surface-slipping mobile microrobots. Proc. Natl Acad. Sci. 117 (7), 3469–3477.CrossRefGoogle ScholarPubMed
Ahmed, D., Lu, M., Nourhani, A., Lammert, P.E., Stratton, Z., Muddana, H.S., Crespi, V.H. & Huang, T.J. 2015 Selectively manipulable acoustic-powered microswimmers. Sci. Rep. 5 (1), 9744.CrossRefGoogle ScholarPubMed
Ahmed, D., Ozcelik, A., Bojanala, N., Nama, N., Upadhyay, A., Chen, Y., Hanna-Rose, W. & Huang, T.J. 2016 Rotational manipulation of single cells and organisms using acoustic waves. Nat. Commun. 7 (1), 11085.CrossRefGoogle ScholarPubMed
Backus, G. 1958 A class of self-sustaining dissipative spherical dynamos. Ann. Phys. 4, 372–447.CrossRefGoogle Scholar
Backus, G. 1986 Poloidal and toroidal fields in geomagnetic field modelling. Rev. Geophys. 24 (1), 75–109.CrossRefGoogle Scholar
Bertin, N., Spelman, T.A., Stephan, O., Gredy, L., Bouriau, M., Lauga, E. & Marmottant, P. 2015 Propulsion of bubble-based acoustic microswimmers. Phys. Rev. Appl. 4 (6), 064012.CrossRefGoogle Scholar
Birkin, P.R., Offin, D.G. & Leighton, T.G. 2016 An activated fluid stream – new techniques for cold water cleaning. Ultrason. Sonochem. 29, 612–618.CrossRefGoogle ScholarPubMed
Bolanos-Jimenez, R., Rossi, M., Fernandez Rivas, D., Kähler, C.J. & Marin, A. 2017 Streaming flow by oscillating bubbles: quantitative diagnostics via particle tracking velocimetry. J. Fluid Mech. 820, 529–548.CrossRefGoogle Scholar
Boyce, W.E. & DiPrima, R.C. 2001 Elementary Differential Equations and Boundary Value Problems. Wiley.Google Scholar
Bullard, E.C. & Gellman, H. 1954 Homogeneous dynamics and terrestrial magnetism. Phil. Trans. R. Soc. Lond. A 247, 213–278.Google Scholar
Chadwick, P. & Trowbridge, E.A. 1967 Elastic wave fields generated by scalar wave functions. Proc. Camb. Phil. Soc. 63 (4), 1177–1187.CrossRefGoogle Scholar
Chandrasekhar, S. 1961 Hydrodynamics and Hydromagnetic Stability. Clarendon Press.Google Scholar
Chang, C.T., Bostwick, J.B., Daniel, S. & Steen, P.H. 2015 Dynamics of sessile drops. Part 2. Experiment. J. Fluid Mech. 768, 442–467.CrossRefGoogle Scholar
Cleve, S., Guedra, M., Mauger, C., Inserra, C. & Blanc-Benon, P. 2019 Microstreaming induced by acoustically trapped, non-spherically oscillating microbubbles. J. Fluid Mech. 875, 597–621.CrossRefGoogle Scholar
Davidson, B.J. & Riley, N. 1971 Cavitation microstreaming. J. Sound Vib. 15 (2), 217–233.CrossRefGoogle Scholar
Dijkink, R.J., van der Dennen, J.P., Ohl, C.D. & Prosperetti, A. 2006 The acoustic scallop: a bubble powered actuator. J. Micromech. Microeng. 16 (8), 1653–1659.CrossRefGoogle Scholar
Doinikov, A.A., Cleve, S., Regnault, G., Mauger, C. & Inserra, C. 2019 a Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. I. Case of modes 0 and m . Phys. Rev. E 100 (3), 033104.CrossRefGoogle Scholar
Doinikov, A.A., Cleve, S., Regnault, G., Mauger, C. & Inserra, C. 2019 b Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. II. Case of modes 1 and m . Phys. Rev. E 100 (3), 033105.CrossRefGoogle ScholarPubMed
Elder, S. 1959 Cavitation microstreaming. J. Acoust. Soc. Am. 31 (1), 54–64.CrossRefGoogle Scholar
Elsasser, W.M. 1946 Induction effects in terrestrial magnetism. Phys. Rev. 69 (3–4), 106–116.CrossRefGoogle Scholar
Fan, C.H., Liu, H.L., Ting, C.Y., Lee, Y.H., Huang, C.Y., Ma, Y.J., Wei, K.C., Yen, T.C. & Yeh, C.K. 2014 Submicron-bubble-enhanced focused ultrasound for blood–brain barrier disruption and improved CNS drug delivery. PLoS ONE 9 (5), e96327.CrossRefGoogle ScholarPubMed
Fauconnier, M., Bera, J.C. & Inserra, C. 2020 Nonspherical modes nondegeneracy of a tethered bubble. Phys. Rev. E 102 (3), 033108.CrossRefGoogle ScholarPubMed
Fauconnier, M., Mauger, C., Bera, J.C. & Inserra, C. 2022 Nonspherical dynamics and microstreaming of a wall-attached microbubble. J. Fluid Mech. 935, A22.CrossRefGoogle Scholar
Guédra, M., Inserra, C., Mauger, C. & Gilles, B. 2016 Experimental evidence of nonlinear mode coupling between spherical and nonspherical oscillations of microbubbles. Phys. Rev. E 94, 053115.Google ScholarPubMed
Inserra, C., Regnault, G., Cleve, S., Mauger, C. & Doinikov, A.A. 2020 a Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. III. Case of self-interacting modes n–n . Phys. Rev. E 101 (1), 013111.CrossRefGoogle ScholarPubMed
Inserra, C., Regnault, G., Cleve, S., Mauger, C. & Doinikov, A.A. 2020 b Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. IV. Case of modes n and m. Phys. Rev. E 102 (4), 043103.CrossRefGoogle ScholarPubMed
Kim, W., Kim, T.-H., Choi, J. & Kim, H.-Y. 2009 Mechanism of particle removal by megasonic waves. Appl. Phys. Lett. 94 (8), 081908.CrossRefGoogle Scholar
Lajoinie, G., De Cock, I., Coussios, C.C., Lentacker, I., Le Gac, S., Stride, E. & Versluis, M. 2016 In vitro methods to study bubble–cell interactions: fundamentals and therapeutic applications. Biomicrofluidics 10 (1), 011501.CrossRefGoogle ScholarPubMed
Lamb, H. 1881 On the oscillations of a viscous spheroid. Proc. Lond. Math. Soc. 13, 51–66.CrossRefGoogle Scholar
Lamb, H. 1916 Hydrodynamics. Cambridge University Press.Google Scholar
Landau, L.D. & Lifshitz, E.M. 1987 Fluid Mechanics. Pergamon Press.Google Scholar
Lee, K.H., Lee, J.H., Won, J.M., Rhee, K. & Chung, S.K. 2012 Micromanipulation using cavitational microstreaming generated by acoustically oscillating twin bubbles. Sensors Actuators A: Phys. 188, 442–449.CrossRefGoogle Scholar
Li, P., Collis, J.F., Brumley, D.R., Schneiders, L. & Sader, J.E. 2023 Structure of the streaming flow generated by a sphere in a fluid undergoing rectilinear oscillation. J. Fluid Mech. 974, A37.CrossRefGoogle Scholar
Longuet-Higgins, M.S. 1998 Viscous streaming from an oscillating spherical bubble. Proc. R. Soc. Lond. A 454 (1970), 725–742.CrossRefGoogle Scholar
Maksimov, A. 2020 Splitting of the surface modes for bubble oscillations near a boundary. Phys. Fluids 32 (10), 102104.CrossRefGoogle Scholar
Maksimov, A.O. 2007 Viscous streaming from surface waves on the wall of acoustically-driven gas bubbles. Eur. J. Mech. B/Fluids 26 (1), 28–42.CrossRefGoogle Scholar
Marin, A., Rossi, M., Rallabandi, B., Wang, C., Hilgenfeldt, S. & Kähler, C.J. 2015 Three-dimensional phenomena in microbubble acoustic streaming. Phys. Rev. Appl. 3 (4), 041001.CrossRefGoogle Scholar
Marmottant, P. & Hilgenfeldt, S. 2003 Controlled vesicle deformation and lysis by single oscillating bubbles. Nature 423 (6936), 153–156.CrossRefGoogle ScholarPubMed
Marmottant, P., Biben, T. & Hilgenfeldt, S. 2008 Deformation and rupture of lipid vesicles in the strong shear flow generated by ultrasound-driven microbubbles. Proc. R. Soc.Lond. A 464 (2095), 1781–1800.Google Scholar
Marmottant, P., Versluis, M., de Jong, N., Hilgenfeldt, S. & Lohse, D. 2006 High-speed imaging of an ultrasound-driven bubble in contact with a wall: ‘Narcissus’ effect and resolved acoustic streaming. Exp. Fluids 41 (2), 147–153.CrossRefGoogle Scholar
Mason, T.J. 1999 Sonochemistry, vol. 2. Oxford University Press. ISBN: 9780198503712.CrossRefGoogle Scholar
Mekki-Berrada, F., Thibault, P. & Marmottant, P. 2016 Acoustic pulsation of a microbubble confined between elastic walls. Phys. Fluids 28 (3), 032004.CrossRefGoogle Scholar
Mie, G. 1908 Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen. Ann. Phys. 25, 377–445.CrossRefGoogle Scholar
Mohanty, S., Zhang, J., McNeill, J.M., Kuenen, T., Linde, F.P., Rouwkema, J. & Misra, S. 2021 Acoustically-actuated bubble-powered rotational micro-propellers. Sensors Actuators: B Chem. 347, 130589.CrossRefGoogle Scholar
Nyborg, W.L. 1958 Acoustic streaming near a boundary. J. Acoust. Soc. Am. 30 (4), 329–339.CrossRefGoogle Scholar
Ohl, C.D. & Wolfrum, B. 2003 Detachment and sonoporation of adherent HeLa-cells by shock wave-induced cavitation. Biochimica et Biophysica Acta – General Sub. 1624 (1–3), 131–138.CrossRefGoogle ScholarPubMed
Padmavati, B.S. & Amaranath, T. 2002 A note on decomposition of solenoidal fields. Appl. Maths Lett. 15 (7), 803–805.CrossRefGoogle Scholar
Pereno, V. 2018 Layered acoustofluidic resonators for the simultaneous optical and acoustic characterisation of cavitation dynamics, microstreaming, and biological effects. Biomicrofluidics 12 (3), 034109.CrossRefGoogle ScholarPubMed
Prosperetti, A. 1977 Viscous effects on perturbed spherical flows. Q. Appl. Maths 34 (4), 339–352.CrossRefGoogle Scholar
Rallabandi, B., Wang, C. & Hilgenfeldt, S. 2014 Two-dimensional streaming flows driven by sessile semicylindrical microbubbles. J. Fluid Mech. 739, 57–71.CrossRefGoogle Scholar
Regnault, G., Mauger, C., Blanc-Benon, P., Doinikov, A.A. & Inserra, C. 2021 Signatures of microstreaming patterns induced by non-spherically oscillating bubbles. J. Acoust. Soc. Am. 150 (2), 1188–1197.CrossRefGoogle ScholarPubMed
Reuter, F., Lauterborn, S., Mettin, R. & Lauterborn, W. 2017 Membrane cleaning with ultrasonically driven bubbles. Ultrason. Sonochem. 37, 542–560.CrossRefGoogle ScholarPubMed
Saint-Michel, B. & Garbin, V. 2020 Acoustic bubble dynamics in a yield-stress fluid. Soft Matt. 16 (46), 10405–10418.CrossRefGoogle Scholar
Shaw, S.J. 2006 Translation and oscillation of a bubble under axisymmetric deformation. Phys. Fluids 18 (7), 2006.CrossRefGoogle Scholar
Spelman, T.A. & Lauga, E. 2017 Arbitrary axisymmetric steady streaming: flow, force and propulsion. J. Eng. Maths 105 (1), 31–65.CrossRefGoogle Scholar
Tho, P., Manasseh, R. & Ooi, A. 2007 Cavitation microstreaming patterns in single and multiple bubble systems. J. Fluid Mech. 576, 191–233.CrossRefGoogle Scholar
Varshalovich, D.A., Moskalev, A.N. & Khersonskii, V.K. 1988 Quantum Theory of Angular Momentum. World Scientific.CrossRefGoogle Scholar
Volk, A., Rossi, M., Rallabandi, B., Kähler, C.J., Hilgenfeldt, S. & Marin, A. 2020 Size-dependent particle migration and trapping in three-dimensional microbubble streaming flows. Phys. Rev. Fluids 5 (11), 114201.CrossRefGoogle Scholar
Wang, C., Rallabandi, B. & Hilgenfeldt, S. 2013 Frequency dependence and frequency control of microbubble streaming flows. Phys. Fluids 25 (2), 022002.CrossRefGoogle Scholar
Wu, J. & Nyborg, W.L. 2008 Ultrasound, cavitation bubbles and their interaction with cells. Adv. Drug Deliv. Rev. 60 (10), 1103–1116.CrossRefGoogle ScholarPubMed
Zhou, Y., Dai, L. & Jiao, N. 2022 Review of bubble applications in microrobotics: propulsion, manipulation, and assembly. Micromachines 13 (7), 1068.CrossRefGoogle ScholarPubMed
Figure 0

Figure 1.Coordinate system used in calculations.

Figure 1

Figure 2.Numerical shape deformations of a bubble experiencing (a) the zonal harmonic $Y_{5}^{0}$, (b) the tesseral harmonic $Y_{5}^{3}$, and (c) the sectoral harmonic $Y_{5}^{5}$. Every isometric view of the spherical harmonics is displayed in top (upper) and side (lower) views.

Figure 2

Figure 3.Microstreaming patterns associated with the three classes of spherical harmonics, i.e. (a) zonal, (b) tesseral and (c) sectoral, experienced by a wall-attached bubble and observed with microscope in a top view configuration. Images are reproduced from Fauconnier et al. (2022).

Figure 3

Figure 4.Presentation of the microstreaming pattern observed for a bubble (equilibrium radius 73.8 µm) exhibiting zonal oscillations on the mode $(n=3, m=0$): (a) numerical contours of the bubble interface in the top, side and three-dimensional views; (b) top, (c) side and (d) isometric views of the resulting flow.

Figure 4

Figure 5.Presentation of the microstreaming pattern observed for a bubble (equilibrium radius 73.8 µm) exhibiting tesseral oscillations on the mode $(n=3, m=1$): (a) numerical contours of the bubble interface in the top, side and three-dimensional views; (b) top, (c) side and (d) isometric views of the resulting flow.

Figure 5

Figure 6.Presentation of the microstreaming pattern observed for a bubble (equilibrium radius 73.8 µm) exhibiting tesseral oscillations on the mode $(n=3, m=2$): (a) numerical contours of the bubble interface in the top, side and three-dimensional views; (b) top, (c) side and (d) isometric views of the resulting flow.

Figure 6

Figure 7.Presentation of the microstreaming pattern observed for a bubble (equilibrium radius 73.8 µm) exhibiting sectoral oscillations on the mode $(n=3, m=3$): (a) numerical contours of the bubble interface in the top, side and three-dimensional views; (b) top, (c) side and (d) isometric views of the resulting flow.

Figure 7

Figure 8.Presentation of the microstreaming pattern observed for various sectoral oscillations, with the degree $n$ ranging from 2 to 5. Every line corresponds to a given degree $n$. (a,d,g,j) The projected side view microstreaming pattern. (b,e,h,k) The deformation of the bubble interface exhibiting the sectoral oscillation. (c,f,i,l) The top view microstreaming pattern.

Figure 8

Figure 9.Evolution of the microstreaming pattern for the sectoral oscillation $n=m=3$ for varying values of the dimensionless thickness of the viscous boundary layer. The liquid dynamic viscosities and resulting dimensionless thickness of the viscous boundary layer are (a) $\eta =10^{-6}$ Pa s, $\gamma =0.06$; (b) $\eta =4\times10^{-6}$ Pa s, $\gamma =0.12$, and (c) $\eta =16\times10^{-6}$ Pa s, $\gamma =0.24$. Only the top and side views are displayed.

Figure 9

Figure 10.Investigation of the orientation of the flow at the zenith of the bubble ($\theta =0$). (a) Three-dimensional view of the bubble interface oscillating on asymmetric modes with $n\in [2,5]$. (b) Evolution of the value of the colatitude $\theta _{max}$ at which the displacement of the bubble interface is maximum, as a function of the order $m$ for various degrees $n$. (c) Evolution of the amplitude of the radial velocity at the zenith ($\theta =0$) of the bubble, at the distance $r=3R_{0}$ as a function of the order $m$ for various degrees $n$. A positive value corresponds to the fountain-like streaming, while a negative value is typical of the anti-fountain behaviour.

Figure 10

Figure 11.Evolution of the microstreaming pattern in the case of a quasi-stationary wave, for the sectoral oscillation $n=m=3$. Three ratios of the amplitudes of the two modes $s^{(n,-m)}/s^{(n,m)}$ are investigated for (a) $s^{(n,-m)}/s^{(n,m)}=0.99$, (b) $s^{(n,-m)}/s^{(n,m)}=0.93$ and (c) $s^{(n,-m)}/s^{(n,m)}=0.66$. Only the top and isometric views are displayed.