# Acoustic microstreaming produced by two interacting gas bubbles undergoing axisymmetric shape oscillations

Alexander A. Doinikov, Gabriel Regnault, Cyril Mauger, Philippe

Blanc-Benon, Claude Inserra

► To cite this version:

Alexander A. Doinikov, Gabriel Regnault, Cyril Mauger, Philippe Blanc-Benon, Claude Inserra. Acoustic microstreaming produced by two interacting gas bubbles undergoing axisymmetric shape oscillations. Journal of Fluid Mechanics, 2022, 931, pp.A19. 10.1017/jfm.2021.926. hal-03775117

HAL Id: hal-03775117

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

Submitted on 10 Dec 2024

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 two interacting gas bubbles undergoing axisymmetric shape oscillations

*Accepted Manuscript Version*

Alexander A. Doinikov,¹ Gabriel Regnault,² Cyril Mauger,² Philippe Blanc-Benon² and
Claude Inserra³

¹Univ Lyon, Collegium de Lyon, Ecole Centrale de Lyon, INSA Lyon, CNRS, LMFA UMR 5509, F-69002 Lyon, France

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

³Univ Lyon, Université Claude Bernard Lyon 1, Centre Léon Bérard, INSERM, UMR 1032, LabTAU, F-69003 Lyon, France
An analytical theory is developed that describes acoustic microstreaming produced by two interacting bubbles. The bubbles are assumed to undergo axisymmetric oscillation modes, which can include radial oscillations, translation and shape modes. Analytical solutions are derived in terms of complex amplitudes of oscillation modes, which means that the modal 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 radii to the viscous penetration depth and the distance between the bubbles. The interaction between the bubbles is considered both when the linear velocity field is calculated and when the second-order velocity field of acoustic microstreaming is calculated. Capabilities of the analytical theory are illustrated by computational examples.

**Key words:** acoustic streaming, bubble dynamics, shape modes
# 1. Introduction

Much theoretical and experimental work has been done on the study of steady vortical flows, called acoustic microstreaming, that are produced by an oscillating bubble. A detailed discussion of this work is provided in a series of recent publications, which are devoted to acoustic microstreaming generated by nonspherical oscillations of a single gas bubble in an unbounded liquid (Doinikov et al. 2019a,b; Inserra et al. 2020a,b). Recently, an excellent review of earlier studies on acoustic microstreaming has also been performed by Mobadersany & Sarkar (2019). It is pointed out in their review that the continued interest in this physical phenomenon is caused by its implication in numerous biomedical effects, such as haemolysis, sonoporation, targeted drug delivery and bone healing, as well as in microfluidic technologies, such as micromixing and ultrasonic cleaning.

The case of two bubbles is the next natural step after the case of a single bubble in the theoretical modelling of acoustic microstreaming. The study of this case allows one to understand how the acoustic interaction between bubbles affects the behaviour of the acoustic microstreaming. This case also paves the way to more practical situations where the streaming flow is induced by a bubble cluster with a large number of cavitation bubbles.

Experimental data on acoustic microstreaming induced by two interacting bubbles are still scarce. The main difficulties in studying a multiple-bubble system are its three-dimensional nature and the difficulty of precisely controlling bubble locations in space. For these reasons, experimental studies mainly focus on substrate-attached bubbles (Tho et al. 2007; Bolanos-Jimenez et al. 2017), pancake-like bubbles constrained in microchannels (Mekki-Berrada et al. 2016) or arrays of armoured bubbles printed on the base of a microfluidic channel (Bertin et al. 2017). In the presence of a boundary, secondary flows appear in the vicinity of the bubble interface, which considerably complicate the analysis of the flow trajectories in three dimensions. Therefore, experimental studies usually consider the case of a streaming flow resulting from purely volume oscillations of bubbles, the interaction between volume and translational oscillations, or the interaction with an ellipsoidal mode in the case of armoured bubbles. Investigation of streaming interactions in the case of axisymmetric shape modes for free, being far from boundaries, bubbles would require a certain type of optical or acoustical trapping (Garbin et al. 2007). Such a trapping has recently been performed in the case
of two interacting bubbles that underwent axisymmetric shape oscillations (Regnault *et al.* 2020),
but no associated streaming pattern has been shown.

A theoretical study on acoustic microstreaming generated by two bubbles has been
carried out by Wang & Chen (2013). However, their analytical model is based on considerable
constraints. It assumes that the bubbles only undergo radial and translational oscillations. Shape
modes are not considered. That is not mentioned explicitly but their solutions for the linear
velocity field are only accurate up to leading terms in the inverse interbubble distance, which
means that the distance between the equilibrium centres of the bubbles should be large compared
to the equilibrium bubble radii. The acoustic interaction between the bubbles is only considered
in the linear approximation. This means that, when the linear velocity field is calculated,
boundary conditions at the bubble surfaces allow for the effect of the scattered wave from one
bubble on the oscillation of the other bubble. However, when the time-averaged liquid velocity
(velocity of acoustic streaming) is calculated, boundary conditions at the bubble surfaces ignore
the effect of second-order streaming flows produced by the neighbouring bubble. Furthermore,
their derivation assumes that the viscous boundary layer is small compared to the bubble radii,
i.e. viscous effects are considered essential only within a thin boundary layer next to the bubble
interfaces and the liquid motion is assumed inviscid beyond this layer. A more accurate
theoretical study on acoustic microstreaming induced by an interacting bubble pair has been
performed by Doinikov & Bouakaz (2016). They did not impose restrictions on the thickness of
the viscous boundary layer but all the other restrictions mentioned above are kept in their model,
including the assumption that the bubbles only undergo radial and translational oscillations.

In the present study, we propose an analytical theory that describes acoustic
microstreaming in the case of two interacting bubbles undergoing oscillation modes of arbitrary
order, which means radial oscillations, translation and nonspherical oscillations known as shape
modes. Our derivation provides analytical solutions in terms of complex amplitudes of
oscillation modes, which means that the modal amplitudes are assumed to be known (for
example, obtained from experimental microbubble dynamics) and serve as input data when the
velocity field of acoustic microstreaming is calculated. It should be mentioned that, if necessary,
the modal amplitudes can be evaluated theoretically, for example, in terms of the amplitude of
the driving acoustic pressure (Plesset 1954; Francescutto & Nabergoj 1978; Hall & Seminara
1980; Shaw 2006; Guédra & Inserra 2018). No restrictions are imposed on the ratio of the bubble
radii to the viscous penetration depth and the distance between the bubbles. Therefore, the obtained analytical solutions are valid for any values of liquid viscosity and any separation distances between the bubbles. The interaction between the bubbles is considered both when the linear velocity field is calculated and when the second-order velocity field of acoustic microstreaming is calculated.

The paper is organised as follows. In §2, the core of the theoretical derivation is described. In §3, capabilities of the analytical theory are illustrated by computational examples. In Appendices A to E, details of the analytical derivation are provided in order to make the understanding of mathematical aspects easier for the reader.

## 2. Theory

We consider two gas bubbles that are immersed in a viscous incompressible liquid and undergo axisymmetric oscillations. The distance between the equilibrium centres of the bubbles is denoted by $d$. We use two spherical coordinate systems in our calculations, $(r_1, \theta_1, \epsilon_1)$ and $(r_2, \theta_2, \epsilon_2)$, which are originated at the equilibrium centres of bubbles 1 and 2, respectively, and the direction $\theta_1 = \theta_2 = 0$ corresponds to the $z$ axis; see figure 1.

<figure><img src="image_2.png" /><figcaption>Figure 1. Coordinate systems used in calculations.</figcaption></figure>

### 2.1. First-order scattered field
Before we will proceed to calculations, a comment should be made. In the case of a single bubble undergoing shape modes, a scattered wave generated by a certain mode has the same angular dependence as the generating mode. For example, a shape mode depending on $P_n(\cos \cos \theta)$, where $P_n$ is the Legendre polynomial of degree $n$, generates a scattered wave depending on $P_n(\cos \cos \theta)$. In the case of two bubbles, the structure of the total scattered wave is different. Due to multiple rescattering between the bubbles, every oscillation mode generates a scattered wave with $P_m(\cos \cos \theta)$ of all degrees, i.e. $0 \le m < \infty$. This effect can be seen in solutions obtained by Doinikov & Bouakaz (2015). Although Doinikov & Bouakaz (2015) assume that two bubbles only undergo pulsation and translation, the scattered wave includes terms with $P_m(\cos \cos \theta)$ of all degrees and the amplitudes of these terms are related to one another. This fact is taken into account in the derivation below.

We assume that each bubble undergoes shape modes, which can have different orders and different frequencies because some of them are excited parametrically. The modes can be divided into groups in accordance with their frequency. These groups can be considered separately because the interaction of modes having the same frequency can only contribute to acoustic microstreaming (Doinikov et al. 2019a,b; Inserra et al. 2020a,b). We consider modes with a frequency $\omega$, which can be equal to the driving frequency or any other sub- or ultra-harmonic of the driving frequency if the modes are excited parametrically.

We assume that the $j$th bubble undergoes $N_j$ modes with a frequency $\omega$. The modes have numbers $M_1^{(j)}, M_2^{(j)}, \dots, M_{N_j}^{(j)}$. In this case, the $\omega$-dependent deformation of the surface of the $j$th bubble can be represented as

$$r_s^{(j)} = e^{-i\omega t} \sum_{n=M_1^{(j)}}^{M_{N_j}^{(j)}} s_n^{(j)} P_n(\mu_j), \qquad (2.1)$$

where $\mu_j = \cos \cos \theta_j$ and $s_n^{(j)}$ is the complex amplitude of the $n$th mode of the $j$th bubble, which is assumed to be known (for example, obtained from experimental bubble dynamics).

The total first-order liquid velocity generated by both bubbles is given by

$$v = v^{(1)} + v^{(2)}, \qquad (2.2)$$
where $v^{(j)}$ is the first-order liquid velocity generated by the jth bubble.

$v^{(j)}$ can be represented as

$$v^{(j)} = \nabla\varphi^{(j)} + \nabla \times \psi^{(j)}, \qquad (2.3)$$

where $\varphi^{(j)}$ and $\psi^{(j)}$ are the scalar and the vector velocity potentials, respectively.

In view of axial symmetry, $\varphi^{(j)}$ and $\psi^{(j)}$ are written as (Doinikov et al. 2019a)

$$\varphi^{(j)} = e^{-i\omega t} \sum_{n=0}^{\infty} a_n^{(j)} \left(\frac{R_{j0}}{r_j}\right)^{n+1} P_n(\mu_j), \qquad (2.4)$$

$$\psi^{(j)} = \psi^{(j)}(r_j, \theta_j, t)e_{\epsilon j} = e^{-i\omega t} e_{\epsilon j} \sum_{n=1}^{\infty} b_n^{(j)} h_n^{(1)}(k_v r_j) P_n^1(\mu_j), \qquad (2.5)$$

where $R_{j0}$ is the equilibrium radius of the jth bubble, $P_n^1$ is the associated Legendre polynomial of the first order and degree $n$, $h_n^{(1)}$ is the spherical Hankel function of the first kind, $k_v = \frac{(1+i)}{\delta}$ is the viscous wavenumber, $\delta = \sqrt{\frac{2\nu}{\omega}}$ is the viscous penetration depth, $\nu = \frac{\eta}{\rho}$ is the kinematic liquid viscosity, $\eta$ is the dynamic liquid viscosity, $\rho$ is the equilibrium liquid density and $e_{\epsilon j}$ is the unit azimuth vector of the jth bubble. Note that axial symmetry allows us to set $\epsilon_1 = \epsilon_2$ and $e_{\epsilon 1} = e_{\epsilon 2}$.

Substitution of (2.4) and (2.5) into (2.3) yields

$$v_r^{(j)} = - \frac{e^{-i\omega t}}{r_j} \sum_{n=0}^{\infty} (n+1) \left[ a_n^{(j)} \left( \frac{R_{j0}}{r_j} \right)^{n+1} + n b_n^{(j)} h_n^{(1)}(k_v r_j) \right] P_n(\mu_j), \qquad (2.6)$$

$$v_{\theta}^{(j)} = \frac{e^{-i\omega t}}{r_j} \sum_{n=1}^{\infty} \left\{ a_n^{(j)} \left( \frac{R_{j0}}{r_j} \right)^{n+1} - b_n^{(j)} \left[ h_n^{(1)}(k_v r_j) + k_v r_j h_n^{(1)}(k_v r_j) \right] \right\} P_n^1(\mu_j), \qquad (2.7)$$

where the prime denotes the derivative with respect to the argument in brackets. Note that when calculating (2.6) and (2.7), we have used mathematical properties of $P_n(\mu)$ and $P_n^1(\mu)$ (Abramowitz & Stegun 1972).

To find $a_n^{(j)}$ and $b_n^{(j)}$, called linear scattering coefficients, boundary conditions at the bubble surfaces are applied. Since this calculation is cumbersome, it is performed in Appendix A.
As an example, figure 2 shows the magnitudes of $a_n^{(j)}$ and $b_n^{(j)}$ calculated for bubbles with $R_{10} = R_{20} = 50$ µm at the interbubble distance $d = 300$ µm and the oscillation frequency $f = 30$ kHz. The surrounding liquid is water. The bubbles are assumed to undergo the radial oscillation with an amplitude of 5 µm. As one can see, the scattering coefficients with $n > 0$ are not zero despite the fact that both bubbles only pulsate, but the magnitude of the scattering coefficients quickly decreases with increasing $n$. It is worth noticing that the coefficients $a_n^{(1)}$ are six orders of magnitudes larger than the coefficients $b_n^{(1)}$. This difference is explained by the weight of the radial functions $(\frac{R_{j0}}{r_j})^{n+1}$ and Hankel function $h_n^{(1)}(k_v r_j)$ that are respectively associated to the coefficients $a_n^{(1)}$ and $b_n^{(1)}$ in (2.6) and (2.7).

<figure><img src="image_3.png" /><figcaption>Figure 2. Magnitude of the linear scattering coefficients for bubble 1. Due to the equivalence of the calculation parameters, the linear scattering coefficients for bubble 2 have the same magnitude.</figcaption></figure>

## 2.2. Acoustic streaming

As shown by Doinikov et al. (2019a), the Eulerian velocity of acoustic streaming can be written as

$$v_E = \nabla \times \Psi, \qquad (2.8)$$
where $\Psi$ obeys the vorticity equation derived by Westervelt (1953) for the case of solenoidal first-order motion,

$$ \Delta^2 \Psi = -\frac{1}{\nu} \nabla \times \langle \mathbf{v} \cdot \nabla \mathbf{v} \rangle, \qquad (2.9) $$

with $\langle \rangle$ denoting the time average.

Equation (2.9) is shown by Doinikov et al. (2019a) to be transformed to

$$ \Delta^2 \Psi = \frac{1}{2\nu} \mathrm{Re} \{ \nabla \times [(\Delta \Psi^*) \times \mathbf{v}] \}, \qquad (2.10) $$

where Re means “the real part of” and the asterisk denotes the complex conjugate. Substituting $\mathbf{v} = \mathbf{v}^{(1)} + \mathbf{v}^{(2)}$ and $\Psi = \Psi^{(1)} + \Psi^{(2)}$, one obtains

$$ \Delta^2 \Psi = \frac{1}{2\nu} \mathrm{Re} \{ \nabla \times [\Delta(\Psi^{(1)*}) \times \mathbf{v}^{(1)}] + \nabla \times [\Delta(\Psi^{(2)*}) \times \mathbf{v}^{(2)}] \\ + \nabla \times [\Delta(\Psi^{(1)*}) \times \mathbf{v}^{(2)}] + \nabla \times [\Delta(\Psi^{(2)*}) \times \mathbf{v}^{(1)}] \}. \qquad (2.11) $$

Equation (2.11) can be represented in the coordinates $(r_1, \theta_1)$ or $(r_2, \theta_2)$. It is shown in Appendix B that $\Psi$ in the coordinates $(r_1, \theta_1)$ is written as

$$ \Psi(r_1, \theta_1) = \Psi(r_1, \theta_1) e_{\epsilon_1}, \qquad (2.12) $$

where $\Psi(r_1, \theta_1)$ obeys the following equation:

$$ \left( \Delta_{r\theta}^{(1)} - \frac{1}{r_1^2 \sin^2 \theta_1} \right)^2 \Psi(r_1, \theta_1) = \frac{1}{2\nu R_{10}^4} \left\{ \begin{aligned} & \sum_{n=0}^{\infty} P_n(\mu_1) P_m^1(\mu_1) [K_{1nm}^{(1)}(r_1) + L_{1nm}^{(1)}(r_1) + M_{1nm}^{(1)}(r_1)] \\ & + \sqrt{1 - \mu_1^2} \sum_{n,m=1}^{\infty} P_n^1(\mu_1) P_m^{1*}(\mu_1) [K_{2nm}^{(1)}(r_1) + K_{2mn}^{(1)}(r_1) + L_{2nm}^{(1)}(r_1) + L_{2mn}^{(1)}(r_1) + M_{2nm}^{(1)}(r_1) + M_{2mn}^{(1)}(r_1)] \end{aligned} \right\}. \qquad (2.13) $$
Here, the operator $\Delta_{r\theta}^{(1)}$ is defined by (B5) and the $r_1$-dependent functions are calculated by (B9), (B10), (B13), (B14), (B20) and (B21).

Equation (2.13) is valid for $r_1 < d$ since this restriction is imposed on coordinate transformations (Appendix A) that were used to derive (2.13). Therefore, the solution of (2.13) describes acoustic streaming in the domain $D_1$ shown in figure 3.

As shown in Appendix B, $\Psi$ in the coordinates $(r_2, \theta_2)$ is written by

$$ \Psi(r_2, \theta_2) = \Psi(r_2, \theta_2)e_{\varepsilon_2}, \qquad (2.14) $$

where $\Psi(r_2, \theta_2)$ is calculated by

$$ \left( \Delta_{r\theta}^{(2)} - \frac{1}{r_2^2 \sin^2 \theta_2} \right)^2 \Psi(r_2, \theta_2) = \frac{1}{2\nu R_{20}^4} \left\{ \begin{aligned} & \sum_{n=0}^{\infty} P_n(\mu_2) P_m^1(\mu_2) [K_{1nm}^{(2)}(r_2) + L_{1nm}^{(2)}(r_2) + M_{1nm}^{(2)}(r_2)] \\ & + \sqrt{1 - \mu_2^2} \sum_{n,m=1}^{\infty} P_n^1(\mu_2) P_m^{1/2}(\mu_2) [K_{2nm}^{(2)}(r_2) + K_{2mn}^{(2)}(r_2) + L_{2nm}^{(2)}(r_2) + L_{2mn}^{(2)}(r_2) + M_{2nm}^{(2)}(r_2) + M_{2mn}^{(2)}(r_2)] \end{aligned} \right\}. \qquad (2.15) $$

Here, the operator $\Delta_{r\theta}^{(2)}$ is defined by (B5) and the $r_2$-dependent functions are calculated by (B26), (B27), (B30), (B31), (B34) and (B35).

Equation (2.15) is valid for $r_2 < d$ and hence its solution describes acoustic streaming in the domain $D_2$ shown in figure 3. Equations (2.13) and (2.15) must give the same results in the space between the bubbles so the domains $D_1$ and $D_2$ can be combined to get acoustic streaming around both bubbles as shown in figure 3.

It should be emphasised that the above restrictions on the solutions of (2.13) and (2.15) are not related to the ratio of the bubble radii to the viscous penetration depth or to the inter-bubble distance. The restrictions shown in figure 3 result from geometrical transformations of coordinates that are used in our derivation and hence are of purely geometrical nature. They are necessary to provide the convergence of infinite sums used in the solutions.
<figure><img src="image_4.png" /><figcaption>Figure 3. Regions where (2.13) and (2.15) are valid.</figcaption></figure>

Equations (2.13) and (2.15) are solved in Appendix C. The result is given by

