# Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. III. Case of self-interacting modes n-n

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. III. Case of self-interacting modes n-n. Physical Review E, 2020, 101, pp.013111. 10.1103/PhysRevE.101.013111. hal-03234699

HAL Id: hal-03234699

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

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. III.
Case of self-interacting modes *n*-n.

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: January 10, 2020)

This paper is the continuation 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 an acoustically driven 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. After solving these general equations for some particular cases of modal interactions (cases 0-*n*, 1-1 and 1-*m*), in this paper the general equations are solved analytically for the case that acoustic microstreaming results from the self-interaction of an arbitrary surface mode *n* ≥ 1. 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 acoustic microstreaming. No restrictions are imposed on the ratio of the bubble radius to the viscous penetration depth. The self-interaction results in specific streaming patterns: a large-scale cross pattern and small recirculation zones in the vicinity of the bubble interface. Particularly the spatial organization of the recirculation zones is unique for a given surface mode and therefore appears as a signature of the *n*-*n* interaction. Experimental streaming patterns related to this interaction are obtained and good agreement is observed with the theoretical model.

# I. INTRODUCTION

Applications of bubble-induced phenomena mainly concern medical therapy and material surface cleaning. In both cases, the immediate candidate usually considered for damage of substrate material or biological cell is the liquid jet impact and shock wave induced by bubble collapses [1, 2]. However the existence of additional mechanisms that result in controlled and weak damages have recently received growing attention. Such mechanisms are induced by highly nonlinear oscillatory behavior of bubbles, and bubble-mediated liquid flows. Applications of soft interaction of bubbles with surrounding media concern the bubble penetration within cracks between substrate and contamination [3], the permeabilization of biological cells [4] and the blood-brain barrier opening [5], to name a few. Among the different mechanical effects of steady-oscillating microbubbles, microstreaming is commonly considered as an important one due to the associated shear stress applied to surrounding solid or elastic walls. Microstreaming is a slow mean flow induced by a fast oscillating body [6]. In opposition to acoustic streaming, caused by the attenuation of an acoustic wave in the fluid, microstreaming is driven by streaming inside the oscillatory boundary layer around the bubble, 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, including the radial mode, the translational one, and the surface modes. Most of the available theoretical studies restrict their analysis to the simplest analytical cases. Davidson and Riley [7] are the first to consider the case of a translating bubble, where microstreaming is produced by the dipole (translation) mode alone. In the present study, we shall refer to this study as the microstreaming induced by the self-interacting translation mode, or 1-1 interaction. Wu and Du [8] and Longuet-Higgins [9] considered the case of a translating and pulsating bubble, where streaming is produced by the dipole mode and the monopole (pulsation) mode. In more recent studies, attempts have been made to consider shape modes of higher order [10, 11, 16]. Doinikov and Bouakaz [10] and Spelman and Lauga [11] developed theories that include modes of all orders but their calculations assume that all the modes oscillate at the same frequency, which means that the parametric generation of shape modes is disregarded. This makes these theories difficult to apply to real experimental conditions, such as these reported by Guédra et al. [12, 13] and Cleve et al. [14, 15], where predominant oscillatory modes are generated parametrically and hence oscillate at half the driving frequency. Only Maksimov [16] derived an asymptotic solution for microstreaming generated by a parametrically excited shape mode of order *n*, assuming that the amplitude of the shape mode is much greater than the one of the pulsation mode so that the contribution of the former to the microstreaming is dominant. Maksimov's analysis corresponds to the self-interacting nonspherical mode, or *n*-*n* interaction, in the

<footnote>* claude.inserra@inserm.fr</footnote>
present model. However, Maksimov's model is developed to describe liquid flows induced by millimeter-sized bubble driven at a frequency of a few kilohertz [17]. This solution is therefore not valid for microfluidic and biomedical applications, where one has to deal with micron-sized bubbles driven at high kHz and MHz frequencies [18]. It should also be emphasized that the above-mentioned theoretical studies assume that the bubble radius $R_0$ is much greater than the viscous penetration depth $\delta$. This means that viscous effects are considered to be sufficiently small and to verify the condition $\delta/R_0 \ll 1$. As a consequence, the above-mentioned theoretical works are performed approximately, 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 [19, 20], we developed a general theory for modeling the velocity field of acoustic microstreaming produced by axisymmetric nonspherical oscillations of an acoustically-driven gas bubble. A first study [19] 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 a second study 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 present study we derive the exact analytical solution of the Lagrangian streaming velocity induced by the self-interacting mode $n$, without restriction on the bubble radius to the viscous penetration depth. In Section II the solution for the Lagrangian streaming velocity is derived, compared to previous theories available in the literature, and numerical examples are provided. In Section III experimental streaming patterns exhibiting the specific signature of $n-n$ interaction are analyzed.

## 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

