# Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. IV. Case of modes n and m

Claude Inserra, Gabriel Regnault, Sarah Cleve, Cyril Mauger, Alexander A.

Doinikov

► To cite this version:

Claude Inserra, Gabriel Regnault, Sarah Cleve, Cyril Mauger, Alexander A. Doinikov. Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. IV. Case of modes n and m. Physical Review E, 2020, 102, pp.043103. 10.1103/PhysRevE.102.043103. hal-03234690

HAL Id: hal-03234690

https://hal.science/hal-03234690v1

Submitted on 6 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. IV.
## Case of modes *n* and *m*

Claude Inserra,¹, * Gabriel Regnault,² Sarah Cleve,² Cyril Mauger,² and Alexander A. Doinikov²

¹Univ Lyon, Université Claude Bernard Lyon 1, Centre Léon Bérard, INSERM, UMR 1032, LabTAU, F-69003, Lyon, France

²Univ Lyon, École Centrale de Lyon, INSA de Lyon, CNRS, LMFA UMR 5509, F-69134 Écully, France

(Dated: September 21, 2020)

This paper is the conclusion of work done in our previous papers [A.A. Doinikov et al., Phys. Rev. E 100, 033104, 033105 (2019)]. The overall aim of the study is to develop a theory for modeling the velocity field of acoustic microstreaming produced by nonspherical oscillations of a gas bubble. In our previous papers, general equations have been derived to describe the velocity field of acoustic microstreaming produced by modes *m* and *n* of bubble oscillations. Particular cases of mode interaction have been derived, such as the 0-*n*, 1-1, 1-*m* and *n*-*n* interaction. Here, the general case of interaction between modes *n* and *m*, *n* > *m*, is solved analytically. Solutions are expressed in terms of complex mode amplitudes, meaning that the mode amplitudes are assumed to be known and serve as input data for the calculation of the velocity field of microstreaming. No restrictions are imposed on the ratio of the bubble radius to the viscous penetration depth. The *n*-*m* interaction results in specific streaming patterns: at large distance from the bubble interface the pattern exhibits 2|*n* − *m*| lobes while 2min(*m*, *n*) lobes exist in the bubble vicinity. The spatial organization of the recirculation zones is unique for the interaction of two distinct nonspherical modes and therefore appears as a signature of the *n*-*m* interaction.

## I. INTRODUCTION

Microstreaming is a mean flow induced by an oscillating body submerged in a fluid. In the case of an oscillating bubble [1], the flow is driven by streaming inside the oscillatory boundary layer around the bubble interface, the so-called Stokes layer. Nonlinear second-order effects are responsible for the extension of the streaming patterns much further than the Stokes layer. The flow pattern is directly related to the bubble oscillation modes. When a bubble oscillates purely spherically, the spherosymmetry of the first-order acoustic velocity field prevents the generation of vorticity [2]. In order to describe bubble-induced streaming flows, it is thus necessary to consider supplementary oscillation modes, including the translational one, and the surface modes. The interaction between the radial oscillations and the translational motion of the bubble was first considered by Wu and Du [3] and Longuet-Higgins [4]. The physical interest in such an interaction lies in the understanding of bubbles oscillating in the vicinity of a wall, where reflected waves naturally induce the translational motion of the bubble center at the same frequency as the radial oscillations. It is worth noting that theoretical studies on acoustic streaming produced by a bubble on a wall [5] are restricted to the case of interaction induced by the monopole and dipole modes. The main theoretical challenge is the derivation of bubble-induced shear stress exerted on an elastic wall for medical applications such as the permeabilization of biological cells [21]. Recently, experimental observations have highlighted a wide variety

of microstreaming patterns induced specifically by shape modes [7, 8], with particular interest for targeted and localized drug delivery. The physical mechanisms underlying the nonspherical instability of the bubble interface is the parametric excitation, meaning that the predominant surface modes oscillate at half the driving frequency $\omega$ [9]. As so, if a parametrically-triggered surface mode oscillates at $\omega/2$ in combination with the monopole oscillation (at frequency $\omega$), then the only contribution to microstreaming comes from the self-interaction of the nonspherical mode, specifically discussed by Maksimov [10] and later Spelman and Lauga [11]. Moreover, a nonspherical mode would interact with the monopole oscillation only if it oscillates at the frequency $\omega$ and is therefore excited on its second parametric resonance, as experimentally described by Cleve et al. [7]. As a consequence, in order to generate any kind of mode interactions, a secondary nonspherical mode is required in addition to the parametrically-excited one. Such triggering is possible at sufficient acoustic pressures due to the nonlinear energy coupling between modes. Such coupling has been theoretically investigated by Doinikov [12] and Shaw [13] and experimentally evidenced by Guedra et al. [9]. When a mode *n* is parametrically excited, secondary translational *m* = 1 and nonspherical modes *m* $\neq$ *n* are expected. Particularly, even modes can only excite other even modes through the nonlinear coupling while odd modes can excite all other modes [13]. This property may lead to the generation of a broad spectrum of nonspherical modes. For the generation of the second-order mean flow, these modes will interact only if they oscillate at the same frequency. The case of interactions between modes *n* and *m* is thus of peculiar importance in order to capture the whole picture of single bubble-induced microstream-

<footnote>* claude.inserra@inserm.fr</footnote>
ing. At this point, it is necessary to indicate that all previous theoretical investigations on bubble-induced microstreaming assumed that the bubble radius $R_0$ is much greater than the viscous penetration depth $\delta$. This hypothesis results in an approximate solution for the fluid velocity field only up to leading terms with respect to the ratio $\delta/R_0$. It makes these theories invalid when this ratio is not small compared to unity, as it is the case with high-viscosity liquids and/or micron-sized bubbles. In our previous studies [2, 14, 15], we have developed a general theory for modeling the velocity field of acoustic microstreaming produced by axisymmetric nonspherical oscillations of an acoustically-driven gas bubble. The first study [2] was dedicated to the case that acoustic microstreaming is generated by the interaction of the radial mode (mode 0) with a mode of arbitrary order $m \ge 1$. In the second study [14], the contribution of the translational mode was investigated through the self-interacting mode 1-1, and the interaction of translation with any arbitrary nonspherical mode $m > 1$. In the third study [15], we investigated the self-interacting contribution of a nonspherical mode $n$.

In the present study, we derive an exact analytical solution for the Lagrangian streaming velocity induced by the interactions of two modes $n$ and $m \neq n$, without any restrictions on the ratio of the bubble radius to the viscous penetration depth. In Sec. II, the solution for the Lagrangian streaming velocity is derived. In Sec. III, numerical examples of the microstreaming induced by the interaction of two modes $n$ and $m$ are provided and compared to previous theories available in the literature.

## II. THEORY

We consider a gas bubble undergoing axisymmetric oscillations, which include the radial pulsation (mode 0), translation (mode 1), and shape modes of order $m > 1$. The liquid motion produced by the bubble oscillations is described in the spherical coordinates $r$ and $\theta$ whose origin is at the equilibrium center of the bubble, and the axis $z$ is the axis of symmetry. The geometry of the problem is depicted in Fig. 1.