$$ \Psi(r_j, \theta_j) = \sum_{l=1}^{\infty} \Psi_l^{(j)}(r_j) P_l^1(\mu_j), \quad j=1,2, \qquad (2.16) $$

where $\Psi_l^{(j)}(r_j)$ is calculated by (C8). Substituting (2.16) into (2.8), one obtains the components of the Eulerian streaming velocity,

$$ v_{Er}(r_j, \theta_j) = -\frac{1}{r_j} \sum_{l=1}^{\infty} l(l+1) \Psi_l^{(j)}(r_j) P_l(\mu_j), \qquad (2.17) $$

$$ v_{E\theta}(r_j, \theta_j) = -\frac{1}{r_j} \sum_{l=1}^{\infty} \left[ \Psi_l^{(j)}(r_j) + r_j \Psi_l^{(j)}(r_j) \right] P_l^1(\mu_j), \qquad (2.18) $$

where $\Psi_l^{(j)}(r_j)$ is calculated by (C17).

The components of the Lagrangian streaming velocity are defined by
$$v_{Lr}(r_j, \theta_j) = v_{Er}(r_j, \theta_j) + v_{Sr}(r_j, \theta_j), \qquad (2.19)$$

$$v_{L\theta}(r_j, \theta_j) = v_{E\theta}(r_j, \theta_j) + v_{S\theta}(r_j, \theta_j), \qquad (2.20)$$

where $v_{Sr}(r_j, \theta_j)$ and $v_{S\theta}(r_j, \theta_j)$ denote the components of the Stokes drift velocity (Longuet-Higgins 1998), which are calculated in Appendix E and given by (E15) and (E16).

Combining the solutions in the coordinates $(r_1, \theta_1)$ and $(r_2, \theta_2)$ provides acoustic microstreaming in the area shown in figure 3. As a consequence, when the distance between the bubbles is reduced, the boundary of the spatial domain where our second-order solutions are valid shifts to the bubble interfaces. It is also worth mentioning that the domain $D = \{r_1 < d \cup r_2 < d\}$ is a theoretical upper limit for the spatial domain where (2.19) and (2.20) are valid. Since in the process of numerical calculations one has to truncate the infinite sums in (2.17) and (2.18), one can only approach this upper limit by increasing the number of terms in the truncated sums, which naturally leads to more time-consuming computations.

### 3. Numerical simulations

Preparatory to the numerical analysis of specific cases of microstreaming induced by an interacting bubble pair, attention is paid to the convergence of the numerical modelling.

#### 3.1. Convergence analysis

In the present modelling, the interaction between bubbles is considered both when the linear velocity field is calculated and when the second-order velocity field of acoustic microstreaming is calculated. Physically, the interaction is realised by re-scattering of sound between the bubbles. One can see that in the course of the calculation of the linear velocity, multiple scattering effects are included into the derivation of the linear scattering coefficients $a_n^{(j)}$ and $b_n^{(j)}$, which are given by the system of equations (A41) and (A42) for $j = 1$ and $n \ge 1$. The infinite summation appearing in (A41) and (A42) must be truncated at some finite value of a given number $n$ of the linear scattering coefficients, while the coefficients appearing in this system of equations, which are given by equations (A43) – (A48) also involve infinite sums,
which determine the number of re-scattering events that are allowed for. These sums must also be truncated to some finite value $m_{lim}$. The numerical case under consideration is the interaction between two bubbles of identical size $R_{10} = R_{20} = 50\mu m$ separated by the inter-bubble distance $d = 150\mu m$ that oscillate radially in phase, with the amplitude 5 $\mu m$ at the frequency 30 kHz. In order to estimate the number of the linear scattering coefficients required to accurately describe multiple scattering effects, figure 4(a) shows the evolution of the total energy $\sum_{i=0}^{n} \|a_i^{(1)}\|^2$ carried by the linear scattering coefficients $a_i^{(1)}$ up to the value $n$. The energy reaches a plateau at $n = 5$, providing the limiting number of first-order interactions that should be taken into account. The linear scattering coefficient of the highest amplitude $a_1^{(1)}$, in the case $n \neq 0$, is plotted in figure 4(b) as a function of the limit value $m_{lim}$. Setting $m_{lim} = 5$ in our case is enough to ensure the convergence of the linear scattering coefficients.

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

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

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

<figure><img src="image_8.png" /></figure>
Figure 4. Convergence analysis of the numerical modelling for the case of two spherically-oscillating bubbles of identical radii ($R_{10} = R_{20} = 50\mu m$) driven by a 30 kHz ultrasound wave in water. (a) Convergence of the energy spectrum of the linear scattering coefficients as a function of the number $n$ of scattered modes. (b) Convergence of the highest scattering coefficient $a_1^{(1)}$ as a function of the finite upper limit $m_{lim}$ of the sums appearing in (A41) - (A48). (c),(d) Convergence of the radial $V_r$ and tangential $V_\theta$ components of the Lagrangian streaming velocity as a function of the finite upper limit $l_{lim}$ of the sums appearing in (2.17) and (2.18).

In the second-order streaming velocities, the velocity field is also defined by an infinite sum of the angular functions $P_n(\mu)$ and $P_n^1(\mu)$, such as in the system of equations (2.17) and (2.18) for the components of the Eulerian streaming velocity. The numerical truncation of these infinite sums to a finite value $l_{lim}$ is quantified on the radial and tangential components of the Lagrangian streaming velocity in figure 4(c,d). The velocity components are computed at the location ($r_1 = R_{01} + 5\mu m, \theta_1 = \frac{\pi}{2}$) in the reference frame of bubble 1. The velocity components have converged for $l_{lim} \ge 6$. In the following, numerical results are only considered for which the convergence has been achieved.

## 3.2. *Effect of the inter-bubble distance*

A single bubble that pulsates radially in an unbounded fluid does not produce any microstreaming (Doinikov *et al.* 2019a). Due to the presence of a second bubble, multiple sound scattering leads to a non-zero microstreaming flow. Figure 5(a,b,c) shows streamline patterns, for different inter-bubble distances, around two bubbles with identical equilibrium radii $R_{10} = R_{20} = 50\mu m$, oscillating radially in phase with the amplitude 5µm at the frequency 30 kHz. The bubbles are indicated by full circles and the flow pattern is only displayed within the region of validity $D = \{r_1 < d \cup r_2 < d\}$.

The streamline patterns consist of large lobes located above and below each bubble. The form and the size of the lobes are changed with increasing inter-bubble distance, namely the
centre of rotation is shifted away from the bubbles and the extent of the vortices increases. Along the horizontal direction $\theta_1 = 0$, the streaming flow is directed towards the bubbles in the inter-bubble space. Two pairs of recirculation vortices are observed in the distance from the bubbles along the horizontal axis. The size of these vortices decreases when the inter-bubble distance increases. Figure 5(d) depicts the evolution of the magnitude of the Lagrangian streaming velocity as a function of the inter-bubble distance at the locations ($r_1 = 1.3R_{01}$, $\theta_1 = \frac{\pi}{2}$) and ($r_2 = 1.5R_{01}$, $\theta_2 = \pi$). The velocity magnitude rapidly decreases to values of the order of 10 µm/s. Such small velocities can not be observed experimentally. It is worth noticing that these results are in good qualitative agreement with the experimental observations of Mekki-Berrada et al. (2016). Though their study is devoted to the analysis of streaming flows induced by pancake-like bubbles, the same kind of long-range streaming patterns is observed around the interacting bubble pair.

<figure><img src="image_9.png" /><figcaption>Figure 5. Evolution of the streaming flow as a function of the inter-bubble distance $d$ in the case of two spherically oscillating bubbles (with the same parameters as in figure 4). (a) $d = 150$µm. (b) $d = 300$µm. (c) $d = 750$µm. (d) Evolution of the magnitude of the Lagrangian streaming</figcaption></figure>
velocity as a function of the inter-bubble distance. The velocity is measured at the locations ($r_1 = 1.3R_{01}, \theta_1 = \frac{\pi}{2}$) and ($r_2 = 1.5R_{01}, \theta_2 = \pi$).

### 3.3. Contribution of the present model

In order to illustrate the contribution of the present model, we perform a comparison between the microstreaming pattern that is obtained by the superimposition of two streaming flows calculated by the single bubble theory and the pattern predicated by the present model. The microstreaming pattern induced by each bubble in the non-interacting case is calculated by the theoretical approach of Inserra et al. (2020b). Numerical results are shown in figure 6 for two identical bubbles with $R_{10} = R_{20} = 50\mu m$ at the inter-bubble distance $d = 200\mu m$, exhibiting radial ($s_0^{(1,2)} = 5\mu m$) and mode 2 ($s_2^{(1,2)} = 10\mu m$) oscillations at the driving frequency 30 kHz in water. As one can see, the microstreaming patterns in figures 6(a) and 6(b) look markedly different.

<figure><img src="image_10.png" /><figcaption>Figure 6. Comparison of (a) the microstreaming pattern obtained by the superimposition of two flows given by the single bubble theory and (b) the pattern predicted by the present model.</figcaption></figure>

In order to reinforce the validity of these results, a computational fluid dynamics simulation has been performed using StarCCM+ software. The CFD simulation is based on the interface motion of a deformable spheroid that sets in motion the fluid in its surroundings. The interface motion is described by (2.1) for each bubble. The interface deformation is achieved by the morphing module implemented in StarCCM+. This method was first validated successfully
when investigating the microstreaming flow induced by a single axisymmetric bubble, whose
theoretical modelling is known (Doinikov *et al.* 2019a,b; Inserra *et al.* 2020a,b). In the present
study, the numerical modelling has been extended to two bubbles undergoing axisymmetric
shape oscillations. The computational domain consists of a quarter of a millimetre sphere (liquid
medium) containing two quarters of sub-millimetre spheroids (i.e. two bubbles), see Appendix F.
The axes of each quarter of the sphere are aligned together. The Navier-Stokes and continuity
equations are resolved through a laminar, incompressible and unsteady model using a pressure
boundary applied at the far field and a slip condition at the interfaces. The smallest mesh size at
the interface is 1 µm. To compare the fluid field given by the CFD simulation to that given by
our analytical model, the velocity field is averaged over one period. Figure 7 displays the
microstreaming pattern obtained by the CFD simulation. The microstreaming pattern is very
similar to the analytic one shown in figure 6(b). The location and the size of the lobes are well
captured and the flow direction is identical to that obtained analytically. Figure 7(b) shows a
partial 3D view of the streamlines induced by the bubble oscillations. This figure highlights the
generation of vortex rings around the bubbles, which correspond to the lobes shown in the
different cross-view.

It is worth noting that the computational time of the CFD simulation is of the order of 24
hours, while the calculation of streamlines by the analytical model takes about 2 hours.

<figure><img src="image_11.png" /><figcaption>Figure 7. Microstreaming flow obtained by CFD simulation. (a) Cross-section view of the flow. (b) 3D view of the streamlines illustrating the vortex ring generation around the bubbles.</figcaption></figure>
We have also made a comparison with the theory developed by Doinikov & Bouakaz (2016). Recall that this theory, first, assumes the distance between the bubbles to be large compared to the bubble radii and, second, ignores the interaction between the bubbles in second-order solutions. Figure 8(a) demonstrates a streamline pattern predicted by this theory, while figure 8(b) shows a pattern given by the present model. The simulations were made for bubbles with $R_{10} = R_{20} = 50$ µm at $d = 300$ µm, assuming that both bubbles undergo radial oscillations with the amplitude 5 µm at the frequency 30 kHz in water. As one can see, figure 8(a) predicts an opposite direction of the streaming flow, which is likely to be a consequence of the fact that the theory of Doinikov & Bouakaz (2016) only allows for the interaction of the bubbles in the first-order solutions.

<figure><img src="image_12.png" /><figcaption>Figure 8. Comparison of (a) the streamline pattern given by the theory of Doinikov & Bouakaz (2016) and (b) the pattern predicted by the present model.</figcaption></figure>

### 3.4. Effect of viscosity

Figures 9 – 12 illustrate the effect of the liquid viscosity on the behaviour of the microstreaming. The simulations were made for bubbles with $R_{10} = R_{20} = 50$ µm at the frequency 30 kHz, three values of the inter-bubble distance $d$ and three values of the viscosity $\eta = 0.001, 0.01$ and $0.1$ Pa s.

Figure 9 shows streamline patterns in the case that both bubbles undergo radial oscillations with the amplitude 5 µm. As one can see, if $d$ is not large, with increasing $\eta$, the
form of vortices is changed noticeably. The vortices become increasingly detached from the bubble surfaces and their extent increases. Figure 10 shows how the magnitude of the Lagrangian streaming velocity $|v_L|$ depends on the liquid viscosity. $|v_L|$ is depicted as a function of the distance from the bubble surface at two values of the inter-bubble distance along two directions. Direction 1 is counted from the surface of bubble 2 at $\theta_2 = 0$, while direction 2 is counted from the surface of bubble 1 at $\theta_1 = -\pi/2$. As one can see, increasing $\eta$ leads to increasing $|v_L|$ and the peak value of $|v_L|$ is shifted away from the bubble surface. It is interesting to note that along direction 2, the increase of $|v_L|$ at $\eta = 0.01$ Pa s is higher than at $\eta = 0.1$ but the region where $|v_L|$ is relatively high remains wider at $\eta = 0.1$.

<figure><img src="image_13.png" /></figure>
Figure 9. Streamline patterns produced by two radially oscillating bubbles of equal radii at different inter-bubble distances and different values of the liquid viscosity. (a),(d),(g) $\eta = 0.001$ Pa s; (b),(e),(h) $\eta = 0.01$ Pa s; (c),(f),(i) $\eta = 0.1$ Pa s.

Figure 11 shows streamline patterns in the case that both bubbles undergo radial and mode 2 oscillations. The amplitude of both modes is 5 µm and the phase shift between modes 0 and 2 is zero. As one can see, with increasing $\eta$, the vortices become increasingly detached from the bubble surfaces. At $d = 600$ µm and $\eta = 0.1$ Pa s, a complicated structure of streamlines is observed in the inter-bubble space. The dependence of the magnitude of the Lagrangian streaming velocity on the liquid viscosity is illustrated by figure 12, which is analogous to figure 10. As one can see, along direction 1, $|v_L|$ at $\eta = 0.01$ Pa s is higher than at $\eta = 0.001$ Pa s but further increasing the viscosity to 0.1 Pa s results in $|v_L|$ that is lower than at $\eta = 0.001$ Pa s. Other effects are similar to those in figure 10.
<figure><img src="image_14.png" /><figcaption>Figure 10. The magnitude of the Lagrangian streaming velocity $|v_L|$, which is generated by two radially oscillating identical bubbles, as a function of the distance from the bubble surface along two directions. Direction 1 is counted from the surface of bubble 2 at $\theta_2 = 0$, while direction 2 is counted from the surface of bubble 1 at $\theta_1 = -\pi/2$. The other parameters are as in figure 9.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 11. Streamline patterns produced by bubbles of equal radii, undergoing radial and mode 2 oscillations, at different inter-bubble distances and different values of the liquid viscosity.</figcaption></figure>

(a),(d),(g) $\eta = 0.001$ Pa s; (b),(e),(h) $\eta = 0.01$ Pa s; (c),(f),(i) $\eta = 0.1$ Pa s.
<figure><img src="image_16.png" /><figcaption>Figure 12. The magnitude of the Lagrangian streaming velocity $|v_L|$, which is generated by two identical bubbles undergoing radial and mode 2 oscillations, as a function of the distance from the bubble surface along two directions. Direction 1 is counted from the surface of bubble 2 at $\theta_2 = 0$, while direction 2 is counted from the surface of bubble 1 at $\theta_1 = -\pi/2$. The other parameters are as in figure 11.</figcaption></figure>

### 3.5. *Comparison with existing experimental observations*

Streaming flows induced in the vicinity of an interacting bubble pair have been reported by Mekki-Berrada *et al.* (2016) for flattened bubbles undergoing nonspherical oscillations in a microchannel. They consider the case of two bubbles of radii $R_{10} = R_{20} = 38\mu m$, separated by $d = 200\mu m$, undergoing spherical, translational and mode 2 oscillations. The associated mode dynamics, meaning the values of the complex amplitudes $s_{0,1,2}^{(1,2)}$, is provided. The resulting streaming pattern, which is shown in figure 13(a), consists of two anti-fountain vortices on both
sides of the bubble pair and two pairs of smaller vortices that are located between the bubbles. Along the line joining the bubble centres, the streaming flow is directed towards the bubbles. Experimental streaming velocities are found to be of the order of mm/s. The measured input data for $s_{0,1,2}^{(1,2)}$ provided by Mekki-Berrada *et al.* (2016) allow one to perform a qualitative comparison with our model. Although our theory implies that bubbles are initially of spherical shape and located in an unbounded fluid, one can expect that their dynamics should exhibit similar qualitative features as in the case investigated by Mekki-Berrada *et al.* (2016). The theoretical streaming pattern is shown in figure 13(b) and is in good agreement with the experimental one: both the large lobes and the small vortices between the bubbles are observed. Also, our theory provides velocity magnitudes of the same order as in the experimental results.

<figure><img src="image_17.png" /><figcaption>Figure 13. Comparison between (a) experimental observations of the steaming flow induced by a bubble pair undergoing radial, translational and mode 2 oscillations (reproduced from Mekki-Berrada *et al.* (2016), see the copyright form) and (b) the theoretical results given by the present model.</figcaption></figure>

4. Conclusion
Analytical equations have been derived that describe acoustic microstreaming induced by two interacting bubbles undergoing arbitrary axisymmetric shape oscillations. The developed theory imposes no restrictions on the bubble sizes and the liquid viscosity. Numerical simulations have been performed to illustrate the influence of the interaction of the bubbles on the microstreaming. In particular, it was shown that interacting bubbles undergoing only radial oscillations could generate microstreaming in contrast to a single radially oscillating bubble. Examples of streamlines for bubbles undergoing shape modes were presented and successfully compared to computational fluid dynamics simulations and to experimental results available in the literature.

It should be mentioned that the approach applied in our study can be generalised to the case of a bubble cluster. Theoretically, that is possible but will lead to very complicated and cumbersome calculations. Another possibility is to apply the so-called two-particle approximation, according to which, under certain conditions, pairwise interactions, calculated by a two-particle theory, can be used to model the behavior of a multi-particle system.

Another point to be mentioned is the effect of Bjerknes forces experienced by bubbles. It should be borne in mind in this connection that this effect is a nonlinear effect of the second order, just as acoustic microstreaming. Therefore, theoretically, they are considered independent, which means that within the framework of the second-order approximation, Bjerknes forces and acoustic microstreaming are assumed not to affect each other.

## Acknowledgements

This work was supported by the LabEx CeLyA of the University of Lyon (Grants No. ANR-10-LABX-0060 and No. ANR-11-IDEX-0007). A. A. D. gratefully acknowledges support from Collegium de Lyon.

## Declaration of interest

The authors report no conflict of interest.

## Appendix A. Calculation of linear scattering coefficients
The boundary condition for the liquid velocity requires that the normal component of $\mathbf{v}$ at $r_j = R_{j0}$ be equal to the normal component of the surface velocity of the jth bubble,

$$v_r|_{r_j=R_{j0}} = \frac{dr_s^{(j)}}{dt} = -i\omega e^{-i\omega t} \sum_{n=M_{N_j}^{(j)}}^{M_{N_j}^{(j)}} s_n^{(j)} P_n(\mu_j) \quad (A1)$$

The second boundary condition requires that the tangential stress vanish at the bubble surfaces because the gas viscosity is much lower than the liquid viscosity,

$$\sigma_{r\theta} = \eta \left( \frac{1}{r_j} \frac{\partial v_r}{\partial \theta_j} + \frac{\partial v_{\theta}}{\partial r_j} - \frac{v_{\theta}}{r_j} \right) = 0 \quad \text{at} \quad r_j = R_{j0} \quad (A2)$$