<figure><img src="image_2.png" /><figcaption>FIG. 1. Geometry of the system under study. (a) Three-dimensional representation of the zonal spherical harmonics (here mode 4 is shown), where $z$ is the axis of axial symmetry. (b) Axial symmetry allows using polar coordinates $(r,\theta)$ to parametrize the bubble interface $r_s$.</figcaption></figure>

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 equated to the normal component of the fluid particle velocity at the interface. Secondly, 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. (7) that describes the time-independent velocity field of acoustic microstreaming produced by the bubble oscillations. The solution of Eq. (7) 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 Appendix A.

### A. Linear solutions

The bubble oscillation is decomposed over $N$ axisymmetric surface modes, corresponding to the basis of Legendre polynomials [21]. The oscillation frequencies of the modes may differ between each other 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 tangential, $v_{1\theta}$, components of the first-order liquid velocity [19]:

$$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 [22]. 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. Their expression is given by [19]

$$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 (4)$$

$$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 (5)$$

## B. Acoustic microstreaming produced by the mode n alone

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 linear solution of first-order fluid velocity, and averaging it over time. Time averaging leads to the result that nonzero contributions to acoustic streaming can come either from pairs of modes that oscillate at the same frequency or from the interaction of a mode with itself. We already provided exact analytical solutions for the cases of microstreaming induced by the interaction of the radial mode [19] or the translation mode $n = 1$ [20] with any arbitrary surface mode $n$. Here the case of the self-interaction of a surface mode is considered. According to the theory developed in Part I [19], in the self-interacting case, the Eulerian streaming velocity is represented by

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

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

$$D^2 \langle \psi_2^{nn} \rangle = \frac{n+1}{2\nu r^2} \left( \frac{R_0}{r} \right)^{n+1} P_n(\mu) P_n^1(\mu) \mathrm{Re} \{ \alpha_1(x_n) \} - \frac{\sqrt{1-\mu^2}}{\nu r^2} P_n^1(\mu) P_n^{1\prime}(\mu) \mathrm{Re} \{ \alpha_2(x_n) \}, \quad (7)$$

where $D$ is the linear operator given by

$$\begin{aligned} 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) \\ &\quad + \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] \end{aligned} \quad (8)$$

and the functions $\alpha_{1,2}(x_n)$ are

$$\alpha_1(x_n) = k_n^2 a_n b_n^* \left[ (n+1)h_n^{(1)}(x_n) - x_n h_n^{(1)\prime}(x_n) \right]^*, \quad (9)$$

$$\alpha_2(x_n) = k_n^2 h_n^{(1)*}(x_n) \left[ a_n b_n^* \left( \frac{R_0}{r} \right)^{n+1} - b_n b_n^* x_n h_n^{(1)\prime}(x_n) \right], \quad (10)$$

where the asterisk denotes the complex conjugate. In order to find an exact solution for the equation ruling the vector potential $\psi_2^{nn}$, one may simplify the angular dependency in Equation (7) which is decomposed over two terms involving Legendre polynomials and associated Legendre polynomials. This is done through the transformation

$$\sqrt{1 - \mu^2} P_n^1(\mu) P_n^{1\prime}(\mu) = n^2 P_n(\mu) P_n^1(\mu) - \sum_{i=1}^{n-1} (2i + 1) P_i(\mu) P_i^1(\mu), \quad (11)$$

obtained when using common identities of associated Legendre functions [22]. Substituting Eq. (11) into Equation (7) and using the set of equations (4, 5) for the coefficients $a_n$ and $b_n$, one obtains

$$D^2 \langle \psi_2^{nn}(x_n, \mu) \rangle = \frac{|b_n|^2 k_n^4}{\nu} \mathrm{Re} \left\{ \sum_{k=1}^{n} \left[ \delta_{kn} G_1(x_n) + (1 - \delta_{kn})(2k+1)G_2(x_n) \right] P_k(\mu) P_k^1(\mu) \right\} \quad (12)$$

where $\delta_{kn}$ is the Kronecker delta, and the functions $G_1(x_n)$ and $G_2(x_n)$ are defined by

$$G_1(x_n) = \frac{n^2}{x_n} h_n^{(1)*}(x_n) h_n^{(1)\prime}(x_n) - \frac{\bar{x}_n^2 h_n^{(1)\prime\prime}(\bar{x}_n) - (n^2 + n - 2)h_n^{(1)}(\bar{x}_n)}{4(n+2)x_n^2} \left( \frac{\bar{x}_n}{x_n} \right)^{n+1} \times \left[ (n^2 - 2n - 1)h_n^{(1)}(x_n) + (n+1)x_n h_n^{(1)\prime}(x_n) \right]^* \quad (13)$$
$$G_2(x_n) = \left[ \frac{\bar{x}_n^2 h_n^{(1)\prime\prime}(\bar{x}_n) - (n^2 + n - 2)h_n^{(1)}(\bar{x}_n)}{2(n+2)} \left( \frac{\bar{x}_n}{x_n} \right)^{n+1} - x_n h_n^{(1)\prime}(x_n) \right] \frac{h_n^{(1)*}(x_n)}{x_n^2}. \quad (14)$$

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