Our derivation assumes that the amplitudes of the bubble oscillation modes are small compared to the equilibrium bubble radius. This assumption allows us to linearize the equations of liquid motion (Navier-Stokes equations) and to find their solutions, assuming that the amplitudes of the bubble oscillation modes are given quantities. These solutions give us the linear velocity field produced by the bubble in the liquid. This calculation is performed in Sec. II A and leads to expressions for the linear radial and tangential velocity components given by in Eqs. (2) and (3). To calculate the first-order liquid velocity, boundary conditions at the bubble surface are applied. Firstly, the normal component of the bubble surface velocity is set equal to the normal component of the fluid particle velocity at the interface. Secondly,

<figure><img src="image_2.png" /><figcaption>FIG. 1. Geometry of the system under study. (a) Three-dimensional representation of the bubble interface, where $z$ is the axis of axial symmetry. Two axisymmetric modes $n = 5$ and $m = 4$, with equal amplitudes, are depicted. (b) Axial symmetry allows using polar coordinates $(r, \theta)$ to parametrize the bubble interface $r_s$.</figcaption></figure>

the tangential stress is supposed to vanish on the bubble surface when considering a non-contaminated interface. In the next step, the equations of liquid motion are written with accuracy up to terms of the second order of smallness with respect to the linear solutions and averaged over time. This operation leads to Eq. (10) which describes the time-independent velocity field of acoustic microstreaming produced by the bubble oscillations. The solution of Eq. (10) is derived in Sec. II B. To derive the velocity field of acoustic streaming, the boundary conditions of vanishing normal velocity and tangential stress of the Lagrangian streaming are applied at the mean position of the interface. The use of these boundary conditions is described in the Appendix.

### A. Linear solutions

The bubble oscillation is decomposed over $N$ axisymmetric surface modes, corresponding to the basis of Legendre polynomials [16]. The oscillation frequencies of the modes may differ due to the parametric behavior of nonspherical bubble dynamics. The bubble surface is hence represented by

$$r_s = R_0 + \sum_{n=0}^{N} s_n e^{-i\omega_n t} P_n(\mu), \quad (1)$$

where $R_0$ is the bubble radius at rest, $s_n$ is the complex amplitude of the $n$th mode, $\omega_n$ is the angular frequency of the $n$th mode, $\mu = \cos\theta$, and $P_n$ is the Legendre polynomial of order $n$. It is assumed that $|s_n|/R_0 \ll 1$. The values of $s_n$ and $\omega_n$ are considered as known quantities, possibly experimentally measured, and serve as input data in the proposed model.

The linearized equations of an incompressible viscous liquid allow us to determine the radial, $v_{1r}$, and tangen-
tial, $v_{1\theta}$, components of the first-order liquid velocity [2]:

$$v_{1r} = -\frac{1}{R_0} \sum_{n=0}^{N} (n+1) e^{-i\omega_n t} \left[ a_n \left( \frac{\bar{x}_n}{x_n} \right)^{n+2} + nb_n \frac{\bar{x}_n}{x_n} h_n^{(1)}(x_n) \right] P_n(\mu), \quad (2)$$

$$v_{1\theta} = \frac{1}{R_0} \sum_{n=0}^{N} e^{-i\omega_n t} \left[ a_n \left( \frac{\bar{x}_n}{x_n} \right)^{n+2} - b_n \frac{\bar{x}_n}{x_n} \{ h_n^{(1)}(x_n) + x_n h_n^{(1)\prime}(x_n) \} \right] P_n^1(\mu), \quad (3)$$

where $x_n = k_n r$, $k_n = (1 + i)/\delta_n$, $\delta_n = \sqrt{2\nu/\omega_n}$, $\nu$ is the kinematic liquid viscosity, $\bar{x}_n = k_n R_0$, $h_n^{(1)}$ is the spherical Hankel function of the first kind of order $n$, $h_n^{(1)\prime}(x_n) = dh_n^{(1)}(x_n)/dx_n$ and $P_n^1$ is the associated Legendre polynomial of the first order and of degree $n$, as defined in [17]. The terms $a_n$ and $b_n$ are linear scattering coefficients of, respectively, the potential and vortical parts of the scattered wave from the bubble. They are deduced from the boundary conditions at the bubble interface:

$$v_{1r}|_{r=R_0} = \frac{dr_s}{dt} \quad (4)$$

$$\sigma_{r\theta} = \eta \left( \frac{1}{r} \frac{\partial v_{1r}}{\partial \theta} + \frac{\partial v_{1\theta}}{\partial r} - \frac{v_{1\theta}}{r} \right) = 0 \quad \text{at } r = R_0, \quad (5)$$

where $\eta$ is the dynamic liquid viscosity. Their expression is given by [2]

$$a_0 = iR_0\omega_0 s_0 \quad (6)$$

$$a_n = \frac{iR_0\omega_n s_n[\bar{x}_n^2 h_n^{1\prime\prime}(\bar{x}_n) - (n^2 + n - 2)h_n^{(1)}(\bar{x}_n)]}{(n+1)[\bar{x}_n^2 h_n^{1\prime\prime}(\bar{x}_n) + (n^2 + 3n + 2)h_n^{(1)}(\bar{x}_n)]}, \quad \text{for } n \ge 1, \quad (7)$$

$$b_n = \frac{2iR_0(n+2)\omega_n s_n}{(n+1)[\bar{x}_n^2 h_n^{1\prime\prime}(\bar{x}_n) + (n^2 + 3n + 2)h_n^{(1)}(\bar{x}_n)]}, \quad \text{for } n \ge 1, \quad (8)$$

and the coefficient $b_n$ is undefined for $n = 0$.

## B. Acoustic microstreaming produced by modes m and n

The derivation of the equations of acoustic streaming relies on taking the nonlinear incompressible Navier-Stokes equations up to second-order terms with respect to the first-order fluid velocity, and averaging them over time. Time averaging leads to the result that nonzero contributions to acoustic streaming can come either from pairs of different modes that oscillate at the same frequency or from the interaction of a mode with itself. We have already provided exact analytical solutions for the cases of microstreaming induced by the interaction of the

radial mode [2] or the translation mode $n = 1$ [14] with any arbitrary surface mode $n$, as well as the case of the self-interaction of a surface mode [15]. Here, the general case of interacting modes $m$ and $n$ is considered. According to the theory developed in Part I [2], for the case of interacting modes $m$ and $n$, the Eulerian streaming velocity is represented by

$$\langle \mathbf{v}_2^{nm} \rangle = \nabla \times [\langle \psi_2^{nm}(r, \theta) \rangle \mathbf{e}_\epsilon], \quad (9)$$

where $\langle \rangle$ denotes the time average, $\mathbf{e}_\epsilon$ is the unit azimuthal vector, and $\langle \psi_2^{nm}(r, \theta) \rangle$ is the amplitude of the vector potential of the streaming velocity that is calculated from Eq. (32) of Part I:

$$D^2 \langle \psi_2^{nm} \rangle = \frac{n+1}{2\nu r^2} P_n(\mu) P_m^1(\mu) \text{Re} \{ G_1(x_n) \} \\ + \frac{m+1}{2\nu r^2} P_m(\mu) P_n^1(\mu) \text{Re} \{ G_2(x_n) \} \quad (10) \\ - \frac{1}{2\nu r^2} \sqrt{1-\mu^2} [P n^1(\mu) P m^1(\mu)]' \text{Re} \{ G_3(x_n) \},$$