It is worthwhile to explain the origin of (A1) and (A2). These equations represent a generally accepted approximation that allows one not to consider the motion of gas inside the bubbles and thus to simplify calculations. The approximation is based on the fact that the viscosity of the gas is small compared to the viscosity of the ambient liquid. This fact allows one to assume the slippage of the liquid over the bubble surface, which means that the boundary condition for the tangential liquid velocity can be omitted so that the boundary condition for the normal liquid velocity, (A1), is only applied. Next, the exact boundary condition for tangential stress implies the equality of the tangential stresses, $\sigma_{r\theta}$, on the outer and inner sides of the bubble surface. However, since $\sigma_{r\theta}$ is proportional to viscosity and the gas viscosity is much smaller than the liquid viscosity, the tangential stress in the gas can be considered as negligible (zero) with respect to the tangential stress in the liquid, which results in (A2). Finally, the boundary condition for normal stress is used to calculate the amplitudes of oscillation modes in terms of the amplitude of an imposed acoustic pressure. However, we assume that we know the mode amplitudes because we measure them experimentally. Therefore, the boundary condition for normal stress is not necessary in our derivation.

In order to satisfy (A1) and (A2), we need to express $\varphi^{(2)}$ and $\psi^{(2)}$ in the coordinates of bubble 1 and vice versa.
For $\varphi^{(1)}$ and $\varphi^{(2)}$, the coordinate transformation can be done by the following identities
(Varshalovich et al. 1988; Doinikov & Bouakaz 2015):

$$ \frac{P_n(\mu_1)}{r_1^{n+1}} = \frac{1}{d^{n+1}} \sum_{m=0}^{\infty} (-1)^m C_{nm} \left(\frac{r_2}{d}\right)^m P_m(\mu_2), \qquad (A3) $$

$$ \frac{P_n(\mu_2)}{r_2^{n+1}} = \frac{(-1)^n}{d^{n+1}} \sum_{m=0}^{\infty} C_{nm} \left(\frac{r_1}{d}\right)^m P_m(\mu_1), \qquad (A4) $$

where $C_{nm} = (n+m)!/(n!m!)$ and $d$ is the distance between the centers of the bubbles.
Equations (A3) and (A4) are valid for $r_{1,2} < d$.

Substituting (A3) into (2.4) at $j=1$ and (A4) into (2.4) at $j=2$, one obtains

$$ \varphi^{(1)}(r_2, \theta_2) = e^{-i\omega t} \sum_{n,m=0}^{\infty} a_n^{(1)} \xi_1^{n+1} (-1)^m C_{nm} \left(\frac{r_2}{d}\right)^m P_m(\mu_2), \qquad (A5) $$

$$ \varphi^{(2)}(r_1, \theta_1) = e^{-i\omega t} \sum_{n,m=0}^{\infty} (-1)^n a_n^{(2)} \xi_2^{n+1} C_{nm} \left(\frac{r_1}{d}\right)^m P_m(\mu_1), \qquad (A6) $$

where $\xi_j = R_{j0}/d$. Equation (A5) gives $\varphi^{(1)}$ in the coordinates $(r_2, \theta_2)$ near the surface of bubble 2, while (A6) gives $\varphi^{(2)}$ in the coordinates $(r_1, \theta_1)$ near the surface of bubble 1, i.e., these equations can be used in the corresponding boundary conditions.

For the coordinate transformation of $\Psi^{(1)}$ and $\Psi^{(2)}$, the following identities are used
(Varshalovich et al. 1988):

$$ h_n^{(1)}(k_v r_1) P_n^1(\mu_1) = \frac{\sqrt{n(n+1)}}{(2n+1)i^n} \sum_{l_1, l_2=0}^{\infty} \frac{i^{l_1-l_2} (-1)^{l_2} (2l_1+1)(2l_2+1)}{\sqrt{l_1(l_1+1)}} C_{l_1 0 l_2 0}^{n0} C_{l_1 1 l_2 0}^{n1} h_2^{(1)}(k_v d) j_{l_1}(k_v r_2) P_{l_1}^1(\mu_2), \qquad (A7) $$

$$ h_n^{(1)}(k_v r_2) P_n^1(\mu_2) = \frac{\sqrt{n(n+1)}}{(2n+1)i^n} \sum_{l_1, l_2=0}^{\infty} \frac{i^{l_1-l_2} (2l_1+1)(2l_2+1)}{\sqrt{l_1(l_1+1)}} C_{l_1 0 l_2 0}^{n0} C_{l_1 1 l_2 0}^{n1} h_2^{(1)}(k_v d) j_{l_1}(k_v r_1) P_{l_1}^1(\mu_1), \qquad (A8) $$
where $C_{l_1 m_1 l_2 m_2}^{LM}$ are the Clebsch-Gordan coefficients (Varshalovich et al. 1988; Zwillinger 2003) and $j_n$ is the spherical Bessel function. Equations (A7) and (A8) follow from equation (34) of § 5.17 of the book by Varshalovich et al. (1988). They are valid for $r_{1,2} < d$.

Substituting (A7) into (2.5) at $j=1$ and (A8) into (2.5) at $j=2$, one obtains

$$ \begin{aligned} \Psi^{(1)}(r_2, \theta_2) &= \Psi^{(1)}(r_2, \theta_2) e_{\varepsilon 2} \\ &= e^{-i\omega t} e_{\varepsilon 2} \sum_{n,m,l=0}^{\infty} \frac{\sqrt{n(n+1)}}{(2n+1)i^n} b_n^{(1)} \frac{i^{m-l} (-1)^l (2m+1)(2l+1)}{\sqrt{m(m+1)}} C_{m0l0}^{n0} C_{m1l0}^{n1} h_l^{(1)}(k_v d) j_m(k_v r_2) P_m^1(\mu_2) \end{aligned} \quad (A9) $$

$$ \begin{aligned} \Psi^{(2)}(r_1, \theta_1) &= \Psi^{(2)}(r_1, \theta_1) e_{\varepsilon 1} \\ &= e^{-i\omega t} e_{\varepsilon 1} \sum_{n,m,l=0}^{\infty} \frac{\sqrt{n(n+1)}}{(2n+1)i^n} b_n^{(2)} \frac{i^{m-l} (2m+1)(2l+1)}{\sqrt{m(m+1)}} C_{m0l0}^{n0} C_{m1l0}^{n1} h_l^{(1)}(k_v d) j_m(k_v r_1) P_m^1(\mu_1) \end{aligned} \quad (A10) $$

Here, we have also used the fact that $e_{\varepsilon 1} = e_{\varepsilon 2}$. Just as (A5) and (A6), (A9) and (A10) are valid near bubbles 2 and 1, respectively, and hence can be used in the boundary conditions at the bubble surfaces.

We first apply (A1) to bubble 1. In this case, $v_r^{(1)}$ is calculated by (2.6) while $v_r^{(2)}$ is calculated by (A6) and (A10) to give

$$ v_r^{(2)}(r_1, \theta_1) = e^{-i\omega t} \sum_{n,m=0}^{\infty} P_m(\mu_1) \left\{ \frac{(-1)^n m}{d} C_{nm} \xi_2^{n+1} a_n^{(2)} \left(\frac{r_1}{d}\right)^{m-1} \right. \\ \left. - \frac{j_m(k_v r_1)}{r_1} \frac{\sqrt{n(n+1)}}{(2n+1)i^n} (2m+1) i^m \sqrt{m(m+1)} b_n^{(2)} \sum_{l=0}^{\infty} i^{-l} (2l+1) C_{m0l0}^{n0} C_{m1l0}^{n1} h_l^{(1)}(k_v d) \right\} \quad (A11) $$

Setting $j=1$ in (A1), substituting (2.6) with $j=1$ and (A11) and then equating terms at the same Legendre polynomials, we obtain

$$ (n+1)a_n^{(1)} + n(n+1)h_n^{(1)}(k_v R_{10})b_n^{(1)} - \sum_{m=0}^{\infty} \{(-1)^m n C_{mn} \xi_1^n \xi_2^{m+1} a_m^{(2)} $$
$$
\left.
\begin{aligned}
& -i^n (2n+1) \sqrt{n(n+1)} j_n(k_v R_{10}) \frac{\sqrt{m(m+1)}}{(2m+1)i^m} b_m^{(2)} \sum_{l=0}^{\infty} i^{-l} (2l+1) C_{n0l0}^{m0} C_{nl10}^{m1} h_l^{(1)}(k_v d) \\
& \qquad = i\omega R_{10} \sum_{m=M_1^{(1)}}^{M_{N_1}^{(1)}} \delta_{nm} s_m^{(1)}, \quad n \ge 0
\end{aligned}
\right\} , \qquad (A12)
$$

where $\delta_{nm}$ is the Kronecker delta.

When (A1) is applied to bubble 2, $v_r^{(2)}$ is calculated by (2.6) while $v_r^{(1)}$ is calculated by (A5) and (A9) to give

$$
v_r^{(1)}(r_2, \theta_2) = e^{-i\omega t} \sum_{n,m=0}^{\infty} P_m(\mu_2) \left\{ \frac{(-1)^m m}{d} C_{nm} \xi_1^{n+1} a_n^{(1)} \left(\frac{r_2}{d}\right)^{m-1} \right. \\
\left. - \frac{j_m(k_v r_2)}{r_2} \frac{\sqrt{n(n+1)}}{(2n+1)i^n} (2m+1)i^m \sqrt{m(m+1)} b_n^{(1)} \sum_{l=0}^{\infty} (-1)^l i^{-l} (2l+1) C_{m0l0}^{n0} C_{ml10}^{n1} h_l^{(1)}(k_v d) \right\} . \quad (A13)
$$

Setting $j = 2$ in (A1), substituting (2.6) with $j = 2$ and (A13) and then equating terms at
the same Legendre polynomials, we obtain

$$
\begin{align}
& (n+1)a_n^{(2)} + n(n+1)h_n^{(1)}(k_v R_{20})b_n^{(2)} - \sum_{m=0}^{\infty} \{(-1)^n n C_{mn} \xi_1^{m+1} \xi_2^n a_m^{(1)} \nonumber \\
& \quad -i^n (2n+1)\sqrt{n(n+1)}j_n(k_v R_{20})\frac{\sqrt{m(m+1)}}{(2m+1)i^m} b_m^{(1)}\sum_{l=0}^{\infty} (-1)^l i^{-l} (2l+1)C_{n0l0}^{m0}C_{nl10}^{m1}h_l^{(1)}(k_v d) \} \nonumber \\
& = i\omega R_{20}\sum_{m=M_1^{(2)}}^{M_{N_2}^{(2)}} \delta_{nm} s_m^{(2)}, \quad n \ge 0 \tag{A14}
\end{align}
$$