$$\langle \psi_2^{nn}(x, \mu) \rangle = \frac{|b_n|^2}{\nu} \mathrm{Re} \left\{ \sum_{k=1}^{n} F_k(x_n) P_k(\mu) P_k^1(\mu) \right\}, \quad (15)$$

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

$$\langle v_{2r}^{nn} \rangle = \frac{|b_n|^2}{\nu r} \mathrm{Re} \left\{ \sum_{k=1}^{n} F_k(x_n) \left[ (P_k^1(\mu))^2 - k(k+1)(P_k(\mu))^2 \right] \right\}, \quad (16)$$

$$\langle v_{2\theta}^{nn} \rangle = -\frac{|b_n|^2}{\nu r} \mathrm{Re} \left\{ \sum_{k=1}^{n} P_k(\mu) P_k^1(\mu) [F_k(x_n) + x F'_k(x_n)] \right\}, \quad (17)$$

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

In the process of calculating Eqs. (16) and (17), we have also obtained the following expressions for the components of the Stokes drift velocity [9] (see Appendix A):

$$v_{Sr}^{nn} = -\frac{(n+1)|b_n|^2}{2\nu R_0} \times \left[ (P_n^1(\mu))^2 - n(n+1)(P_n(\mu))^2 \right] \mathrm{Re}\{S_1(x_n)\}, \quad (18)$$

$$v_{S\theta}^{nn} = \frac{(n+1)|b_n|^2}{2\nu R_0} P_n(\mu) P_n^1(\mu) \mathrm{Re}\{S_2(x_n)\}, \quad (19)$$

where the functions $S_1(x_n)$ and $S_2(x_n)$ are defined by Eqs. (A33) and (A34). The summation of the Eulerian streaming velocity and the Stokes drift velocity provides the Lagrangian streaming velocity,

$$v_L^{nn} = \langle v_2^{nn} \rangle + v_S^{nn}, \quad (20)$$

for which the radial and tangential components can be calculated using Eqs. (16)-(19). It has been verified that the Lagrangian streaming velocity given by Eq. (20) corresponds to the one derived in our previous work [20] for the particular case $n = 1$.

### C. Numerical examples

Numerical simulations were performed at the following values of physical parameters: $\rho = 1000 \text{ kg/m}^3$, $\eta = 0.001 \text{ Pa} \cdot \text{s}$, $f = 50 \text{ kHz}$, and $R_0 = 50 \text{ }\mu\text{m}$. Figure 2 exemplifies Lagrangian streamline patterns produced by the self-interaction of modes 1 to 4, respectively. As one

can see, the main vortices present a two-scale pattern. In the bubble vicinity, the vortices have a form of lobes. According to the presented examples, the number of lobes is equal to $4n$ for the $n-n$ interaction. Far from the bubble, the vortices have a cross-like shape. This cross-like shape is similar to the one induced by the purely translating bubble, case 1-1 (Figure 2(a)). In the present analysis, a specific pattern is revealed and appears as a signature of the case of self-interaction. Indeed, when considering the interaction $0-n$ between the radial mode and a nonspherical mode [19], the streamline patterns looked like lobes whose number equals $2n$. When considering the interaction between the translational mode and any nonspherical mode $n$ [20], streamlines form lobes whose number is equal to $2(n-1)$.

### D. Comparison with previous theories

The present model is here compared to existing theoretical models, which take into account self-interaction of axisymmetric nonspherical oscillations. Such cases have been derived previously by Maksimov [16] and Spelman and Lauga [11]. Before comparing our results to the above-mentioned studies, it is worth noting that in both of them the physical situation is different from ours. Both studies assume that the bubble is fixed while the liquid oscillates around it, whereas we assume that the bubble is moving while the liquid at infinity is at rest. This means that, in their case, at infinity the amplitude of the first-order liquid velocity tends to a nonzero constant. In this sense, their derivations follow the same mathematical approach as the one proposed by Longuet-Higgins [9] who described the streaming velocity induced by purely translating bubble and coupled radial-translational oscillation. In addition, both theories assume that the viscous penetration length scale is small compared to the body size. Mathematically it means that their solutions are derived approximatively, in powers of the small parameter $\epsilon = \delta/R_0$, within the viscous boundary layer. A solution outside the boundary layer is derived, and the matching between the inner and outer solution provides the overall streaming velocity. The aim of the comparison we propose is hence to provide an idea of the difference between our solution and the results of previous studies available in the literature.

First our model is compared to the theory of Spelman and Lauga who derived a theory that includes axisymmetric modes of arbitrary orders, oscillating at the same frequency. Their 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)} \right) \left( \int_{\mu}^{1} P_k(x) dx \right), \quad (21)$$

where $T_k$ and $S_k$ are coefficients resulting from the matching between the inner and outer solution of the Eulerian streaming. In particular, their calculation de-
<figure><img src="image_3.png" /></figure>

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

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