where $D$ is the linear operator given by

$$D = \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 \frac{\partial}{\partial r} \right) + \frac{1}{r^2 \sin \theta} \frac{\partial}{\partial \theta} \left( \sin \theta \frac{\partial}{\partial \theta} \right) - \frac{1}{r^2 \sin^2 \theta} \\ = \frac{k_n^2}{x_n^2} \frac{\partial}{\partial x_n} \left( x_n^2 \frac{\partial}{\partial x_n} \right) \\ + \frac{k_n^2}{x_n^2} \left[ (1 - \mu^2) \frac{\partial^2}{\partial \mu^2} - 2\mu \frac{\partial}{\partial \mu} - \frac{1}{1 - \mu^2} \right] \quad (11)$$

and the functions $G_{1,2,3}(x_n)$ are defined by

$$G_1(x_n) = k_n^2 a_n b_m^* \left(\frac{R_0}{r}\right)^{n+1} \left[(n+1)h_m^{(1)}(x_n) - x_n h_m^{(1)\prime}(x_n)\right]^* \\ - nk_n^2 b_n b_m^* \left[x_n h_n^{(1)\prime}(x_n) h_m^{(1)*}(x_n) + x_n^* h_n^{(1)}(x_n) h_m^{(1)\prime\ast}(x_n)\right], \quad (12)$$

$$G_2(x_n) = k_n^2 a_m b_n^* \left(\frac{R_0}{r}\right)^{m+1} \left[(m+1)h_m^{(1)}(x_n) - x_n h_m^{(1)\prime}(x_n)\right]^* \\ - m k_n^2 b_m b_n^* \left[x_n h_m^{(1)\prime}(x_n) h_m^{(1)*}(x_n) + x_n^* h_m^{(1)}(x_n) h_m^{(1)\prime\ast}(x_n)\right], \quad (13)$$

$$G_3(x_n) = b_m^* k_n^2 h_m^{(1)*}(x_n) \left[ a_n \left( \frac{R_0}{r} \right)^{n+1} \right. \\ \left. - b_n \left[ h_n^{(1)}(x_n) + x_n h_n^{(1)\prime}(x_n) \right] \right] \\ + b_n^* k_n^2 h_n^{(1)*}(x_n) \left[ a_m \left( \frac{R_0}{r} \right)^{m+1} \right. \\ \left. - b_m \left[ h_m^{(1)}(x_n) + x_n h_m^{(1)\prime}(x_n) \right] \right], \quad (14)$$

where the asterisk denotes the complex conjugate. In order to find an exact solution for the equation ruling
the vector potential $\psi_2^{nm}$, one may simplify the angular dependency in Eq. (10), which is decomposed over three terms involving Legendre polynomials and associated Legendre polynomials. Such a decomposition is performed along the orthogonal basis of the $P_k^i(\mu)$ functions by using the properties of overlapping integrals of three associated Legendre polynomials [18]:

$$I(l_1, m_1, l_2, m_2, l_3, m_3) = \int_{-1}^{1} P_{l_1}^{m_1}(x) P_{l_2}^{m_2}(x) P_{l_3}^{m_3}(x) dx, \quad (15)$$

whose exact analytical expression exists in the special case $m_3 = m_1 + m_2$:

$$I = 2(-1)^{m_3} \sqrt{\frac{(l_1 + m_1)!(l_2 + m_2)!(l_3 + m_3)!}{(l_1 - m_1)!(l_2 - m_2)!(l_3 - m_3)!}} \times \begin{pmatrix} l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l_3 \\ m_1 & m_2 & -m_3 \end{pmatrix}, \quad (16)$$

where we have introduced the Wigner 3-j symbols. The decomposition of the first two angular terms in Eq. (10) is straightforward and leads to

$$P_n(\mu) P_m^1(\mu) = \sum_{k=n-m}^{n+m} \alpha_{knm} P_k^1(\mu), \quad (17)$$

$$P_m(\mu) P_n^1(\mu) = \sum_{k=n-m}^{n+m} \alpha_{kmn} P_k^1(\mu), \quad (18)$$

assuming $n > m$ in the investigated case and by introducing the coefficient

$$\alpha_{knm} = -\frac{(2k+1)}{k(k+1)} \sqrt{(m+1)m(k+1)k} \times \begin{pmatrix} n & m & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} n & m & k \\ 0 & 1 & -1 \end{pmatrix}. \quad (19)$$

It is worth noting that the limits of summation $n - m \le k \le n + m$ in the Eqs. (17,18) come from the existence property of the Wigner 3-j symbols. Based on the same mathematical framework and by using common identities of associated Legendre functions [17], the third angular term appearing in Eq. (10) is decomposed as

$$\sqrt{1 - \mu^2} [P_n^1(\mu) P_m^1(\mu)]' = \sum_{k=n-m}^{n+m} \gamma_{knm} P_k^1(\mu), \quad (20)$$

where the coefficient $\gamma_{knm}$ is defined by

$$\gamma_{knm} = \frac{n(n+1)}{2}\alpha_{knm} + \frac{m(m+1)}{2}\alpha_{kmn} - \frac{1}{2}\beta_{knm} - \frac{1}{2}\beta_{kmn}, \quad (21)$$

$$\beta_{knm} = \frac{(2k+1)}{k(k+1)} \sqrt{\frac{(n+1)n(k+1)k}{(m-2)!} \begin{pmatrix} n & k & m \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} n & k & m \\ 1 & 1 & -2 \end{pmatrix}}. \quad (22)$$

In the projection of the third angular term on the $P_k^i(\mu)$ basis, the condition $n > m \ge 2$ is imposed by the existence property of Wigner 3-j symbols. By noticing that the particular cases ($m = 0$, any $n$), ($m = 1$, any $n$) and ($m, n = m$) have been previously derived [2, 14, 15], then all interacting cases will be considered. Substituting Eqs. (17,18,20) into Eq. (10), one obtains

$$D^2 \langle \psi_2^{nm} \rangle = \text{Re} \left\{ \sum_{k=n-m}^{n+m} H_k(x_n) P_k^1(\mu) \right\}, \quad (23)$$

$$H_k(x_n) = \frac{n+1}{2\nu r^2} \alpha_{knm} G_1(x_n) + \frac{m+1}{2\nu r^2} \alpha_{kmn} G_2(x_n) - \frac{1}{2\nu r^2} \gamma_{knm} G_3(x_n). \quad (24)$$

Solving Eq. (23) requires cumbersome calculations that are provided in the Appendix. As a result, the solution of Eq. (23) is found to be

$$\langle \psi_2^{nm}(x, \mu) \rangle = \text{Re} \left\{ \sum_{k=n-m}^{n+m} F_k(x_n) P_k^1(\mu) \right\}, \quad (25)$$

where the functions $F_k(x_n)$ are defined by Eq. (A6). Equation (25) leads to the following expressions for the components of the Eulerian streaming velocity:

$$\langle v_{2r}^{nm} \rangle = -\frac{1}{r} \text{Re} \left\{ \sum_{k=n-m}^{n+m} k(k+1) F_k(x_n) P_k(\mu) \right\}, \quad (26)$$