To apply the boundary condition given by (A2), we need to calculate $v_{\theta}^{(1)}$ in the coordinates $(r_2, \theta_2)$ and $v_{\theta}^{(2)}$ in the coordinates $(r_1, \theta_1)$. By using (A5) and (A9) for calculating $v_{\theta}^{(1)}$ and (A6) and (A10) for calculating $v_{\theta}^{(2)}$, we obtain
$$
v_{\theta}^{(1)}(r_2, \theta_2) = \frac{e^{-i\omega t}}{r_2} \sum_{n,m=0}^{\infty} P_m^1(\mu_2) \left\{ (-1)^m C_{nm} \xi_1^{n+1} a_n^{(1)} \left(\frac{r_2}{d}\right)^m - \left[j_m(k_v r_2) + k_v r_2 j_m'(k_v r_2)\right] \frac{\sqrt{n(n+1)}}{(2n+1)i^n} \frac{(2m+1)i^m}{\sqrt{m(m+1)}} b_n^{(1)} \sum_{l=0}^{\infty} (-1)^l i^{-l} (2l+1) C_{m0l0}^{n0} C_{m1l0}^{n1} h_l^{(1)}(k_v d) \right\}, \quad (A15)
$$

$$
v_{\theta}^{(2)}(r_1, \theta_1) = \frac{e^{-i\omega t}}{r_1} \sum_{n,m=0}^{\infty} P_m^1(\mu_1) \left\{ (-1)^n C_{nm} \xi_2^{n+1} a_n^{(2)} \left(\frac{r_1}{d}\right)^m - \left[j_m(k_v r_1) + k_v r_1 j_m'(k_v r_1)\right] \frac{\sqrt{n(n+1)}}{(2n+1)i^n} \frac{(2m+1)i^n}{\sqrt{m(m+1)}} b_n^{(2)} \sum_{l=0}^{\infty} i^{-l} (2l+1) C_{m0l0}^{n0} C_{m1l0}^{n1} h_l^{(1)}(k_v d) \right\}. \quad (A16)
$$

Let us apply (A2) to bubble 1. Setting $j = 1$, substituting (2.6) and (2.7) with $j = 1$, using (A11) and (A16) and then equating terms at the same associated Legendre polynomials, we obtain

$$
\begin{equation}
\begin{aligned}
& 2(n+2)a_n^{(1)} + [k_v^2 R_{10}^2 h_n^{(1)\prime\prime}(k_v R_{10}) + (n^2+n-2)h_n^{(1)}(k_v R_{10})]b_n^{(1)} \\
& \quad - \sum_{m=0}^{\infty} \left\{ (-1)^m 2(n-1)C_{mn} \xi_1^n \xi_2^{m+1} a_m^{(2)} \right. \\
& \qquad \left. - [k_v^2 R_{10}^2 j_n^{\prime\prime}(k_v R_{10}) + (n^2+n-2)j_n(k_v R_{10})] \frac{\sqrt{m(m+1)}}{(2m+1)i^m} \frac{(2n+1)i^n}{\sqrt{n(n+1)}} b_m^{(2)} \right. \\
& \qquad \left. \times \sum_{l=0}^{\infty} i^{-l} (2l+1) C_{n0l0}^{m0} C_{n1l0}^{m1} h_l^{(1)}(k_v d) \right\} = 0, \quad n \ge 1.
\end{aligned}
\tag{A17}
\end{equation}
$$

For bubble 2, we set $j = 2$ in (A2), substitute (2.6) and (2.7) with $j = 2$, use (A13) and (A15) and then equate terms at the same associated Legendre polynomials. As a result, we obtain

$$
\begin{equation}
\begin{aligned}
& 2(n+2)a_n^{(2)} + [k_v^2 R_{20}^2 h_n^{(1)\prime\prime}(k_v R_{20}) + (n^2+n-2)h_n^{(1)}(k_v R_{20})]b_n^{(2)} - \sum_{m=0}^{\infty} \left\{ (-1)^n 2(n-1)C_{mn} \xi_1^{m+1} \xi_2^n a_m^{(1)} \right. \\
& \qquad \left. - [k_v^2 R_{20}^2 j_n^{\prime\prime}(k_v R_{20}) + (n^2+n-2)j_n(k_v R_{20})] \frac{\sqrt{m(m+1)}}{(2m+1)i^m} \frac{(2n+1)i^n}{\sqrt{n(n+1)}} b_m^{(1)} \right. \\
& \qquad \left. \times \sum_{l=0}^{\infty} (-1)^l i^{-l} (2l+1) C_{n0l0}^{m0} C_{n1l0}^{m1} h_l^{(1)}(k_v d) \right\} = 0, \quad n \ge 1.
\end{aligned}
\tag{A18}
\end{equation}
$$
From (A12) and (A14) at $n = 0$, it follows that

$$a_0^{(1)} = i\omega R_{10} \sum_{m=M_1^{(1)}}^{M_1^{(1)}} \delta_{0m} s_m^{(1)}, \qquad (A19)$$

$$a_0^{(2)} = i\omega R_{20} \sum_{m=M_1^{(2)}}^{M_2^{(2)}} \delta_{0m} s_m^{(2)}, \qquad (A20)$$

which means that $a_0^{(j)}$ is nonzero if mode 0 exists and hence $s_0^{(j)}$ is nonzero.

For $n \ge 1$, (A12), (A14), (A17) and (A18) are combined in the following system:

$$a_n^{(1)} + f_n^{(1)} b_n^{(1)} - \sum_{m=1}^{\infty} \{n\alpha_{2nm}a_m^{(2)} - \beta_{2nm}b_m^{(2)}\} = B_{1n}, \qquad (A21)$$

$$a_n^{(2)} + f_n^{(2)} b_n^{(2)} - \sum_{m=1}^{\infty} \{n\alpha_{1nm}a_m^{(1)} - \beta_{1nm}b_m^{(1)}\} = B_{2n}, \qquad (A22)$$

$$a_n^{(1)} + g_n^{(1)} b_n^{(1)} - \sum_{m=1}^{\infty} \left\{ \frac{n^2 - 1}{n + 2} \alpha_{2nm} a_m^{(2)} - \gamma_{2nm} b_m^{(2)} \right\} = B_{3n}, \qquad (A23)$$

$$a_n^{(2)} + g_n^{(2)} b_n^{(2)} - \sum_{m=1}^{\infty} \left\{ \frac{n^2 - 1}{n + 2} \alpha_{1nm} a_m^{(1)} - \gamma_{1nm} b_m^{(1)} \right\} = B_{4n}, \qquad (A24)$$

where the following designations are used:

$$f_n^{(j)} = n h_n^{(1)} (k_v R_{j0}), \qquad (A25)$$

$$g_n^{(j)} = \frac{k_v^2 R_{j0}^2 h_n^{(1)\prime\prime} (k_v R_{j0}) + (n^2 + n - 2)h_n^{(1)} (k_v R_{j0})}{2(n+2)}, \qquad (A26)$$

$$\alpha_{1nm} = \frac{(-1)^n C_{mn} \xi_1^{m+1} \xi_2^n}{n+1}, \qquad (A27)$$

$$\alpha_{2nm} = \frac{(-1)^m C_{mn} \xi_1^n \xi_2^{m+1}}{n+1}, \qquad (A28)$$

$$\beta_{1nm} = nj_n(k_v R_{20})G_{nm}, \qquad (A29)$$
$$ \beta_{2nm} = nj_n(k_v R_{10})H_{nm}, \quad (A30) $$

$$ \gamma_{1nm} = \frac{[k_v^2 R_{20}^2 j_n''(k_v R_{20}) + (n^2 + n - 2)j_n(k_v R_{20})]G_{nm}}{2(n+2)}, \quad (A31) $$

$$ \gamma_{2nm} = \frac{[k_v^2 R_{10}^2 j_n''(k_v R_{10}) + (n^2 + n - 2)j_n(k_v R_{10})]H_{nm}}{2(n+2)}, \quad (A32) $$

$$ G_{nm} = \frac{i^{n-m}(2n+1)\sqrt{m(m+1)}}{(2m+1)\sqrt{n(n+1)}} \sum_{l=|n-m|}^{n+m} (-i)^{-l}(2l+1)C_{n0l0}^{m0} C_{n1l0}^{m1} h_l^{(1)}(k_v d), \quad (A33) $$

$$ H_{nm} = \frac{i^{n-m}(2n+1)\sqrt{m(m+1)}}{(2m+1)\sqrt{n(n+1)}} \sum_{l=|n-m|}^{n+m} i^{-l}(2l+1)C_{n0l0}^{m0} C_{n1l0}^{m1} h_l^{(1)}(k_v d), \quad (A34) $$

$$ B_{1n} = \frac{n}{n+1} \xi_1^n \xi_2 a_0^{(2)} + \frac{i\omega R_{10}}{n+1} \sum_{m=M_1^{(1)}}^{M_1^{(1)}} \delta_{nm} s_m^{(1)}, \quad (A35) $$

$$ B_{2n} = \frac{(-1)^n n}{n+1} \xi_1^n \xi_2 a_0^{(1)} + \frac{i\omega R_{20}}{n+1} \sum_{m=M_1^{(2)}}^{M_2^{(2)}} \delta_{nm} s_m^{(2)}, \quad (A36) $$

$$ B_{3n} = \frac{n-1}{n+2} \xi_1^n \xi_2 a_0^{(2)}, \quad (A37) $$

$$ B_{4n} = \frac{(-1)^n (n-1)}{n+2} \xi_1^n \xi_2 a_0^{(1)}. \quad (A38) $$

In (A33) and (A34), the triangle inequality has been used, which requires that the indices of $C_{l_1 m_1 l_2 m_2}^{LM}$ satisfy the following conditions: $l_1 + l_2 - L \ge 0$, $l_1 - l_2 + L \ge 0$, $-l_1 + l_2 + L \ge 0$ (Zwillinger 2003).

From (A22) and (A24), it follows that

$$ a_n^{(2)} = \frac{g_n^{(2)} B_{2n} - f_n^{(2)} B_{4n}}{g_n^{(2)} - f_n^{(2)}} - \frac{(n^2 - 1) f_n^{(2)} - n(n+2) g_n^{(2)}}{(n+2)(g_n^{(2)} - f_n^{(2)})} \sum_{m=1}^{\infty} \alpha_{1nm} a_m^{(1)} \\ - \frac{1}{g_n^{(2)} - f_n^{(2)}} \sum_{m=1}^{\infty} (g_n^{(2)} \beta_{1nm} - f_n^{(2)} \gamma_{1nm}) b_m^{(1)}, \quad (A39) $$
$$b_n^{(2)} = \frac{B_{2n} - B_{4n}}{f_n^{(2)} - g_n^{(2)}} + \frac{2n+1}{(n+2)(f_n^{(2)} - g_n^{(2)})} \sum_{m=1}^{\infty} \alpha_{1nm} a_m^{(1)} - \frac{1}{f_n^{(2)} - g_n^{(2)}} \sum_{m=1}^{\infty} (\beta_{1nm} - \gamma_{1nm}) b_m^{(1)} \quad (A40)$$

Given $a_n^{(1)}$ and $b_n^{(1)}$, (A39) and (A40) allow one to calculate $a_n^{(2)}$ and $b_n^{(2)}$.

From (A21) and (A23), by using (A39) and (A40), one obtains

$$\sum_{k=1}^{\infty} A_{1nk} a_k^{(1)} + \sum_{k=1}^{\infty} A_{2nk} b_k^{(1)} = D_{1n}, \quad (A41)$$

$$\sum_{k=1}^{\infty} A_{3nk} a_k^{(1)} + \sum_{k=1}^{\infty} A_{4nk} b_k^{(1)} = D_{2n}, \quad (A42)$$

where

$$A_{1nk} = \delta_{nk} (g_n^{(1)} - f_n^{(1)}) - \sum_{m=1}^{\infty} \frac{(2m+1)(g_n^{(1)} \beta_{2nm} - f_n^{(1)} \gamma_{2nm}) \alpha_{1mk}}{(m+2)(g_m^{(2)} - f_m^{(2)})} \\ - \frac{(n^2-1)f_n^{(1)} - n(n+2)g_n^{(1)}}{n+2} \sum_{m=1}^{\infty} \frac{[(m^2-1)f_m^{(2)} - m(m+2)g_m^{(2)}] \alpha_{2nm} \alpha_{1mk}}{(m+2)(g_m^{(2)} - f_m^{(2)})}, \quad (A43)$$

$$A_{2nk} = \sum_{m=1}^{\infty} \frac{(\beta_{1mk} - \gamma_{1mk})(g_n^{(1)}\beta_{2nm} - f_n^{(1)}\gamma_{2nm})}{g_m^{(2)} - f_m^{(2)}} \\ + \frac{n(n+2)g_n^{(1)} - (n^2-1)f_n^{(1)}}{n+2} \sum_{m=1}^{\infty} \frac{\alpha_{2nm}(g_m^{(2)}\beta_{1mk} - f_m^{(2)}\gamma_{1mk})}{g_m^{(2)} - f_m^{(2)}}, \quad (A44)$$

$$D_{1n} = g_n^{(1)} B_{1n} - f_n^{(1)} B_{3n} + \frac{n(n+2)g_n^{(1)} - (n^2-1)f_n^{(1)}}{n+2} \sum_{m=1}^{\infty} \frac{(g_m^{(2)} B_{2m} - f_m^{(2)} B_{4m}) \alpha_{2nm}}{g_m^{(2)} - f_m^{(2)}} \\ + \sum_{m=1}^{\infty} \frac{(g_n^{(1)} \beta_{2nm} - f_n^{(1)} \gamma_{2nm})(B_{2m} - B_{4m})}{g_m^{(2)} - f_m^{(2)}}, \quad (A45)$$

$$A_{3nk} = \frac{2n+1}{n+2} \sum_{m=1}^{\infty} \frac{[(m^2-1)f_m^{(2)} - m(m+2)g_m^{(2)}] \alpha_{2nm} \alpha_{1mk}}{(m+2)(g_m^{(2)} - f_m^{(2)})} - \sum_{m=1}^{\infty} \frac{(2m+1)(\beta_{2nm} - \gamma_{2nm}) \alpha_{1mk}}{(m+2)(g_m^{(2)} - f_m^{(2)})}, \quad (A46)$$

$$A_{4nk} = (f_n^{(1)} - g_n^{(1)}) \delta_{nk} + \frac{2n+1}{n+2} \sum_{m=1}^{\infty} \frac{(g_m^{(2)} \beta_{1mk} - f_m^{(2)} \gamma_{1mk}) \alpha_{2nm}}{g_m^{(2)} - f_m^{(2)}} + \sum_{m=1}^{\infty} \frac{(\beta_{1mk} - \gamma_{1mk})(\beta_{2nm} - \gamma_{2nm})}{g_m^{(2)} - f_m^{(2)}}, \quad (A47)$$
$$D_{2n} = B_{1n} - B_{3n} + \frac{2n+1}{n+2} \sum_{m=1}^{\infty} \frac{(g_m^{(2)} B_{2m} - f_m^{(2)} B_{4m}) \alpha_{2nm}}{g_m^{(2)} - f_m^{(2)}} + \sum_{m=1}^{\infty} \frac{(\beta_{2nm} - \gamma_{2nm}) (B_{2m} - B_{4m})}{g_m^{(2)} - f_m^{(2)}}. \quad (A48)$$

To find $a_n^{(1)}$ and $b_n^{(1)}$, we have to solve (A41) and (A42). The scattering coefficients decrease for $n \to \infty$; see figure 2. Therefore, we can truncate (A41) and (A42). For example, we can take $N_a$ coefficients of $a_n^{(1)}$ with $n$ from 1 to $N_a$ and $N_b$ coefficients of $b_n^{(1)}$ with $n$ from 1 to $N_b$. As a result, we will obtain the following truncated equations:

$$\sum_{k=1}^{N_a} A_{1nk} a_k^{(1)} + \sum_{k=1}^{N_b} A_{2nk} b_k^{(1)} = D_{1n}, \quad n=1,2,\dots,N_a, \qquad (A49)$$

$$\sum_{k=1}^{N_a} A_{3nk} a_k^{(1)} + \sum_{k=1}^{N_b} A_{4nk} b_k^{(1)} = D_{2n}, \quad n=1,2,\dots,N_b. \qquad (A50)$$

As one can see, we obtain $N_a$ equations from (A49) and $N_b$ equations from (A50), which give in total a system of $N_a+N_b$ equations in $N_a$ unknowns $a_n^{(1)}$ and $N_b$ unknowns $b_n^{(1)}$.

With the help of (A33) and (A34), the coordinate transformations for $\psi^{(j)}$, $v_r^{(j)}$, and $v_\theta^{(j)}$ are written in a shorter form,

$$\psi^{(1)}(r_2, \theta_2) = e^{-i\omega t} \sum_{n,m=1}^{\infty} b_m^{(1)} G_{nm} j_n(k_v r_2) P_n^1(\mu_2), \qquad (A51)$$

$$v_r^{(1)}(r_2, \theta_2) = \frac{e^{-i\omega t}}{r_2} \sum_{\substack{n=1 \\ m=0}}^{\infty} n \left[ a_m^{(1)} (-1)^n C_{mn} \xi_1^{m+1} \left(\frac{r_2}{d}\right)^n - b_m^{(1)} (n+1) G_{nm} j_n(k_v r_2) \right] P_n(\mu_2), \quad (A52)$$

$$v_{\theta}^{(1)}(r_2, \theta_2) = \frac{e^{-i\omega t}}{r_2} \sum_{\substack{n=1 \\ m=0}}^{\infty} \left\{ a_m^{(1)} (-1)^n C_{mn} \xi_1^{m+1} \left(\frac{r_2}{d}\right)^n - b_m^{(1)} G_{nm} \left[ j_n(k_v r_2) + k_v r_2 j_n'(k_v r_2) \right] \right\} P_n^1(\mu_2), \quad (A53)$$

$$\psi^{(2)}(r_1, \theta_1) = e^{-i\omega t} \sum_{n,m=1}^{\infty} b_m^{(2)} H_{nm} j_n(k_v r_1) P_n^1(\mu_1), \qquad (A54)$$
$$v_r^{(2)}(r_1, \theta_1) = \frac{e^{-i\omega t}}{r_1} \sum_{\substack{n=1 \\ m=0}}^{\infty} n \left[ a_m^{(2)}(-1)^m C_{mn} \xi_2^{m+1} \left(\frac{r_1}{d}\right)^n - b_m^{(2)} (n+1) H_{nm} j_n(k_v r_1) \right] P_n(\mu_1), \quad (A55)$$

$$v_{\theta}^{(2)}(r_1, \theta_1) = \frac{e^{-i\omega t}}{r_1} \sum_{\substack{n=1 \\ m=0}}^{\infty} \left\{ a_m^{(2)}(-1)^m C_{mn} \xi_2^{m+1} \left(\frac{r_1}{d}\right)^n - b_m^{(2)} H_{nm} \left[ j_n(k_v r_1) + k_v r_1 j'_n(k_v r_1) \right] \right\} P_n^1(\mu_1) \quad . \quad (A56)$$

These equations are used in the calculation of acoustic streaming.

## Appendix B. Representation of the equation for $\Psi$ in different coordinates

Equation (2.11) can be written as

$$\Delta^2 \Psi = T_1 + T_2 + T_{12}, \quad (B1)$$

where

$$T_j = \frac{1}{2\nu} \mathrm{Re} \{ \nabla \times [(\Delta \boldsymbol{\psi}^{(j)*}) \times \boldsymbol{v}^{(j)}] \} \quad j=1,2, \quad (B2)$$

$$T_{12} = \frac{1}{2\nu} \mathrm{Re} \{ \nabla \times [(\Delta \boldsymbol{\psi}^{(1)*}) \times \boldsymbol{v}^{(2)}] + \nabla \times [(\Delta \boldsymbol{\psi}^{(2)*}) \times \boldsymbol{v}^{(1)}] \} \quad (B3)$$

By using the coordinate transformations given in Appendix A, (B1) – (B3) can be expressed in the coordinates $(r_1, \theta_1)$ or $(r_2, \theta_2)$. Both representations are needed to calculate acoustic streaming around two bubbles.

Let us first express (B1) – (B3) in the coordinates $(r_1, \theta_1)$.

As shown by Doinikov et al. (2019a), (B2) with $j=1$ can be represented as

$$T_1(r_1, \theta_1) = \frac{e_{\epsilon 1}}{2\nu r_1} \mathrm{Re} \left\{ \left[ \frac{\partial}{\partial r_1} \left[ r_1 v_r^{(1)} \left( \Delta_{r\theta}^{(1)} \psi^{(1)*} - \frac{\psi^{(1)*}}{r_1^2 \sin^2 \theta_1} \right) \right] + \frac{\partial}{\partial \theta_1} \left[ v_{\theta}^{(1)} \left( \Delta_{r\theta}^{(1)} \psi^{(1)*} - \frac{\psi^{(1)*}}{r_1^2 \sin^2 \theta_1} \right) \right] \right\}, \quad (B4)$$

where the operator $\Delta_{r\theta}^{(j)}$ is defined by

$$\Delta_{r\theta}^{(j)} = \frac{1}{r_j^2} \frac{\partial}{\partial r_j} \left( r_j^2 \frac{\partial}{\partial r_j} \right) + \frac{1}{r_j^2 \sin \theta_j} \frac{\partial}{\partial \theta_j} \left( \sin \theta_j \frac{\partial}{\partial \theta_j} \right). \quad (B5)$$
$\Psi^{(j)}$ obeys the equation $\Delta\Psi^{(j)} = -k_v^2\Psi^{(j)}$, which gives

$$ \Delta_{r\theta}^{(j)} \psi^{(j)} - \frac{\psi^{(j)}}{r_j^2 \sin^2 \theta_j} = -k_v^2 \psi^{(j)}. \quad (B6) $$

Considering (B6), (B4) reduces to

$$ T_1(r_1, \theta_1) = T_1(r_1, \theta_1) e_{\varepsilon 1} = \frac{e_{\varepsilon 1}}{2\nu r_1} \text{Re} \left\{ k_v^2 \left[ \frac{\partial}{\partial r_1} (r_1 v_r^{(1)} \psi^{(1)*}) + \frac{\partial}{\partial \theta_1} (v_\theta^{(1)} \psi^{(1)*}) \right] \right\}. \quad (B7) $$

Substitution of (2.5) – (2.7) into (B7) yields

$$ T_1(r_1, \theta_1) = \frac{1}{2\nu R_{10}^4} \left[ \sum_{n=0}^{\infty} P_n(\mu_1) P_m^1(\mu_1) K_{1nm}^{(1)}(r_1) + \sqrt{1-\mu_1^2} \sum_{n,m=1}^{\infty} \left[ P_n^1(\mu_1) P_m^1(\mu_1) \right]' K_{2nm}^{(1)}(r_1) \right], \quad (B8) $$

where

$$ K_{1nm}^{(1)}(r_1) = \text{Re} \left\{ b_m^{(1)*} \frac{(n+1)k_v^2 R_{10}^3}{r_1} \left\{ a_n^{(1)} \left( \frac{R_{10}}{r_1} \right)^{n+1} \left[ \frac{(n+1)R_{10}}{r_1} h_m^{(1)*}(k_v r_1) - k_v^* R_{10} h_m^{(1)/*}(k_v r_1) \right] \right. \right. \\ \left. \left. -nb_n^{(1)} \left[ k_v R_{10} h_n^{(1)/*}(k_v r_1) h_m^{(1)*}(k_v r_1) + k_v^* R_{10} h_n^{(1)}(k_v r_1) h_m^{(1)/*}(k_v r_1) \right] \right\} \right\}, \quad (B9) $$

$$ K_{2nm}^{(1)}(r_1) = \text{Re} \left\{ b_m^{(1)*} \frac{k_v^2 R_{10}^4 h_m^{(1)*}(k_v r_1)}{r_1^2} \left[ b_n^{(1)}[h_n^{(1)}(k_v r_1) + k_v r_1 h_n^{(1)}/(k_v r_1)] - a_n^{(1)} \left(\frac{R_{10}}{r_1}\right)^{n+1} \right] \right\}. \quad (B10) $$

In a similar way, $T_2$ in the coordinates $(r_1, \theta_1)$ is written as

$$ T_2(r_1, \theta_1) = T_2(r_1, \theta_1) e_{\varepsilon 1} = \frac{e_{\varepsilon 1}}{2\nu r_1} \text{Re} \left\{ k_v^2 \left[ \frac{\partial}{\partial r_1} (r_1 v_r^{(2)} \psi^{(2)*}) + \frac{\partial}{\partial \theta_1} (v_\theta^{(2)} \psi^{(2)*}) \right] \right\}, \quad (B11) $$

where $\psi^{(2)}$, $v_r^{(2)}$, and $v_\theta^{(2)}$ in the coordinates $(r_1, \theta_1)$ are given by (A54) – (A56). Substitution of these equations into (B11) yields

$$ T_2(r_1, \theta_1) = \frac{1}{2\nu R_{10}^4} \left[ \sum_{n,m=1}^{\infty} P_n(\mu_1) P_m^1(\mu_1) L_{1nm}^{(1)}(r_1) + \sqrt{1-\mu_1^2} \sum_{n,m=1}^{\infty} \left[ P_n^1(\mu_1) P_m^1(\mu_1) \right]' L_{2nm}^{(1)}(r_1) \right], \quad (B12) $$
where

$$
L_{1nm}^{(1)}(r_1) = \operatorname{Re} \left\{ \frac{nk_v^2 R_{10}^3}{r_1} \sum_{k=1}^{\infty} (b_k^{(2)} H_{mk})^* \left\{ \left(\frac{r_1}{d}\right)^n \left[ k_v^* R_{10} j_m'(k_v^* r_1) + \frac{nR_{10}}{r_1} j_m(k_v^* r_1) \right] \sum_{k=0}^{\infty} a_k^{(2)} (-1)^k C_{kn} \xi_2^{k+1} \right. \right. \\
\left. \left. -(n+1) \left[ k_v R_{10} j_n'(k_v r_1) j_m(k_v^* r_1) + k_v^* R_{10} j_n(k_v r_1) j_m'(k_v^* r_1) \right] \sum_{k=1}^{\infty} b_k^{(2)} H_{nk} \right\} \right\}, \quad (B13)
$$

$$
L_{2nm}^{(1)}(r_1) = \operatorname{Re} \left\{ \frac{k_v^2 R_{10}^4 j_m(k_v^* r_1)}{r_1^2} \sum_{k=1}^{\infty} (b_k^{(2)} H_{mk})^* \right. \\
\left. \times \left[ [j_n(k_v r_1) + k_v r_1 j_n'(k_v r_1)] \sum_{k=1}^{\infty} b_k^{(2)} H_{nk} - \left(\frac{r_1}{d}\right)^n \sum_{k=0}^{\infty} a_k^{(2)} (-1)^k C_{kn} \xi_2^{k+1} \right] \right\}. \quad (B14)
$$

Before representing (B3) in the coordinates ( $^{r_1}, \theta_1$ ), let us transform it to a more
convenient form.

As $\boldsymbol{\psi}^{(j)}$ obeys the equation $\Delta\boldsymbol{\psi}^{(j)} = -k_v^2\boldsymbol{\psi}^{(j)}$, (B3) can be written as

$$
\mathbf{T}_{12} = \frac{1}{2\nu} \mathrm{Re} \left\{ k_v^2 \nabla \times (\boldsymbol{\psi}^{(1)*} \times \boldsymbol{v}^{(2)}) + k_v^2 \nabla \times (\boldsymbol{\psi}^{(2)*} \times \boldsymbol{v}^{(1)}) \right\}. \quad (B15)
$$

With the help of the identity

$$
\nabla \times (\mathbf{a} \times \mathbf{b}) = (\mathbf{b} \cdot \nabla)\mathbf{a} - (\mathbf{a} \cdot \nabla)\mathbf{b} + \mathbf{a}(\nabla \cdot \mathbf{b}) - \mathbf{b}(\nabla \cdot \mathbf{a}), \quad (B16)
$$

taking into account that $\nabla \cdot \mathbf{v}^{(j)} = 0$ and $\nabla \cdot \boldsymbol{\psi}^{(j)} = 0$, (B15) is transformed to

$$
\mathbf{T}_{12} = \frac{1}{2\nu} \mathrm{Re} \left\{ k_v^2 \left[ (\mathbf{v}^{(2)} \cdot \nabla) \boldsymbol{\psi}^{(1)*} - (\boldsymbol{\psi}^{(1)*} \cdot \nabla) \mathbf{v}^{(2)} + (\mathbf{v}^{(1)} \cdot \nabla) \boldsymbol{\psi}^{(2)*} - (\boldsymbol{\psi}^{(2)*} \cdot \nabla) \mathbf{v}^{(1)} \right] \right\}. \quad (B17)
$$

Expressed in the coordinates ( $^{r_1}, \theta_1$ ), (B17) gives

$$
\mathbf{T}_{12}(r_1, \theta_1) = \mathbf{T}_{12}(r_1, \theta_1) e_{\epsilon 1} = \frac{e_{\epsilon 1}}{2\nu} \operatorname{Re} \left\{ k_v^2 \left[ \begin{aligned}[t]
& v_r^{(1)} \frac{\partial \psi^{(2)*}}{\partial r_1} + v_r^{(2)} \frac{\partial \psi^{(1)*}}{\partial r_1} + \frac{v_\theta^{(1)}}{r_1} \frac{\partial \psi^{(2)*}}{\partial \theta_1} + \frac{v_\theta^{(2)}}{r_1} \frac{\partial \psi^{(1)*}}{\partial \theta_1} \\
& - \frac{\psi^{(1)*} (v_r^{(2)} \sin \theta_1 + v_\theta^{(2)} \cos \theta_1) + \psi^{(2)*} (v_r^{(1)} \sin \theta_1 + v_\theta^{(1)} \cos \theta_1)}{r_1 \sin \theta_1}
\end{aligned} \right] \right\}, \quad (B18)
$$
where $\Psi^{(1)}$, $v_r^{(1)}$ and $v_\theta^{(1)}$ are calculated by (2.5) – (2.7) while $\Psi^{(2)}$, $v_r^{(2)}$ and $v_\theta^{(2)}$ are calculated by (A54) – (A56). Substitution of these equations into (B18) yields

$$ T_{12}(r_1, \theta_1) = \frac{1}{2\nu R_{10}^4} \left[ \sum_{n=0}^{\infty} P_n(\mu_1) P_m^1(\mu_1) M_{1nm}^{(1)}(r_1) + \sqrt{1-\mu_1^2} \sum_{n,m=1}^{\infty} [P_n^1(\mu_1) P_m^1(\mu_1)]' M_{2nm}^{(1)}(r_1) \right], \quad (B19) $$

where

$$ M_{1nm}^{(1)}(r_1) = \operatorname{Re} \left\{ \frac{k_v^2 R_{10}^4}{r_1^2} \left[ (n+1) \left(\frac{R_{10}}{r_1}\right)^{n+1} \left[ (n+1) j_m(k_v^* r_1) - k_v^* r_1 j_m'(k_v^* r_1) \right] a_n^{(1)} \sum_{k=1}^{\infty} b_k^{(2)*} H_{mk}^* \right. \right. \\ 
+ n \left(\frac{r_1}{d}\right)^n \left[ n h_m^{(1)*}(k_v r_1) + k_v^* r_1 h_m^{(1)/*}(k_v r_1) \right] b_m^{(1)*} \sum_{k=0}^{\infty} a_k^{(2)} (-1)^k C_{kn} \xi_2^{k+1} \\ 
- n(n+1) \left[ k_v r_1 j_m(k_v^* r_1) h_n^{(1)}/(k_v r_1) + k_v^* r_1 j_m'(k_v^* r_1) h_n^{(1)}(k_v r_1) \right] b_n^{(1)} \sum_{k=1}^{\infty} b_k^{(2)*} H_{mk}^* \\ 
\left. \left. - n(n+1) \left[ k_v^* r_1 j_n(k_v r_1) h_m^{(1)/*}(k_v r_1) + k_v r_1 j_n'(k_v r_1) h_m^{(1)/*}(k_v r_1) \right] b_m^{(1)*} \sum_{k=1}^{\infty} b_k^{(2)} H_{nk} \right] \right\}, \quad (B20) $$

$$ M_{2nm}^{(1)}(r_1) = \operatorname{Re} \left\{ \frac{k_v^2 R_{10}^4}{r_1^2} \left[ j_m(k_v r_1) \left(\frac{R_{10}}{r_1}\right)^{n+1} a_n^{(1)*} \sum_{k=1}^{\infty} b_k^{(2)} H_{mk} + h_n^{(1)}(k_v r_1) \left(\frac{r_1}{d}\right)^m b_n^{(1)} \sum_{k=0}^{\infty} a_k^{(2)*} (-1)^k C_{km} \xi_2^{k+1} \right. \right. \\ 
+ \left. \left[ k_v r_1 j_m'(k_v r_1) h_n^{(1)*}(k_v r_1) - k_v^* r_1 j_m(k_v r_1) h_n^{(1)/*}(k_v r_1) \right] b_n^{(1)*} \sum_{k=1}^{\infty} b_k^{(2)} H_{mk} \right] \right\}. \quad (B21) $$

It follows from the above derivation that $\Psi$ can be represented as

$$ \Psi(r_1, \theta_1) = \Psi(r_1, \theta_1) e_{\epsilon 1}, \quad (B22) $$

where $\Psi(r_1, \theta_1)$ obeys the following equation:

$$ \left( \Delta_{r\theta}^{(1)} - \frac{1}{r_1^2 \sin^2 \theta_1} \right)^2 \Psi(r_1, \theta_1) = \frac{1}{2\nu R_{10}^4} \left\{ \sum_{n=0}^{\infty} P_n(\mu_1) P_m^1(\mu_1) \left[ K_{1nm}^{(1)}(r_1) + L_{1nm}^{(1)}(r_1) + M_{1nm}^{(1)}(r_1) \right] \right\} $$
$$ + \sqrt{1 - \mu_1^2} \sum_{n,m=1}^{\infty} P_n^1(\mu_1) P_m^{1\dagger}(\mu_1) [K_{2nm}^{(1)}(r_1) + K_{2mn}^{(1)}(r_1) + L_{2nm}^{(1)}(r_1) + L_{2mn}^{(1)}(r_1) + M_{2nm}^{(1)}(r_1) + M_{2mn}^{(1)}(r_1)] \}. \quad (B23) $$

To express (B1) – (B3) in the coordinates $(r_2, \theta_2)$, the following manipulations are performed.

Interchanging indices 1 and 2 in (B11), one obtains

$$ T_1(r_2, \theta_2) = T_1(r_2, \theta_2) e_{\varepsilon 2} = \frac{e_{\varepsilon 2}}{2\nu r_2} \text{Re} \left\{ k_v^2 \left[ \frac{\partial}{\partial r_2} (r_2 v_r^{(1)} \psi^{(1)*}) + \frac{\partial}{\partial \theta_2} (v_\theta^{(1)} \psi^{(1)*}) \right] \right\}. \quad (B24) $$

Substitution of (A51) – (A53) into (B24) yields

$$ T_1(r_2, \theta_2) = \frac{1}{2\nu R_{20}^4} \left[ \sum_{n,m=1}^{\infty} P_n(\mu_2) P_m^1(\mu_2) K_{1nm}^{(2)}(r_2) + \sqrt{1 - \mu_2^2} \sum_{n,m=1}^{\infty} \left[ P_n^1(\mu_2) P_m^1(\mu_2) \right]' K_{2nm}^{(2)}(r_2) \right], \quad (B25) $$

where

$$ K_{1nm}^{(2)}(r_2) = \text{Re} \left\{ \frac{n k_v^2 R_{20}^3}{r_2} \sum_{k=1}^{\infty} (b_k^{(1)} G_{mk})^* \left\{ (-1)^n \left(\frac{r_2}{d}\right)^n \left[ k_v^* R_{20} j_m^\prime(k_v^* r_2) + \frac{n R_{20}}{r_2} j_m(k_v^* r_2) \right] \sum_{k=0}^{\infty} a_k^{(1)} C_{kn} \xi_1^{k+1} \right. \right. \\ \left. \left. -(n+1) \left[ k_v R_{20} j_n^\prime(k_v r_2) j_m(k_v^* r_2) + k_v^* R_{20} j_n(k_v r_2) j_m^\prime(k_v^* r_2) \right] \sum_{k=1}^{\infty} b_k^{(1)} G_{nk} \right\} \right\}, \quad (B26) $$

$$ K_{2nm}^{(2)}(r_2) = \text{Re} \left\{ \frac{k_v^2 R_{20}^4 j_m(k_v^* r_2)}{r_2^2} \sum_{k=1}^{\infty} (b_k^{(1)} G_{mk})^* \times \left[ \left[ j_n(k_v r_2) + k_v r_2 j_n^\prime(k_v r_2) \right] \sum_{k=1}^{\infty} b_k^{(1)} G_{nk} - (-1)^n \left(\frac{r_2}{d}\right)^n \sum_{k=0}^{\infty} a_k^{(1)} C_{kn} \xi_1^{k+1} \right] \right\}. \quad (B27) $$

Replacing index 1 in (B7) with index 2, one obtains

$$ T_2(r_2, \theta_2) = T_2(r_2, \theta_2) e_{\varepsilon 2} = \frac{e_{\varepsilon 2}}{2\nu r_2} \text{Re} \left\{ k_v^2 \left[ \frac{\partial}{\partial r_2} (r_2 v_r^{(2)} \psi^{(2)*}) + \frac{\partial}{\partial \theta_2} (v_\theta^{(2)} \psi^{(2)*}) \right] \right\}. \quad (B28) $$

Substitution of (2.5) – (2.7) into (B28) yields
$$ T_2(r_2, \theta_2) = \frac{1}{2\nu R_{20}^4} \left[ \sum_{\substack{n=0 \\ m=1}}^{\infty} P_n(\mu_2) P_m^1(\mu_2) L_{1nm}^{(2)}(r_2) + \sqrt{1-\mu_2^2} \sum_{n,m=1}^{\infty} \left[ P_n^1(\mu_2) P_m^1(\mu_2) \right]' L_{2nm}^{(2)}(r_2) \right], \quad (B29) $$

where

$$ L_{1nm}^{(2)}(r_2) = \operatorname{Re} \left\{ b_m^{(2)*} \frac{(n+1)k_v^2 R_{20}^4}{r_2} \left\{ a_n^{(2)} \left( \frac{R_{20}}{r_2} \right)^{n+1} \left[ \frac{(n+1)R_{20}}{r_2} h_m^{(1)*}(k_v r_2) - k_v^* R_{20} h_m^{(1)/*}(k_v r_2) \right] \right. \right. \\ \left. \left. -nb_n^{(2)} \left[ k_v R_{20} h_n^{(1)}/(k_v r_2) h_m^{(1)}/(k_v r_2) + k_v^* R_{20} h_n^{(1)}/(k_v r_2) h_m^{(1)}/(k_v r_2) \right] \right\} \right\}, \quad (B30) $$

$$ L_{2nm}^{(2)}(r_2) = \operatorname{Re} \left\{ b_m^{(2)*} \frac{k_v^2 R_{20}^4 h_m^{(1)*}(k_v r_2)}{r_2^2} \left[ b_n^{(2)} \left[ h_n^{(1)}(k_v r_2) + k_v r_2 h_n^{(1)}/(k_v r_2) \right] - a_n^{(2)} \left( \frac{R_{20}}{r_2} \right)^{n+1} \right] \right\}. \quad (B31) $$

Interchanging indices 1 and 2 in (B18), one obtains

$$ T_{12}(r_2, \theta_2) = T_{12}(r_2, \theta_2) e_{\varepsilon 2} = \frac{e_{\varepsilon 2}}{2\nu} \operatorname{Re} \left\{ k_v^2 \left[ v_r^{(2)} \frac{\partial \psi^{(1)*}}{\partial r_2} + v_r^{(1)} \frac{\partial \psi^{(2)*}}{\partial r_2} + \frac{v_\theta^{(2)}}{r_2} \frac{\partial \psi^{(1)*}}{\partial \theta_2} + \frac{v_\theta^{(1)}}{r_2} \frac{\partial \psi^{(2)*}}{\partial \theta_2} \right. \right. \\ \left. \left. - \frac{\psi^{(2)*} (v_r^{(1)} \sin \theta_2 + v_\theta^{(1)} \cos \theta_2) + \psi^{(1)*} (v_r^{(2)} \sin \theta_2 + v_\theta^{(2)} \cos \theta_2)}{r_2 \sin \theta_2} \right] \right\}, \quad (B32) $$

where $\psi^{(2)}$, $v_r^{(2)}$ and $v_\theta^{(2)}$ are calculated by (2.5) – (2.7) while $\psi^{(1)}$, $v_r^{(1)}$ and $v_\theta^{(1)}$ are calculated by (A51) – (A53). Substitution of these equations into (B32) yields

$$ T_{12}(r_2, \theta_2) = \frac{1}{2\nu R_{20}^4} \left[ \sum_{\substack{n=0 \\ m=1}}^{\infty} P_n(\mu_2) P_m^1(\mu_2) M_{1nm}^{(2)}(r_2) + \sqrt{1-\mu_2^2} \sum_{n,m=1}^{\infty} \left[ P_n^1(\mu_2) P_m^1(\mu_2) \right]' M_{2nm}^{(2)}(r_2) \right], \quad (B33) $$

where

$$ M_{1nm}^{(2)}(r_2) = \operatorname{Re} \left\{ \frac{k_v^2 R_{20}^4}{r_2^2} \left[ (n+1) \left( \frac{R_{20}}{r_2} \right)^{n+1} \left[ (n+1) j_m(k_v^* r_2) - k_v^* r_2 j_m'(k_v^* r_2) \right] a_n^{(2)} \sum_{k=1}^{\infty} b_k^{(1)*} G_{mk}^* \right. \right. \\ \left. \left. + (-1)^n n \left( \frac{r_2}{d} \right)^n \left[ n h_m^{(1)*}(k_v r_2) + k_v^* r_2 h_m^{(1)/*}(k_v r_2) \right] b_m^{(2)*} \sum_{k=0}^{\infty} a_k^{(1)} C_{kn} \xi_1^{k+1} \right] \right\} $$
$$-n(n+1)\left[k_v r_2 j_m(k_v^* r_2) h_n^{(1)\dagger}(k_v r_2) + k_v^* r_2 j'_m(k_v^* r_2) h_n^{(1)}(k_v r_2)\right] b_n^{(2)\dagger} \sum_{k=1}^{\infty} b_k^{(1)\dagger} G_{mk}^* \\ -n(n+1)\left[k_v^* r_2 j_n(k_v r_2) h_m^{(1)\dagger\dagger}(k_v r_2) + k_v r_2 j'_n(k_v r_2) h_m^{(1)\dagger\dagger}(k_v r_2)\right] b_m^{(2)\dagger\dagger} \sum_{k=1}^{\infty} b_k^{(1)\dagger\dagger} G_{nk}^*\right\}, \quad (B34)$$

$$M_{2nm}^{(2)}(r_2) = \operatorname{Re} \left\{ \frac{k_v^2 R_{20}^4}{r_2^2} \left[ j_m(k_v r_2) \left( \frac{R_{20}}{r_2} \right)^{n+1} a_n^{(2)*} \sum_{k=1}^{\infty} b_k^{(1)} G_{mk} + (-1)^m h_n^{(1)}(k_v r_2) \left( \frac{r_2}{d} \right)^m b_n^{(2)} \sum_{k=0}^{\infty} a_k^{(1)*} C_{km} \xi_1^{k+1} \right. \right. \\ \left. \left. + \left[ k_v r_2 j'_m(k_v r_2) h_n^{(1)*}(k_v r_2) - k_v^* r_2 j_m(k_v r_2) h_n^{(1)\dagger\dagger}(k_v r_2) \right] b_n^{(2)*} \sum_{k=1}^{\infty} b_k^{(1)} G_{mk} \right] \right\}. \quad (B35)$$

Representing $\Psi$ as

$$\Psi(r_2, \theta_2) = \Psi(r_2, \theta_2) e_{\varepsilon_2} \quad (B36)$$

and using (B24) – (B35), one obtains

$$\left( \Delta_{r\theta}^{(2)} - \frac{1}{r_2^2 \sin^2 \theta_2} \right)^2 \Psi(r_2, \theta_2) = \frac{1}{2\nu R_{20}^4} \left\{ \begin{aligned} & \sum_{n=0}^{\infty} P_n(\mu_2) P_m^1(\mu_2) \left[ K_{1nm}^{(2)}(r_2) + L_{1nm}^{(2)}(r_2) + M_{1nm}^{(2)}(r_2) \right] \\ & + \sqrt{1 - \mu_2^2} \sum_{n,m=1}^{\infty} P_n^1(\mu_2) P_m^1(\mu_2) \left[ K_{2nm}^{(2)}(r_2) + K_{2mn}^{(2)}(r_2) + L_{2nm}^{(2)}(r_2) + L_{2mn}^{(2)}(r_2) + M_{2nm}^{(2)}(r_2) + M_{2mn}^{(2)}(r_2) \right] \end{aligned} \right\}. \quad (B37)$$

## Appendix C. Solution of (2.13) and (2.15)

Equations (2.13) and (2.15) can be represented as

$$\left( \Delta_{r\theta}^{(j)} - \frac{1}{r_j^2 \sin^2 \theta_j} \right)^2 \Psi(r_j, \theta_j) = E^{(j)}(r_j, \mu_j), \quad j=1,2, \quad (C1)$$

where

$$E^{(j)}(r_j, \mu_j) = \frac{1}{2\nu R_{j0}^4} \left\{ \sum_{\substack{n=0 \\ m=1}}^{\infty} P_n(\mu_j) P_m^1(\mu_j) \left[ K_{1nm}^{(j)}(r_j) + L_{1nm}^{(j)}(r_j) + M_{1nm}^{(j)}(r_j) \right] \right\}$$
$$ + \sqrt{1 - \mu_j^2} \sum_{n,m=1}^{\infty} P_n^{l'}(\mu_j) P_m^{l'}(\mu_j) [K_{2nm}^{(j)}(r_j) + K_{2mn}^{(j)}(r_j) + L_{2nm}^{(j)}(r_j) + L_{2mn}^{(j)}(r_j) + M_{2nm}^{(j)}(r_j) + M_{2mn}^{(j)}(r_j)] . \quad (C2) $$

The function $E_l^{(j)}(r_j, \mu_j)$ can be expanded as follows:

$$ E_l^{(j)}(r_j, \mu_j) = \sum_{l'=1}^{\infty} E_{l'}^{(j)}(r_{j'}) P_l^{l'}(\mu_j), \quad (C3) $$

where the expansion coefficients $E_l^{(j)}(r_j)$ are calculated by

$$ E_l^{(j)}(r_j) = \frac{2l+1}{2l(l+1)} \int_{-1}^{1} E^{(j)}(r_j, \mu_j) P_l^{l'}(\mu_j) d\mu_j . \quad (C4) $$

With the help of (D5), (D9) and (D10), the calculation of $E_l^{(j)}(r_j)$ yields

$$ \begin{aligned} E_l^{(j)}(r_j) = & \frac{2l+1}{2\nu R_{j0}^4 \sqrt{l(l+1)}} \left\{ \sum_{n=0}^{\infty} \sum_{\substack{m=|n-l| \\ m \ge 1}}^{n+l} \frac{\sqrt{m(m+1)} C_{n0l0}^{m0} C_{n0l1}^{m1}}{2m+1} \left\{ K_{1nm}^{(j)}(r_j) + L_{1nm}^{(j)}(r_j) + M_{1nm}^{(j)}(r_j) \right. \right. \\ & \left. + n^2 \left[ K_{2nm}^{(j)}(r_j) + K_{2mn}^{(j)}(r_j) + L_{2nm}^{(j)}(r_j) + L_{2mn}^{(j)}(r_j) + M_{2nm}^{(j)}(r_j) + M_{2mn}^{(j)}(r_j) \right] \right\} \\ & - \sum_{n=2}^{\infty} \sum_{\substack{m=1 \\ k=|n-2m-l|}}^{[n/2]} \sum_{k \ge 1}^{n-2m+l} \frac{(2n-4m+1)\sqrt{k(k+1)} C_{(n-2m)0l0}^{k0} C_{(n-2m)0l1}^{k1}}{2k+1} \\ & \left. \times \left[ K_{2nk}^{(j)}(r_j) + K_{2kn}^{(j)}(r_j) + L_{2nk}^{(j)}(r_j) + L_{2kn}^{(j)}(r_j) + M_{2nk}^{(j)}(r_j) + M_{2kn}^{(j)}(r_j) \right] \right\}, \end{aligned} \quad (C5) $$

where $[]$ means the integer part of an expression in brackets.

Considering (C3), a solution to (C1) is sought as

$$ \Psi(r_j, \theta_j) = \sum_{l=1}^{\infty} \Psi_l^{(j)}(r_j) P_l^{l'}(\mu_j) . \quad (C6) $$

Substituting (C3) and (C6) into (C1) gives the following equation for $\Psi_l^{(j)}(r_j)$:
$$ \Psi_l^{(j)IV} + \frac{4}{r_j} \Psi_l^{(j)III} - \frac{2l(l+1)}{r_j^2} \Psi_l^{(j)II} + \frac{l(l+1)(l^2+l-2)}{r_j^4} \Psi_l^{(j)} = E_l^{(j)}(r_j), \quad (C7) $$

where the prime denotes the derivative with respect to $r_j$.

Equation (C7) is solved by the method of variation of constants (Boyce & DiPrima 2001). As a result, one obtains

$$ \Psi_l^{(j)}(r_j) = \frac{C_{1l}^{(j)}(r_j)}{r_j^{l-1}} + \frac{C_{2l}^{(j)}(r_j)}{r_j^{l+1}} + r_j' C_{3l}^{(j)}(r_j) + r_j^{l+2} C_{4l}^{(j)}(r_j), \quad (C8) $$

where the functions $C_{1l}^{(j)}(r_j) - C_{4l}^{(j)}(r_j)$ are calculated by

$$ C_{1l}^{(j)}(r_j) = C_{1l0}^{(j)} + \frac{1}{2(2l-1)(2l+1)} \int_{R_{j0}}^{r_j} s^{l+2} E_l^{(j)}(s) ds, \quad (C9) $$

$$ C_{2l}^{(j)}(r_j) = C_{2l0}^{(j)} - \frac{1}{2(2l+1)(2l+3)} \int_{R_{j0}}^{r_j} s^{l+4} E_l^{(j)}(s) ds, \quad (C10) $$

$$ C_{3l}^{(j)}(r_j) = C_{3l0}^{(j)} - \frac{1}{2(2l-1)(2l+1)} \int_{R_{j0}}^{r_j} s^{3-l} E_l^{(j)}(s) ds, \quad (C11) $$

$$ C_{4l}^{(j)}(r_j) = C_{4l0}^{(j)} + \frac{1}{2(2l+1)(2l+3)} \int_{R_{j0}}^{r_j} s^{1-l} E_l^{(j)}(s) ds, \quad (C12) $$

Here, $C_{1l0}^{(j)} - C_{4l0}^{(j)}$ are constants to be determined by boundary conditions.

We apply the boundary conditions at the bubble surfaces, which require that the normal velocity and the tangential stress of the Lagrangian streaming vanish at the equilibrium bubble surfaces,

$$ v_{Lr}(r_j, \theta_j) = 0 \quad \text{at} \quad r_j = R_{j0}, \quad (C13) $$

$$ \tau_L(r_j, \theta_j) = \eta \left[ \frac{1}{r_j} \frac{\partial v_{Lr}(r_j, \theta_j)}{\partial \theta_j} + \frac{\partial v_{L\theta}(r_j, \theta_j)}{\partial r_j} - \frac{v_{L\theta}(r_j, \theta_j)}{r_j} \right] = 0 \quad \text{at} \quad r_j = R_{j0}. \quad (C14) $$
Here, $v_{Lr}(r_j, \theta_j) = v_{Er}(r_j, \theta_j) + v_{Sr}(r_j, \theta_j)$ and $v_{L\theta}(r_j, \theta_j) = v_{E\theta}(r_j, \theta_j) + v_{S\theta}(r_j, \theta_j)$ are the components of the Lagrangian streaming velocity, $v_{Er}(r_j, \theta_j)$ and $v_{E\theta}(r_j, \theta_j)$ are the components of the Eulerian streaming velocity $v_E$, $v_{Sr}(r_j, \theta_j)$ and $v_{S\theta}(r_j, \theta_j)$ are the components of the Stokes drift velocity $v_S$ and $\tau_L(r_j, \theta_j)$ is the Lagrangian tangential stress.

The components of $v_S$ are calculated in Appendix E and given by (E15) and (E16). The components of $v_E$ are calculated by substituting (C6) into Eq. (2.8), which gives

$$v_{Er}(r_j, \theta_j) = -\frac{1}{r_j} \sum_{l=1}^{\infty} l(l+1) \Psi_l^{(j)}(r_j) P_l(\mu_j), \quad (C15)$$

$$v_{E\theta}(r_j, \theta_j) = -\frac{1}{r_j} \sum_{l=1}^{\infty} \left[ \Psi_l^{(j)}(r_j) + r_j \Psi_l^{(j)\prime}(r_j) \right] P_l^1(\mu_j), \quad (C16)$$

where, as follows from (C8) and the properties of the functions $C_{1l}^{(j)}(r_j) - C_{4l}^{(j)}(r_j)$, $\Psi_l^{(j)\prime}(r_j)$ is calculated by

$$\Psi_l^{(j)\prime}(r_j) = -\frac{(l-1)C_{1l}^{(j)}(r_j)}{r_j^l} - \frac{(l+1)C_{2l}^{(j)}(r_j)}{r_j^{l+2}} + l r_j^{l-1} C_{3l}^{(j)}(r_j) + (l+2) r_j^{l+1} C_{4l}^{(j)}(r_j). \quad (C17)$$

Substituting (C15), (C16), (E15) and (E16) into (C13) and (C14), and using (C8), (C17) and (C18),

$$\Psi_l^{(j)\prime\prime}(r_j) = \frac{l(l-1)C_{1l}^{(j)}(r_j)}{r_j^{l+1}} + \frac{(l+1)(l+2)C_{2l}^{(j)}(r_j)}{r_j^{l+3}} + l(l-1)r_j^{l-2}C_{3l}^{(j)}(r_j) + (l+1)(l+2)r_j^l C_{4l}^{(j)}(r_j), \quad (C18)$$

one obtains equations for $C_{1l0}^{(j)} - C_{4l0}^{(j)}$,

$$C_{1l0}^{(j)} + \frac{1}{R_{j0}^2} C_{2l0}^{(j)} + R_{j0}^{2l-1} C_{3l0}^{(j)} + R_{j0}^{2l+1} C_{4l0}^{(j)} = X_{1l}^{(j)}, \quad (C19)$$
$$ (l^2 - 1)C_{1l0}^{(j)} + \frac{l(l+2)}{R_{j0}^2}C_{2l0}^{(j)} + (l^2 - 1)R_{j0}^{2l-1}C_{3l0}^{(j)} + l(l+2)R_{j0}^{2l+1}C_{4l0}^{(j)} = X_{2l}^{(j)}, \quad (C20) $$

where

$$ X_{1l}^{(j)} = \frac{R_{j0}^{l}}{l(l+1)} V_{Srl}^{(j)}(R_{j0}), \quad (C21) $$

$$ X_{2l}^{(j)} = \frac{R_{j0}^{l}}{2} [V_{Srl}^{(j)}(R_{j0}) - V_{S\theta l}^{(j)}(R_{j0}) + R_{j0}V_{S\theta l}^{(j)\prime}(R_{j0})], \quad (C22) $$

Since $j = 1,2$, we have four equations. We need four more equations to find all the constants. The wanting equations can be obtained as follows.

Let us consider the solution in the coordinates $(r_1, \theta_1)$. Let us take a point with the coordinates $\theta_1 = 0$ and $r_1 = d - R_{20}$. This point is located on the surface of bubble 2 and hence $v_{Lr}(r_1, \theta_1)$ must vanish at this point. This requirement leads to the following equation:

$$ C_{1l0}^{(1)} + \frac{C_{2l0}^{(1)}}{(d - R_{20})^2} + (d - R_{20})^{2l-1} C_{3l0}^{(1)} + (d - R_{20})^{2l+1} C_{4l0}^{(1)} = X_{3l}^{(1)}, \quad (C23) $$

where

$$ X_{3l}^{(1)} = \frac{(d - R_{20})^l}{l(l+1)} V_{Srl}^{(1)}(d - R_{20}) + \frac{1}{2(2l-1)(2l+1)} \int_{R_{10}}^{d-R_{20}} \left[ (d - R_{20})^{2l-1} s^{3-l} - s^{l+2} \right] E_l^{(1)}(s) ds \\ + \frac{1}{2(2l+1)(2l+3)} \int_{R_{10}}^{d-R_{20}} \left[ \frac{s^{l+4}}{(d - R_{20})^2} - \frac{(d - R_{20})^{2l+1}}{s^{l-1}} \right] E_l^{(1)}(s) ds. \quad (C24) $$

For the solution in the coordinates $(r_2, \theta_2)$, we take a point with the coordinates $\theta_2 = \pi$ and $r_2 = d - R_{10}$, which is located on the surface of bubble 1. The requirement that $v_{Lr}(r_2, \theta_2) = 0$ at this points, yields

$$ C_{1l0}^{(2)} + \frac{C_{2l0}^{(2)}}{(d - R_{10})^2} + (d - R_{10})^{2l-1} C_{3l0}^{(2)} + (d - R_{10})^{2l+1} C_{4l0}^{(2)} = X_{3l}^{(2)}, \quad (C25) $$
where

$$X_{3l}^{(2)} = \frac{(d - R_{10})^l}{l(l+1)} V_{Srl}^{(2)}(d - R_{10}) + \frac{1}{2(2l-1)(2l+1)} \int_{R_{20}}^{d-R_{10}} \left[ (d - R_{10})^{2l-1} s^{3-l} - s^{l+2} \right] E_l^{(2)}(s) ds \\ + \frac{1}{2(2l+1)(2l+3)} \int_{R_{20}}^{d-R_{10}} \left[ \frac{s^{l+4}}{(d - R_{10})^2} - \frac{(d - R_{10})^{2l+1}}{s^{l-1}} \right] E_l^{(2)}(s) ds \quad . \quad (C26)$$

Let now us consider the Lagrangian tangential stress $\tau_L(r_j, \theta_j)$, which is defined by

$$\tau_L(r_j, \theta_j) = \eta \left[ \frac{\partial v_{Lr}(r_j, \theta_j)}{\partial \theta_j} + \frac{\partial v_{L\theta}(r_j, \theta_j)}{\partial r_j} - \frac{v_{L\theta}(r_j, \theta_j)}{r_j} \right] \\ = \eta \sum_{l=1}^{\infty} \left[ \frac{2(1-l^2)}{r_j^{l+1}} C_{1l}^{(j)}(r_j) - \frac{2l(l+2)}{r_j^{l+3}} C_{2l}^{(j)}(r_j) + 2(1-l^2)r_j^{l-2}C_{3l}^{(j)}(r_j) - 2l(l+2)r_j^l C_{4l}^{(j)}(r_j) \right. \\ \left. + \frac{V_{Srl}^{(j)}(r_j) - V_{S\theta l}^{(j)}(r_j)}{r_j} + V_{S\theta l}^{(j)\prime}(r_j) \right] P_l^1(\mu_j) \quad . \quad (C27)$$

As one can see, $\tau_L(r_1, \theta_1) = 0$ at $\theta_1 = 0$ as $P_l^1(\mu_1) = 0$ at $\theta_1 = 0$. However, $\tau_L(r_2, \theta_2)$ must also be equal to zero in a small neighborhood of the point with the coordinates $\theta_1 = 0$ and $r_1 = d - R_{20}$ as this area is located on the surface of bubble 2. To satisfy this condition, the expression in the square brackets in the sum of (C27) must vanish at $r_1 = d - R_{20}$. Similar considerations are also valid for $\tau_L(r_2, \theta_2)$ at the point $(\theta_2 = \pi, r_2 = d - R_{10})$. As a result, we obtain two last equations,

$$(l^2 - 1)C_{1l0}^{(j)} + \frac{l(l+2)}{(d - R_{(3-j)0})^2} C_{2l0}^{(j)} + (l^2 - 1)(d - R_{(3-j)0})^{2l-1} C_{3l0}^{(j)} + l(l+2)(d - R_{(3-j)0})^{2l+1} C_{4l0}^{(j)} = X_{4l0}^{(j)} \quad , \quad (C28)$$

where
$$
\begin{aligned}
X_{4l}^{(j)} = & \frac{(d - R_{(3-j)0})^l}{2} [V_{Srl}^{(j)}(d - R_{(3-j)0}) - V_{S\theta l}^{(j)}(d - R_{(3-j)0}) + (d - R_{(3-j)0})V_{S\theta l}^{(j)l}(d - R_{(3-j)0})] \\
& + \frac{l^2 - 1}{2(2l-1)(2l+1)} \int_{R_{j0}}^{d-R_{(3-j)0}} [(d - R_{(3-j)0})^{2l-1} s^{3-l} - s^{l+2}] E_l^{(j)}(s) ds \\
& + \frac{l(l+2)}{2(2l+1)(2l+3)} \int_{R_{j0}}^{d-R_{(3-j)0}} \left[ \frac{s^{l+4}}{(d - R_{(3-j)0})^2} - \frac{(d - R_{(3-j)0})^{2l+1}}{s^{l-1}} \right] E_l^{(j)}(s) ds.
\end{aligned}
\quad (C29)
$$

Combining (C19), (C20), (C23), (C25) and (C28), one finds

$$
C_{1l0}^{(j)} = \frac{(d - R_{(3-j)0})^{2l-1} [X_{2l}^{(j)} - l(l+2)X_{1l}^{(j)}] - R_{j0}^{2l-1} [X_{4l}^{(j)} - l(l+2)X_{3l}^{(j)}]}{(2l+1)[R_{j0}^{2l-1} - (d - R_{(3-j)0})^{2l-1}]}, \quad (C30)
$$

$$
C_{2l0}^{(j)} = \frac{R_{j0}^2 (d - R_{(3-j)0})^2 \{(d - R_{(3-j)0})^{2l+1} [(l^2 - 1)X_{1l}^{(j)} - X_{2l}^{(j)}] - R_{j0}^{2l+1} [(l^2 - 1)X_{3l}^{(j)} - X_{4l}^{(j)}]\}}{(2l+1)[R_{j0}^{2l+3} - (d - R_{(3-j)0})^{2l+3}]}, \quad (C31)
$$

$$
C_{3l0}^{(j)} = \frac{l(l+2)(X_{1l}^{(j)} - X_{3l}^{(j)}) - X_{2l}^{(j)} + X_{4l}^{(j)}}{(2l+1)[R_{j0}^{2l-1} - (d - R_{(3-j)0})^{2l-1}]}, \quad (C32)
$$

$$
C_{4l0}^{(j)} = \frac{R_{j0}^2 [X_{2l}^{(j)} - (l^2 - 1)X_{1l}^{(j)}] - (d - R_{(3-j)0})^2 [X_{4l}^{(j)} - (l^2 - 1)X_{3l}^{(j)}]}{(2l+1)[R_{j0}^{2l+3} - (d - R_{(3-j)0})^{2l+3}]}. \quad (C33)
$$

**Appendix D. Mathematical formulae used in calculations**

Here, auxiliary mathematical formulas are provided that are used in the principal
calculations. In equations that follow, the prime denotes the derivative with respect to an
argument in brackets, [ ] means the integer part of an expression in brackets, $C_{l_1 m_1 l_2 m_2}^{LM}$ are the
Clebsch-Gordan coefficients, and $\delta_{nm}$ is the Kronecker delta.

We use the following formulas for $P_n(\mu)$ and $P_n^1(\mu)$:
$$ \sqrt{1 - \mu^2} P_n^{1/}(\mu) = n(n+1)P_n(\mu) - \mu P_n'(\mu), \qquad (D1) $$

$$ \mu P_n'(\mu) = n P_n(\mu) + P_{n-1}(\mu), \qquad (D2) $$

$$ (1 - \mu^2) P_n'(\mu) = \frac{n(n+1)}{2n+1} [P_{n-1}(\mu) - P_{n+1}(\mu)], \qquad (D3) $$

$$ P_n'(\mu) = \sum_{k=1}^{[(n+1)/2]} (2n - 4k + 3) P_{n-2k+1}(\mu). \qquad (D4) $$

With the help of these formulas, the following identities are derived:

$$ \sqrt{1 - \mu^2} P_n^{1/}(\mu) P_m^{1/}(\mu) = n^2 P_n(\mu) P_m^{1/}(\mu) - \sum_{k=1}^{[n/2]} (2n - 4k + 1) P_{n-2k}(\mu) P_m^{1/}(\mu), \qquad (D5) $$

$$ P_n^1(\mu) P_m^1(\mu) = \frac{m(m+1)}{2m+1} \sum_{k=1}^{[(n+1)/2]} (2n-4k+3) [P_{m-1}(\mu) P_{n-2k+1}(\mu) - P_{m+1}(\mu) P_{n-2k+1}(\mu)] . \qquad (D6) $$

To calculate integrals with products of $P_n(\mu)$ and $P_m^1(\mu)$, we use the following identity:

$$ P_{l_1}^{m_1}(\mu) P_{l_2}^{m_2}(\mu) = \sqrt{\frac{(l_1+m_1)!(l_2+m_2)!}{(l_1-m_1)!(l_2-m_2)!}} \sum_{l=l_1-l_2}^{l_1+l_2} \sqrt{\frac{(l-m_1-m_2)!}{(l+m_1+m_2)!}} C_{l_1 0 l_2 0}^{l_1 0} C_{l_1 m_1 l_2 m_2}^{l(m_1+m_2)} P_l^{m_1+m_2}(\mu) . \qquad (D7) $$

This identity follows from equation (9) of § 5.6 of the book by Varshalovich *et al.* (1988).

From (D7) and the orthogonality condition for $P_n^1(\mu)$,

$$ \int_{-1}^{1} P_{n}^{1}(\mu) P_{m}^{1}(\mu) d\mu = \frac{2n(n+1)}{2n+1} \delta_{nm}, \qquad (D8) $$

one obtains the following integrals:

$$ \int_{-1}^{1} P_{n}(\mu) P_{m}^{1}(\mu) P_{l}^{1}(\mu) d\mu = 2\sqrt{l(l+1)} \sum_{\substack{m=n-l \\ m \ge 1}}^{n+l} \frac{\sqrt{m(m+1)} C_{n0l0}^{m0} C_{n0l1}^{m1}}{2m+1}, \qquad (D9) $$

$$ \int_{-1}^{1} P_{n-2k}(\mu) P_{m}^{1}(\mu) P_{l}^{1}(\mu) d\mu = 2\sqrt{l(l+1)} \sum_{\substack{m=n-2k-l \\ m \ge 1}}^{n-2k+l} \frac{\sqrt{m(m+1)} C_{(n-2k)0l0}^{m0} C_{(n-2k)0l1}^{m1}}{2m+1}. \qquad (D10) $$
From (D7) and the orthogonality condition for $P_n(\mu)$,

$$ \int_{-1}^{1} P_n(\mu) P_m(\mu) d\mu = \frac{2}{2n+1} \delta_{nm}, \qquad (D11) $$

one obtains

$$ \int_{-1}^{1} P_l(\mu) P_n(\mu) P_m(\mu) d\mu = \sum_{m=|l-n|}^{l+n} \frac{2}{2m+1} (C_{l0n0}^{m0})^2, \qquad (D12) $$

$$ \int_{-1}^{1} P_l(\mu) P_{m-1}(\mu) P_{n-2k+1}(\mu) d\mu = \sum_{m-1=|l-(n-2k+1)|}^{l+n-2k+1} \frac{2}{2m-1} (C_{l0(n-2k+1)0}^{(m-1)0})^2, \qquad (D13) $$

$$ \int_{-1}^{1} P_l(\mu) P_{m+1}(\mu) P_{n-2k+1}(\mu) d\mu = \sum_{m+1=|l-(n-2k+1)|}^{l+n-2k+1} \frac{2}{2m+3} (C_{l0(n-2k+1)0}^{(m+1)0})^2. \qquad (D14) $$

## Appendix E. Calculation of the Stokes drift velocity

The Stokes drift velocity is calculated by (Doinikov et al. 2019a)

$$ \mathbf{v}_S = \frac{1}{2\omega} \mathrm{Re} \{i(\mathbf{v} \cdot \nabla)\mathbf{v}^*\}. \qquad (E1) $$

Substitution of $\mathbf{v} = \mathbf{v}^{(1)} + \mathbf{v}^{(2)}$ into (E1) yields

$$ \mathbf{v}_S = \frac{1}{2\omega} \mathrm{Re} \left\{ i(\mathbf{v}^{(1)} \cdot \nabla)\mathbf{v}^{(1)*} + i(\mathbf{v}^{(2)} \cdot \nabla)\mathbf{v}^{(2)*} + i(\mathbf{v}^{(1)} \cdot \nabla)\mathbf{v}^{(2)*} + i(\mathbf{v}^{(2)} \cdot \nabla)\mathbf{v}^{(1)*} \right\}. \qquad (E2) $$

To represent $\mathbf{v}_S$ in the coordinates $(r_1, \theta_1)$, $\nabla$ is taken in the coordinates $(r_1, \theta_1)$, $\mathbf{v}^{(1)}$ is calculated by (2.6) and (2.7) and $\mathbf{v}^{(2)}$ is calculated by (A55) and (A56). As a result, we obtain the following equations for the components of $\mathbf{v}_S$:

$$ \mathbf{v}_{Sr}(r_1, \theta_1) = \frac{1}{2\omega} \mathrm{Re} \left\{ i(\mathbf{v}_r^{(1)} + \mathbf{v}_r^{(2)}) \frac{\partial(\mathbf{v}_r^{(1)*} + \mathbf{v}_r^{(2)*})}{\partial r_1} + \frac{i(\mathbf{v}_\theta^{(1)} + \mathbf{v}_\theta^{(2)})}{r_1} \frac{\partial(\mathbf{v}_r^{(1)*} + \mathbf{v}_r^{(2)*})}{\partial \theta_1} \right\} $$
$$ = \frac{1}{2\omega r_1^3} \left[ \sum_{n,m=0}^{\infty} S_{1nm}^{(1)}(r_1) P_n(\mu_1) P_m(\mu_1) + \sum_{n,m=1}^{\infty} S_{2nm}^{(1)}(r_1) P_n^1(\mu_1) P_m^1(\mu_1) \right], \quad (E3) $$

$$ \begin{aligned} v_{S\theta}(r_1, \theta_1) &= \frac{1}{2\omega} \operatorname{Re} \left\{ i(v_r^{(1)} + v_r^{(2)}) \frac{\partial(v_\theta^{(1)*} + v_\theta^{(2)*})}{\partial r_1} + \frac{i(v_\theta^{(1)} + v_\theta^{(2)})}{r_1} \frac{\partial(v_\theta^{(1)*} + v_\theta^{(2)*})}{\partial \theta_1} \right. \\ & \qquad \left. + \frac{i(v_r^{(1)*} + v_r^{(2)*})(v_\theta^{(1)} + v_\theta^{(2)})}{r_1} \right\} \\ &= \frac{1}{2\omega r_1^3} \left[ \sum_{\substack{n=0 \\ m=1}}^{\infty} Q_{1nm}^{(1)}(r_1) P_n(\mu_1) P_m^1(\mu_1) + \sqrt{1-\mu_1^2} \sum_{n,m=1}^{\infty} Q_{2nm}^{(1)}(r_1) P_n^1(\mu_1) P_m^{1/2}(\mu_1) \right], \quad (E4) \end{aligned} $$

where

$$ \begin{aligned} S_{1nm}^{(1)}(r_1) &= \operatorname{Im} \left\{ \left[ (n+1)a_n^{(1)} \left( \frac{R_{10}}{r_1} \right)^{n+1} + n(n+1)b_n^{(1)} h_n^{(1)}(k_v r_1) \right. \right. \\ & \qquad \left. \left. - n \left( \frac{r_1}{d} \right)^n \sum_{k=0}^{\infty} a_k^{(2)} (-1)^k C_{kn} \xi_2^{k+1} + n(n+1) j_n(k_v r_1) \sum_{k=1}^{\infty} b_k^{(2)} H_{nk} \right] \right. \\ & \qquad \left. \times \left[ (m+1)(m+2)a_m^{(1)*} \left( \frac{R_{10}}{r_1} \right)^{m+1} + m(m+1)b_m^{(1)*} \left[ h_m^{(1)}(k_v r_1) - k_v r_1 h_m^{(1)/(2)}(k_v r_1) \right]^* \right. \right. \\ & \qquad \left. \left. + m(m-1) \left( \frac{r_1}{d} \right)^m \sum_{k=0}^{\infty} a_k^{(2)*} (-1)^k C_{km} \xi_2^{k+1} + m(m+1) \left[ j_m(k_v r_1) - k_v r_1 j_m'(k_v r_1) \right]^* \sum_{k=1}^{\infty} (b_k^{(2)} H_{mk})^* \right] \right\}, \quad (E5) \\ & n, m \ge 0, \end{aligned} $$

$$ S_{2nm}^{(1)}(r_1) = \operatorname{Im} \left\{ \left[ (n+1)a_n^{(1)*} \left( \frac{R_{10}}{r_1} \right)^{n+1} + n(n+1)b_n^{(1)*} h_n^{(1)*}(k_v r_1) \right] \right\} $$
$$
\begin{aligned}
& -n \left(\frac{r_1}{d}\right)^n \sum_{k=0}^{\infty} a_k^{(2)*} (-1)^k C_{kn} \xi_2^{k+1} + n(n+1) j_n (k_v^* r_1) \sum_{k=1}^{\infty} (b_k^{(2)} H_{nk})^* \\
& \times \left[ a_m^{(1)} \left(\frac{R_{10}}{r_1}\right)^{m+1} - b_m^{(1)} \left[ h_m^{(1)}(k_v r_1) + k_v r_1 h_m^{(1)\prime}(k_v r_1) \right] \right. \\
& \left. + \left(\frac{r_1}{d}\right)^m \sum_{k=0}^{\infty} a_k^{(2)} (-1)^k C_{km} \xi_2^{k+1} - \left[ j_m(k_v r_1) + k_v r_1 j_m'(k_v r_1) \right] \sum_{k=1}^{\infty} b_k^{(2)} H_{mk} \right] \Biggr\} , \quad n, m \ge 1
\end{aligned}
\quad (E6)
$$

$$
\begin{aligned}
Q_{1nm}^{(1)}(r_1) = \operatorname{Im} \left\{ \left[ (n+1)a_n^{(1)}\left(\frac{R_{10}}{r_1}\right)^{n+1} + n(n+1)b_n^{(1)}h_n^{(1)}(k_v r_1) \right. \right. \\
\left. - n\left(\frac{r_1}{d}\right)^n \sum_{k=0}^{\infty} a_k^{(2)}(-1)^k C_{kn}\xi_2^{k+1} + n(n+1)j_n(k_v r_1)\sum_{k=1}^{\infty} b_k^{(2)}H_{nk} \right] \\
\left. \times \left[ -(m+3)a_m^{(1)*}\left(\frac{R_{10}}{r_1}\right)^{m+1} + b_m^{(1)*}\left[2h_m^{(1)}(k_v r_1) - (k_v r_1)^2 h_m^{(1)\prime\prime}(k_v r_1)\right]^* \right. \right. \\
\left. + (m-2)\left(\frac{r_1}{d}\right)^m \sum_{k=0}^{\infty} a_k^{(2)*}(-1)^k C_{km}\xi_2^{k+1} + \left[2j_m(k_v r_1) - (k_v r_1)^2 j_m''(k_v r_1)\right]^* \sum_{k=1}^{\infty} (b_k^{(2)}H_{mk})^* \right] \Biggr\} , \quad n \ge 0, m \ge 1 \\
\left. \right\}
\end{aligned}
\quad (E7)
$$

$$
\begin{aligned}
Q_{2nm}^{(1)}(r_1) = \operatorname{Im} \left\{ \left\{ a_n^{(1)} \left( \frac{R_{10}}{r_1} \right)^{n+1} - b_n^{(1)} \left[ h_n^{(1)}(k_v r_1) + k_v r_1 h_n^{(1)\prime}(k_v r_1) \right] \right. \right. \\
\left. + \left( \frac{r_1}{d} \right)^n \sum_{k=0}^{\infty} a_k^{(2)} (-1)^k C_{kn} \xi_2^{k+1} - \left[ j_n(k_v r_1) + k_v r_1 j_n'(k_v r_1) \right] \sum_{k=1}^{\infty} b_k^{(2)} H_{nk} \right\} \\
\left. \times \left\{ a_m^{(1)*} \left( \frac{R_{10}}{r_1} \right)^{m+1} - b_m^{(1)*} \left[ h_m^{(1)}(k_v r_1) + k_v r_1 h_m^{(1)\prime}(k_v r_1) \right]^* \right. \right.
\end{aligned}
$$
$$ + \left(\frac{r_1}{d}\right)^m \sum_{k=0}^{\infty} a_k^{(2)*} (-1)^k C_{km} \xi_2^{k+1} - \left[ j_m(k_v r_1) + k_v r_1 j_m(k_v r_1) \right]^* \sum_{k=1}^{\infty} (b_k^{(2)} H_{mk})^* \right\} , \quad n,m \ge 1 \quad (E8) $$

To represent (E2) in the coordinates $(r_2, \theta_2)$, $\nabla$ is taken in the coordinates $(r_2, \theta_2)$, $v^{(1)}$ is calculated by (A52) and (A53) and $v^{(2)}$ is calculated by (2.6) and (2.7). As a result, we obtain

$$ \begin{aligned} v_{Sr}(r_2, \theta_2) &= \frac{1}{2\omega} \operatorname{Re} \left\{ i(v_r^{(1)} + v_r^{(2)}) \frac{\partial(v_r^{(1)*} + v_r^{(2)*})}{\partial r_2} + \frac{i(v_\theta^{(1)} + v_\theta^{(2)})}{r_2} \frac{\partial(v_r^{(1)*} + v_r^{(2)*})}{\partial \theta_2} \right\} \\ &= \frac{1}{2\omega r_2^3} \left[ \sum_{n,m=0}^{\infty} S_{1nm}^{(2)}(r_2) P_n(\mu_2) P_m(\mu_2) + \sum_{n,m=1}^{\infty} S_{2nm}^{(2)}(r_2) P_n^1(\mu_2) P_m^1(\mu_2) \right], \end{aligned} \quad (E9) $$

$$ \begin{aligned} v_{S\theta}(r_2, \theta_2) &= \frac{1}{2\omega} \operatorname{Re} \left\{ i(v_r^{(1)} + v_r^{(2)}) \frac{\partial(v_\theta^{(1)*} + v_\theta^{(2)*})}{\partial r_2} + \frac{i(v_\theta^{(1)} + v_\theta^{(2)})}{r_2} \frac{\partial(v_\theta^{(1)*} + v_\theta^{(2)*})}{\partial \theta_2} \right. \\ & \qquad \left. + \frac{i(v_r^{(1)*} + v_r^{(2)*})(v_\theta^{(1)} + v_\theta^{(2)})}{r_2} \right\} \\ &= \frac{1}{2\omega r_2^3} \left[ \sum_{\substack{n=0 \\ m=1}}^{\infty} Q_{1nm}^{(2)}(r_2) P_n(\mu_2) P_m^1(\mu_2) + \sqrt{1-\mu_2^2} \sum_{n,m=1}^{\infty} Q_{2nm}^{(2)}(r_2) P_n^1(\mu_2) P_m^{1/2}(\mu_2) \right], \end{aligned} \quad (E10) $$

where

$$ \begin{aligned} S_{1nm}^{(2)}(r_2) = \operatorname{Im} \left\{ \left[ (n+1)a_n^{(2)} \left( \frac{R_{20}}{r_2} \right)^{n+1} + n(n+1)b_n^{(2)} h_n^{(1)}(k_v r_2) \right. \right. \\ \left. \left. \qquad -(-1)^n n \left( \frac{r_2}{d} \right)^n \sum_{k=0}^{\infty} a_k^{(1)} C_{kn} \xi_1^{k+1} + n(n+1) j_n(k_v r_2) \sum_{k=1}^{\infty} b_k^{(1)} G_{nk} \right] \right. \\ \left. \qquad \times \left[ (m+1)(m+2) a_m^{(2)*} \left( \frac{R_{20}}{r_2} \right)^{m+1} + m(m+1) b_m^{(2)*} \left[ h_m^{(1)}(k_v r_2) - k_v r_2 h_m^{(1/2)}(k_v r_2) \right] \right]^* \right\} \end{aligned} $$
$$ +(-1)^m m(m-1) \left(\frac{r_2}{d}\right)^m \sum_{k=0}^{\infty} a_k^{(1)*} C_{km} \xi_1^{k+1} + m(m+1) \left[j_m(k_v r_2) - k_v r_2 j_m'(k_v r_2)\right]^* \sum_{k=1}^{\infty} (b_k^{(1)} G_{mk})^* \right] , \qquad (E11) $$

$$ S_{2nm}^{(2)}(r_2) = \operatorname{Im} \left\{ \left[ (n+1)a_n^{(2)*} \left( \frac{R_{20}}{r_2} \right)^{n+1} + n(n+1)b_n^{(2)*} h_n^{(1)*} (k_v r_2) \right. \right. \\ \left. \left. -(-1)^n n \left( \frac{r_2}{d} \right)^n \sum_{k=0}^{\infty} a_k^{(1)*} C_{kn} \xi_1^{k+1} + n(n+1) j_n(k_v r_2) \sum_{k=1}^{\infty} (b_k^{(1)} G_{nk})^* \right] \right. \\ \left. \times \left[ a_m^{(2)} \left( \frac{R_{20}}{r_2} \right)^{m+1} - b_m^{(2)} \left[ h_m^{(1)}(k_v r_2) + k_v r_2 h_m^{(1)\prime} (k_v r_2) \right] \right. \right. \\ \left. \left. +(-1)^m \left( \frac{r_2}{d} \right)^m \sum_{k=0}^{\infty} a_k^{(1)} C_{km} \xi_1^{k+1} - \left[ j_m(k_v r_2) + k_v r_2 j_m'(k_v r_2) \right] \sum_{k=1}^{\infty} b_k^{(1)} G_{mk} \right] \right\}, \quad n, m \ge 1, \qquad (E12) $$

$$ Q_{1nm}^{(2)}(r_2) = \operatorname{Im} \left\{ \left[ (n+1)a_n^{(2)} \left( \frac{R_{20}}{r_2} \right)^{n+1} + n(n+1)b_n^{(2)} h_n^{(1)}(k_v r_2) \right. \right. \\ \left. \left. -(-1)^n n \left( \frac{r_2}{d} \right)^n \sum_{k=0}^{\infty} a_k^{(1)} C_{kn} \xi_1^{k+1} + n(n+1) j_n(k_v r_2) \sum_{k=1}^{\infty} b_k^{(1)} G_{nk} \right] \right. \\ \left. \times \left[ -(m+3)a_m^{(2)*} \left( \frac{R_{20}}{r_2} \right)^{m+1} + b_m^{(2)*} \left[ 2h_m^{(1)}(k_v r_2) - (k_v r_2)^2 h_m^{(1)\prime\prime} (k_v r_2) \right]^* \right. \right. \\ \left. \left. +(-1)^m (m-2) \left( \frac{r_2}{d} \right)^m \sum_{k=0}^{\infty} a_k^{(1)*} C_{km} \xi_1^{k+1} + \left[ 2j_m(k_v r_2) - (k_v r_2)^2 j_m''(k_v r_2) \right]^* \sum_{k=1}^{\infty} (b_k^{(1)} G_{mk})^* \right] \right\}, \quad n \ge 0, m \ge 1, \qquad (E13) $$
$$
\begin{aligned}
Q_{2nm}^{(2)}(r_2) = \operatorname{Im} & \left\{ \left\{ a_n^{(2)} \left( \frac{R_{20}}{r_2} \right)^{n+1} - b_n^{(2)} \left[ h_n^{(1)}(k_v r_2) + k_v r_2 h_n^{(1)/(k_v r_2)} \right] \right. \right. \\
& \left. \left. + (-1)^n \left( \frac{r_2}{d} \right)^n \sum_{k=0}^{\infty} a_k^{(1)} C_{kn} \xi_1^{k+1} - \left[ j_n(k_v r_2) + k_v r_2 j'_n(k_v r_2) \right] \sum_{k=1}^{\infty} b_k^{(1)} G_{nk} \right\} \\
& \left. \times \left\{ a_m^{(2)*} \left( \frac{R_{20}}{r_2} \right)^{m+1} - b_m^{(2)*} \left[ h_m^{(1)}(k_v r_2) + k_v r_2 h_m^{(1)/(k_v r_2)} \right]^* \right. \right. \\
& \left. \left. + (-1)^m \left( \frac{r_2}{d} \right)^m \sum_{k=0}^{\infty} a_k^{(1)*} C_{km} \xi_1^{k+1} - \left[ j_m(k_v r_2) + k_v r_2 j'_m(k_v r_2) \right]^* \sum_{k=1}^{\infty} (b_k^{(1)} G_{mk})^* \right\} \right\} , \quad n,m \ge 1 \\
& \qquad , \quad (E14)
\end{aligned}
$$

Equations (E3), (E4), (E9) and (E10) can be represented as

$$
v_{Sr}(r_j, \theta_j) = \sum_{l=0}^{\infty} V_{Srl}^{(j)}(r_j) P_l(\mu_j), \quad (E15)
$$

$$
v_{S\theta}(r_j, \theta_j) = \sum_{l=1}^{\infty} V_{S\theta l}^{(j)}(r_j) P_l^1(\mu_j), \quad (E16)
$$

where $V_{Srl}^{(j)}(r_j)$ is calculated by

$$
V_{Srl}^{(j)}(r_j) = \frac{2l+1}{4\omega r_j^3} \int_{-1}^{1} \left[ \sum_{n,m=0}^{\infty} S_{1nm}^{(j)}(r_j) P_n(\mu_j) P_m(\mu_j) + \sum_{n,m=1}^{\infty} S_{2nm}^{(j)}(r_j) P_n^1(\mu_j) P_m^1(\mu_j) \right] P_l(\mu_j) d\mu_j, \quad (E17)
$$

while $V_{S\theta l}^{(j)}(r_j)$ is calculated by

$$
V_{S\theta l}^{(j)}(r_j) = \frac{2l+1}{4l(l+1)\omega r_j^3} \\
\times \int_{-1}^{1} \left\{ \sum_{n=0}^{\infty} Q_{1nm}^{(j)}(r_j) P_n(\mu_j) P_m^1(\mu_j) + \sqrt{1-\mu_j^2} \sum_{n,m=1}^{\infty} Q_{2nm}^{(j)}(r_j) P_n^1(\mu_j) P_m^{1/(\mu_j)} \right\} P_l^1(\mu_j) d\mu_j . \quad (E18)
$$
With the help of (D5), (D6), (D9), (D10), (D12) – (D14), the integration of (E17) and (E18) yields final expressions for $V_{Srl}^{(j)}(r_j)$ and $V_{Sθl}^{(j)}(r_j)$,

$$ V_{Srl}^{(j)}(r_j) = \frac{2l+1}{2\omega r_j^3} \left\{ \sum_{n=0}^{\infty} \sum_{m=|l-n|}^{l+n} \frac{C_{l0n0}^{m0}}{2m+1} S_{1nm}^{(j)}(r_j) \right. \\ \left. + \sum_{n=1}^{\infty} \sum_{m=1}^{[(n+1)/2]} \sum_{k=|l-(n-2m+1)|}^{l+n-2m+1} \frac{(2n-4m+3)(C_{l0(n-2m+1)0}^{k0})^2}{2k+1} \left[ \frac{(k+1)(k+2)}{2k+3} S_{2n(k+1)}^{(j)}(r_j) - \frac{k(k-1)}{2k-1} S_{2n(k-1)}^{(j)}(r_j) \right] \right\}, \quad (E19) $$

$$ V_{S\theta l}^{(j)}(r_j) = \frac{2l+1}{2\sqrt{l(l+1)\omega r_j^3}} \left\{ \sum_{n=0}^{\infty} \sum_{\substack{m=|n-l| \\ m \ge 1}}^{n+l} \frac{\sqrt{m(m+1)} C_{n0l0}^{m0} C_{n0l1}^{m1}}{2m+1} [Q_{1nm}^{(j)}(r_j) - n^2 Q_{2nm}^{(j)}(r_j)] \right. \\ \left. + \sum_{n=2}^{\infty} \sum_{\substack{m=1 \\ k \ge 1}}^{[n/2]} \sum_{\substack{k=n-2m-l \\ k \ge 1}}^{n-2m+l} \frac{(2n-4m+1)\sqrt{k(k+1)} C_{(n-2m)0l0}^{k0} C_{(n-2m)0l1}^{k1}}{2k+1} Q_{2nk}^{(j)}(r_j) \right\}. \quad (E20) $$

When calculating the Lagrangian streaming velocity, the term with $l=0$ in (E15) should be omitted as this term is a function of $r_j$ alone and hence does not contribute to acoustic streaming.

In our calculations, the derivative $V_{S\theta l}^{(j)\prime}(r_j)$ is used; see (C22) and (C29). It is calculated by

$$ V_{S\theta l}^{(j)\prime}(r_j) = -\frac{3V_{S\theta l}^{(j)}(r_j)}{r_j} + \frac{2l+1}{2\sqrt{l(l+1)\omega r_j^3}} \left\{ \sum_{n=0}^{\infty} \sum_{m=|n-l|}^{n+l} \frac{\sqrt{m(m+1)} C_{n0l0}^{m0} C_{n0l1}^{m1}}{2m+1} [Q_{1nm}^{(j)\prime}(r_j) - n^2 Q_{2nm}^{(j)\prime}(r_j)] \right. \\ \left. + \sum_{n=2}^{\infty} \sum_{\substack{m=1 \\ k \ge 1}}^{[n/2]} \sum_{\substack{k=n-2m-l \\ k \ge 1}}^{n-2m+l} \frac{(2n-4m+1)\sqrt{k(k+1)} C_{(n-2m)0l0}^{k0} C_{(n-2m)0l1}^{k1}}{2k+1} Q_{2nk}^{(j)\prime}(r_j) \right\}, \quad (E21) $$

where
$$
\begin{aligned}
Q_{1nm}^{(1)\prime}(r_1) = \operatorname{Im} & \left\{ \left[ -\frac{(n+1)^2}{R_{10}} a_n^{(1)} \left( \frac{R_{10}}{r_1} \right)^{n+2} + n(n+1)b_n^{(1)} k_v h_n^{(1)\prime} (k_v r_1) \right. \right. \\
& \left. -\frac{n^2}{d} \left(\frac{r_1}{d}\right)^{n-1} \sum_{k=0}^{\infty} a_k^{(2)} (-1)^k C_{km} \xi_2^{k+1} + n(n+1)k_v j_n'(k_v r_1) \sum_{k=1}^{\infty} b_k^{(2)} H_{nk} \right] \\
& \left. \times \left[ -(m+3)a_m^{(1)*} \left( \frac{R_{10}}{r_1} \right)^{m+1} + b_m^{(1)*} \left[ 2h_m^{(1)}(k_v r_1) - (k_v r_1)^2 h_m^{(1)\prime\prime} (k_v r_1) \right]^* \right. \right. \\
& \left. + (m-2) \left(\frac{r_1}{d}\right)^m \sum_{k=0}^{\infty} a_k^{(2)*} (-1)^k C_{km} \xi_2^{k+1} + \left[ 2j_m(k_v r_1) - (k_v r_1)^2 j_m^{\prime\prime} (k_v r_1) \right]^* \sum_{k=1}^{\infty} (b_k^{(2)} H_{mk})^* \right] \\
& \left. + \left[ (n+1)a_n^{(1)} \left( \frac{R_{10}}{r_1} \right)^{n+1} + n(n+1)b_n^{(1)} h_n^{(1)} (k_v r_1) \right. \right. \\
& \left. - n \left(\frac{r_1}{d}\right)^n \sum_{k=0}^{\infty} a_k^{(2)} (-1)^k C_{kn} \xi_2^{k+1} + n(n+1)j_n(k_v r_1) \sum_{k=1}^{\infty} b_k^{(2)} H_{nk} \right] \right\} \\
& \times \left[ \frac{(m+1)(m+3)}{R_{10}} a_m^{(1)*} \left( \frac{R_{10}}{r_1} \right)^{m+2} + b_m^{(1)*} k_v^* \left[ 2h_m^{(1)\prime}(k_v r_1) - 2k_v r_1 h_m^{(1)\prime\prime}(k_v r_1) - (k_v r_1)^2 h_m^{(1)\prime\prime\prime}(k_v r_1) \right]^* \right. \\
& \left. + \frac{m(m-2)}{d} \left(\frac{r_1}{d}\right)^{m-1} \sum_{k=0}^{\infty} a_k^{(2)*} (-1)^k C_{km} \xi_2^{k+1} , \right. \\
& \left. + k_v^* \left[ 2j_m'(k_v r_1) - 2k_v r_1 j_m''(k_v r_1) - (k_v r_1)^2 j_m'''(k_v r_1) \right]^* \sum_{k=1}^{\infty} (b_k^{(2)} H_{mk})^* \right] \right\}, \quad n \ge 0, m \ge 1, \quad (E22)
\end{aligned}
$$

$$
Q_{2nm}^{(1)\prime}(r_1) = \operatorname{Im} \left\{ \left\{ -\frac{n+1}{R_{10}} a_n^{(1)} \left( \frac{R_{10}}{r_1} \right)^{n+2} - b_n^{(1)} k_v \left[ 2h_n^{(1)\prime}(k_v r_1) + k_v r_1 h_n^{(1)\prime\prime}(k_v r_1) \right] \right\} \right\}
$$
$$
\begin{aligned}
& + \frac{n}{d} \left(\frac{r_1}{d}\right)^{n-1} \sum_{k=0}^{\infty} a_k^{(2)} (-1)^k C_{km} \xi_2^{k+1} - k_v \left[2j'_n(k_v r_1) + k_v r_1 j''_n(k_v r_1)\right] \sum_{k=1}^{\infty} b_k^{(2)} H_{nk} \} \\
& \times \left\{ a_m^{(1)*} \left(\frac{R_{10}}{r_1}\right)^{m+1} - b_m^{(1)*} \left[h_m^{(1)}(k_v r_1) + k_v r_1 h_m^{(1)/(k_v r_1)}\right]^* \right. \\
& \left. + \left(\frac{r_1}{d}\right)^m \sum_{k=0}^{\infty} a_k^{(2)*} (-1)^k C_{km} \xi_2^{k+1} - \left[j_m(k_v r_1) + k_v r_1 j'_m(k_v r_1)\right]^* \sum_{k=1}^{\infty} (b_k^{(2)} H_{mk})^* \right\} \\
& + \left\{ a_n^{(1)} \left(\frac{R_{10}}{r_1}\right)^{n+1} - b_n^{(1)} \left[h_n^{(1)}(k_v r_1) + k_v r_1 h_n^{(1)/(k_v r_1)}\right] \right. \\
& \left. + \left(\frac{r_1}{d}\right)^n \sum_{k=0}^{\infty} a_k^{(2)} (-1)^k C_{kn} \xi_2^{k+1} - \left[j_n(k_v r_1) + k_v r_1 j'_n(k_v r_1)\right] \sum_{k=1}^{\infty} b_k^{(2)} H_{nk} \right\} \\
& \times \left\{ -\frac{m+1}{R_{10}} a_m^{(1)*} \left(\frac{R_{10}}{r_1}\right)^{m+2} - b_m^{(1)*} k_v^* \left[2h_m^{(1)/(k_v r_1)} + k_v r_1 h_m^{(1)\prime/(k_v r_1)}\right]^* \right. \\
& \left. + \frac{m}{d} \left(\frac{r_1}{d}\right)^{m-1} \sum_{k=0}^{\infty} a_k^{(2)*} (-1)^k C_{km} \xi_2^{k+1} \right. \\
& \left. - k_v^* \left[2j_m(k_v r_1) + k_v r_1 j''_m(k_v r_1)\right]^* \sum_{k=1}^{\infty} (b_k^{(2)} H_{mk})^* \right\}, \quad n, m \ge 1,
\end{aligned}
\quad (E23)
$$

$$
\begin{aligned}
Q_{1nm}^{(2)\prime}(r_2) = \operatorname{Im} \left\{ \left[ -\frac{(n+1)^2}{R_{20}} a_n^{(2)} \left(\frac{R_{20}}{r_2}\right)^{n+2} + n(n+1)b_n^{(2)} k_v h_n^{(1)\prime} (k_v r_2) \right. \right. \\
\left. \left. - \frac{(-1)^n n^2}{d} \left(\frac{r_2}{d}\right)^{n-1} \sum_{k=0}^{\infty} a_k^{(1)} C_{kn} \xi_1^{k+1} + n(n+1)k_v j'_n(k_v r_2) \sum_{k=1}^{\infty} b_k^{(1)} G_{nk} \right] \right\}
\end{aligned}
$$
$$
\begin{aligned}
& \times \left[ -(m+3)a_m^{(2)*} \left( \frac{R_{20}}{r_2} \right)^{m+1} + b_m^{(2)*} \left[ 2h_m^{(1)}(k_v r_2) - (k_v r_2)^2 h_m^{(1)\parallel}(k_v r_2) \right]^* \right. \\
& \left. \qquad + (-1)^m (m-2) \left( \frac{r_2}{d} \right)^m \sum_{k=0}^{\infty} a_k^{(1)*} C_{km} \xi_1^{k+1} + \left[ 2j_m(k_v r_2) - (k_v r_2)^2 j_m^{\parallel}(k_v r_2) \right]^* \sum_{k=1}^{\infty} (b_k^{(1)} G_{mk})^* \right] \\
& \qquad + \left[ (n+1)a_n^{(2)} \left( \frac{R_{20}}{r_2} \right)^{n+1} + n(n+1)b_n^{(2)} h_n^{(1)}(k_v r_2) \right. \\
& \qquad \left. \qquad - (-1)^n n \left( \frac{r_2}{d} \right)^n \sum_{k=0}^{\infty} a_k^{(1)} C_{kn} \xi_1^{k+1} + n(n+1) j_n(k_v r_2) \sum_{k=1}^{\infty} b_k^{(1)} G_{nk} \right] \\
& \times \left[ \frac{(m+1)(m+3)}{R_{20}} a_m^{(2)*} \left( \frac{R_{20}}{r_2} \right)^{m+2} + b_m^{(2)*} k_v^* \left[ 2h_m^{(1)\parallel}(k_v r_2) - 2k_v r_2 h_m^{(1)\parallel}(k_v r_2) - (k_v r_2)^2 h_m^{(1)\parallel\parallel}(k_v r_2) \right]^* \right. \\
& \qquad \left. \qquad + \frac{(-1)^m m(m-2)}{d} \left( \frac{r_2}{d} \right)^{m-1} \sum_{k=0}^{\infty} a_k^{(1)*} C_{km} \xi_1^{k+1}, \right. \\
& \qquad \left. \qquad + k_v^* \left[ 2j_m'(k_v r_2) - 2k_v r_2 j_m''(k_v r_2) - (k_v r_2)^2 j_m'''(k_v r_2) \right]^* \sum_{k=1}^{\infty} (b_k^{(1)} G_{mk})^* \right] \right\},
\end{aligned}
\quad (E24)
$$

$$
\begin{aligned}
Q_{2nm}^{(2)\prime}(r_2) = \operatorname{Im} & \left\{ \left\{ -\frac{n+1}{R_{20}} a_n^{(2)} \left( \frac{R_{20}}{r_2} \right)^{n+2} - b_n^{(2)} k_v \left[ 2h_n^{(1)\prime}(k_v r_2) + k_v r_2 h_n^{(1)\prime\prime}(k_v r_2) \right] \right. \right. \\
& \left. \left. + \frac{(-1)^n n}{d} \left( \frac{r_2}{d} \right)^{n-1} \sum_{k=0}^{\infty} a_k^{(1)} C_{kn} \xi_1^{k+1} - k_v \left[ 2j_n'(k_v r_2) + k_v r_2 j_n''(k_v r_2) \right] \sum_{k=1}^{\infty} b_k^{(1)} G_{nk} \right\} \right\} \\
& \times \left\{ a_m^{(2)*} \left( \frac{R_{20}}{r_2} \right)^{m+1} - b_m^{(2)*} \left[ h_m^{(1)}(k_v r_2) + k_v r_2 h_m^{(1)\prime}(k_v r_2) \right]^* \right\}
\end{aligned}
$$
$$
\begin{equation}
\begin{aligned}
& \left( -(-1)^m \left( \frac{r_2}{d} \right)^m \sum_{k=0}^{\infty} a_k^{(1)*} C_{km} \xi_1^{k+1} - \left[ j_m(k_v r_2) + k_v r_2 j_m'(k_v r_2) \right]^* \sum_{k=1}^{\infty} (b_k^{(1)} G_{mk})^* \right) \\
& + \left\{ a_n^{(2)} \left( \frac{R_{20}}{r_2} \right)^{n+1} - b_n^{(2)} \left[ h_n^{(1)}(k_v r_2) + k_v r_2 h_n^{(1)\prime}(k_v r_2) \right] \right. \\
& \left. + (-1)^n \left( \frac{r_2}{d} \right)^n \sum_{k=0}^{\infty} a_k^{(1)} C_{kn} \xi_1^{k+1} - \left[ j_n(k_v r_2) + k_v r_2 j_n'(k_v r_2) \right] \sum_{k=1}^{\infty} b_k^{(1)} G_{nk} \right\} \\
& \times \left\{ -\frac{m+1}{R_{20}} a_m^{(2)*} \left( \frac{R_{20}}{r_2} \right)^{m+2} - b_m^{(2)*} k_v^* \left[ 2h_m^{(1)\prime}(k_v r_2) + k_v r_2 h_m^{(1)\prime\prime}(k_v r_2) \right]^* \right. \\
& \left. + \frac{(-1)^m m}{d} \left( \frac{r_2}{d} \right)^{m-1} \sum_{k=0}^{\infty} a_k^{(1)*} C_{km} \xi_1^{k+1} \right. \\
& \left. - k_v^* \left[ 2j_m'(k_v r_2) + k_v r_2 j_m''(k_v r_2) \right]^* \sum_{k=1}^{\infty} (b_k^{(1)} G_{mk})^* \right\}, \quad n,m \ge 1
\end{aligned}
\tag{E25}
\end{equation}
$$

Appendix F. Geometry and mesh of StarCCM+ simulations

The computational domain (figure 14) consists of a quarter of a millimetre sphere (liquid medium) containing two quarters of sub-millimetre spheres (i.e. two bubbles). The axisymmetric assumption is used to model only a quarter of the entire fluid problem. The axes of each quarter of the sphere are aligned together. The mesh size is sampled from 1 µm at the bubble interface to 2 µm at the edge of a sphere of radius 500 µm around the bubble centres. The initial condition is chosen in such a way that each bubble interface is deformed as follows:

$$
r_s^{(j)} = R_{j0} + s_0^{(j)} \cos(2\pi f_0 t) + s_2^{(j)} P_2(\cos\theta) \cos(2\pi f_0 t) \quad (F1)
$$

With $R_{j0} = 50\mu m$, $s_0^{(j)} = 10\mu m$, $s_2^{(j)} = 10\mu m$ and $f_0 = 30kHz$. The initially deformed bubble interface is visible in figure 14(b).
<figure><img src="image_18.png" /><figcaption>Figure 14. (a) View of the geometry. A quarter of fluid volume containing two bubbles at the separation distance *d* is considered. (b) View of the mesh resolution in a plane containing the bubble centres. The inlet shows the mesh sizes in the vicinity of the bubble interface.</figcaption></figure>

REFERENCES

Abramowitz, M. & Stegun, I. A. 1972 *Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables*. US Department of Commerce.

Bertin, N., Spelman, T. A., Combriat, T., Hue, H., Stéphan, O., Lauga, E. & Marmottant, P. 2017
Bubble-based acoustic micropropulsors: active surface and mixers. *Lab on a Chip* **17**,
1515–1528.

Bolanos-Jimenez, R., Rossi, M., Fernandez Rivas, D., Kähler, C. J. & Marin, A. 2017
Bubble-based acoustic micropropulsors: active surface and mixers. *J. Fluid Mech.* **820**,
529–548.

Boyce, W. E. & DiPrima, R. C. 2001 *Elementary Differential Equations and Boundary Value Problems*. Wiley.

Doinikov, A. A. & Bouakaz, A. 2015 Theoretical model for coupled radial and translational motion of two bubbles at arbitrary separation distances. *Phys. Rev. E* **92**, 043001.

Doinikov, A. A. & Bouakaz, A. 2016 Microstreaming generated by two acoustically induced gas bubbles. *J. Fluid Mech.* **796**, 318–339.

Doinikov, A. A., Cleve, S., Regnault, G., Mauger, C. & Inserra, C. 2019a Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. I. Case of modes 0 and m. *Phys. Rev. E* **100**, 033104.
Doinikov, A. A., Cleve, S., Regnault, G., Mauger, C. & Inserra, C. 2019b Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. II. Case of modes 1 and *m*. *Phys. Rev. E* **100**, 033105.

Francescutto, A. & Nabergoj R. 1978 Pulsation amplitude threshold for surface waves on oscillating bubbles. *Acustica* **41**, 215–220.

Garbin, V., Cojoc, D., Ferrari, E., Di Fabrizio, D., Overvelde, M. L. J., van der Meer, S. M., de Jong, N., Lohse, D. & Versluis, M. 2007 Changes in microbubble dynamics near a boundary revealed by combined optical micromanipulation and high-speed imaging. *Appl. Phys. Lett.* **90**, 114103.

Guédra, M. & Inserra, C. 2018 Bubble shape oscillations of finite amplitude. *J. Fluid Mech.* **857**, 681–703.

Hall, P. & Seminara, G. 1980 Nonlinear oscillations of non-spherical cavitation bubbles in acoustic fields. *J. Fluid Mech.* **101**, 423–444.

Inserra, C., Regnault, G., Cleve, S., Mauger, C. & Doinikov, A. A. 2020a Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. III. Case of self-interacting modes *n-n*. *Phys. Rev. E* **101**, 013111.

Inserra, C., Regnault, G., Cleve, S., Mauger, C. & Doinikov, A. A. 2020b Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. IV. Case of modes *n* and *m*. *Phys. Rev. E* **102**, 043103.

Longuet-Higgins, M. S. 1998 Viscous streaming from an oscillating spherical bubble. Proc. R. Soc. London, Ser. A **454**, 725–742.

Mekki-Berrada, F., Combriat, T., Thibault, P. & Marmottant, P. 2016 Interactions enhance the acoustic streaming around flattened microfluidic bubbles. *J. Fluid Mech.* **797**, 851–873.

Mobadersany, N. & Sarkar, K. 2019 Acoustic microstreaming near a plane wall due to a pulsating free or coated bubble: velocity, vorticity and closed streamlines. *J. Fluid Mech.* **875**, 781–806.

Plesset, M. S. 1954 On the stability of fluid flows with spherical symmetry. *J. Appl. Phys.* **25**, 96–98.

Regnault, G., Mauger, C., Blanc-Benon, P. & Inserra, C. 2020 Secondary radiation force between two closely spaced acoustic bubbles. *Phys. Rev. E* **102**, 031101.
Shaw, S. J. 2006 Translation and oscillation of a bubble under axisymmetric deformation. *Phys. Fluids* **18**, 072104.

Tho, P., Manasseh, R. & Ooi, A. 2007 Cavitation microstreaming patterns in single and multiple bubble systems. *J. Fluid Mech.* **576**, 191–233.

Varshalovich, D. A., Moskalev, A. N. & Khersonskii, V. K. 1988 *Quantum Theory of Angular Momentum*. World Scientific.

Wang, C. & Cheng, J. 2013 Cavitation microstreaming generated by a bubble pair in an ultrasound field. *J. Acoust. Soc. Am.* **134**, 1675–1682.

Westervelt, P. J. 1953 The theory of steady rotational flow generated by a sound field. *J. Acoust. Soc. Am.* **25**, 60–67.

Zwillinger, D. 2003 *Standard Mathematical Tables and Formulae*. CRC Press.