<figure><img src="image_6.png" /><figcaption>FIG. 2. Numerical examples of streamline patterns produced by the self-interaction mode in cases: (a) 1-1, (b) 2-2, (c) 3-3, (d) 4-4.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 3. Comparison of the present model to the theory of Spelman and Lauga [11] for the self-interaction of mode 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/5$, as a function of the normalized distance $r/R_0$.</figcaption></figure>

pends on the investigated mode pair contribution and phase shift between modes. It is shown that, if the modal amplitudes $s_n$ 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$. This case naturally includes the case of a self-interacting mode. It is thus necessary to derive expressions to the upper order in the expansion, up to the third order of the ratio $\delta/R_0$ in their analysis. Figure 3 compares the streamline pattern and Lagrangian streaming velocities between the theory of Spelman and Lauga and the present model, in case of the 2-2 interaction (see Fig. 4 in [11]). The velocity components are normalized by the factor $\omega_n|s_n|^2/R_0$. The parameters are the same as in Fig. 2. The streamline patterns present a qualitative agreement, with a smaller radial extension of the vortices in the vicinity of the bubble in the present model. Concerning the radial and tangential components of the Lagrangian velocity (Fig. 3(c,d)), both models show a quantitative agreement, particularly far from the bubble interface. The present model also predicts a smaller radial velocity amplitude, confirming this smaller extension of the streamline pattern.

We secondly compare our model to the theory developed by Maksimov [16]. The Maksimov model considers acoustic streaming from a bubble undergoing a parametrically-excited surface mode, hence oscillating at $\omega/2$, where $\omega$ is the angular frequency of the driving ultrasound wave. This case corresponds exactly to the self-interacting mode $n-n$ discussed here, as no interaction arises from the radial and the nonspherical modes at the second order of approximation. Maksimov's model provides the following expressions for the radial and tangential components of the Lagrangian velocity in the $n-n$ interaction:

$$ u_r = \frac{(\omega/2)|s_n^2|}{2\pi R_0} \frac{2n+1}{4\pi} \frac{(n+2)(2n+1)(n+3)}{(n+1)(4n3)} \left(\frac{R_0}{r}\right)^{2n-1} \times \left[1 - \left(\frac{R_0}{r}\right)^2\right] \frac{\sin((2n+1)\theta)}{\sin\theta}, \quad (22) $$

$$ u_{\theta} = \frac{(\omega/2)|s_n^2|}{2\pi R_0} \frac{2n+1}{4\pi} \frac{(n+2)(2n-1)(n+3)}{(n+1)(4n3)} \left(\frac{R_0}{r}\right)^{2n+1} \times \left[1 - \frac{2n+1}{2n-1} \left(\frac{R_0}{r}\right)^2\right] \frac{\cos((2n+1)\theta)}{\sin\theta}. \quad (23) $$
<figure><img src="image_8.png" /><figcaption>FIG. 4. Comparison of the present model to the theory of Maksimov [16] for the self-interaction of nonspherical mode $n = 10$. Evolution of the normalized (a) radial and (b) tangential components of the Lagrangian streaming velocity given by both theories and calculated for the angle $\theta = \pi/5$, as a function of the normalized distance $r/R_0$.</figcaption></figure>

In comparison to our derivation, the Lagrangian velocities obtained in Eqs. (22) and (23) contain simplified dependences over the angle $\theta$. Such simplification arises because Maksimov considered high-order distorsion modes $n \gg 1$. This angular dependence nevertheless describes a 4n-lobes pattern as the ones shown in Fig. 2. In addition, Maksimov hypothesizes that the streaming is dominated by the surface wave contribution in the region of the liquid from the bubble wall out to distances equivalent to several wavelengths of that surface wave. This explains why the streamline pattern derived by Maksimov does not contain the cross-like shape at larger distances from the bubble interface. Figure 4 compares the dependence on $r$ given by Eqs. (22) and (23) to the results of our theory for the particular angle $\theta = \pi/5$. The velocity components are normalized by the factor $\omega|s_n|^2/R_0$. The parameters are the same as in Fig. 2. As one can see, the theory of Maksimov predicts a greater velocity amplitude. However, it should be emphasized that two different physical situations are compared.

### III. EXPERIMENTAL RESULTS

In order to illustrate experimentally these theoretical findings, one needs to observe the microstreaming around a free (far from boundaries) non-spherically oscillating bubble. In addition, both the bubble interface dynamics and the surrounding fluid motion must be captured at the same time. This experimental challenge has been recently overcome, and is explained in more detail in previous works [15, 19]. A short summary of the here relevant points is given in the following. The experiments rely on the nucleation, trapping and excitation of a unique bubble into a 31 kHz standing wave field established within a 8 cm cubic acoustic levitation chamber. In order to control the bubble nonspherical oscillations, the bubble