$$\langle v_{2\theta}^{nm} \rangle = -\frac{1}{r} \text{Re} \left\{ \sum_{k=n-m}^{n+m} [F_k(x_n) + x_n F'_k(x_n)] P_k^1(\mu) \right\}, \quad (27)$$

where the functions $F'_k(x_n)$ are defined by Eq. (A15).

In the process of calculating Eqs. (26) and (27), we have also obtained the following expressions for the components of the Stokes drift velocity [4] (see the Appendix):

$$v_{Sr}^{nm} = \text{Re} \left\{ \sum_{k=n-m}^{n+m} T_k(x_n) P_k(\mu) \right\}, \quad (28)$$

$$v_{S\theta}^{nm} = \text{Re} \left\{ \sum_{k=n-m}^{n+m} U_k(x_n) P_k^1(\mu) \right\}, \quad (29)$$

where the functions $T_k(x_n)$ and $U_k(x_n)$ are defined by Eqs. (A36) and (A40), respectively. The sum of the Eulerian streaming velocity and the Stokes drift velocity provides the Lagrangian streaming velocity,

$$\mathbf{v}_{L}^{nm} = \langle \mathbf{v}_{2}^{nm} \rangle + \mathbf{v}_{S}^{nm}, \quad (30)$$

for which the radial and tangential components can be calculated using Eqs. (26)-(29). It has been verified that the Lagrangian streaming velocity given by Eq. (30) corresponds to the one derived in our previous work [15] for the particular case $n = m, m \ge 2$. A MATLAB code for the calculation of the Lagrangian velocity is provided as Supplementary Material [19].
<figure><img src="image_3.png" /><figcaption>FIG. 2. Numerical examples of streamline patterns produced by the interaction of modes $n$ and $m$ in cases: (a) ($n = 3, m = 2$), (b) ($n = 4, m = 2$), (c) ($n = 4, m = 3$), (d) ($n = 5, m = 2$), (e) ($n = 5, m = 3$), (f) ($n = 5, m = 4$). Inserts show the streaming close to the bubble interface in the near-bubble region limited by the dashed line.</figcaption></figure>

### III. NUMERICAL EXAMPLES

#### A. Streamlines given by the present model

Numerical simulations were performed at the following values of physical parameters: liquid density $\rho = 1000 \text{ kg/m}^3$, dynamic liquid viscosity $\eta = 0.001 \text{ Pa} \cdot \text{s}$, driving frequency $f = 50 \text{ kHz}$, and equilibrium radius $R_0 = 50 \text{ µm}$. Figure 2 exemplifies Lagrangian streamline patterns produced by the interaction of two modes $n \in \{3, 4, 5\}$ and $m \in \{2, 3, 4\}$ with the conditions $n \neq m, n > m \geq 2$. As one can see, the main vortices form a two-scale pattern. In the vicinity of the bubble surface, the vortices have a form of lobes. According to the presented examples, the number of lobes is equal to $2\min(m, n)$ for the $n-m$ interaction. For instance, the patterns resulting from the $n-2$ interactions always reveal a 4-lobes picture in the bubble vicinity, clearly visible in Fig. 2(a,b) and noticeable in Fig. 2(d). In the latter figure, each lobe is decomposed into two counter-rotating vortices. Far from the bubble, the vortices have a lobe shape in which the number of lobes equals $2|n-m|$. In the present analysis, a specific pattern is revealed and appears as a signature of the case of the interaction between modes $n$ and $m$. We recall that, when considering

the interaction 0-$n$ between the radial mode and a nonspherical mode [2], the streamline patterns looked like lobes whose number equals $2n$. When considering the interaction between the translational mode and any nonspherical mode $n$ [14], streamlines form lobes whose number is equal to $2(n-1)$. When considering self-interacting case $n-n$, streamlines exhibit $4n$ lobes in the bubble vicinity with a cross-like shape far from the bubble [15]. It is worth noting that the presence of the vortices close to the bubble interface may be difficult to be observed experimentally. Thus the observable streaming pattern will look as that far from the bubble interface, consisting of $2(n-m)$ lobes. These patterns will look similar to the ones obtained through the interaction between the radial mode and any nonspherical mode $n$, whose lobe number equals $2n$. For instance, the pattern induced by the interaction 5-2 will be similar to the one induced by the interaction 0-3. This highlights the importance of assessing experimentally the bubble interface dynamics when discussing the induced streaming flows. Fig. 3(a,b) shows the radial and tangential components of the Lagrangian streaming velocity as a function of the distance from the bubble surface $r/R_0$ for various values of the phase shift $\Delta\phi$ between modes $n=3$ and $m=2$, calculated at the angle $\theta = \pi/4$. The velocity components are
<figure><img src="image_4.png" /><figcaption>FIG. 3. Dependency of the Lagrangian streaming velocity components on the distance from the bubble interface, for the 3 - 2 interaction: (a, b) at various values of the phase shift $\Delta\phi$ between modes, and (c, d) at various dynamic viscosity $\eta$ for the surrounding medium. The velocity components are normalized by the factor $\omega_n|s_n||s_m|/R_0$.</figcaption></figure>

normalized by the factor $\omega_n|s_n||s_m|/R_0$. For both velocity components, the amplitude of the streaming velocity increases considerably as the phase shift varies from 0 to $\pi/2$. Another observation is that the velocity amplitude decays within a short distance from the bubble surface. As the present model captures the microstreaming flow whatever the bubble size or the viscous penetration depth, Fig. 3(c,d) represents the dependence of the radial and tangential components of the Lagrangian streaming velocity for increasing values of the liquid viscosity and hence increasing ratio of the viscous penetration depth to the bubble radius. Numerical simulations are performed for the in-phase ($\Delta\phi = 0$) 3 - 2 interaction, at the angular position $\theta = \pi/4$. Results evidence that two orders of magnitude of values of liquid viscosity are required in order to significantly modify the Lagrangian streaming velocity. The full knowledge of the streaming pattern induced by a nonspherically-oscillating microbubble is of particular interest for applications such as acoustic cleaning [20], biological cell permeabilization leading to the sonoporation process [21], or any bubble-induced drug delivery technique. For instance, in therapeutic applications, it has been evidenced that bubble-induced mi-

crostreaming is responsible for cell deformability [22] or acts as a transport mechanism for the release of lipids or nanoparticles from microbubbles [23]. Characterizing the velocity field of the fluid flow, particularly in the vicinity of the bubble interface, matters when close bubble-cell interactions occur, such as in the case of capillary-confined microbubbles or in sonoporation experiments involving a large number of biological cells and a dense bubble cloud. Here partitioning of the streaming pattern is observed between confined lobes in the near field of the bubble interface that possess higher velocities relatively to the far-field lobes. The near-field high-velocity lobes act as recirculation vortices which could enhance the deformation of biological particles, as well as exerting shear stress on them, possibly until their rupture. Potential applications of near-interface high-speed vortices concern ultrasound-induced thrombolysis [24], where one must guarantee both destruction of blood clots and limited size of the clot fragments in order to avoid secondary vessel occlusion further downstream the area of the treated vessel. By trapping the fragments into near-bubble re-circulation vortices, continuous size reduction of these fragments may occur if they are subjected to sufficiently
<figure><img src="image_5.png" /><figcaption>FIG. 4. Comparison of the present model to the theory of Spelman and Lauga [11] for the interaction of modes $n = 3$ and $m = 2$. Upper line: Streamline pattern according to (a) the theory of Spelman and Lauga and (b) the present model. Evolution of the normalized (c) radial and (d) tangential components of the Lagrangian streaming velocity given by both theories and calculated for the angle $\theta = \pi/4$, as a function of the normalized distance $r/R_0$. The velocity components are normalized by the factor $\omega_n|s_n||s_m|/R_0$.</figcaption></figure>

intense shear stresses. In addition, a strong mixing efficiency of stably-oscillating bubbles has been shown to hasten enzymatic fibrinolysis, increasing the delivery of plasminogens into blood clots [25]. We therefore expect that the present model will help in designing customized microbubbles or ultrasound-emission sequences in order to enhance specific streaming flows amongst the large variety of obtained patterns.

## B. Comparison with previous theories

The present model is compared to the theory of Spelman and Lauga [11], which considers axisymmetric modes of arbitrary orders, oscillating at the same frequency. To the best of our knowledge, this model was the only one describing the microstreaming induced by arbitrary nonspherical modes. We emphasize the fact that, contrary to the Spelman model [11], we do not assume the viscous penetration length scale being small in comparison to the bubble radius. This means that, while their solutions are derived approximatively, in powers of

the small parameter $\epsilon = \delta/R_0$, our solution is available whatever the liquid viscosity and the bubble size. The Spelman model provides the following expression for the external leading order Lagrangian streaming:

$$ \langle \psi_L(r, \mu) \rangle = \sum_{k=1}^{\infty} \left( T_k r^{-k} + S_k r^{-(k-2)} - \sum_{n=0}^{\infty} \sum_{m=1}^{\infty} Y_{knm} \right. \\ \left. \times r^{-(n+m+3)} \right) \times \left( \int_{\mu}^{1} P_k(x) dx \right), \quad (31) $$

where the coefficients $T_k$, $S_k$, and $Y_{knm}$ result from the matching at the viscous boundary layer between the inner and outer solution of the Eulerian streaming. In particular, their calculation depends on the investigated mode pair contribution and the phase shift between modes. It is shown that, if the modal amplitudes ($s_n, s_m$) are in-phase or $\pi$ out-of-phase with each other, then the contribution of these mode pairs to the steady streaming is identical to zero at the first order of the expansion over $\delta/R_0$. It is thus necessary to derive expressions to upper order in the expansion, up to the third order of the ratio $\delta/R_0$ in their analysis. For the sake of simplicity, we here
compare the streamline patterns obtained by the present model and the Spelman model in the case of in-phase interacting modes $n$ and $m$. We emphasize that no examples of $n-m$ interacting cases are presented in Ref. [11]. Fig. 4 compares the streamline patterns and Lagrangian streaming velocities given by the theory of Spelman and Lauga and the present model, in case of the 3 - 2 interaction. The velocity components are normalized by the factor $\omega_n|s_n||s_m|/R_0$. The parameters are as in Fig. 2. The streamline patterns show a good qualitative agreement, except the lack of vortices in the Spelman model in the vicinity of the bubble interface. This trend has been noticed for other $n-m$ interaction cases. Concerning the radial and tangential components of the Lagrangian velocity (Fig. 4(c, d)), both models show a quantitative agreement far from the bubble interface.

## IV. CONCLUSIONS

In the present paper, the general theory developed in the previous study [2] has been applied to the case that acoustic microstreaming results from the interaction between two axisymmetric modes $n$ and $m$ such that $n > m \ge 2$. Analytical solutions are derived in terms of complex amplitudes of oscillation modes, which means that the mode amplitudes are assumed to be known and serve as input data when the velocity field of acoustic microstreaming is calculated. No restrictions are imposed

on the ratio of the bubble radius to the viscous penetration depth. The $n-m$ interaction results in a specific streamline pattern. This pattern exhibits a $2|n-m|$ lobes-like shape far from the bubble in addition to small recirculation zones in the vicinity of the bubble interface. The number of lobes in the bubble vicinity equals $2\min(m, n)$. In summary, this series of theoretical studies has evidenced the signature of each case of interacting modes. The interaction of the radial oscillation and a nonspherical mode $n$ results in a flower-type pattern with $2n$ lobes. The $1-n$ interaction generates a $2(n-1)$ lobes pattern, while the self-interacting case $n-n$ produces a cross-like shape with $4n$ lobes in the bubble vicinity. Finally, the $n-m$ ($n > m \ge 2$) interaction produces a $2|n-m|$ lobes-like shape with $2\min(m, n)$ vortices in the bubble vicinity. These features highlight the importance of determining the interface dynamics and the modal content when analyzing microstreaming flows.

## ACKNOWLEDGMENTS

This work was supported by the LabEx CeLyA of the University of Lyon (ANR-10-LABX-0060 / ANR-11-IDEX-0007). A.A.D. gratefully acknowledges the financial support from Institut National des Sciences Appliquées de Lyon (INSA de Lyon) and the Collegium de Lyon.

## Appendix A: Solution of Eq. (23)

The right-hand side of Eq. (23) suggests that a solution can be sought in the following form:

$$ \langle \psi_2^{nm}(x, \mu) \rangle = \mathrm{Re} \left\{ \sum_{k=n-m}^{n+m} F_k(x) P_k^1(\mu) \right\}, \quad (A1) $$

where $x = x_n = k_n r$ is introduced for the sake of simplicity, and $F_k(x)$ is a function to be found. The action of the operator $D^2$ on the proposed solution (A1) results in

$$ D^2 \langle \psi_2^{nm}(x, \mu) \rangle = \mathrm{Re} \left\{ \sum_{k=n-m}^{n+m} \frac{k_n^4}{x^4} \left\{ \begin{aligned} & x^4 F_k^{(IV)}(x) + 4x^3 F_k'''(x) - 2k(k+1)x^2 F_k''(x) \\ & + k(k+1)[k(k+1)-2]F_k(x) \end{aligned} \right\} P_k^1(\mu) \right\}. \quad (A2) $$

After substituting Eq. (A2) into Eq. (23), the identification for each term over the functions $P_k^1$ that forms a set of linearly independent functions provides the following equation for the $F_k(x)$ functions:

$$ \frac{d^4 F_k}{dx^4} + \frac{4}{x} \frac{d^3 F_k}{dx^3} - \frac{2k(k+1)}{x^2} \frac{d^2 F_k}{dx^2} + \frac{k(k+1)[k(k+1)-2]}{x^4} F_k(x) = \frac{1}{k_n^4} H_k(x), \quad (A3) $$

where the function $H_k(x)$ is defined by Eq. (24). Equation (A3) can be solved by the method of variation of parameters, also known as the Lagrange method [26]. According to this method, we first need to find solutions to a homogeneous equation that corresponds to the left-hand side of Eq. (A3). The solutions are sought in the form $x^\lambda$. Substitution of $x^\lambda$ into the homogeneous form of Eq. (A3) leads to a polynomial of fourth order in $\lambda$,

$$ \lambda(\lambda - 1)(\lambda - 2)(\lambda + 1) - 2k(k+1)\lambda(\lambda - 1) + k(k+1)(k(k+1) - 2) = 0. \quad (A4) $$
The roots of this polynomial are $-(k-1), -(k+1), k, k+2$. Therefore, the general solution of the homogeneous form of Eq. (A3) is written as

$$F_k(x) = C_{k1}x^{k+2} + C_{k2}x^k + \frac{C_{k3}}{x^{k-1}} + \frac{C_{k4}}{x^{k+1}}, \qquad (A5)$$

where $C_{ki}$ are constants. According to the Lagrange method, the solution of the inhomogeneous equation (A3) is obtained by setting the coefficients $C_{ki}$ to be functions of $x$, such as

$$F_k(x) = C_{k1}(x)x^{k+2} + C_{k2}(x)x^k + \frac{C_{k3}(x)}{x^{k-1}} + \frac{C_{k4}(x)}{x^{k+1}}, \qquad (A6)$$

where the $C_{ki}$ functions should obey the following system of equations:

$$\begin{aligned} C'_{k1}y_1 + C'_{k2}y_2 + C'_{k3}y_3 + C'_{k4}y_4 &= 0 \\ C'_{k1}y'_1 + C'_{k2}y'_2 + C'_{k3}y'_3 + C'_{k4}y'_4 &= 0 \\ C'_{k1}y''_1 + C'_{k2}y''_2 + C'_{k3}y''_3 + C'_{k4}y''_4 &= 0 \\ C'_{k1}y'''_1 + C'_{k2}y'''_2 + C'_{k3}y'''_3 + C'_{k4}y'''_4 &= H_k(x)/k_n^4 \end{aligned} \qquad (A7)$$

Here the prime denotes the derivative with respect to $x$ and the functions $y_i$ are given by

$$y_1 = x^{k+2}, \quad y_2 = x^k, \quad y_3 = x^{-(k-1)}, \quad y_4 = x^{-(k+1)}. \qquad (A8)$$

Equations (A7) are a system of algebraic equations with the unknowns $C'_{ki}$. Solving this system and integrating the solutions over $x$, one obtains

$$C_{k1}(x) = C_{k10} + \frac{1}{2(2k+1)(2k+3)} \int_{x_{n0}}^{x} \frac{s^{1-k}}{k_n^4} H_k(s) ds, \qquad (A9)$$

$$C_{k2}(x) = C_{k20} - \frac{1}{2(2k-1)(2k+1)} \int_{x_{n0}}^{x} \frac{s^{3-k}}{k_n^4} H_k(s) ds, \qquad (A10)$$

$$C_{k3}(x) = C_{k30} + \frac{1}{2(2k-1)(2k+1)} \int_{x_{n0}}^{x} \frac{s^{2+k}}{k_n^4} H_k(s) ds, \qquad (A11)$$

$$C_{k4}(x) = C_{k40} - \frac{1}{2(2k+1)(2k+3)} \int_{x_{n0}}^{x} \frac{s^{4+k}}{k_n^4} H_k(s) ds, \qquad (A12)$$

where the constants $C_{ki0}$, $i = 1, 2, 3, 4$, have to be determined according to the boundary conditions. To apply the boundary conditions, we first calculate the components of the Eulerian streaming velocity by using Eq. (A1):

$$\begin{aligned} \langle v_{2r}^{nm} \rangle &= -\frac{1}{r} \frac{\partial}{\partial \mu} \{ \langle \psi_2^{nm} \rangle \sqrt{1 - \mu^2} \} \qquad (A13) \\ &= -\frac{1}{r} \mathrm{Re} \left\{ \sum_{k=n-m}^{n+m} k(k+1) F_k(x) P_k(\mu) \right\}, \end{aligned}$$

$$\begin{aligned} \langle v_{2\theta}^{nm} \rangle &= -\frac{1}{r} \frac{\partial}{\partial x} \{x \langle \psi_2^{nm} \rangle\} \qquad (A14) \\ &= -\frac{1}{r} \mathrm{Re} \left\{ \sum_{k=n-m}^{n+m} [F_k(x) + x F'_k(x)] P_k^1(\mu) \right\}, \end{aligned}$$

where the first derivative of the $F_k$ function is written as

$$F'_{k}(x) = (k + 2)C_{k1}(x)x^{k+1} + kC_{k2}(x)x^{k-1} - (k - 1)C_{k3}(x)x^{-k} - (k + 1)C_{k4}(x)x^{-k-2}. \qquad (A15)$$
The condition of zero streaming velocity at infinity requires that $F_k(x)/r \to 0$ for $r \to \infty$, which leads to

$$C_{k10} = -\frac{1}{2(2k+1)(2k+3)} \int_{x_{n0}}^{\infty} \frac{s^{1-k}}{k_n^4} H_k(s) ds, \qquad (A16)$$

$$C_{k20} = \frac{1}{2(2k-1)(2k+1)} \int_{x_{n0}}^{\infty} \frac{s^{3-k}}{k_n^4} H_k(s) ds. \qquad (A17)$$

In order to calculate the coefficients $C_{k30}$ and $C_{k40}$, boundary conditions at the bubble surface have to be applied. Equations (A13) and (A14) give the components of the Eulerian streaming velocity. To apply the boundary conditions at the bubble surface, we need to know the Lagrangian streaming velocity, which is defined by

$$v_L = \langle v_2 \rangle + v_S, \qquad (A18)$$

where $v_S$ is the Stokes drift velocity given by [4]

$$v_S = \left\langle \left( \int v_1 dt \cdot \nabla \right) v_1 \right\rangle_{nm}, \qquad (A19)$$

where $v_1$ is the first-order liquid velocity. For the case of the interaction between modes $n$ and $m$, Eq. (A19) reduces to

$$v_S^{nm} = \frac{1}{2\omega_n} \mathrm{Re} \{ \iota (v_{1n} + v_{1m} \cdot \nabla) (v_{1n} + v_{1m})^* \}, \qquad (A20)$$

where $v_{1i}$ is the first-order liquid velocity produced by the mode $i$. As a result, the components of the Stokes drift velocity associated to the $m-n$ interaction are given by

$$v_{Sr}^{nm} = \frac{1}{2\omega_n} \mathrm{Re} \left\{ \iota (v_{nr} + v_{mr}) \frac{\partial(v_{nr} + v_{mr})^*}{\partial r} + \frac{\iota}{r} (v_{n\theta} + v_{m\theta}) \frac{\partial(v_{nr} + v_{mr})^*}{\partial \theta} \right\}, \qquad (A21)$$

$$v_{S\theta}^{nm} = \frac{1}{2\omega_n} \mathrm{Re} \left\{ \iota (v_{nr} + v_{mr}) \frac{\partial(v_{n\theta} + v_{m\theta})^*}{\partial r} + \frac{\iota}{r} (v_{n\theta} + v_{m\theta}) \left[ v_{nr} + v_{mr} + \frac{\partial(v_{n\theta} + v_{m\theta})}{\partial \theta} \right]^* \right\}. \qquad (A22)$$

The components of the first-order liquid velocity induced by the oscillation of any arbitrary axisymmetric mode $n$ have already been determined in our previous study [2]:

$$v_{1nr} = -\frac{1}{R_0}(n+1)e^{-i\omega_n t} \left[ a_n \left(\frac{R_0}{r}\right)^{n+2} + n b_n \frac{R_0}{r} h_n^{(1)}(x) \right] P_n(\mu), \qquad (A23)$$

$$v_{1n\theta} = \frac{1}{R_0} e^{-i\omega_n t} \left[ a_n \left(\frac{R_0}{r}\right)^{n+2} - b_n \frac{R_0}{r} [h_n^{(1)}(x) + xh_n^{(1)'}(x)] \right] P_n^1(\mu). \qquad (A24)$$

Substituting Eqs. (A23) and (A24) into Eqs. (A21) and (A22) leads to the calculation of the Stokes drift velocity components:

$$\begin{aligned} v_{Sr}^{nm} = & -\frac{1}{2\nu R_0} [P_n^1(\mu)P_m^1(\mu) - n(n+1)P_n(\mu)P_m(\mu)] \mathrm{Re}\{S_{11}(x)\} \\ & -\frac{1}{2\nu R_0} [P_n^1(\mu)P_m^1(\mu) - m(m+1)P_n(\mu)P_m(\mu)] \mathrm{Re}\{S_{12}(x)\} \\ & +\frac{1}{2\nu R_0} P_n^1(\mu)P_m^1(\mu) \mathrm{Re}\{S_{13}(x)\}, \end{aligned} \qquad (A25)$$

$$\begin{aligned} v_{S\theta}^{nm} = & \frac{(n+1)}{2\nu R_0} P_n(\mu) P_m^1(\mu) \mathrm{Re}\{S_{21}(x)\} + \frac{(m+1)}{2\nu R_0} P_m(\mu) P_n^1(\mu) \mathrm{Re}\{S_{22}(x)\} \\ & + \frac{1}{2\nu R_0} \sqrt{1 - \mu^2} [P_n^1(\mu) P_m^{1'}(\mu) - P_m^1(\mu) P_n^{1'}(\mu)] \mathrm{Re}\{S_{23}(x)\}, \end{aligned} \qquad (A26)$$
where we introduce the functions

$$ S_{11}(x) = -\frac{m+1}{x^2}a_n a_m^* \left(\frac{\bar{x}}{x}\right)^{n+m+3} + m(m+1)\left(\frac{\bar{x}}{x^2}\right) b_n b_m^* h_n^{(1)\prime}(x) h_m^{(1)*}(x) \\ + a_m b_n^* \left(\frac{\bar{x}}{x}\right)^{m+2} (m+1) \left[\frac{m+1}{x^2}h_n^{(1)}(x) + \frac{1}{x}h_n^{(1)\prime}(x)\right]^* \quad (A27) $$

$$ S_{12}(x) = -\frac{n+1}{x^2}a_m a_n^* \left(\frac{\bar{x}}{x}\right)^{n+m+3} + n(n+1)\left(\frac{\bar{x}}{x^2}\right) b_m b_n^* h_m^{(1)\prime}(x) h_n^{(1)*}(x) \\ + a_n b_m^* \left(\frac{\bar{x}}{x}\right)^{n+2} (n+1) \left[\frac{n+1}{x^2}h_m^{(1)}(x) + \frac{1}{x}h_m^{(1)\prime}(x)\right]^* \quad (A28) $$

$$ S_{13}(x) = \frac{\bar{x}}{x^3} b_n b_m^* h_n^{(1)}(x) h_m^{(1)*}(x) [n(n+1) - m(m+1)] + a_n b_m^* \left(\frac{\bar{x}}{x}\right)^{n+2} h_m^{(1)*}(x) \frac{1}{x^2} [m(m+1) - n(n+1)] \\ + a_m b_n^* \left(\frac{\bar{x}}{x}\right)^{m+2} h_n^{(1)*}(x) \frac{1}{x^2} [n(n+1) - m(m+1)] \quad (A29) $$

$$ S_{21}(x) = \frac{\bar{x}}{x^3} \left[ a_n \left(\frac{\bar{x}}{x}\right)^{n+1} + nb_n h_n^{(1)}(x) \right] \left[ -(m+3)a_m \left(\frac{\bar{x}}{x}\right)^{m+1} + b_m \left\{ 2h_m^{(1)} - x^2 h_m^{(1)\prime\prime}(x) \right\} \right]^* \quad (A30) $$

$$ S_{22}(x) = S_{21}(x)|_{n \leftrightarrow m} \quad (A31) $$

$$ S_{23}(x) = \frac{\bar{x}}{x^3} \left[ a_n \left(\frac{\bar{x}}{x}\right)^{n+1} - b_n \left\{ h_n^{(1)}(x) + xh_n^{(1)\prime}(x) \right\} \right] \left[ a_m \left(\frac{\bar{x}}{x}\right)^{m+1} - b_m \left\{ h_m^{(1)}(x) + xh_m^{(1)\prime}(x) \right\} \right]^* \quad (A32) $$

Now the boundary conditions at the bubble surface can be applied. It is known that, for arbitrary surface periodic deformations, both normal velocity component and tangential stress of the Lagrangian streaming must vanish at the mean position of the interface. These two conditions are written in the following form:

$$ v_{Lr}^{nn} = \langle v_{2r}^{nn} \rangle + v_{Sr}^{nn} = 0 \quad \text{at} \ r = R_0, \quad (A33) $$

$$ \frac{1}{r} \frac{\partial v_{Lr}^{nn}}{\partial \theta} + \frac{\partial v_{L\theta}^{nn}}{\partial r} - \frac{v_{L\theta}^{nn}}{r} = 0 \quad \text{at} \ r = R_0. \quad (A34) $$

In order to apply the first boundary condition (Eq. (A33)), it is required to expand the radial component of the Stokes drift velocity over the same set of orthogonal $P_k(\mu)$ functions as the radial component of the Eulerian velocity. This is performed using the properties of double and triple overlapping integrals of Legendre polynomials, leading to

$$ v_{Sr}^{nm} = \operatorname{Re} \left\{ \sum_{k=n-m}^{n+m} T_k(x) P_k(\mu) \right\}, \quad (A35) $$

where the function $T_k(x)$ is defined by

$$ T_k(x) = -\frac{1}{2\nu R_0} [a_{knm} - n(n+1)b_{knm}] S_{11}(x) - \frac{1}{2\nu R_0} [a_{knm} - m(m+1)b_{knm}] S_{12}(x) + \frac{1}{2\nu R_0} a_{knm} S_{13}(x), \quad (A36) $$

and the following coefficients are introduced

$$ a_{knm} = -(2k+1)\sqrt{(n+1)n(m+1)m} \begin{pmatrix} k & m & n \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} k & m & n \\ 0 & 1 & -1 \end{pmatrix}, \quad (A37) $$

$$ b_{knm} = (2k+1) \begin{pmatrix} k & n & m \\ 0 & 0 & 0 \end{pmatrix}^2. \quad (A38) $$

In order to apply the second boundary condition (Eq. (A34)), it is required to expand the tangential component of the Stokes drift velocity over the same set of orthogonal $P_k^1(\mu)$ functions as the tangential component of the Eulerian velocity. This is performed using the properties of double and triple overlapping integrals of Legendre polynomials, leading to

$$ v_{S\theta}^{nm} = \operatorname{Re} \left\{ \sum_{k=n-m}^{n+m} U_k(x) P_k^1(\mu) \right\}, \quad (A39) $$

where the function $U_k(x)$ is defined by

$$ U_k(x) = \frac{n+1}{2\nu R_0} \alpha_{knm} S_{21}(x) + \frac{m+1}{2\nu R_0} \alpha_{kmn} S_{22}(x) + \frac{1}{2\nu R_0} S_{23}(x) \left[ \frac{m(m+1)}{2} \alpha_{kmn} - \frac{n(n+1)}{2} \alpha_{knm} - \frac{1}{2} \beta_{knm} + \frac{1}{2} \beta_{kmn} \right]. \quad (A40) $$
Substituting Eqs. (A13), (A14), (A35) and (A39) into the conditions (A33) and (A34), after a cumbersome but straightforward calculation, one obtains the coefficients

$$C_{k30} = -C_{k20}x_{n0}^{2k-1} + \frac{R_0}{2(2k+1)}x_{n0}^{k-1}\left[\frac{k+3}{k+1}T_k(x_{n0}) + U_k(x_{n0}) - x_{n0}U'_k(x_{n0})\right], \quad (A41)$$

$$C_{k40} = -C_{k10}x_{n0}^{2k+3} - \frac{R_0}{2(2k+1)}x_{n0}^{k+1}\left[\frac{k-2}{k}T_k(x_{n0}) + U_k(x_{n0}) - x_{n0}U'_k(x_{n0})\right]. \quad (A42)$$

To sum up, we have shown that the solution of Eq. (23) is given by Eq. (A1), which in turns leads to the calculation of the components of the Eulerian streaming velocity given by Eqs. (A13) and (A14). We have calculated all the quantities that appear in the above equations. In the course of this calculation, we have also calculated the Stokes drift velocity, which, when being added to the Eulerian streaming velocity, gives the Lagrangian streaming velocity.

[1] J. Kolb and W.L. Nyborg, "Small-scale acoustic streaming effects in liquid", The Journal of the Acoustical Society of America 26(5), 1237-1242 (1956).

[2] A.A. Doinikov, S. Cleve, G. Regnault, C. Mauger and C. Inserra "Acoustic microstreaming produced by non-spherical oscillations of a gas bubble. I. Case of modes 0 and m", Physical Review E 100, 033104 (2019).

[3] J. Wu and G. Du, "Streaming generated by a bubble in an ultrasound field", The Journal of the Acoustical Society of America 101, 1899-1907 (1997).

[4] M.S. Longuet-Higgins, "Viscous streaming from an oscillating spherical bubble", Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 454, 725-742 (1998).

[5] A.A. Doinikov and A. Bouakaz "Effect of a distant rigid wall on microstreaming generated by an acoustically driven gas bubble", Journal of Fluid Mechanics 742, 425 (2014).

[6] I. Lentacker, I. De Cock, R. Deckers, S. De Smedt and C. Moonen, "Understanding ultrasound induced sonoporation: definitions and underlying mechanisms", Advanced Drug Delivery Reviews 29, 550-562 (2016).

[7] S. Cleve, M. Guédra, C. Mauger, C. Inserra and P. Blanc-Benon, "Microstreaming around acoustically trapped, non-spherically oscillating microbubbles", Journal of Fluid Mechanics 875, 597-621 (2019).

[8] F. Mekki-Berrada, T. Combriat, P. Thibault and P. Marrottant "Interactions enhance the acoustic streaming around flattened microfluidic bubbles", Journal of Fluid Mechanics 797, 851-873 (2016).

[9] M. Guédra, C. Inserra, C. Mauger and B. Gilles, "Experimental evidence of nonlinear mode coupling between spherical and nonspherical oscillations of microbubbles", Physical Review E 94, 053115 (2016).

[10] A.O. Maksimov, "Viscous streaming from surface waves on the wall of acoustically-driven gas bubbles", European Journal of Mechanics-B/Fluids 26(1), 28-42 (2007).

[11] T.A. Spelman and E. Lauga, "Arbitrary axisymmetric steady streaming: Flow, force and propulsion", Journal of Engineering Mathematics 105(1), 31-65 (2017).

[12] A.A. Doinikov "Translational motion of a bubble undergoing shape oscillations", Journal of Fluid Mechanics 501, 1-24 (2004).

[13] S.J. Shaw "Translation and oscillation of a bubble un-der axisymmetric deformation", Physiccs of Fluids 18(7), 072104 (2006).

[14] A.A. Doinikov, S. Cleve, G. Regnault, C. Mauger and C. Inserra "Acoustic microstreaming produced by non-spherical oscillations of a gas bubble. II. Case of modes 1 and m", Physical Review E 100, 033105 (2019).

[15] C. Inserra, G. Regnault, S. Cleve, C. Mauger and A.A. Doinikov "Acoustic microstreaming produced by non-spherical oscillations of a gas bubble. III. Case of self-interacting modes n-n", Physical Review E 101, 013111 (2020).

[16] M.S. Plesset "On the stability of fluid flows with spherical symmetry", Journal of Applied Physics 25, 96 (1954).

[17] M. Abramowitz and I.N. Stegun "Handbook of Mathematical Functions" (Dover Publications, New York, 1965).

[18] S.H. Dong and R. Lemus "The overlap integral of three associated Legendre polynomials", Applied Mathematics Letters 15, 541-546 (2002).

[19] See Supplemental Material at [URL will be inserted by publisher] for a MATLAB code that implements the calculation of the Lagrangian streaming velocity for the n-m case (MainProgram_Case_n.m).

[20] F. Reuter and R. Mettin "Mechanisms of single bubble cleaning", Ultrasonics Sonochemistry 29, 550-562 (2016).

[21] I. Lentacker, I. De Cock, R. Deckers, S.C. De Smedt and C.T.W. Moonen "Understanding ultrasound induced sonoporation: Definitions and underlying mechanisms", Advanced Drug Delivery Reviews 72(15), 49-64 (2014).

[22] Y. Xie, N. Nama, P. Li, Z. Mao, P.H. Huang, C. Zhao, F. Costanzo and T.J. Huang "Probing cell deformability via acoustically actuated bubbles", Small 12(7), 902-910 (2016).

[23] G. Lajoinie, Y. Luan, E. Gelderblom, B. Dollet, F. Mastik, H. Dewitte, I. Lentacker, N. De Jong and M. Versluis "Non-spherical oscillations drive the ultrasound-mediated release from targeted microbubbles", Communication Physics 1, 22 (2018).

[24] K.B. Bader, G. Bouchoux and C.K. Holland "Sonothrombolysis". In: Escoffre JM., Bouakaz A. (eds) Therapeutic Ultrasound. Advances in Experimental Medicine and Biology, vol 880. Springer, Cham. (2016)

[25] J.T. Sutton, K.J. Haworth, G. Pyne-Geithman and C.K. Holland "Ultrasound-mediated drug delivery for cardio-
vascular disease", Expert Opinion on Drug Delivery
10(5), 573-592 (2014).

[26] W.E. Boyce and R.C. DiPrima "Elementary Differential

Equations and Boundary Value Problems" (Wiley, New
York, 2001).