coalescence technique [14] has demonstrated its interest in establishing steady-state oscillation of an axisymmetric surface mode with, in addition, a controlled symmetry axis. By varying the bubble size and acoustic pressure, several surface modes can be obtained on their first or secondary parametric resonance. The establishment of a steady-state regime of nonspherical oscillations is mandatory for the generation of a microstreaming pattern associated to a captured temporal dynamics of the bubble interface. By decomposing the bubble interface over the set of Legendre polynomials, similarly to Eq. (1), the modal amplitudes $s_n$ of each mode, including the radial ($n = 0$), translational ($n = 1$) and nonspherical ($n > 1$) modes, are assessed. These modal coefficients are associated to a given microstreaming pattern, and can then be used as input parameters in the developed theory.

In our experiment, surface modes are parametrically excited. They hence oscillate primarily at half the driving frequency if excited on their first parametric resonance. If the energy transfer between modes (including translational and nonspherical modes $m \neq n$) is weak enough, then the modal content of the bubble interface dynamics contains one (volumic) oscillation at $\omega$ and one (nonspherical) oscillation at $\omega/2$. However, it is known that contributions to microstreaming only come from pairs of modes oscillating at the same frequency or from self-interacting nonspherical modes. In the case discussed here, recalling that radial oscillations do not produce streaming flow by itself [19], then the predominant contribution to the liquid velocity field is the self-interaction of mode $n$. In our experimental configuration, parametrically-excited modes on their first resonance, hence oscillating at $\omega/2$, are modes $n = 2$ or $n = 3$.

Figure 5 presents the experimental and numerical results corresponding to the self-interacting cases 2-2 and 3-3. In Fig. 5(a), a bubble of equilibrium radius 48 µm
## 2-2 interaction

<figure><img src="image_9.png" /><figcaption>FIG. 5. Experimental and numerical investigation of the microstreaming patterns induced by a self-interacting mode. First line: A bubble of equilibrium radius 48 µm is driven by a 31.25 kHz ultrasound field at the acoustic pressure amplitude 20 kPa. (a) Streak photography of the resulting microstreaming pattern (image size 620 × 620 µm). (b) Numerical streaming pattern of the 2-2 interaction. Second line: A bubble of equilibrium radius 64.5 µm is driven by a 31.25 kHz ultrasound field at the acoustic pressure amplitude 12 kPa. (c) Streak photography of the resulting microstreaming pattern (image size 700 × 700 µm). (d) Numerical streaming pattern of the 3-3 interaction.</figcaption></figure>

is driven at the acoustic amplitude 20 kPa. This set of parameters leads to a predominant mode 2 oscillation. High-speed imaging of the bubble temporal dynamics reveals that mode 2 oscillates with an amplitude $s_2 = 7$ µm, such that $s_2/R_0 \sim 0.15$. The experimental streamline pattern reveals a cross-like shape with two small recirculation zones in the vicinity of the bubble interface, each of them composed of two vortices. The corresponding numerical streamline pattern, when applying the experimental parameters $s_2$ in the theoretical modeling, is presented in Fig. 5(b). This pattern also presents 8 lobes with a cross-like shape far from the bubble interface. In particular, the asymmetry in the lobe ordering

is respected. In Fig. 5(c), a bubble of equilibrium radius 64.5 µm is driven at the acoustic amplitude 12 kPa. This set of parameters leads to a predominant mode 3 oscillation. High-speed imaging of the bubble temporal dynamics reveals that mode 3 oscillates with an amplitude $s_3 = 20$ µm, such as $s_3/R_0 \sim 0.3$. The experimental streamline pattern reveals a cross-like shape with four vortices in the vicinity of the bubble interface. The corresponding numerical streamline pattern, when applying the experimental parameters $s_3$ in the theoretical modeling, is presented in Fig. 5(d). This pattern presents 12 lobes and a cross-like shape far from the bubble interface. Clearly, the numerical pattern captures more
details in the recirculation zones than one can resolve experimentally. However, a good agreement is observed in the overall pattern, as well as in the location of the four main vortices in the vicinity of the bubble interface.

## IV. CONCLUSIONS

In the present paper, the general theory developed in a previous study [19] has been applied to the case that acoustic microstreaming results from the self-interaction of an arbitrary surface mode $n \ge 1$. 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 self-interaction of a mode results in a specific streamline pattern. This pattern exhibits a cross-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 $4n$ for the self-interacting case $n-n$. In opposite to the $2n$ lobes obtained in the $0-n$ interaction case and the $2(n-1)$ lobes for the $1-n$ interaction case, this specific pattern appears as a signature for the case of self-interaction. Experimental streamline patterns are presented for parametrically excited surface modes of a bubble oscillating in an unbounded liquid. These experimental patterns exhibit the particular signature of self-interacting nonspherical mode, in agreement with our numerical simulations.

## 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).

## Appendix A: Solution of Eq. 12

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

$$ \langle \psi_2^{nn}(x, \mu) \rangle = \frac{|b_n|^2}{\nu} \mathrm{Re} \left\{ \sum_{k=1}^{n} F_k(x) P_k(\mu) 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^{nn}(x, \mu) \rangle = \frac{|b_n|^2 k_n^4}{\nu} \mathrm{Re} \left\{ \sum_{k=1}^{n} \left[ \sum_{l=k}^{n} \left[ g_l^{(1)} \delta_{kl} + g_l^{(2)} \beta_{lk} + g_l^{(3)} \left( \sum_{m=k}^{l} \beta_{lm} \beta_{mk} \right) \right] \right] P_k P_k^1 \right\}, \quad (A2) $$

where the $g_l$ functions are defined by

$$ g_l^{(1)} = \frac{1}{x^4} \left[ 4x^3 \frac{d^3 F_k}{dx^3} + x^4 \frac{d^4 F_k}{dx^4} \right], \quad (A3) $$

$$ g_l^{(2)} = \frac{1}{x^4} \left[ 2F_k(x) + 2x^2 \frac{d^2 F_k}{dx^2} \right], \quad (A4) $$

$$ g_l^{(3)} = \frac{1}{x^4} F_k(x), \quad (A5) $$

and by introducing the coefficients $\beta_{nk} = 2(2k+1)[1-(1+k)\delta_{kn}]$. After substituting Eq. (A2) into Eq. (12), the identification for each term over the functions $P_k P_k^1$ that forms a set of linearly independent functions provides the following equation for the $F_k(x)$ functions:

$$ \begin{aligned} & \frac{d^4 F_k}{dx^4} + \frac{4}{x} \frac{d^3 F_k}{dx^3} - \frac{4k(2k+1)}{x^2} \frac{d^2 F_k}{dx^2} \\ & \qquad + \frac{4k(2k+1)(2k^2+k-1)}{x^4} F_k \\ &= \delta_{kn} G_1(x) + (1-\delta_{kn})(2k+1)G_2(x) \\ & \qquad - \sum_{l=k+1}^{n} \left[ \frac{4(2k+1)}{x^2} \frac{d^2 F_l}{dx^2} - \frac{4(2k+1)}{x^4} [l(2l+1)] \right. \\ & \qquad \left. + k(2k+1) - (1-\delta_{k+1,l})E_{kl} - 1 \right] F_l, \quad (A6) \end{aligned} $$

where the coefficient $E_{kl} = \sum_{m=k+1}^{l-1}(2m+1)$ has been introduced. At a first sight Eq. (A6) is a set of differential equations that couple the whole set of $F_k$ functions. However, the equation for $k=n$ involves only the $F_n$ term, that can therefore be calculated. The equation for $k=n-1$ involves the two terms $F_n$ and $F_{n-1}$, and the knowledge of the $F_n$ term allows calculating $F_{n-1}$. Following this procedure we can sequentially calculate all the functions $F_k$. At all stages of this process we need to solve the equation

$$ \begin{aligned} & \frac{d^4 F_k}{dx^4} + \frac{4}{x} \frac{d^3 F_k}{dx^3} - \frac{4k(2k+1)}{x^2} \frac{d^2 F_k}{dx^2} \\ & \qquad + \frac{4k(2k+1)(2k^2+k-1)}{x^4} F_k = H_k(x), \quad (A7) \end{aligned} $$

where the right-hand side function $H_k(x)$ is defined by

$$ \begin{aligned} H_k(x) = & \delta_{kn} G_1(x) + (1 - \delta_{kn})(2k+1)G_2(x) \\ & - \sum_{l=k+1}^{n} \left[ \frac{4(2k+1)}{x^2} \frac{d^2 F_l}{dx^2} - \frac{4(2k+1)}{x^4} [l(2l+1)] \right. \\ & \left. + k(2k+1) - (1 - \delta_{k+1,l})E_{kl} - 1 \right] F_l. \quad (A8) \end{aligned} $$

Eq. (A7) can be solved by the method of variation of parameters, also known as the Lagrange method [23].
According to this method, we first need to find solutions to a homogeneous equation that corresponds to the left-hand side of Eq. (A7),

$$ \frac{d^4 F_k}{dx^4} + \frac{4}{x} \frac{d^3 F_k}{dx^3} - \frac{4k(2k+1)}{x^2} \frac{d^2 F_k}{dx^2} + \frac{4k(2k+1)(2k^2+k-1)}{x^4} F_k = 0. \quad (A9) $$

The solutions are sought in the form $x^\lambda$. Substitution of $x^\lambda$ into Eq. (A9) leads to a polynomial of fourth order in $\lambda$,

$$ \lambda(\lambda - 1)(\lambda - 2)(\lambda + 1) - 4k(2k + 1)\lambda(\lambda - 1) + 4k(2k + 1)(2k^2 + k - 1) = 0. \quad (A10) $$

The roots of this polynomial are $-(2k-1)$, $-(2k+1)$, $2k$, $2k+2$. Therefore, the general solution of Eq. (A9) is written as

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

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

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

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). \end{aligned} \quad (A13) $$

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

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

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

$$ C_{k1}(x) = C_{k10} + \frac{1}{2(4k+1)(4k+3)} \int_{\bar{x}_n}^{x} s^{-2k+1} H_k(s) ds, \quad (A15) $$

$$ C_{k2}(x) = C_{k20} - \frac{1}{2(4k-1)(4k+1)} \int_{\bar{x}_n}^{x} s^{-2k+3} H_k(s) ds, \quad (A16) $$

$$ C_{k3}(x) = C_{k30} + \frac{1}{2(4k-1)(4k+1)} \int_{\bar{x}_n}^{x} s^{2k+2} H_k(s) ds, \quad (A17) $$

$$ C_{k4}(x) = C_{k40} - \frac{1}{2(4k+1)(4k+3)} \int_{\bar{x}_n}^{x} s^{2k+4} H_k(s) ds, \quad (A18) $$

where the constants $C_{ki0}$, $i=1-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}^{nn} \rangle &= -\frac{1}{r} \frac{\partial}{\partial \mu} \left\{ \langle \psi_2^{nn} \rangle \sqrt{1 - \mu^2} \right\} && (A19) \\ &= \frac{|b_n|^2}{\nu r} \operatorname{Re} \left\{ \sum_{k=1}^{n} F_k(x) \left[ (P_k^1(\mu))^2 - k(k+1)(P_k(\mu))^2 \right] \right\}, \end{aligned} $$

$$ \begin{aligned} \langle v_{2\theta}^{nn} \rangle &= -\frac{1}{r} \frac{\partial}{\partial x} \{x \langle \psi_2^{nn} \rangle\} && (A20) \\ &= -\frac{|b_n|^2}{\nu r} \operatorname{Re} \left\{ \sum_{k=1}^{n} P_k(\mu) P_k^1(\mu) [F_k(x) + x F'_k(x)] \right\}, \end{aligned} $$

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

$$ \begin{aligned} F'_{k}(x) = 2(k+1)C_{k1}(x)x^{2k+1} + 2kC_{k2}(x)x^{2k-1} && (A21) \\ - (2k-1)C_{k3}(x)x^{-2k} - (2k+1)C_{k4}(x)x^{-2k-2}. && \end{aligned} $$

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(4k+1)(4k+3)} \int_{x_{n0}}^{\infty} s^{-2k+1} H_k(s) ds, \quad (A22) $$

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

In order to calculate the coefficients $C_{k30}$ and $C_{k40}$, boundary conditions at the bubble surface have to be applied. Equations (A19) and (A20) 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, \quad (A24) $$

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

$$ v_S = \left\langle \left( \int v_1 dt \cdot \nabla \right) v_1 \right\rangle_{nn}, \quad (A25) $$

where $v_1$ is the first-order liquid velocity. For the case of a self-interacting mode $n$, Eq. (A25) reduces to

$$ v_S^{nn} = \frac{1}{2\omega_n} \operatorname{Re} \{ \iota (v_{1n} \cdot \nabla) v_{1n}^* \}, \quad (A26) $$

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

$$ v_{Sr}^{nn} = \frac{1}{2\omega_n} \operatorname{Re} \left\{ \iota v_{1nr} \frac{\partial v_{1nr}^*}{\partial r} + \frac{\iota v_{1n\theta}}{r} \frac{\partial v_{1nr}^*}{\partial \theta} \right\}, \quad (A27) $$
$$v_{S\theta}^{nn} = \frac{1}{2\omega_n} \mathrm{Re} \left\{ \imath v_{1nr} \left( \frac{\partial v_{1n\theta}}{\partial r} - \frac{v_{1n\theta}}{r} \right)^* + \frac{\imath v_{1n\theta}}{r} \frac{\partial v_{1n\theta}^*}{\partial \theta} \right\}. \quad (A28)$$

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 [19]:

$$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), \quad (A29)$$

$$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} \left[ h_n^{(1)}(x) + x h_n^{(1)\prime}(x) \cdot \right] \right] P_n^1(\mu) \quad (A30)$$

Substituting Eqs. (A29) and (A30) into Eqs. (A27) and (A28) leads to the calculation of the Stokes drift velocity components:

$$v_{Sr}^{nn} = -\frac{(n+1)|b_n|^2}{2\nu R_0} \left[ (P_n^1(\mu))^2 - n(n+1)(P_n(\mu))^2 \right] \mathrm{Re}\{S_1(x)\}, \quad (A31)$$

$$v_{S\theta}^{nn} = \frac{(n+1)|b_n|^2}{2\nu R_0} P_n(\mu) P_n^1(\mu) \mathrm{Re}\{S_2(x)\}, \quad (A32)$$

where we introduce the functions:

$$S_1(x) = \frac{n\bar{x}_n}{x^2} h_n^{(1)*}(x) h_n^{(1)\prime}(x) + \frac{[\bar{x}_n^2 h_n^{(1)\prime\prime}(\bar{x}_n) - (n^2 + n - 2)h_n^{(1)}(\bar{x}_n)]}{2(n+2)} \left(\frac{\bar{x}_n}{x}\right)^{n+2} \left[ \frac{n+1}{x^2} h_n^{(1)}(x) + \frac{1}{x} h_n^{(1)\prime}(x) \right]^*, \quad (A33)$$

$$S_2(x) = \frac{n\bar{x}_n}{x} h_n^{(1)*}(x) h_n^{(1)\prime\prime}(x) + \frac{[\bar{x}_n^2 h_n^{(1)\prime\prime}(\bar{x}_n) - (n^2 + n - 2)h_n^{(1)}(\bar{x}_n)]^*}{2(n+2)} \left(\frac{\bar{x}_n}{x}\right)^{n+2} \left[ h_n^{(1)\prime\prime}(x) - \frac{(n+1)(n+2)}{x^2} h_n^{(1)}(x) \right]. \quad (A34)$$

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 (A35)$$

$$\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 (A36)$$

Substituting Eqs. (A19), (A20), (A31) and (A32) into the conditions (A35) and (A36), after a cumbersome but straightforward calculation, one obtains the coefficients:

$$C_{k30} = -C_{k20} \bar{x}_n^{4k-1} - (1 - \delta_{kn}) \frac{2k+1}{4k+1} \bar{x}_n^{2k-1} \sum_{l=k+1}^{n} \left[ C_{l10} \bar{x}_n^{2l+2} + C_{l20} \bar{x}_n^{2l} + \frac{C_{l30}}{\bar{x}_n^{2l-1}} + \frac{C_{l40}}{\bar{x}_n^{2l+1}} \right] \\ - \frac{(n+1)\bar{x}_n^{2k-1}}{4(4k+1)} \{ \delta_{kn} [\bar{x}_n S_2'(\bar{x}_n) - S_2(\bar{x}_n) - 2n(2n+3)S_1(\bar{x}_n)] \\ - (1 - \delta_{kn})2(2k+1)S_1(\bar{x}_n) \}, \quad (A37)$$

$$C_{k40} = -C_{k10} \bar{x}_n^{4k+3} + (1 - \delta_{kn}) \frac{2k+1}{4k+1} \bar{x}_n^{2k+1} \sum_{l=k+1}^{n} \left[ C_{l10} \bar{x}_n^{2l+2} + C_{l20} \bar{x}_n^{2l} + \frac{C_{l30}}{\bar{x}_n^{2l-1}} + \frac{C_{l40}}{\bar{x}_n^{2l+1}} \right] \\ + \frac{(n+1)\bar{x}_n^{2k+1}}{4(4k+1)} \{ \delta_{kn} [\bar{x}_n S_2'(\bar{x}_n) - S_2(\bar{x}_n) - 2(n-1)(2n+1)S_1(\bar{x}_n)] \\ - (1 - \delta_{kn})2(2k+1)S_1(\bar{x}_n) \}. \quad (A38)$$
To sum up, we have shown that the solution of Eq. (12) is given by Eq. (A1), which in turns leads to the calculation of the components of the Eulerian streaming velocity given by Eqs. (A19) and (A20). 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] C.D. Ohl, M. Arora, R. Ikink, N. De Jong, M. Versluis, M. Delius and D. Lohse, "Sonoporation from jetting cavitation bubbles", Biophysical Journal **91**(11), 4285-4295 (2006).

[2] P. Prentice, A. Cuschieri, K. Dholakia, M. Prausnitz and P. Campbell, "Membrane disruption by optically controlled microbubble cavitation", Nature Physics **1**(2), 107 (2005).

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

[4] 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).

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

[6] 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).

[7] B.J. Davidson and N. Riley, "Cavitation microstreaming", Journal of Sound and Vibration **15**(2), 217-233 (1971).

[8] 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).

[9] 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).

[10] A.A. Doinikov and A. Bouakaz, "Acoustic microstreaming around a gas bubble", The Journal of the Acoustical Society of America **127**(2), 703-709 (2010).

[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] 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).

[13] M. Guédra, S. Cleve, C. Mauger, P. Blanc-Benon and C. Inserra, "Dynamics of nonspherical microbubble oscillations above instability threshold", Physical Review E **96**, 063104 (2017).

[14] S. Cleve, M. Guédra, C. Inserra, C. Mauger and P. Blanc-Benon, "Surface modes with controlled axisymmetry triggered by bubble coalescence in a high-amplitude acoustic field", Physical Review E **98**, 033115 (2018).

[15] 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).

[16] 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).

[17] T.G. Leighton, "From seas to surgeries, from babbling brooks to baby scans: The acoustics of gas bubbles in liquids", International Journal of Modern Physics B **18**(25), 3267-3314 (2004).

[18] P. Marmottant, J. Raven, H. Gardeniers, J. Bomer and S. Hilgenfeldt "Microfluidics with ultrasound-driven bubbles", Journal of Fluid Mechanics **568**, 109-118 (2006).

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

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

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

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

[23] W.E. Boyce and R.C. DiPrima "Elementary Differential Equations and Boundary Value Problems" (Wiley, New York, 2001).