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 non-spherical oscillations of a single gas bubble in an unbounded liquid (Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a,Reference Doinikov, Cleve, Regnault, Mauger and Inserrab; Inserra et al. Reference Inserra, Regnault, Cleve, Mauger and Doinikov2020a,Reference Inserra, Regnault, Cleve, Mauger and Doinikovb). Recently, an excellent review of earlier studies on acoustic microstreaming has also been performed by Mobadersany & Sarkar (Reference Mobadersany and Sarkar2019). 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, Manasseh & Ooi Reference Tho, Manasseh and Ooi2007; Bolanos-Jimenez et al. Reference Bolanos-Jimenez, Rossi, Fernandez Rivas, Kähler and Marin2017), pancake-like bubbles constrained in microchannels (Mekki-Berrada et al. Reference Mekki-Berrada, Combriat, Thibault and Marmottant2016) or arrays of armoured bubbles printed on the base of a microfluidic channel (Bertin et al. Reference Bertin, Spelman, Combriat, Hue, Stéphan, Lauga and Marmottant2017). 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. Reference Garbin, Cojoc, Ferrari, Di Fabrizio, Overvelde, van der Meer, de Jong, Lohse and Versluis2007). Such a trapping has recently been performed in the case of two interacting bubbles that underwent axisymmetric shape oscillations (Regnault et al. Reference Regnault, Mauger, Blanc-Benon and Inserra2020), 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 (Reference Wang and Cheng2013). 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 inter-bubble distance, which means that the distance between the equilibrium centres of the bubbles should be large compared with 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 with 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 (Reference Doinikov and Bouakaz2016). 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 non-spherical 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 the 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 Reference Plesset1954; Francescutto & Nabergoj Reference Francescutto and Nabergoj1978; Hall & Seminara Reference Hall and Seminara1980; Shaw Reference Shaw2006; Guédra & Inserra Reference Guédra and Inserra2018). 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},{\varepsilon _1})$ and 
$({r_2},{\theta _2},{\varepsilon _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 1. Coordinate systems used in calculations.
2.1. First-order scattered field
Before we 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 \theta )$, where 
${P_n}$ is the Legendre polynomial of degree n, generates a scattered wave depending on 
${P_n}(\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 \theta )$ of all degrees, i.e. 
${P_m}(\cos \theta )$. This effect can be seen in solutions obtained by Doinikov & Bouakaz (Reference Doinikov and Bouakaz2015). Although Doinikov & Bouakaz (Reference Doinikov and Bouakaz2015) assume that two bubbles only undergo pulsation and translation, the scattered wave includes terms with 
${P_m}(\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. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a,Reference Doinikov, Cleve, Regnault, Mauger and Inserrab; Inserra et al. Reference Inserra, Regnault, Cleve, Mauger and Doinikov2020a,Reference Inserra, Regnault, Cleve, Mauger and Doinikovb). 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 jth bubble undergoes 
${N_j}$ modes with a frequency 
$\omega $. The modes have numbers 
$M_1^{(j)},M_2^{(j)},M_{{N_j}}^{(j)}$. In this case, the 
$\omega $-dependent deformation of the surface of the jth bubble can be represented as
\begin{equation}r_s^{(j)} = {\textrm{e}^{ - \textrm{i}\omega t}}\mathop \sum \limits_{n = M_1^{(j)}}^{M_{{N_j}}^{(j)}} s_n^{(j)}{P_n}({\mu _j}),\end{equation}
where 
${\mu _j} = \cos {\theta _j}$ and 
$s_n^{(j)}$ is the complex amplitude of the nth mode of the jth bubble, which is assumed to be known (for example, obtained from the experimental bubble dynamics).
The total first-order liquid velocity generated by both bubbles is given by
\begin{equation}\boldsymbol{v} = {\boldsymbol{v}^{(1)}} + {\boldsymbol{v}^{(2)}},\end{equation}
where 
${\boldsymbol{v}^{(j)}}$ is the first-order liquid velocity generated by the jth bubble. The value of 
${\boldsymbol{v}^{(j)}}$ can be represented as
\begin{equation}{\boldsymbol{v}^{(j)}} = \boldsymbol{\nabla }{\varphi ^{(j)}} + \boldsymbol{\nabla } \times {\boldsymbol{\psi }^{(j)}},\end{equation}
where 
${\varphi ^{(j)}}$ and 
${\boldsymbol{\psi }^{(j)}}$ are the scalar and the vector velocity potentials, respectively.
In view of axial symmetry, 
${\varphi ^{(j)}}$ and 
${\boldsymbol{\psi }^{(j)}}$ are written as (Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a)
\begin{gather}{\varphi ^{(j)}} = {\textrm{e}^{ - \textrm{i}\omega t}}\mathop \sum \limits_{n = 0}^\infty a_n^{(j)}{\left( {\frac{{{R_{j0}}}}{{{r_j}}}} \right)^{n + 1}}{P_n}({\mu _j}),\end{gather}
\begin{gather}{\boldsymbol{\psi }^{(j)}} = {\psi ^{(j)}}({r_j},{\theta _j},t){\boldsymbol{e}_{\boldsymbol{\varepsilon j}}} = {\textrm{e}^{ - \textrm{i}\omega t}}{\boldsymbol{e}_{\boldsymbol{\varepsilon j}}}\mathop \sum \limits_{n = 1}^\infty b_n^{(j)}h_n^{(1)}({k_v}{r_j})P_n^1({\mu _j}),\end{gather}
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} = (1 + i)/\delta $ is the viscous wavenumber, 
$\delta = \sqrt {2\nu /\omega } $ is the viscous penetration depth, 
$\nu = \eta /\rho $ is the kinematic liquid viscosity, 
$\eta $ is the dynamic liquid viscosity, 
$\rho $ is the equilibrium liquid density and 
${\boldsymbol{e}_{\boldsymbol{\varepsilon j}}}$ is the unit azimuth vector of the jth bubble. Note that axial symmetry allows us to set 
${\varepsilon _1} = {\varepsilon _2}$ and 
${\boldsymbol{e}_{\boldsymbol{\varepsilon }1}} = {\boldsymbol{e}_{\boldsymbol{\varepsilon }2}}$.
Substitution of (2.4) and (2.5) into (2.3) yields
\begin{gather}v_r^{(j)} ={-} \frac{{{\textrm{e}^{ - \textrm{i}\omega t}}}}{{{r_j}}}\sum\limits_{n = 0}^\infty {(n + 1)} \left[ {a_n^{(j)}{{\left( {\frac{{{R_{j0}}}}{{{r_j}}}} \right)}^{n + 1}} + nb_n^{(j)}h_n^{(1)}({k_v}{r_j})} \right]{P_n}({\mu _j}),\end{gather}
\begin{gather}v_\theta ^{(j)} = \frac{{{\textrm{e}^{ - \textrm{i}\omega t}}}}{{{r_j}}}\sum\limits_{n = 1}^\infty {\left\{ {a_n^{(j)}{{\left( {\frac{{{R_{j0}}}}{{{r_j}}}} \right)}^{n + 1}} - b_n^{(j)}[h_n^{(1)}({k_v}{r_j}) + {k_v}{r_j}h_n^{(1)/}({k_v}{r_j})]} \right\}P_n^1({\mu _j})} ,\end{gather}
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 Reference Abramowitz and Stegun1972).
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\;\mathrm{\mu }\textrm{m}$ at the inter-bubble distance 
$d = 300\;\mathrm{\mu }\textrm{m}$ and the oscillation frequency 
$f = 30\;\textrm{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 
${({R_{j0}}/{r_j})^{n + 1}}$ and Hankel function 
$h_n^{(1)}({k_v}{r_j})$ that are respectively associated with the coefficients 
$a_n^{(1)}$ and 
$b_n^{(1)}$ in (2.6) and (2.7).
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.
2.2. Acoustic streaming
As shown by Doinikov et al. (Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a), the Eulerian velocity of acoustic streaming can be written as
\begin{equation}{\boldsymbol{v}_E} = \boldsymbol{\nabla } \times \boldsymbol{\varPsi },\end{equation}
where 
$\boldsymbol{\varPsi }$ obeys the vorticity equation derived by Westervelt (Reference Westervelt1953) for the case of solenoidal first-order motion,
\begin{equation}{\varDelta ^2}\boldsymbol{\varPsi } ={-} \frac{1}{\nu }\boldsymbol{\nabla } \times \langle \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}\rangle ,\end{equation}
with 
$\left\langle {} \right\rangle$ denoting the time average.
Equation (2.9) is shown by Doinikov et al. (Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a) to be transformed to
\begin{equation}{\varDelta ^2}\boldsymbol{\varPsi } = \frac{1}{{2\nu }}Re \{ \boldsymbol{\nabla } \times [(\varDelta {\boldsymbol{\psi }^ \ast }) \times \boldsymbol{v}]\} ,\end{equation}
where Re means “the real part of” and the asterisk denotes the complex conjugate. Substituting 
$\boldsymbol{v} = {\boldsymbol{v}^{(1)}} + {\boldsymbol{v}^{(2)}}$ and 
$\boldsymbol{\psi } = {\boldsymbol{\psi }^{(1)}} + {\boldsymbol{\psi }^{(2)}}$, one obtains
\begin{equation}\begin{array}{ccccc} {\varDelta ^2}\boldsymbol{\varPsi } & = \dfrac{1}{{2\nu }}Re \{ \boldsymbol{\nabla } \times [\varDelta ({\boldsymbol{\psi }^{(1) \ast }}) \times {\boldsymbol{v}^{(1)}}] + \boldsymbol{\nabla } \times [\varDelta ({\boldsymbol{\psi }^{(2) \ast }}) \times {\boldsymbol{v}^{(2)}}]\\ & \quad + \boldsymbol{\nabla } \times [\varDelta ({\boldsymbol{\psi }^{(1) \ast }}) \times {\boldsymbol{v}^{(2)}}] + \boldsymbol{\nabla } \times [\varDelta ({\boldsymbol{\psi }^{(2) \ast }}) \times {\boldsymbol{v}^{(1)}}]\} . \end{array}\end{equation} 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 
$\boldsymbol{\varPsi }$ in the coordinates 
$({r_1},{\theta _1})$ is written as
\begin{equation}\boldsymbol{\varPsi }({r_1},{\theta _1}) = \varPsi ({r_1},{\theta _1}){\boldsymbol{e}_{\varepsilon 1}},\end{equation}
where 
$\varPsi ({r_1},{\theta _1})$ obeys the following equation:
\begin{align}
& {\left( {\Delta_{r\theta }^{(1)} - \dfrac{1}{{r_1^2{{\sin
}^2}{\theta_1}}}} \right)^2}\varPsi ({r_1},{\theta _1}) =
\dfrac{1}{{2\nu R_{10}^4}}\left\{\displaystyle\sum_{\substack{n = 0\\ m = 1}}^\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})]}\right.\nonumber\\ &
\quad \left. + \sqrt {1 - \mu_1^2} \displaystyle\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})]}\right\}.
\end{align}
Here, the operator 
$\varDelta _{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, 
$\boldsymbol{\varPsi }$ in the coordinates 
$({r_2},{\theta _2})$ is written by
\begin{equation}\boldsymbol{\varPsi }({r_2},{\theta _2}) = \varPsi ({r_2},{\theta _2}){\boldsymbol{e}_{\varepsilon 2}},\end{equation}
where 
$\varPsi ({r_2},{\theta _2})$ is calculated by
\begin{align}&
{\left( {\Delta_{r\theta }^{(2)} - \dfrac{1}{{r_2^2{{\sin
}^2}{\theta_2}}}} \right)^2}\varPsi ({r_2},{\theta _2}) =
\dfrac{1}{{2\nu R_{20}^4}}\left\{
\displaystyle\sum_{\substack{n = 0\\ m = 1}}^\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})]}\right.\nonumber\\ &
\quad \left. { + \sqrt {1 - \mu_2^2} \displaystyle\sum_{n,m =
1}^\infty
{P_n^1({\mu_2})P_m^{1/}({\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})]} } \right\}.
\end{align}
Here, the operator 
$\varDelta _{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.
Equations (2.13) and (2.15) are solved in Appendix C. The result is given by
\begin{equation}\varPsi ({r_j},{\theta _j}) = \sum\limits_{l = 1}^\infty {\varPsi _l^{(j)}({r_j})P_l^1({\mu _j})} ,\quad j = 1,2,\end{equation}
where 
$\varPsi _l^{(j)}({r_j})$ is calculated by (C8). Substituting (2.16) into (2.8), one obtains the components of the Eulerian streaming velocity,
\begin{gather}{v_{Er}}({r_j},{\theta _j}) ={-} \frac{1}{{{r_j}}}\sum\limits_{l = 1}^\infty {l(l + 1)\varPsi _l^{(j)}({r_j}){P_l}({\mu _j})} ,\end{gather}
\begin{gather}{v_{E\theta }}({r_j},{\theta _j}) ={-} \frac{1}{{{r_j}}}\sum\limits_{l = 1}^\infty {[\varPsi _l^{(j)}({r_j}) + {r_j}\varPsi _l^{(j)/}({r_j})]P_l^1({\mu _j})} ,\end{gather}
where 
$\varPsi _l^{(j)/}({r_j})$ is calculated by (C17).
The components of the Lagrangian streaming velocity are defined by
\begin{gather}{v_{Lr}}({r_j},{\theta _j}) = {v_{Er}}({r_j},{\theta _j}) + {v_{Sr}}({r_j},{\theta _j}),\end{gather}
\begin{gather}{v_{L\theta }}({r_j},{\theta _j}) = {v_{E\theta }}({r_j},{\theta _j}) + {v_{S\theta }}({r_j},{\theta _j}),\end{gather}
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 Reference Longuet-Higgins1998), 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 (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 (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\;\mathrm{\mu }\textrm{m}$ separated by the inter-bubble distance 
$d = 150\;\mathrm{\mu }\textrm{m}$ that oscillate radially in phase, with the amplitude 5 μ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\nolimits_{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_i^{(1)}$, in the case 
$n \ne 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 4. Convergence analysis of the numerical modelling for the case of two spherically oscillating bubbles of identical radii 
$({R_{10}} = {R_{20}} = 50\;\mathrm{\mu }\textrm{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)}$ 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 (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\;\mathrm{\mu }\textrm{m,}{\theta _1} = {\rm \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. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a). Due to the presence of a second bubble, multiple sound scattering leads to a non-zero microstreaming flow. Figures 5(a)–5(c) shows streamline patterns, for different inter-bubble distances, around two bubbles with identical equilibrium radii 
${R_{10}} = {R_{20}} = 50\;\mathrm{\mu }\textrm{m}$, oscillating radially in phase with the amplitude 
$5\;\mathrm{\mu }\textrm{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} \}$.
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\;\mathrm{\mu }\textrm{m}$, (b) 
$d = 300\;\mathrm{\mu }\textrm{m}$, (c) 
$d = 750\;\mathrm{\mu }\textrm{m}$. (d) Evolution of the magnitude of the Lagrangian streaming velocity as a function of the inter-bubble distance. The velocity is measured at the locations 
$({r_1} = 1.3{R_{01}},{\theta _1} = {\rm \pi}/2)$ and 
$({r_2} = 1.5{R_{01}},{\theta _2} = {\rm \pi})$.
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.3{R_{01}},{\theta _1} = {\rm \pi}/2)$ and 
$({r_2} = 1.5{R_{01}},{\theta _2} = {\rm \pi})$. The velocity magnitude rapidly decreases to values of the order of 10 μm s−1. Such small velocities cannot be observed experimentally. It is worth noticing that these results are in good qualitative agreement with the experimental observations of Mekki-Berrada et al. (Reference Mekki-Berrada, Combriat, Thibault and Marmottant2016). Although their study is devoted to the analysis of streaming flows induced by pancake-like bubbles, the same kind of long-range streaming pattern is observed around the interacting bubble pair.
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. (Reference Inserra, Regnault, Cleve, Mauger and Doinikov2020b). Numerical results are shown in figure 6 for two identical bubbles with 
${R_{10}} = {R_{20}} = 50\;\mathrm{\mu }\textrm{m}$ at the inter-bubble distance 
$d = 200\;\mathrm{\mu }\textrm{m}$, exhibiting radial 
$(s_0^{(1,2)} = 5\;\mathrm{\mu }\textrm{m)}$ and mode 2 
$(s_2^{(1,2)} = 10\;\mathrm{\mu }\textrm{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. Note, for example, that the direction of streamlines in the areas above and below the bubbles is different in figures 6(a) and 6(b).
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.
In order to reinforce the validity of these results, a computational fluid dynamics simulation has been performed using StarCCM+ software. The computational fluid dynamics (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. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a,Reference Doinikov, Cleve, Regnault, Mauger and Inserrab; Inserra et al. Reference Inserra, Regnault, Cleve, Mauger and Doinikov2020a,Reference Inserra, Regnault, Cleve, Mauger and Doinikovb). 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 with 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 three-dimensional 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.
Figure 7. Microstreaming flow obtained by CFD simulation. (a) Cross-section view of the flow. (b) Three-dimensional view of the streamlines illustrating the vortex ring generation around the bubbles.
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 approximately 2 hours.
We have also made a comparison with the theory developed by Doinikov & Bouakaz (Reference Doinikov and Bouakaz2016). Recall that this theory, first, assumes the distance between the bubbles to be large compared with 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\;\mathrm{\mu }\textrm{m}$ at 
$d = 300\;\mathrm{\mu }\textrm{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 (Reference Doinikov and Bouakaz2016) only allows for the interaction of the bubbles in the first-order solutions.
Figure 8. Comparison of (a) the streamline pattern given by the theory of Doinikov & Bouakaz (Reference Doinikov and Bouakaz2016) and (b) the pattern predicted by the present model.
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\;\mathrm{\mu }\textrm{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. 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\;\textrm{Pa}\;\textrm{s}$; (b,e,h) 
$\eta = 0.01\;\textrm{Pa}\;\textrm{s}$; (c,f,i) 
$\eta = 0.1\;\textrm{Pa}\;\textrm{s}$.
Figure 10. The magnitude of the Lagrangian streaming velocity 
$|{\boldsymbol{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} ={-} {\rm \pi}/2$. The other parameters are as in figure 9.
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: (a,d,g) 
$\eta = 0.001\;\textrm{Pa}\;\textrm{s}$; (b,e,h) 
$\eta = 0.01\;\textrm{Pa}\;\textrm{s}$; (c,f,i) 
$\eta = 0.1\;\textrm{Pa}\;\textrm{s}$.
Figure 12. The magnitude of the Lagrangian streaming velocity 
$|{\boldsymbol{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} ={-} {\rm \pi}/2$. The other parameters are as in figure 11.
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 
$|{\boldsymbol{v}_L}|$ depends on the liquid viscosity. 
$|{\boldsymbol{v}_L}|$ 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} ={-} {\rm \pi}/2$. As one can see, increasing 
$\eta $ leads to increasing 
$|{\boldsymbol{v}_L}|$ and the peak value of 
$|{\boldsymbol{v}_L}|$ is shifted away from the bubble surface. It is interesting to note that, along direction 2, the increase of 
$|{\boldsymbol{v}_L}|$ at 
$\eta = 0.01\;\textrm{Pa}\;\textrm{s}$ is higher than at 
$\eta = 0.1$ but the region where 
$|{\boldsymbol{v}_L}|$ is relatively high remains wider at 
$\eta = 0.1$.
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\;\mathrm{\mu }\textrm{m}$ and 
$\eta = 0.1\;\textrm{Pa}\;\textrm{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, 
$|{\boldsymbol{v}_L}|$ at 
$\eta = 0.01\;\textrm{Pa}\;\textrm{s}$ is higher than at 
$\eta = 0.001\;\textrm{Pa}\;\textrm{s}$ but further increasing the viscosity to 0.1 Pa s results in a 
$|{\boldsymbol{v}_L}|$ that is lower than at 
$\eta = 0.001\;\textrm{Pa}\;\textrm{s}$. Other effects are similar to those in figure 10.
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. (Reference Mekki-Berrada, Combriat, Thibault and Marmottant2016) for flattened bubbles undergoing non-spherical oscillations in a microchannel. They consider the case of two bubbles of radii 
${R_{10}} = {R_{20}} = 38\;\mathrm{\mu }\textrm{m}$, separated by 
$d = 200\;\mathrm{\mu }\textrm{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 millimetres per second. The measured input data for 
$s_{0,1,2}^{(1,2)}$ provided by Mekki-Berrada et al. (Reference Mekki-Berrada, Combriat, Thibault and Marmottant2016) 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. (Reference Mekki-Berrada, Combriat, Thibault and Marmottant2016). 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 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. (Reference Mekki-Berrada, Combriat, Thibault and Marmottant2016), see the copyright form) and (b) the theoretical results given by the present model.
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 with 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 behaviour 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.
Funding
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 
$\boldsymbol{v}$ at 
${r_j} = {R_{j0}}$ be equal to the normal component of the surface velocity of the jth bubble,
\begin{equation}{v_r}{|_{{r_j} = {R_{j0}}}} = \frac{{\textrm{d}r_s^{(j)}}}{{\textrm{d}t}} ={-} \textrm{i}\omega \,{\textrm{e}^{ - \textrm{i}\omega t}}\sum\limits_{n = M_1^{(j)}}^{M_{{N_j}}^{(j)}} {s_n^{(j)}{P_n}({\mu _j})} .\end{equation}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,
\begin{equation}{\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\;\textrm{at}\;{r_j} = {R_{j0}}.\end{equation}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 with 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 
${\boldsymbol{\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. Reference Varshalovich, Moskalev and Khersonskii1988; Doinikov & Bouakaz Reference Doinikov and Bouakaz2015):
\begin{gather}\frac{{{P_n}({\mu _1})}}{{r_1^{n + 1}}} = \frac{1}{{{d^{n + 1}}}}\sum\limits_{m = 0}^\infty {{{( - 1)}^m}{C_{nm}}{{\left( {\frac{{{r_2}}}{d}} \right)}^m}} {P_m}({\mu _2}),\end{gather}
\begin{gather}\frac{{{P_n}({\mu _2})}}{{r_2^{n + 1}}} = \frac{{{{( - 1)}^n}}}{{{d^{n + 1}}}}\sum\limits_{m = 0}^\infty {{C_{nm}}{{\left( {\frac{{{r_1}}}{d}} \right)}^m}} {P_m}({\mu _1}),\end{gather}
where 
${C_{nm}} = (n + m)!/(n!m!)$ and d is the distance between the centres 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
\begin{gather}{\varphi ^{(1)}}({r_2},{\theta _2}) = {\textrm{e}^{ - \textrm{i}\omega t}}\sum\limits_{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})} ,\end{gather}
\begin{gather}{\varphi ^{(2)}}({r_1},{\theta _1}) = {\textrm{e}^{ - \textrm{i}\omega t}}\sum\limits_{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})} ,\end{gather}
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 
${\boldsymbol{\psi }^{(1)}}$ and 
${\boldsymbol{\psi }^{(2)}}$, the following identities are used (Varshalovich et al. Reference Varshalovich, Moskalev and Khersonskii1988):
\begin{align} &h_n^{(1)}({k_v}{r_1})P_n^1({\mu _1})\nonumber\\ &\quad = \frac{{\sqrt {n(n + 1)} }}{{(2n + 1){i^n}}}\sum\limits_{{l_1},{l_2} = 0}^\infty {\frac{{{i^{{l_1} - {l_2}}}{{( - 1)}^{{l_2}}}(2{l_1} + 1)(2{l_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_{{l_2}}^{(1)}({k_v}d){j_{{l_1}}}({k_v}{r_2})P_{{l_1}}^1({\mu _2})} ,\end{align}
\begin{align} &h_n^{(1)}({k_v}{r_2})P_n^1({\mu _2})\nonumber\\ &\quad = \frac{{\sqrt {n(n + 1)} }}{{(2n + 1){i^n}}}\sum\limits_{{l_1},{l_2} = 0}^\infty {\frac{{{i^{{l_1} - {l_2}}}(2{l_1} + 1)(2{l_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_{{l_2}}^{(1)}({k_v}d){j_{{l_1}}}({k_v}{r_1})P_{{l_1}}^1({\mu _1})} ,\end{align}
where 
$C_{{l_1}{m_1}{l_2}{m_2}}^{LM}$ are the Clebsch–Gordan coefficients (Varshalovich et al. Reference Varshalovich, Moskalev and Khersonskii1988; Zwillinger Reference Zwillinger2003) and 
${j_n}$ is the spherical Bessel function. Equations (A7) and (A8) follow from (34) of § 5.17 of the book by Varshalovich et al. (Reference Varshalovich, Moskalev and Khersonskii1988). 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{align}& {\boldsymbol{\psi
}^{(1)}}({r_2},{\theta _2}) = {\psi ^{(1)}}({r_2},{\theta
_2}){\boldsymbol{e}_{\varepsilon 2}}\nonumber\\ & \quad =
{\textrm{e}^{ - \textrm{i}\omega
t}}{\boldsymbol{e}_{\varepsilon 2}}\sum\limits_{n,m,l =
0}^\infty \dfrac{{\sqrt {n(n + 1)} }}{{(2n +
1){i^n}}}b_n^{(1)}\dfrac{{{i^{m - l}}{{( - 1)}^l}(2m +
1)(2l + 1)}}{{\sqrt {m(m + 1)}
}}\nonumber\\&\qquad \times C_{m0l0}^{n0}C_{m1l0}^{n1}h_l^{(1)}({k_v}d){j_m}({k_v}{r_2})P_m^1({\mu
_2}),\end{align}
\begin{align} & {\boldsymbol{\psi
}^{(2)}}({r_1},{\theta _1}) = {\psi ^{(2)}}({r_1},{\theta
_1}){\boldsymbol{e}_{\varepsilon 1}}\nonumber\\ & \quad =
{\textrm{e}^{ - \textrm{i}\omega
t}}{\boldsymbol{e}_{\varepsilon 1}}\sum\limits_{n,m,l =
0}^\infty \dfrac{{\sqrt {n(n + 1)} }}{{(2n +
1){i^n}}}b_n^{(2)}\dfrac{{{i^{m - l}}(2m + 1)(2l +
1)}}{{\sqrt {m(m + 1)}
}}\nonumber\\ &\qquad \times C_{m0l0}^{n0}C_{m1l0}^{n1}h_l^{(1)}({k_v}d){j_m}({k_v}{r_1})P_m^1({\mu
_1}). \end{align}
Here, we have also used the fact that 
${\boldsymbol{e}_{\varepsilon 1}} = {\boldsymbol{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
\begin{align} &
v_r^{(2)}({r_1},{\theta _1}) = {\textrm{e}^{ -
\textrm{i}\omega t}}\sum\limits_{n,m = 0}^\infty
{{P_m}({\mu _1})\left\{ {\dfrac{{{{( -
1)}^n}m}}{d}{C_{nm}}\xi_2^{n + 1}a_n^{(2)}} \right.}
{\left( {\dfrac{{{r_1}}}{d}} \right)^{m - 1}}\nonumber\\ & \quad -
\dfrac{{{j_m}({k_v}{r_1})}}{{{r_1}}}\left. {\dfrac{{\sqrt
{n(n + 1)} }}{{(2n + 1){i^n}}}(2m + 1){i^m}\sqrt {m(m + 1)}
b_n^{(2)}\sum\limits_{l = 0}^\infty {{i^{ - l}}(2l +
1)C_{m0l0}^{n0}C_{m1l0}^{n1}h_l^{(1)}({k_v}d)} } \right\}.
\end{align}
Setting j = 1 in (A1), substituting (2.6) with j = 1 and (A11) and then equating terms at the same Legendre polynomials, we obtain
\begin{align} &
(n + 1)a_n^{(1)} + n(n +
1)h_n^{(1)}({k_v}{R_{10}})b_n^{(1)} - \sum\limits_{m =
0}^\infty \left\{\vphantom{\dfrac{{\sqrt {m(m + 1)}
}}{{(2m + 1){i^m}}}}{{{( - 1)}^m}n{C_{mn}}\xi_1^n\xi_2^{m +
1}a_m^{(2)}} \right.\nonumber\\ & \qquad - {i^n}(2n + 1)\sqrt {n(n + 1)}
{j_n}({k_v}{R_{10}})\left. {\dfrac{{\sqrt {m(m + 1)}
}}{{(2m + 1){i^m}}}b_m^{(2)}\sum\limits_{l = 0}^\infty
{{i^{ - l}}(2l +
1)C_{n0l0}^{m0}C_{n1l0}^{m1}h_l^{(1)}({k_v}d)} } \right\}\nonumber\\
& \quad = \textrm{i}\omega {R_{10}}\sum\limits_{m =
M_1^{(1)}}^{M_{{N_1}}^{(1)}} {{\delta _{nm}}s_m^{(1)}}
\textrm{,}\quad n \ge 0,
\end{align}
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
\begin{align} &
v_r^{(1)}({r_2},{\theta _2}) = {\textrm{e}^{ -
\textrm{i}\omega t}}\sum\limits_{n,m = 0}^\infty
{{P_m}({\mu _2})\left\{ {\dfrac{{{{( -
1)}^m}m}}{d}{C_{nm}}\xi_1^{n + 1}a_n^{(1)}} \right.}
{\left( {\dfrac{{{r_2}}}{d}} \right)^{m - 1}}\nonumber\\ & \quad
\left. { -
\dfrac{{{j_m}({k_v}{r_2})}}{{{r_2}}}\dfrac{{\sqrt {n(n +
1)} }}{{(2n + 1){i^n}}}(2m + 1){i^m}\sqrt {m(m + 1)}
b_n^{(1)}\sum\limits_{l = 0}^\infty {{{( - 1)}^l}{i^{ -
l}}(2l + 1)C_{m0l0}^{n0}C_{m1l0}^{n1}h_l^{(1)}({k_v}d)} }
\right\}. \end{align}
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\limits_{m =
0}^\infty \left\{\vphantom{\dfrac{{\sqrt {m(m +
1)} }}{{(2m + 1){i^m}}}}{{{( - 1)}^n}n{C_{mn}}\xi_1^{m +
1}\xi_2^na_m^{(1)}}\nonumber \right.\\ & \qquad \left. - {i^n}(2n +
1)\sqrt {n(n + 1)} {j_n}({k_v}{R_{20}})ccb_m^{(1)}\sum\limits_{l = 0}^\infty
{{{( - 1)}^l}{i^{ - l}}(2l +
1)C_{n0l0}^{m0}C_{n1l0}^{m1}h_l^{(1)}({k_v}d)} \right\}\nonumber\\
& \quad = \textrm{i}\omega {R_{20}}\sum\limits_{m =
M_1^{(2)}}^{M_{{N_2}}^{(2)}} {{\delta _{nm}}s_m^{(2)}}
\textrm{,}\quad n \ge 0.
\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
\begin{align} & v_\theta
^{(1)}({r_2},{\theta _2}) = \dfrac{{{\textrm{e}^{ -
\textrm{i}\omega t}}}}{{{r_2}}}\sum\limits_{n,m = 0}^\infty
{P_m^1({\mu _2})\left\{\vphantom{\left. \quad\sum\limits_{l = 0}^\infty {{{( - 1)}^l}{i^{
- l}}(2l + 1)C_{m0l0}^{n0}C_{m1l0}^{n1}h_l^{(1)}({k_v}d)}
\right\}}\vphantom{\dfrac{{\sqrt {m(m + 1)}
}}{{(2m + 1){i^m}}}} {{{( - 1)}^m}{C_{nm}}\xi_1^{n +
1}a_n^{(1)}{{\left( {\dfrac{{{r_2}}}{d}} \right)}^m}}
\right.} \nonumber\\ & \quad \left. - [{j_m}({k_v}{r_2})
+ {k_v}{r_2}j_m^/({k_v}{r_2})]\dfrac{{\sqrt {n(n + 1)}
}}{{(2n + 1){i^n}}}\dfrac{{(2m + 1){i^m}}}{{\sqrt {m(m +
1)} }}b_n^{(1)}\right.\nonumber\\
&\left. \quad\times\sum\limits_{l = 0}^\infty {{{( - 1)}^l}{i^{
- l}}(2l + 1)C_{m0l0}^{n0}C_{m1l0}^{n1}h_l^{(1)}({k_v}d)}
\right\},
\end{align}
\begin{align}&
v_\theta^{(2)}({r_1},{\theta _1}) = \dfrac{{{\textrm{e}^{ -
\textrm{i}\omega t}}}}{{{r_1}}}\sum\limits_{n,m = 0}^\infty
P_m^1({\mu _1})\left\{{{( - 1)}^n}{C_{nm}}\xi_2^{n + 1}a_n^{(2)}{{\left(
{\dfrac{{{r_1}}}{d}} \right)}^m} \right. \nonumber\\
& \quad \left. - [{j_m}({k_v}{r_1}) +
{k_v}{r_1}j_m^/({k_v}{r_1})]\dfrac{{\sqrt {n(n + 1)}
}}{{(2n + 1){i^n}}}\dfrac{{(2m + 1){i^m}}}{{\sqrt {m(m +
1)} }}b_n^{(2)}\right. \nonumber\\ &\left.\quad \times\sum\limits_{l =
0}^\infty {{i^{ - l}}(2l + 1)C_{m0l0}^{n0}C_{m1l0}^{n1}h_l^{(1)}({k_v}d)} \right\}.
\end{align}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{align}& 2(n + 2)a_n^{(1)} + [k_v^2R_{10}^2h_n^{(1)/{/}}({k_v}{R_{10}}) + ({n^2} + n - 2)h_n^{(1)}({k_v}{R_{10}})]b_n^{(1)}\nonumber\\
&\quad - \sum\limits_{m = 0}^\infty \left\{ {{( - 1)}^m}2(n - 1){C_{mn}}\xi_1^n\xi_2^{m + 1}a_m^{(2)} \right.\nonumber\\
& \quad \left. - [k_v^2R_{10}^2j_n^{/{/}}({k_v}{R_{10}}) + ({n^2} + n - 2){j_n}({k_v}{R_{10}})]\dfrac{{\sqrt {m(m + 1)} }}{{(2m + 1){i^m}}}\dfrac{{(2n + 1){i^n}}}{{\sqrt {n(n + 1)} }}b_m^{(2)} \right.\nonumber\\
& \quad \quad \left. \times \sum\limits_{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{align}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{align}
& 2(n + 2)a_n^{(2)} + [k_v^2R_{20}^2h_n^{(1)/{/}}({k_v}{R_{20}}) + ({n^2} + n - 2)h_n^{(1)}({k_v}{R_{20}})]b_n^{(2)}\nonumber\\
&\quad - \sum\limits_{m = 0}^\infty \left\{ ( - 1)^{n}2(n - 1){C_{mn}}\xi_1^{m + 1}\xi_2^na_m^{(1)} \right.\nonumber\\
& \quad \left. - [k_v^2R_{20}^2j_n^{/{/}}({k_v}{R_{20}}) + ({n^2} + n - 2){j_n}({k_v}{R_{20}})]
\dfrac{{\sqrt {m(m + 1)} }}{{(2m + 1){i^m}}}\dfrac{{(2n + 1){i^n}}}{{\sqrt {n(n + 1)} }}b_m^{(1)}\right.\nonumber\\
& \quad \quad \left. \times \sum\limits_{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{align}From (A12) and (A14) at n = 0, it follows that
\begin{gather}a_0^{(1)} = \textrm{i}\omega {R_{10}}\sum\limits_{m = M_1^{(1)}}^{M_{{N_1}}^{(1)}} {{\delta _{0m}}s_m^{(1)}} ,\end{gather}
\begin{gather}a_0^{(2)} = \textrm{i}\omega {R_{20}}\sum\limits_{m = M_1^{(2)}}^{M_{{N_2}}^{(2)}} {{\delta _{0m}}s_m^{(2)}} ,\end{gather}
which means that 
$a_0^{(j)}$ is non-zero if mode 0 exists and hence 
$s_0^{(j)}$ is non-zero.
For 
$n \ge 1$, (A12), (A14), (A17) and (A18) are combined in the following system:
\begin{gather}a_n^{(1)} + f_n^{(1)}b_n^{(1)} - \sum\limits_{m = 1}^\infty {\{ n{\alpha _{2nm}}a_m^{(2)} - {\beta _{2nm}}b_m^{(2)}\} } = {B_{1n}},\end{gather}
\begin{gather}a_n^{(2)} + f_n^{(2)}b_n^{(2)} - \sum\limits_{m = 1}^\infty {\{ n{\alpha _{1nm}}a_m^{(1)} - {\beta _{1nm}}b_m^{(1)}\} } = {B_{2n}},\end{gather}
\begin{gather}a_n^{(1)} + g_n^{(1)}b_n^{(1)} - \sum\limits_{m = 1}^\infty {\left\{ {\frac{{{n^2} - 1}}{{n + 2}}{\alpha_{2nm}}a_m^{(2)} - {\gamma_{2nm}}b_m^{(2)}} \right\}} = {B_{3n}},\end{gather}
\begin{gather}a_n^{(2)} + g_n^{(2)}b_n^{(2)} - \sum\limits_{m = 1}^\infty {\left\{ {\frac{{{n^2} - 1}}{{n + 2}}{\alpha_{1nm}}a_m^{(1)} - {\gamma_{1nm}}b_m^{(1)}} \right\}} = {B_{4n}},\end{gather}where the following designations are used:
\begin{gather}f_n^{(j)} = nh_n^{(1)}({k_v}{R_{j0}}),\end{gather}
\begin{gather}g_n^{(j)} = \frac{{k_v^2R_{j0}^2h_n^{(1)/{/}}({k_v}{R_{j0}}) + ({n^2} + n - 2)h_n^{(1)}({k_v}{R_{j0}})}}{{2(n + 2)}},\end{gather}
\begin{gather}{\alpha _{1nm}} = \frac{{{{( - 1)}^n}{C_{mn}}}}{{n + 1}}\xi _1^{m + 1}\xi _2^n,, \end{gather}
\begin{gather}{\alpha _{2nm}} = \frac{{{{( - 1)}^m}{C_{mn}}}}{{n + 1}}\xi _1^n\xi _2^{m + 1},\end{gather}
\begin{gather}{\beta _{1nm}} = n{j_n}({k_v}{R_{20}}){G_{nm}},\end{gather}
\begin{gather}{\beta _{2nm}} = n{j_n}({k_v}{R_{10}}){H_{nm}},\end{gather}
\begin{gather}{\gamma _{1nm}} = \frac{{[k_v^2R_{20}^2j_n^{/{/}}({k_v}{R_{20}}) + ({n^2} + n - 2){j_n}({k_v}{R_{20}})]{G_{nm}}}}{{2(n + 2)}},\end{gather}
\begin{gather}{\gamma _{2nm}} = \frac{{[k_v^2R_{10}^2j_n^{/{/}}({k_v}{R_{10}}) + ({n^2} + n - 2){j_n}({k_v}{R_{10}})]{H_{nm}}}}{{2(n + 2)}},\end{gather}
\begin{gather}{G_{nm}} = \frac{{{i^{n - m}}(2n + 1)\sqrt {m(m + 1)} }}{{(2m + 1)\sqrt {n(n + 1)} }}\sum\limits_{l = |n - m|}^{n + m} {{{( - i)}^{ - l}}(2l + 1)C_{n0l0}^{m0}C_{n1l0}^{m1}h_l^{(1)}({k_v}d)} ,\end{gather}
\begin{gather}{H_{nm}} = \frac{{{i^{n - m}}(2n + 1)\sqrt {m(m + 1)} }}{{(2m + 1)\sqrt {n(n + 1)} }}\sum\limits_{l = |n - m|}^{n + m} {{i^{ - l}}(2l + 1)C_{n0l0}^{m0}C_{n1l0}^{m1}h_l^{(1)}({k_v}d)} ,\end{gather}
\begin{gather}{B_{1n}} = \frac{n}{{n + 1}}\xi _1^n{\xi _2}a_0^{(2)} + \frac{{\textrm{i}\omega {R_{10}}}}{{n + 1}}\sum\limits_{m = M_1^{(1)}}^{M_{{N_1}}^{(1)}} {{\delta _{nm}}s_m^{(1)}} ,\end{gather}
\begin{gather}{B_{2n}} = \frac{{{{( - 1)}^n}n}}{{n + 1}}{\xi _1}\xi _2^na_0^{(1)} + \frac{{\textrm{i}\omega {R_{20}}}}{{n + 1}}\sum\limits_{m = M_1^{(2)}}^{M_{{N_2}}^{(2)}} {{\delta _{nm}}s_m^{(2)}} ,\end{gather}
\begin{gather}{B_{3n}} = \frac{{n - 1}}{{n + 2}}\xi _1^n{\xi _2}a_0^{(2)},\end{gather}
\begin{gather}{B_{4n}} = \frac{{{{( - 1)}^n}(n - 1)}}{{n + 2}}{\xi _1}\xi _2^na_0^{(1)}.\end{gather}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 Reference Zwillinger2003).
From (A22) and (A24), it follows that
\begin{align} &
a_n^{(2)} = \dfrac{{g_n^{(2)}{B_{2n}} -
f_n^{(2)}{B_{4n}}}}{{g_n^{(2)} - f_n^{(2)}}} -
\dfrac{{({n^2} - 1)f_n^{(2)} - n(n + 2)g_n^{(2)}}}{{(n +
2)(g_n^{(2)} - f_n^{(2)})}}\sum\limits_{m = 1}^\infty
{{\alpha _{1nm}}a_m^{(1)}}\nonumber \\ & \quad -
\dfrac{1}{{g_n^{(2)} - f_n^{(2)}}}\sum\limits_{m =
1}^\infty {(g_n^{(2)}{\beta _{1nm}} - f_n^{(2)}{\gamma
_{1nm}})b_m^{(1)}} , \end{align}
\begin{align}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\limits_{m = 1}^\infty {{\alpha
_{1nm}}a_m^{(1)}} - \frac{1}{{f_n^{(2)} -
g_n^{(2)}}}\sum\limits_{m = 1}^\infty {({\beta _{1nm}}
\!-\! {\gamma _{1nm}})b_m^{(1)}}
.\end{align}
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
\begin{gather}\sum\limits_{k = 1}^\infty {{A_{1nk}}a_k^{(1)}} + \sum\limits_{k = 1}^\infty {{A_{2nk}}b_k^{(1)}} = {D_{1n}},\end{gather}
\begin{gather}\sum\limits_{k = 1}^\infty {{A_{3nk}}a_k^{(1)}} + \sum\limits_{k = 1}^\infty {{A_{4nk}}b_k^{(1)}} = {D_{2n}},\end{gather}where
\begin{align} {A_{1nk}} &= {\delta _{nk}}(g_n^{(1)} - f_n^{(1)}) - \sum\limits_{m = 1}^\infty {\dfrac{{(2m + 1)(g_n^{(1)}{\beta _{2nm}} - f_n^{(1)}{\gamma _{2nm}}){\alpha _{1mk}}}}{{(m + 2)(g_m^{(2)} - f_m^{(2)})}}}\nonumber \\ & \quad - \dfrac{{({n^2} - 1)f_n^{(1)} - n(n + 2)g_n^{(1)}}}{{n + 2}}\sum\limits_{m = 1}^\infty {\dfrac{{[({m^2} - 1)f_m^{(2)} - m(m + 2)g_m^{(2)}]{\alpha _{2nm}}{\alpha _{1mk}}}}{{(m + 2)(g_m^{(2)} - f_m^{(2)})}}} , \end{align}
\begin{align} {A_{2nk}} &= \sum\limits_{m = 1}^\infty {\dfrac{{({\beta _{1mk}} - {\gamma _{1mk}})(g_n^{(1)}{\beta _{2nm}} - f_n^{(1)}{\gamma _{2nm}})}}{{g_m^{(2)} - f_m^{(2)}}}} \nonumber\\ & \quad + \dfrac{{n(n + 2)g_n^{(1)} - ({n^2} - 1)f_n^{(1)}}}{{n + 2}}\sum\limits_{m = 1}^\infty {\dfrac{{{\alpha _{2nm}}(g_m^{(2)}{\beta _{1mk}} - f_m^{(2)}{\gamma _{1mk}})}}{{g_m^{(2)} - f_m^{(2)}}}} , \end{align}
\begin{align}{D_{1n}} &= g_n^{(1)}{B_{1n}} - f_n^{(1)}{B_{3n}} + \dfrac{{n(n + 2)g_n^{(1)} - ({n^2} - 1)f_n^{(1)}}}{{n + 2}}\sum\limits_{m = 1}^\infty {\dfrac{{(g_m^{(2)}{B_{2m}} - f_m^{(2)}{B_{4m}}){\alpha _{2nm}}}}{{g_m^{(2)} - f_m^{(2)}}}} \nonumber\\ & \quad + \sum\limits_{m = 1}^\infty {\dfrac{{(g_n^{(1)}{\beta _{2nm}} - f_n^{(1)}{\gamma _{2nm}})({B_{2m}} - {B_{4m}})}}{{g_m^{(2)} - f_m^{(2)}}}} , \end{align}
\begin{align}{A_{3nk}} &= \frac{{2n + 1}}{{n + 2}}\sum\limits_{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)})}}}\nonumber\\ & \quad- \sum\limits_{m = 1}^\infty {\frac{{(2m + 1)({\beta _{2nm}} - {\gamma _{2nm}}){\alpha _{1mk}}}}{{(m + 2)(g_m^{(2)} - f_m^{(2)})}}} ,\end{align}
\begin{align}{A_{4nk}} &= (f_n^{(1)} - g_n^{(1)}){\delta _{nk}} + \frac{{2n + 1}}{{n + 2}}\sum\limits_{m = 1}^\infty {\frac{{(g_m^{(2)}{\beta _{1mk}} - f_m^{(2)}{\gamma _{1mk}}){\alpha _{2nm}}}}{{g_m^{(2)} - f_m^{(2)}}}}\nonumber\\ &\quad + \sum\limits_{m = 1}^\infty {\frac{{({\beta _{1mk}} - {\gamma _{1mk}})({\beta _{2nm}} - {\gamma _{2nm}})}}{{g_m^{(2)} - f_m^{(2)}}}} ,\end{align}
\begin{align}{D_{2n}} &= {B_{1n}} - {B_{3n}} + \frac{{2n + 1}}{{n + 2}}\sum\limits_{m = 1}^\infty {\frac{{(g_m^{(2)}{B_{2m}} - f_m^{(2)}{B_{4m}}){\alpha _{2nm}}}}{{g_m^{(2)} - f_m^{(2)}}}}\nonumber\\ &\quad + \sum\limits_{m = 1}^\infty {\frac{{({\beta _{2nm}} - {\gamma _{2nm}})({B_{2m}} - {B_{4m}})}}{{g_m^{(2)} - f_m^{(2)}}}} .\end{align} 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:
\begin{gather}\sum\limits_{k = 1}^{{N_a}} {{A_{1nk}}a_k^{(1)}} + \sum\limits_{k = 1}^{{N_b}} {{A_{2nk}}b_k^{(1)}} = {D_{1n}},\quad n = 1,2 \ldots ,{N_a},\end{gather}
\begin{gather}\sum\limits_{k = 1}^{{N_a}} {{A_{3nk}}a_k^{(1)}} + \sum\limits_{k = 1}^{{N_b}} {{A_{4nk}}b_k^{(1)}} = {D_{2n}},\quad n = 1,2, \ldots {N_b}.\end{gather}
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,
\begin{gather}{\psi ^{(1)}}({r_2},{\theta _2}) = {\textrm{e}^{ - \textrm{i}\omega t}}\sum\limits_{n,m = 1}^\infty {b_m^{(1)}{G_{nm}}{j_n}({k_v}{r_2})P_n^1({\mu _2})} ,\end{gather}
\begin{gather}v_r^{(1)}({r_2},{\theta _2}) = \frac{{{\textrm{e}^{ - \textrm{i}\omega t}}}}{{{r_2}}}\displaystyle\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})} ,\end{gather}
\begin{gather}
v_\theta^{(1)}({r_2},{\theta _2}) = \frac{{{\textrm{e}^{ -
\textrm{i}\omega t}}}}{{{r_2}}}\displaystyle\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}\right.\nonumber\\
\left. \quad - b_m^{(1)}{G_{nm}}[{j_n}({k_v}{r_2}) +
{k_v}{r_2}j_n^/({k_v}{r_2})] \right\}P_n^1({\mu _2})
,\end{gather}
\begin{gather}{\psi ^{(2)}}({r_1},{\theta _1}) = {\textrm{e}^{ - \textrm{i}\omega t}}\sum\limits_{n,m = 1}^\infty {b_m^{(2)}{H_{nm}}{j_n}({k_v}{r_1})P_n^1({\mu _1})} ,\end{gather}
\begin{gather}v_r^{(2)}({r_1},{\theta _1}) = \frac{{{\textrm{e}^{ - \textrm{i}\omega t}}}}{{{r_1}}}\displaystyle\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})} ,\end{gather}
\begin{gather}
v_\theta^{(2)}({r_1},{\theta _1}) = \frac{{{\textrm{e}^{ - \textrm{i}\omega t}}}}{{{r_1}}}\displaystyle\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}\right.\nonumber\\
\left. \quad - b_m^{(2)}{H_{nm}}[{j_n}({k_v}{r_1}) + {k_v}{r_1}j_n^/({k_v}{r_1})] \right\} P_n^1({\mu _1}) .\end{gather}These equations are used in the calculation of acoustic streaming.
Appendix B. Representation of the equation for 
$\varPsi $ in different coordinates
Equation (2.11) can be written as
\begin{equation}{\varDelta ^2}\boldsymbol{\varPsi } = {\boldsymbol{T}_1} + {\boldsymbol{T}_2} + {\boldsymbol{T}_{12}},\end{equation}where
\begin{gather}{\boldsymbol{T}_j} = \frac{1}{{2\nu }}Re \{ \boldsymbol{\nabla } \times [(\varDelta {\boldsymbol{\psi }^{(j) \ast }}) \times {\boldsymbol{v}^{(j)}}]\} \quad j = 1,2,\end{gather}
\begin{gather}{\boldsymbol{T}_{12}} = \frac{1}{{2\nu }}Re \{ \boldsymbol{\nabla } \times [(\varDelta {\boldsymbol{\psi }^{(1) \ast }}) \times {\boldsymbol{v}^{(2)}}] + \boldsymbol{\nabla } \times [(\varDelta {\boldsymbol{\psi }^{(2) \ast }}) \times {\boldsymbol{v}^{(1)}}]\} .\end{gather}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. (Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a), (B2) with j = 1 can be represented as
\begin{align}
{\boldsymbol{T}_1}({r_1},{\theta _1}) &=
\frac{{{\boldsymbol{e}_{\varepsilon 1}}}}{{2\nu {r_1}}}Re
\left\{ \frac{\partial }{{\partial {r_1}}}\left[
{{r_1}v_r^{(1)}\left( {\Delta_{r\theta }^{(1)}{\psi^{(1)
\ast }} - \frac{{{\psi^{(1) \ast }}}}{{r_1^2{{\sin
}^2}{\theta_1}}}} \right)} \right]\right.\nonumber\\ &\left. \quad + \frac{\partial
}{{\partial {\theta_1}}}\left[ {v_\theta^{(1)}\left(
{\Delta_{r\theta }^{(1)}{\psi^{(1) \ast }} -
\frac{{{\psi^{(1) \ast }}}}{{r_1^2{{\sin }^2}{\theta_1}}}}
\right)} \right] \right\},\end{align}
where the operator 
$\varDelta _{r\theta }^{(j)}$ is defined by
\begin{equation}\varDelta _{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).\end{equation} Here, 
${\boldsymbol{\psi }^{(j)}}$ obeys the equation 
$\varDelta {\boldsymbol{\psi }^{(j)}} ={-} k_v^2{\boldsymbol{\psi }^{(j)}}$, which gives
\begin{equation}\varDelta _{r\theta }^{(j)}{\psi ^{(j)}} - \frac{{{\psi ^{(j)}}}}{{r_j^2{{\sin }^2}{\theta _j}}} ={-} k_v^2{\psi ^{(j)}}.\end{equation}Considering (B6), (B4) reduces to
\begin{equation}{\boldsymbol{T}_1}({r_1},{\theta _1}) = {{T}_1}({r_1},{\theta _1}){\boldsymbol{e}_{\varepsilon 1}} = \frac{{{\boldsymbol{e}_{\varepsilon 1}}}}{{2\nu {r_1}}}Re \left\{ {k_v^2\left[ {\frac{\partial }{{\partial {r_1}}}({r_1}v_r^{(1)}{\psi^{(1) \ast }}) + \frac{\partial }{{\partial {\theta_1}}}(v_\theta^{(1)}{\psi^{(1) \ast }})} \right]} \right\}.\end{equation}Substitution of (2.5)–(2.7) into (B7) yields
\begin{align}{{T}_1}({r_1},{\theta _1}) &= \frac{1}{{2\nu R_{10}^4}}\left[
\displaystyle\sum_{\substack{n = 0\\ m = 1}}^\infty
{{P_n}({\mu_1})P_m^1({\mu_1})K_{1nm}^{(1)}({r_1})}\right.\nonumber\\ &\quad \left.+ \sqrt
{1 - \mu_1^2} \sum\limits_{n,m = 1}^\infty
{{{[P_n^1({\mu_1})P_m^1({\mu_1})]}^/}K_{2nm}^{(1)}({r_1})}\right],\end{align}
where
\begin{align}
K_{1nm}^{(1)}({r_1}) &= Re \left\{ b_m^{(1) \ast
}\dfrac{{(n + 1)k_v^2R_{10}^3}}{{{r_1}}}\left\{
{a_n^{(1)}{{\left( {\dfrac{{{R_{10}}}}{{{r_1}}}}
\right)}^{n + 1}}\left[ {\dfrac{{(n +
1){R_{10}}}}{{{r_1}}}h_m^{(1) \ast }({k_v}{r_1}) - k_v^
\ast {R_{10}}h_m^{(1)/{\ast} }({k_v}{r_1})} \right]}
\right.\right.\nonumber\\ & \quad \left.\left.-
nb_n^{(1)}[{k_v}{R_{10}}h_n^{(1)/}({k_v}{r_1})h_m^{(1) \ast
}({k_v}{r_1}) + k_v^ \ast
{R_{10}}h_n^{(1)}({k_v}{r_1})h_m^{(1)/{\ast}
}({k_v}{r_1})] \vphantom{{{\left( {\dfrac{{{R_{10}}}}{{{r_1}}}}
\right)}^{n + 1}}}\right\} \right\}, \end{align}
\begin{align}K_{2nm}^{(1)}({r_1}) = Re \left\{ {b_m^{(1) \ast }\frac{{k_v^2R_{10}^4h_m^{(1) \ast }({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\}.\end{align} In a similar way, 
${\boldsymbol{T}_2}$ in the coordinates 
$({r_1},{\theta _1})$ is written as
\begin{equation}{\boldsymbol{T}_2}({r_1},{\theta _1}) = {{T}_2}({r_1},{\theta _1}){\boldsymbol{e}_{\varepsilon 1}} = \frac{{{\boldsymbol{e}_{\varepsilon 1}}}}{{2\nu {r_1}}}Re \left\{ {k_v^2\left[ {\frac{\partial }{{\partial {r_1}}}({r_1}v_r^{(2)}{\psi^{(2) \ast }}) + \frac{\partial }{{\partial {\theta_1}}}(v_\theta^{(2)}{\psi^{(2) \ast }})} \right]} \right\},\end{equation}
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
\begin{align}{{T}_2}({r_1},{\theta _1}) &= \frac{1}{{2\nu R_{10}^4}}\left[
\sum\limits_{n,m = 1}^\infty
{{P_n}({\mu_1})P_m^1({\mu_1})L_{1nm}^{(1)}({r_1})}\right.\nonumber\\ &\left.\quad + \sqrt
{1 - {\mu^2}} \sum\limits_{n,m = 1}^\infty {{{[P_n^1({\mu_1})P_m^1({\mu_1})]}^/}L_{2nm}^{(1)}({r_1})}
\right],\end{align}
where
\begin{align}
L_{1nm}^{(1)}({r_1}) &= Re \left\{
{\dfrac{{nk_v^2R_{10}^3}}{{{r_1}}}\sum\limits_{k =
1}^\infty {{{(b_k^{(2)}{H_{mk}})}^ \ast }\left\{ {{{\left(
{\dfrac{{{r_1}}}{d}} \right)}^n}\left[ {k_v^ \ast
{R_{10}}j_m^/(k_v^ \ast {r_1}) +
\dfrac{{n{R_{10}}}}{{{r_1}}}{j_m}(k_v^ \ast {r_1})}
\right]\sum\limits_{k = 0}^\infty {a_k^{(2)}{{( -
1)}^k}{C_{kn}}\xi_2^{k + 1}} } \right.} } \right.\nonumber\\ & \quad
\left. {\left. { - (n +
1)[{k_v}{R_{10}}j_n^/({k_v}{r_1}){j_m}(k_v^ \ast {r_1}) +
k_v^ \ast {R_{10}}{j_n}({k_v}{r_1})j_m^/(k_v^ \ast
{r_1})]\sum\limits_{k = 1}^\infty {b_k^{(2)}{H_{nk}}} }
\right\}} \right\},
\end{align}
\begin{align}
L_{2nm}^{(1)}({r_1}) &= Re \left\{
{\dfrac{{k_v^2R_{10}^4{j_m}(k_v^ \ast
{r_1})}}{{r_1^2}}\sum\limits_{k = 1}^\infty
{{{(b_k^{(2)}{H_{mk}})}^ \ast }} } \right.\nonumber\\ & \quad \times
\left. {\left[ {[{j_n}({k_v}{r_1}) +
{k_v}{r_1}j_n^/({k_v}{r_1})]\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{nk}}} - {{\left( {\dfrac{{{r_1}}}{d}}
\right)}^n}\sum\limits_{k = 0}^\infty {a_k^{(2)}{{( -
1)}^k}{C_{kn}}\xi_2^{k + 1}} } \right]} \right\}.
\end{align}
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 
$\varDelta {\boldsymbol{\psi }^{(j)}} ={-} k_v^2{\boldsymbol{\psi }^{(j)}}$, (B3) can be written as
\begin{equation}{\boldsymbol{T}_{12}} = \frac{1}{{2\nu }}Re \{ k_v^2\boldsymbol{\nabla } \times ({\boldsymbol{\psi }^{(1) \ast }} \times {\boldsymbol{v}^{(2)}}) + k_v^2\boldsymbol{\nabla } \times ({\boldsymbol{\psi }^{(2) \ast }} \times {\boldsymbol{v}^{(1)}})\} .\end{equation}With the help of the identity
\begin{equation}\boldsymbol{\nabla } \times (\boldsymbol{a} \times \boldsymbol{b}) = (\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{a} - (\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{b} + \boldsymbol{a}(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{b}) - \boldsymbol{b}(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{a}),\end{equation}
taking into account that 
$\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{v}^{(j)}} = 0$ and 
$\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{\psi }^{(j)}} = 0$, (B15) is transformed to
\begin{equation}{\boldsymbol{T}_{12}} = \frac{1}{{2\nu }}Re \{ k_v^2[({\boldsymbol{v}^{(2)}}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol{\psi }^{(1) \ast }} - ({\boldsymbol{\psi }^{(1) \ast }}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol{v}^{(2)}} + ({\boldsymbol{v}^{(1)}}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol{\psi }^{(2) \ast }} - ({\boldsymbol{\psi }^{(2) \ast }}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol{v}^{(1)}}]\} .\end{equation} Expressed in the coordinates 
$({r_1},{\theta _1})$, (B17) gives
\begin{align} &
{\boldsymbol{T}_{12}}({r_1},{\theta _1}) = {{T}_{12}}({r_1},{\theta _1}){\boldsymbol{e}_{\varepsilon 1}}\nonumber\\ &\quad = \dfrac{{{\boldsymbol{e}_{\varepsilon 1}}}}{{2\nu }}Re
\left\{ {k_v^2\left[ {v_r^{(1)}\dfrac{{\partial {\psi^{(2)
\ast }}}}{{\partial {r_1}}} + v_r^{(2)}\dfrac{{\partial
{\psi^{(1) \ast }}}}{{\partial {r_1}}} +
\dfrac{{v_\theta^{(1)}}}{{{r_1}}}\dfrac{{\partial
{\psi^{(2) \ast }}}}{{\partial {\theta_1}}} +
\dfrac{{v_\theta^{(2)}}}{{{r_1}}}\dfrac{{\partial
{\psi^{(1) \ast }}}}{{\partial {\theta_1}}}} \right.}
\right.\nonumber\\ & \qquad \left. {\left. { - \dfrac{{{\psi^{(1)
\ast }}(v_r^{(2)}\sin {\theta_1} + v_\theta^{(2)}\cos
{\theta_1}) + {\psi^{(2) \ast }}(v_r^{(1)}\sin {\theta_1} +
v_\theta^{(1)}\cos {\theta_1})}}{{{r_1}\sin {\theta_1}}}}
\right]} \right\},
\end{align}
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
\begin{equation}{{T}_{12}}({r_1},{\theta _1}) = \frac{1}{{2\nu
R_{10}^4}}\left[ {\displaystyle\sum_{\substack{n = 0\\ m =
1}}^\infty
{{P_n}({\mu_1})P_m^1({\mu_1})M_{1nm}^{(1)}({r_1})} + \sqrt
{1 - \mu_1^2} \sum\limits_{n,m = 1}^\infty
{{{[P_n^1({\mu_1})P_m^1({\mu_1})]}^/}M_{2nm}^{(1)}({r_1})}
} \right],\end{equation}
where
\begin{align}
M_{1nm}^{(1)}({r_1}) &= Re \left\{
{\dfrac{{k_v^2R_{10}^4}}{{r_1^2}}\left[ {(n + 1){{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{n + 1}}[(n +
1){j_m}(k_v^ \ast {r_1}) - k_v^ \ast {r_1}j_m^/(k_v^ \ast
{r_1})]a_n^{(1)}\sum\limits_{k = 1}^\infty {b_k^{(2) \ast
}H_{mk}^ \ast } } \right.} \right.\nonumber\\ & \quad + n{\left(
{\dfrac{{{r_1}}}{d}} \right)^n}[nh_m^{(1) \ast
}({k_v}{r_1}) + k_v^ \ast {r_1}h_m^{(1)/{\ast}
}({k_v}{r_1})]b_m^{(1) \ast }\sum\limits_{k = 0}^\infty
{a_k^{(2)}{{( - 1)}^k}{C_{kn}}\xi _2^{k + 1}}\nonumber \\ & \quad -
n(n + 1)[{k_v}{r_1}{j_m}(k_v^ \ast
{r_1})h_n^{(1)/}({k_v}{r_1}) + k_v^ \ast {r_1}j_m^/(k_v^
\ast {r_1})h_n^{(1)}({k_v}{r_1})]b_n^{(1)}\sum\limits_{k =
1}^\infty {b_k^{(2) \ast }H_{mk}^ \ast }\nonumber \\ & \quad \left.
{\left. { - n(n + 1)[k_v^ \ast
{r_1}{j_n}({k_v}{r_1})h_m^{(1)/{\ast} }({k_v}{r_1}) +
{k_v}{r_1}j_n^/({k_v}{r_1})h_m^{(1) \ast
}({k_v}{r_1})]b_m^{(1) \ast }\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{nk}}} } \right]} \right\},
\end{align}
\begin{align}
M_{2nm}^{(1)}({r_1}) &= Re \left\{
{\dfrac{{k_v^2R_{10}^4}}{{r_1^2}}\left[
{{j_m}({k_v}{r_1}){{\left( {\dfrac{{{R_{10}}}}{{{r_1}}}}
\right)}^{n + 1}}a_n^{(1) \ast }\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{mk}} + h_n^{(1)}({k_v}{r_1}){{\left(
{\dfrac{{{r_1}}}{d}} \right)}^m}b_n^{(1)}\sum\limits_{k =
0}^\infty {a_k^{(2) \ast }{{( - 1)}^k}{C_{km}}\xi_2^{k +
1}} } } \right.} \right.\nonumber\\ & \quad \left. {\left. { +
[{k_v}{r_1}j_m^/({k_v}{r_1})h_n^{(1) \ast }({k_v}{r_1}) -
k_v^ \ast {r_1}{j_m}({k_v}{r_1})h_n^{(1)/{\ast}
}({k_v}{r_1})]b_n^{(1) \ast }\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{mk}}} } \right]} \right\}.
\end{align}
It follows from the above derivation that 
$\boldsymbol{\varPsi }$ can be represented as
\begin{equation}\boldsymbol{\varPsi }({r_1},{\theta _1}) = \varPsi ({r_1},{\theta _1}){\boldsymbol{e}_{\varepsilon 1}},\end{equation}
where 
$\varPsi ({r_1},{\theta _1})$ obeys the following equation:
\begin{align} &
{\left( {\Delta_{r\theta }^{(1)} - \dfrac{1}{{r_1^2{{\sin
}^2}{\theta_1}}}} \right)^2}\varPsi ({r_1},{\theta _1}) =
\dfrac{1}{{2\nu R_{10}^4}}\left\{ \displaystyle\sum_{\substack{n = 0\\ m = 1}}^\infty
{P_n}({\mu_1})P_m^1({\mu_1})[K_{1nm}^{(1)}({r_1}) +
L_{1nm}^{(1)}({r_1}) \right.\nonumber\\ & \left.\quad + M_{1nm}^{(1)}({r_1})] + \sqrt {1 - \mu_1^2} \sum\limits_{n,m =
1}^\infty P_n^1({\mu_1})P_m^{1/}({\mu_1})[K_{2nm}^{(1)}({r_1})\right.\nonumber\\ & \left.\quad + 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})] \vphantom{\displaystyle\sum_{\substack{n = 0\\ m = 1}}^\infty}\right\}.
\end{align}
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
\begin{align}{\boldsymbol{T}_1}({r_2},{\theta _2}) = {{T}_1}({r_2},{\theta _2}){\boldsymbol{e}_{\varepsilon 2}} = \frac{{{\boldsymbol{e}_{\varepsilon 2}}}}{{2\nu {r_2}}}Re \left\{ {k_v^2\left[ {\frac{\partial }{{\partial {r_2}}}({r_2}v_r^{(1)}{\psi^{(1) \ast }}) + \frac{\partial }{{\partial {\theta_2}}}(v_\theta^{(1)}{\psi^{(1) \ast }})} \right]} \right\}.\end{align}Substitution of (A51)–(A53) into (B24) yields
\begin{align}{{T}_1}({r_2},{\theta _2}) &= \frac{1}{{2\nu R_{20}^4}}
\left[\sum\limits_{n,m = 1}^\infty {{P_n}({\mu_2})P_m^1({\mu_2})K_{1nm}^{(2)}({r_2})} + \sqrt{1 - \mu_2^2}
\sum\limits_{n,m = 1}^\infty {{{[P_n^1({\mu_2})P_m^1({\mu_2})]}^/}K_{2nm}^{(2)}({r_2})} \right],\end{align}
where
\begin{align} &
K_{1nm}^{(2)}({r_2}) = Re \left\{
{\dfrac{{nk_v^2R_{20}^3}}{{{r_2}}}\sum\limits_{k =
1}^\infty {{{(b_k^{(1)}{G_{mk}})}^ \ast }\left\{ {{{( -
1)}^n}{{\left( {\dfrac{{{r_2}}}{d}} \right)}^n}\left[ {k_v^
\ast {R_{20}}j_m^/(k_v^ \ast {r_2}) +
\dfrac{{n{R_{20}}}}{{{r_2}}}{j_m}(k_v^ \ast {r_2})}
\right]\sum\limits_{k = 0}^\infty
{a_k^{(1)}{C_{kn}}\xi_1^{k + 1}} } \right.} } \right.\nonumber\\ &
\quad \left. {\left. { - (n +
1)[{k_v}{R_{20}}j_n^/({k_v}{r_2}){j_m}(k_v^ \ast {r_2}) +
k_v^ \ast {R_{20}}{j_n}({k_v}{r_2})j_m^/(k_v^ \ast
{r_2})]\sum\limits_{k = 1}^\infty {b_k^{(1)}{G_{nk}}} }
\right\}} \right\},\end{align}
\begin{align}
K_{2nm}^{(2)}({r_2}) &= Re \left\{ {\dfrac{{k_v^2R_{20}^4{j_m}(k_v^ \ast {r_2})}}{{r_2^2}}}
\right.\sum\limits_{k = 1}^\infty {{{(b_k^{(1)}{G_{mk}})}^
\ast }} \nonumber\\ & \quad \left. { \times \left\{
{[{j_n}({k_v}{r_2}) +
{k_v}{r_2}j_n^/({k_v}{r_2})]\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{nk}}} - {{( - 1)}^n}{{\left(
{\dfrac{{{r_2}}}{d}} \right)}^n}\sum\limits_{k = 0}^\infty
{a_k^{(1)}{C_{kn}}\xi_1^{k + 1}} } \right\}} \right\}.
\end{align}
Replacing index 1 in (B7) with index 2, one obtains
\begin{equation}{\boldsymbol{T}_2}({r_2},{\theta _2}) = {{T}_2}({r_2},{\theta _2}){\boldsymbol{e}_{\varepsilon 2}} = \frac{{{\boldsymbol{e}_{\varepsilon 2}}}}{{2\nu {r_2}}}Re \left\{ {k_v^2\left[ {\frac{\partial }{{\partial {r_2}}}({r_2}v_r^{(2)}{\psi^{(2) \ast }}) + \frac{\partial }{{\partial {\theta_2}}}(v_\theta^{(2)}{\psi^{(2) \ast }})} \right]} \right\}.\end{equation}Substitution of (2.5)–(2.7) into (B28) yields
\begin{align}
{{T}_2}({r_2},{\theta _2}) = \frac{1}{{2\nu R_{20}^4}}\left[
{\displaystyle\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\limits_{n,m = 1}^\infty
{{{[P_n^1({\mu_2})P_m^1({\mu_2})]}^/}L_{2nm}^{(2)}({r_2})}
} \right],\end{align}
where
\begin{align}
L_{1nm}^{(2)}({r_2}) &= Re \left\{ {b_m^{(2) \ast
}\dfrac{{(n + 1)k_v^2R_{20}^3}}{{{r_2}}}\left\{
{a_n^{(2)}{{\left( {\dfrac{{{R_{20}}}}{{{r_2}}}}
\right)}^{n + 1}}\left[ {\dfrac{{(n +
1){R_{20}}}}{{{r_2}}}h_m^{(1) \ast }({k_v}{r_2}) - k_v^
\ast {R_{20}}h_m^{(1)/{\ast} }({k_v}{r_2})} \right]}
\right.} \right.\nonumber\\ & \quad { { -
nb_n^{(2)}[{k_v}{R_{20}}h_n^{(1)/}({k_v}{r_2})h_m^{(1) \ast
}({k_v}{r_2}) + k_v^ \ast
{R_{20}}h_n^{(1)}({k_v}{r_2})h_m^{(1)/{\ast}
}({k_v}{r_2})]} \}} \},
\end{align}
\begin{equation}L_{2nm}^{(2)}({r_2}) = Re \left\{ {b_m^{(2) \ast }\frac{{k_v^2R_{20}^4h_m^{(1) \ast }({k_v}{r_2})}}{{r_2^2}}\left[ {b_n^{(2)}[h_n^{(1)}({k_v}{r_2}) + {k_v}{r_2}h_n^{(1)/}({k_v}{r_2})] - a_n^{(2)}{{\left( {\frac{{{R_{20}}}}{{{r_2}}}} \right)}^{n + 1}}} \right]} \right\}.\end{equation}Interchanging indices 1 and 2 in (B18), one obtains
\begin{align}
{\boldsymbol{T}_{12}}({r_2},{\theta _2}) &=
{{T}_{12}}({r_2},{\theta _2}){\boldsymbol{e}_{\varepsilon
2}} = \dfrac{{{\boldsymbol{e}_{\varepsilon 2}}}}{{2\nu }}Re
\left\{ {k_v^2\left[ {v_r^{(2)}\dfrac{{\partial {\psi^{(1)
\ast }}}}{{\partial {r_2}}} + v_r^{(1)}\dfrac{{\partial
{\psi^{(2) \ast }}}}{{\partial {r_2}}} +
\dfrac{{v_\theta^{(2)}}}{{{r_2}}}\dfrac{{\partial
{\psi^{(1) \ast }}}}{{\partial {\theta_2}}} +
\dfrac{{v_\theta^{(1)}}}{{{r_2}}}\dfrac{{\partial
{\psi^{(2) \ast }}}}{{\partial {\theta_2}}}} \right.} \right.\nonumber\\ & \quad \left. {\left. { - \dfrac{{{\psi^{(2)
\ast }}(v_r^{(1)}\sin {\theta_2} + v_\theta^{(1)}\cos
{\theta_2}) + {\psi^{(1) \ast }}(v_r^{(2)}\sin {\theta_2} +
v_\theta^{(2)}\cos {\theta_2})}}{{{r_2}\sin {\theta_2}}}}
\right]} \right\}, \end{align}
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
\begin{align}{{T}_{12}}({r_2},{\theta _2}) = \frac{1}{{2\nu R_{20}^4}}\left[ {\displaystyle\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\limits_{n,m = 1}^\infty {{{[P_n^1({\mu_2})P_m^1({\mu_2})]}^/}M_{2nm}^{(2)}({r_2})} } \right],\end{align}where
\begin{align}
M_{1nm}^{(2)}({r_2}) &= Re \left\{ {\dfrac{{k_v^2R_{20}^4}}{{r_2^2}}\left[ {(n + 1){{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{n + 1}}[(n +
1){j_m}(k_v^ \ast {r_2}) - k_v^ \ast {r_2}j_m^/(k_v^ \ast
{r_2})]a_n^{(2)}\sum\limits_{k = 1}^\infty {b_k^{(1) \ast
}G_{mk}^ \ast } } \right.} \right.\nonumber\\ & \quad + {( -
1)^n}n{\left( {\dfrac{{{r_2}}}{d}} \right)^n}[nh_m^{(1)
\ast }({k_v}{r_2}) + k_v^ \ast {r_2}h_m^{(1)/{\ast}
}({k_v}{r_2})]b_m^{(2) \ast }\sum\limits_{k = 0}^\infty
{a_k^{(1)}{C_{kn}}\xi _1^{k + 1}} \nonumber\\ & \quad - n(n +
1)[{k_v}{r_2}{j_m}(k_v^ \ast {r_2})h_n^{(1)/}({k_v}{r_2}) +
k_v^ \ast {r_2}j_m^/(k_v^ \ast
{r_2})h_n^{(1)}({k_v}{r_2})]b_n^{(2)}\sum\limits_{k =
1}^\infty {b_k^{(1) \ast }G_{mk}^ \ast } \nonumber\\ & \quad \left.
{\left. { - n(n + 1)[k_v^ \ast
{r_2}{j_n}({k_v}{r_2})h_m^{(1)/{\ast} }({k_v}{r_2}) +
{k_v}{r_2}j_n^/({k_v}{r_2})h_m^{(1) \ast
}({k_v}{r_2})]b_m^{(2) \ast }\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{nk}}} } \right]} \right\},
\end{align}
\begin{align}
M_{2nm}^{(2)}({r_2}) &= Re \left\{
{\dfrac{{k_v^2R_{20}^4}}{{r_2^2}}\left[
{{j_m}({k_v}{r_2}){{\left( {\dfrac{{{R_{20}}}}{{{r_2}}}}
\right)}^{n + 1}}a_n^{(2) \ast }\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{mk}}} } \right.} \right. + {( -
1)^m}h_n^{(1)}({k_v}{r_2}){\left( {\dfrac{{{r_2}}}{d}}
\right)^m}b_n^{(2)}\sum\limits_{k = 0}^\infty {a_k^{(1)
\ast }{C_{km}}\xi _1^{k + 1}} \nonumber\\ & \quad \left. {\left. { +
[{k_v}{r_2}j_m^/({k_v}{r_2})h_n^{(1) \ast }({k_v}{r_2}) -
k_v^ \ast {r_2}{j_m}({k_v}{r_2})h_n^{(1)/{\ast}
}({k_v}{r_2})]b_n^{(2) \ast }\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{mk}}} } \right]} \right\}.
\end{align}
Representing 
$\boldsymbol{\varPsi }$ as
\begin{equation}\boldsymbol{\varPsi }({r_2},{\theta _2}) = \varPsi ({r_2},{\theta _2}){\boldsymbol{e}_{\varepsilon 2}},\end{equation}and using (B24)–(B35), one obtains
\begin{align}
&{\left( {\Delta_{r\theta }^{(2)} - \dfrac{1}{{r_2^2{{\sin
}^2}{\theta_2}}}} \right)^2}\varPsi ({r_2},{\theta _2}) = \dfrac{1}{{2\nu R_{20}^4}}\left\{
\displaystyle\sum_{\substack{n = 0\\ m = 1}}^\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})]} \right.\nonumber\\ &
\quad \left. + \sqrt {1 - \mu_2^2} \sum\limits_{n,m =
1}^\infty
{P_n^1({\mu_2})P_m^{1/}({\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})]} \right\}.
\end{align}
Appendix C. Solution of (2.13) and (2.15)
Equations (2.13) and (2.15) can be represented as
\begin{equation}{\left( {\Delta_{r\theta }^{(j)} - \frac{1}{{r_j^2{{\sin }^2}{\theta_j}}}} \right)^2}\varPsi ({r_j},{\theta _j}) = {E^{(j)}}({r_j},{\mu _j}),\quad j = 1,2,\end{equation}where
\begin{align} &
{E^{(j)}}({r_j},{\mu _j}) = \dfrac{1}{{2\nu R_{j0}^4}}\left\{ {\displaystyle\sum_{\substack{n = 0\\ m =
1}}^\infty
{{P_n}({\mu_j})P_m^1({\mu_j})[K_{1nm}^{(j)}({r_j}) +
L_{1nm}^{(j)}({r_j}) + M_{1nm}^{(j)}({r_j})]} } \right.\nonumber\\ &
\quad \left. { + \sqrt {1 - \mu_j^2} \sum\limits_{n,m =
1}^\infty
{P_n^{1/}({\mu_j})P_m^1({\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})]} } \right\}.
\end{align}
The function 
${E^{(j)}}({r_j},{\mu _j})$ can be expanded as follows:
\begin{equation}{E^{(j)}}({r_j},{\mu _j}) = \sum\limits_{l = 1}^\infty {E_l^{(j)}({r_j})P_l^1({\mu _j})} ,\end{equation}
where the expansion coefficients 
$E_l^{(j)}({r_j})$ are calculated by
\begin{equation}E_l^{(j)}({r_j}) = \frac{{2l + 1}}{{2l(l + 1)}}\int_{ - 1}^1 {{E^{(j)}}({r_j},{\mu _j})P_l^1({\mu _j})\,\textrm{d}{\mu _j}} .\end{equation}
With the help of (D5), (D9) and (D10), the calculation of 
$E_l^{(j)}({r_j})$ yields
\begin{align} &
E_l^{(j)}({r_j}) = \dfrac{{2l + 1}}{{2\nu R_{j0}^4\sqrt
{l(l + 1)} }}\left\{ {\displaystyle\sum_{n = 0}^\infty
{\displaystyle\sum_{\substack{m = |n - l|\\ m \ge 1}}^{n + l} {\dfrac{{\sqrt {m(m +
1)} C_{n0l0}^{m0}C_{n0l1}^{m1}}}{{2m + 1}}\{
K_{1nm}^{(j)}({r_j}) + L_{1nm}^{(j)}({r_j}) +
M_{1nm}^{(j)}({r_j})} } } \right.\nonumber\\ & \quad +
{n^2}[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})]\} \nonumber\\ & \quad -
\sum\limits_{n = 2}^\infty {\sum\limits_{m = 1}^{[{n / 2}]}
{\displaystyle\sum_{\substack{k = |n - 2m - l|\\ k \ge 1}}^{n - 2m + l} {\dfrac{{(2n - 4m +
1)\sqrt {k(k + 1)} C_{(n - 2m)0l0}^{k0}C_{(n -
2m)0l1}^{k1}}}{{2k + 1}}} } }\nonumber \\ & \quad { \times
[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})]} \},
\end{align}
where [ ] means the integer part of an expression in brackets.
Considering (C3), a solution to (C1) is sought as
\begin{equation}\varPsi ({r_j},{\theta _j}) = \sum\limits_{l = 1}^\infty {\varPsi _l^{(j)}({r_j})P_l^1({\mu _j})} .\end{equation}
Substituting (C3) and (C6) into (C1) gives the following equation for 
$\varPsi _l^{(j)}({r_j})$:
\begin{equation}\varPsi _l^{(j)IV} + \frac{4}{{{r_j}}}\varPsi _l^{(j)/{/}/} - \frac{{2l(l + 1)}}{{r_j^2}}\varPsi _l^{(j)/{/}} + \frac{{l(l + 1)({l^2} + l - 2)}}{{r_j^4}}\varPsi _l^{(j)} = E_l^{(j)}({r_j}),\end{equation}
where the prime denotes the derivative with respect to 
${r_j}$.
Equation (C7) is solved by the method of variation of constants (Boyce & DiPrima Reference Boyce and DiPrima2001). As a result, one obtains
\begin{equation}\varPsi _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^lC_{3l}^{(j)}({r_j}) + r_j^{l + 2}C_{4l}^{(j)}({r_j}),\end{equation}
where the functions 
$C_{1l}^{(j)}({r_j}) - C_{4l}^{(j)}({r_j})$ are calculated by
\begin{gather}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)\,\textrm{d}s} ,\end{gather}
\begin{gather}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)\,\textrm{d}s} ,\end{gather}
\begin{gather}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)\,\textrm{d}s} ,\end{gather}
\begin{gather}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)\,\textrm{d}s} .\end{gather}
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,
\begin{gather}{v_{Lr}}({r_j},{\theta _j}) = 0\;\textrm{at}\;{r_j} = {R_{j0}},\end{gather}
\begin{gather}{\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\;\textrm{at}\;{r_j} = {R_{j0}}.\end{gather}
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 
${\boldsymbol{v}_E}$, 
${v_{Sr}}({r_j},{\theta _j})$ and 
${v_{S\theta }}({r_j},{\theta _j})$ are the components of the Stokes drift velocity 
${\boldsymbol{v}_S}$ and 
${\tau _L}({r_j},{\theta _j})$ is the Lagrangian tangential stress.
The components of 
${\boldsymbol{v}_S}$ are calculated in Appendix E and given by (E15) and (E16). The components of 
${\boldsymbol{v}_E}$ are calculated by substituting (C6) into (2.8), which gives
\begin{gather}{v_{Er}}({r_j},{\theta _j}) ={-} \frac{1}{{{r_j}}}\sum\limits_{l = 1}^\infty {l(l + 1)\varPsi _l^{(j)}({r_j}){P_l}({\mu _j})} ,\end{gather}
\begin{gather}{v_{E\theta }}({r_j},{\theta _j}) ={-} \frac{1}{{{r_j}}}\sum\limits_{l = 1}^\infty {[\varPsi _l^{(j)}({r_j}) + {r_j}\varPsi _l^{(j)/}({r_j})]P_l^1({\mu _j})} ,\end{gather}
where, as follows from (C8) and the properties of the functions 
$C_{1l}^{(j)}({r_j}) - C_{4l}^{(j)}({r_j})$, 
$\Psi _l^{(j)/}({r_j})$ is calculated by
\begin{align}\varPsi _l^{(j)/}({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}}} + lr_j^{l - 1}C_{3l}^{(j)}({r_j}) + (l + 2)r_j^{l + 1}C_{4l}^{(j)}({r_j}).\end{align}Substituting (C15), (C16), (E15) and (E16) into (C13) and (C14), and using (C8), (C17) and (C18),
\begin{align}
\varPsi _l^{(j)/{/}}({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}}}\nonumber\\ &\quad + l(l - 1)r_j^{l -
2}C_{3l}^{(j)}({r_j}) + (l + 1)(l +
2)r_j^lC_{4l}^{(j)}({r_j}),\end{align}
one obtains equations for 
$C_{1l0}^{(j)} - C_{4l0}^{(j)}$,
\begin{gather}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)},\end{gather}
\begin{gather}({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)},\end{gather}where
\begin{gather}X_{1l}^{(j)} = \frac{{R_{j0}^l}}{{l(l + 1)}}V_{Srl}^{(j)}({R_{j0}}),\end{gather}
\begin{gather}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)/}({R_{j0}})].\end{gather}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:
\begin{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)},\end{equation}where
\begin{align}
X_{3l}^{(1)} &= \dfrac{{{{(d - {R_{20}})}^l}}}{{l(l +
1)}}V_{Srl}^{(1)}(d - {R_{20}})\nonumber\\ & \quad + \dfrac{1}{{2(2l - 1)(2l +
1)}}\int_{{R_{10}}}^{d - {R_{20}}} {[{{(d - {R_{20}})}^{2l
- 1}}{s^{3 - l}} - {s^{l + 2}}]E_l^{(1)}(s)\,\textrm{d}s} \nonumber\\ & \quad + \dfrac{1}{{2(2l + 1)(2l +
3)}}\int_{{R_{10}}}^{d - {R_{20}}} {\left[ {\dfrac{{{s^{l +
4}}}}{{{{(d - {R_{20}})}^2}}} - \dfrac{{{{(d -
{R_{20}})}^{2l + 1}}}}{{{s^{l - 1}}}}}
\right]E_l^{(1)}(s)\,\textrm{d}s} .
\end{align}
For the solution in the coordinates 
$({r_2},{\theta _2})$, we take a point with the coordinates 
${\theta _2} = {\rm \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
\begin{equation}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)},\end{equation}where
\begin{align}
X_{3l}^{(2)} &= \dfrac{{{{(d - {R_{10}})}^l}}}{{l(l +
1)}}V_{Srl}^{(2)}(d - {R_{10}})\nonumber\\ & \quad + \dfrac{1}{{2(2l - 1)(2l +
1)}}\int_{{R_{20}}}^{d - {R_{10}}} {[{{(d - {R_{10}})}^{2l
- 1}}{s^{3 - l}} - {s^{l + 2}}]E_l^{(2)}(s)\,\textrm{d}s} \nonumber\\ & \quad + \dfrac{1}{{2(2l + 1)(2l +
3)}}\int_{{R_{20}}}^{d - {R_{10}}} {\left[ {\dfrac{{{s^{l +
4}}}}{{{{(d - {R_{10}})}^2}}} - \dfrac{{{{(d -
{R_{10}})}^{2l + 1}}}}{{{s^{l - 1}}}}}
\right]E_l^{(2)}(s)\,\textrm{d}s} .
\end{align}
Let now us consider the Lagrangian tangential stress 
${\tau _L}({r_j},{\theta _j})$, which is defined by
\begin{align}
{\tau _L}({r_j},{\theta _j}) &= \eta \left[
{\dfrac{1}{{{r_j}}}\dfrac{{\partial
{v_{Lr}}({r_j},{\theta_j})}}{{\partial {\theta_j}}} +
\dfrac{{\partial {v_{L\theta
}}({r_j},{\theta_j})}}{{\partial {r_j}}} -
\dfrac{{{v_{L\theta }}({r_j},{\theta_j})}}{{{r_j}}}}\right]\nonumber\\
& = \eta \sum\limits_{l = 1}^\infty {\left[
{\dfrac{{2(1 - {l^2})}}{{r_j^{l + 1}}}C_{1l}^{(j)}({r_j}) -
\dfrac{{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})} \right.} \nonumber\\ & \quad \quad
\left. { - 2l(l + 2)r_j^lC_{4l}^{(j)}({r_j}) + \dfrac{{V_{Srl}^{(j)}({r_j}) - V_{S\theta
l}^{(j)}({r_j})}}{{{r_j}}} + V_{S\theta l}^{(j)/}({r_j})}
\right]P_l^1({\mu _j}). \end{align}
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_1},{\theta _1})$ must also be equal to zero in a small neighbourhood 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} = {\rm \pi},\;{r_2} = d - {R_{10}})$. As a result, we obtain two last equations,
\begin{align}&({l^2} - 1)C_{1l0}^{(j)}
+ \frac{{l(l + 2)}}{{{{(d - {R_{(3 - j)0}})}^2}}}C_{2l0}^{(j)}\nonumber\\
&\quad + ({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_{4l}^{(j)},\end{align}
where
\begin{align}
X_{4l}^{(j)} &= \dfrac{{{{(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)/}(d - {R_{(3 - j)0}})]\nonumber\\ & \quad +
\dfrac{{{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)\,\textrm{d}s} \nonumber\\ & \quad +
\dfrac{{l(l + 2)}}{{2(2l + 1)(2l + 3)}}\int_{{R_{j0}}}^{d -
{R_{(3 - j)0}}} {\left[ {\dfrac{{{s^{l + 4}}}}{{{{(d -
{R_{(3 - j)0}})}^2}}} - \dfrac{{{{(d - {R_{(3 - j)0}})}^{2l
+ 1}}}}{{{s^{l - 1}}}}} \right]E_l^{(j)}(s)\,\textrm{d}s} .
\end{align}
Combining (C19), (C20), (C23), (C25) and (C28), one finds
\begin{equation}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}}]}},\end{equation}
\begin{equation}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}}]}},\end{equation}
\begin{equation}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}}]}},\end{equation}
\begin{equation}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}}]}}.\end{equation}Appendix D. Mathematical formulas 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 )$:
\begin{gather}\sqrt {1 - {\mu ^2}} P_n^{1/}(\mu ) = n(n + 1){P_n}(\mu ) - \mu P_n^/(\mu ),\end{gather}
\begin{gather}\mu P_n^/(\mu ) = n{P_n}(\mu ) + P_{n - 1}^/(\mu ),\end{gather}
\begin{gather}(1 - {\mu ^2})P_n^/(\mu ) = \frac{{n(n + 1)}}{{2n + 1}}[{P_{n - 1}}(\mu ) - {P_{n + 1}}(\mu )],\end{gather}
\begin{gather}P_n^/(\mu ) = \sum\limits_{k = 1}^{[(n + 1)/2]} {(2n - 4k + 3)} {P_{n - 2k + 1}}(\mu ).\end{gather}With the help of these formulas, the following identities are derived:
\begin{equation}\sqrt {1 - {\mu ^2}} P_n^{1/}(\mu )P_m^1(\mu ) = {n^2}{P_n}(\mu )P_m^1(\mu ) - \sum\limits_{k = 1}^{[{n / 2}]} {(2n - 4k + 1)} {P_{n - 2k}}(\mu )P_m^1(\mu ),\end{equation}
\begin{equation}P_n^1(\mu )P_m^1(\mu ) = \frac{{m(m + 1)}}{{2m + 1}}\sum\limits_{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 )]} .\end{equation} To calculate integrals with products of 
${P_n}(\mu )$ and 
$P_m^1(\mu )$, we use the following identity:
\begin{align}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\limits_{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}^{l0}C_{{l_1}{m_1}{l_2}{m_2}}^{l({m_1} + {m_2})}P_l^{{m_1} + {m_2}}(\mu )} .\end{align}This identity follows from (9) of § 5.6 of the book by Varshalovich et al. (Reference Varshalovich, Moskalev and Khersonskii1988).
From (D7) and the orthogonality condition for 
$P_n^1(\mu )$,
\begin{equation}\int_{ - 1}^1 {P_n^1(\mu )P_m^1(\mu )\,\textrm{d}\mu } = \frac{{2n(n + 1)}}{{2n + 1}}{\delta _{nm}},\end{equation}one obtains the following integrals:
\begin{gather}\int_{ - 1}^1 {{P_n}(\mu )P_m^1(\mu )P_l^1(\mu )\,\textrm{d}\mu } = 2\sqrt {l(l + 1)} \displaystyle\sum_{\substack{m = |n - l|\\ m \ge 1}}^{n + l} {\frac{{\sqrt {m(m + 1)} C_{n0l0}^{m0}C_{n0l1}^{m1}}}{{2m + 1}}} ,\end{gather}
\begin{gather}\int_{ - 1}^1 {{P_{n - 2k}}(\mu )P_m^1(\mu )P_l^1(\mu )\,\textrm{d}\mu } = 2\sqrt {l(l + 1)} \displaystyle\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}}} .\end{gather} From (D7) and the orthogonality condition for 
${P_n}(\mu )$,
\begin{equation}\int_{ - 1}^1 {{P_n}(\mu ){P_m}(\mu )\,\textrm{d}\mu } = \frac{2}{{2n + 1}}{\delta _{nm}},\end{equation}one obtains
\begin{gather}\int_{ - 1}^1 {{P_l}(\mu ){P_n}(\mu ){P_m}(\mu )\,\textrm{d}\mu } = \sum\limits_{m = |l - n|}^{l + n} {\frac{2}{{2m + 1}}{{(C_{l0n0}^{m0})}^2}} ,\end{gather}
\begin{gather}\int_{ - 1}^1 {{P_l}(\mu ){P_{m - 1}}(\mu ){P_{n - 2k + 1}}(\mu )\,\textrm{d}\mu } = \sum\limits_{m - 1 = |l - (n - 2k + 1)|}^{l + n - 2k + 1} {\frac{2}{{2m - 1}}{{(C_{l0(n - 2k + 1)0}^{(m - 1)0})}^2}} ,\end{gather}
\begin{gather}\int_{ - 1}^1 {{P_l}(\mu ){P_{m + 1}}(\mu ){P_{n - 2k + 1}}(\mu )\,\textrm{d}\mu } = \sum\limits_{m + 1 = |l - (n - 2k + 1)|}^{l + n - 2k + 1} {\frac{2}{{2m + 3}}{{(C_{l0(n - 2k + 1)0}^{(m + 1)0})}^2}} .\end{gather}Appendix E. Calculation of the Stokes drift velocity
The Stokes drift velocity is calculated by (Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a)
\begin{equation}{\boldsymbol{v}_S} = \frac{1}{{2\omega }}Re \{ i(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol{v}^ \ast }\} .\end{equation}
Substitution of 
$\boldsymbol{v} = {\boldsymbol{v}^{(1)}} + {\boldsymbol{v}^{(2)}}$ into (E1) yields
\begin{align}{\boldsymbol{v}_S} = \frac{1}{{2\omega }}Re \{ i({\boldsymbol{v}^{(1)}}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol{v}^{(1) \ast }} + i({\boldsymbol{v}^{(2)}}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol{v}^{(2) \ast }} + i({\boldsymbol{v}^{(1)}}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol{v}^{(2) \ast }} + i({\boldsymbol{v}^{(2)}}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol{v}^{(1) \ast }}\} .\end{align} To represent 
${\boldsymbol{v}_S}$ in the coordinates 
$({r_1},{\theta _1})$, 
$\boldsymbol{\nabla }$ is taken in the coordinates 
$({r_1},{\theta _1})$, 
${\boldsymbol{v}^{(1)}}$ is calculated by (2.6) and (2.7) and 
${\boldsymbol{v}^{(2)}}$ is calculated by (A55) and (A56). As a result, we obtain the following equations for the components of 
${\boldsymbol{v}_S}$:
\begin{align}
{v_{Sr}}({r_1},{\theta _1}) &= \dfrac{1}{{2\omega }}Re
\left\{ {i(v_r^{(1)} + v_r^{(2)})\dfrac{{\partial (v_r^{(1)
\ast } + v_r^{(2) \ast })}}{{\partial {r_1}}} +
\dfrac{{i(v_\theta^{(1)} +
v_\theta^{(2)})}}{{{r_1}}}\dfrac{{\partial (v_r^{(1) \ast }
+ v_r^{(2) \ast })}}{{\partial {\theta_1}}}} \right\}\nonumber\\ &= \dfrac{1}{{2\omega r_1^3}}\left[ {\sum\limits_{n,m
= 0}^\infty
{S_{1nm}^{(1)}({r_1}){P_n}({\mu_1}){P_m}({\mu_1})} +
\sum\limits_{n,m = 1}^\infty
{S_{2nm}^{(1)}({r_1})P_n^1({\mu_1})P_m^1({\mu_1})} }
\right], \end{align}
\begin{align}
{v_{S\theta }}({r_1},{\theta _1}) &= \dfrac{1}{{2\omega }}Re
\left\{ {i(v_r^{(1)} + v_r^{(2)})\dfrac{{\partial
(v_\theta^{(1) \ast } + v_\theta^{(2) \ast })}}{{\partial
{r_1}}} + \dfrac{{i(v_\theta^{(1)} +
v_\theta^{(2)})}}{{{r_1}}}\dfrac{{\partial (v_\theta^{(1)
\ast } + v_\theta^{(2) \ast })}}{{\partial {\theta_1}}}} \right.\nonumber\\ & \quad \left. { + \dfrac{{i(v_r^{(1) \ast } +
v_r^{(2) \ast })(v_\theta^{(1)} + v_\theta^{(2)})}}{{{r_1}}}} \right\}\nonumber\\ & =
\dfrac{1}{{2\omega r_1^3}}\left[
{\displaystyle\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\limits_{n,m = 1}^\infty
{Q_{2nm}^{(1)}({r_1})P_n^1({\mu_1})P_m^{1/}({\mu_1})} }
\right], \end{align}
where
\begin{align}
S_{1nm}^{(1)}({r_1}) &= {\mathop{\rm Im}\nolimits} \left\{
{\left[ {(n + 1)a_n^{(1)}{{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{n + 1}} + n(n +
1)b_n^{(1)}h_n^{(1)}({k_v}{r_1})} \right.} \right.\nonumber\\ &
\quad \left. { - n{{\left( {\dfrac{{{r_1}}}{d}}
\right)}^n}\sum\limits_{k = 0}^\infty {a_k^{(2)}{{( -
1)}^k}{C_{kn}}\xi_2^{k + 1}} + n(n +
1){j_n}({k_v}{r_1})\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{nk}}} } \right]\nonumber\\ & \quad \times \left[ (m +
1)(m + 2)a_m^{(1) \ast }{{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{m + 1}} + m(m +
1)b_m^{(1) \ast }[h_m^{(1)}({k_v}{r_1}) \right.\nonumber\\ &
\quad \left. \left. - {k_v}{r_1}h_m^{(1)/}({k_v}{r_1})]^ \ast + m(m - 1){{\left(
{\dfrac{{{r_1}}}{d}} \right)}^m}\sum\limits_{k = 0}^\infty
{a_k^{(2) \ast }{{( - 1)}^k}{C_{km}}\xi_2^{k + 1}}\right.\right.\nonumber\\ &
\quad \left. \left. + m(m + 1){{[{j_m}({k_v}{r_1}) - {k_v}{r_1}j_m^/({k_v}{r_1})]}^
\ast }\sum\limits_{k = 1}^\infty {{{(b_k^{(2)}{H_{mk}})}^ \ast }}\right] \right\},\quad n,m \ge 0,
\end{align}
\begin{align}
S_{2nm}^{(1)}({r_1}) &= {\mathop{\rm Im}\nolimits} \left\{
{\left[ {(n + 1)a_n^{(1) \ast }{{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{n + 1}} + n(n +
1)b_n^{(1) \ast }h_n^{(1) \ast }({k_v}{r_1})} \right.}
\right.\nonumber\\ & \quad \left. { - n{{\left( {\dfrac{{{r_1}}}{d}}
\right)}^n}\sum\limits_{k = 0}^\infty {a_k^{(2) \ast }{{( -
1)}^k}{C_{kn}}\xi_2^{k + 1}} + n(n + 1){j_n}(k_v^ \ast
{r_1})\sum\limits_{k = 1}^\infty {{{(b_k^{(2)}{H_{nk}})}^
\ast }} } \right]\nonumber\\ & \quad \times \left[
{a_m^{(1)}{{\left( {\dfrac{{{R_{10}}}}{{{r_1}}}}
\right)}^{m + 1}} - b_m^{(1)}[h_m^{(1)}({k_v}{r_1}) +
{k_v}{r_1}h_m^{(1)/}({k_v}{r_1})]} \right.\nonumber\\ & \quad \left.
{\left. { + {{\left( {\dfrac{{{r_1}}}{d}}
\right)}^m}\sum\limits_{k = 0}^\infty {a_k^{(2)}{{( -
1)}^k}{C_{km}}\xi_2^{k + 1}} - [{j_m}({k_v}{r_1}) +
{k_v}{r_1}j_m^/({k_v}{r_1})]\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{mk}}} } \right]} \right\},\quad n,m \ge 1,
\end{align}
\begin{align}
Q_{1nm}^{(1)}({r_1}) &= {\mathop{\rm Im}\nolimits} \left\{
{\left[ {(n + 1)a_n^{(1)}{{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{n + 1}} + n(n +
1)b_n^{(1)}h_n^{(1)}({k_v}{r_1})} \right.} \right.\nonumber\\ &
\quad \left. { - n{{\left( {\dfrac{{{r_1}}}{d}}
\right)}^n}\sum\limits_{k = 0}^\infty {a_k^{(2)}{{( -
1)}^k}{C_{kn}}\xi_2^{k + 1}} + n(n +
1){j_n}({k_v}{r_1})\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{nk}}} } \right]\nonumber\\ & \quad \times \left[-
(m + 3)a_m^{(1) \ast }{{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{m + 1}} + b_m^{(1)
\ast }{[2h_m^{(1)}({k_v}{r_1}) -
{{({k_v}{r_1})}^2}h_m^{(1)/{/}}({k_v}{r_1})]}^ \ast
\right.\nonumber\\ & \quad \left. \left. + (m - 2){{\left(
{\dfrac{{{r_1}}}{d}} \right)}^m}\sum\limits_{k = 0}^\infty
{a_k^{(2) \ast }{{( - 1)}^k}{C_{km}}\xi_2^{k + 1}} +
[2{j_m}({k_v}{r_1})\right.\right.\nonumber\\ & \left. \left. \quad - {{({k_v}{r_1})}^2}j_m^{/{/}}({k_v}{r_1})]^\ast
\sum\limits_{k = 1}^\infty {{{(b_k^{(2)}{H_{mk}})}^ \ast
}} \right]\right\},\quad n \ge 0,\;m \ge 1,
\end{align}
\begin{align}
Q_{2nm}^{(1)}({r_1}) &= {\mathop{\rm Im}\nolimits} \left\{
{\left\{ {a_n^{(1)}{{\left( {\dfrac{{{R_{10}}}}{{{r_1}}}}
\right)}^{n + 1}} - b_n^{(1)}[h_n^{(1)}({k_v}{r_1}) +
{k_v}{r_1}h_n^{(1)/}({k_v}{r_1})]} \right.} \right.\nonumber\\ &
\quad + \left. {{\left( {\dfrac{{{r_1}}}{d}}
\right)}^n}\sum\limits_{k = 0}^\infty {a_k^{(2)}{{( -
1)}^k}{C_{kn}}\xi_2^{k + 1}} - [{j_n}({k_v}{r_1}) +
{k_v}{r_1}j_n^/({k_v}{r_1})]\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{nk}}} \right\}\nonumber\\ & \quad \times \left\{
a_m^{(1) \ast }{{\left( {\dfrac{{{R_{10}}}}{{{r_1}}}}
\right)}^{m + 1}} - b_m^{(1) \ast }{{[h_m^{(1)}({k_v}{r_1})
+ {k_v}{r_1}h_m^{(1)/}({k_v}{r_1})]}^ \ast } \right.\nonumber\\ &
\quad \left. \left. + {{\left( {\dfrac{{{r_1}}}{d}}
\right)}^m}\sum\limits_{k = 0}^\infty {a_k^{(2) \ast }{{( -
1)}^k}{C_{km}}\xi_2^{k + 1}}
\right.\right.\nonumber\\ & \left. \left. \quad - {{[{j_m}({k_v}{r_1}) +
{k_v}{r_1}j_m^/({k_v}{r_1})]}^ \ast }\sum\limits_{k =
1}^\infty {{{(b_k^{(2)}{H_{mk}})}^ \ast }}\right\}
\right\},\quad n,m \ge 1. \end{align}
To represent (E2) in the coordinates 
$({r_2},{\theta _2})$, 
$\boldsymbol{\nabla }$ is taken in the coordinates 
$({r_2},{\theta _2})$, 
${\boldsymbol{v}^{(1)}}$ is calculated by (A52) and (A53) and 
${\boldsymbol{v}^{(2)}}$ is calculated by (2.6) and (2.7). As a result, we obtain
\begin{align} {v_{Sr}}({r_2},{\theta _2}) &= \dfrac{1}{{2\omega }}Re \left\{ {i(v_r^{(1)} + v_r^{(2)})\dfrac{{\partial (v_r^{(1) \ast } + v_r^{(2) \ast })}}{{\partial {r_2}}} + \dfrac{{i(v_\theta^{(1)} + v_\theta^{(2)})}}{{{r_2}}}\dfrac{{\partial (v_r^{(1) \ast } + v_r^{(2) \ast })}}{{\partial {\theta_2}}}} \right\}\nonumber\\ &= \dfrac{1}{{2\omega r_2^3}}\left[ {\sum\limits_{n,m = 0}^\infty {S_{1nm}^{(2)}({r_2}){P_n}({\mu_2}){P_m}({\mu_2})} + \sum\limits_{n,m = 1}^\infty {S_{2nm}^{(2)}({r_2})P_n^1({\mu_2})P_m^1({\mu_2})} } \right],\end{align}
\begin{align} {v_{S\theta
}}({r_2},{\theta _2}) &= \dfrac{1}{{2\omega }}Re \left\{
{i(v_r^{(1)} + v_r^{(2)})\dfrac{{\partial (v_\theta^{(1)
\ast } + v_\theta^{(2) \ast })}}{{\partial {r_2}}} +
\dfrac{{i(v_\theta^{(1)} +
v_\theta^{(2)})}}{{{r_2}}}\dfrac{{\partial (v_\theta^{(1)
\ast } + v_\theta^{(2) \ast })}}{{\partial {\theta_2}}}}
\right.\nonumber\\ & \quad \left. { + \dfrac{{i(v_r^{(1)
\ast } + v_r^{(2) \ast })(v_\theta^{(1)} +
v_\theta^{(2)})}}{{{r_2}}}} \right\}\nonumber\\ &=
\dfrac{1}{{2\omega r_2^3}}\left[
\displaystyle\sum_{\substack{n = 0\\ m = 1}}^\infty
{Q_{1nm}^{(2)}({r_2}){P_n}({\mu_2})P_m^1({\mu_2})}\right.\nonumber\\
&\quad \left. + \sqrt {1 - \mu_2^2} \sum\limits_{n,m =
1}^\infty
{Q_{2nm}^{(2)}({r_2})P_n^1({\mu_2})P_m^{1/}({\mu_2})}
\right],
\end{align}where
\begin{align}
S_{1nm}^{(2)}({r_2}) &= {\mathop{\rm Im}\nolimits} \left\{
{\left[ {(n + 1)a_n^{(2)}{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{n + 1}} + n(n +
1)b_n^{(2)}h_n^{(1)}({k_v}{r_2})} \right.}
\right.\nonumber\\ & \quad \left. { - {{( - 1)}^n}n{{\left(
{\dfrac{{{r_2}}}{d}} \right)}^n}\sum\limits_{k = 0}^\infty
{a_k^{(1)}{C_{kn}}\xi_1^{k + 1}} + n(n +
1){j_n}({k_v}{r_2})\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{nk}}} } \right]\nonumber\\
& \quad \times
\left[ (m + 1)(m + 2)a_m^{(2) \ast }{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{m + 1}}\right. \nonumber\\
&\left. \quad + m(m + 1)b_m^{(2) \ast }{{[h_m^{(1)}({k_v}{r_2}) -
{k_v}{r_2}h_m^{(1)/}({k_v}{r_2})]}^ \ast }
\right.\nonumber\\ & \quad + \left. \left. {{( -
1)}^m}m(m - 1){{\left( {\dfrac{{{r_2}}}{d}}
\right)}^m}\sum\limits_{k = 0}^\infty {a_k^{(1) \ast
}{C_{km}}\xi_1^{k + 1}}\right.\right.\nonumber\\
&\left.\left.\quad + m(m + 1){{[{j_m}({k_v}{r_2}) -
{k_v}{r_2}j_m^/({k_v}{r_2})]}^ \ast }\sum\limits_{k =
1}^\infty {{{(b_k^{(1)}{G_{mk}})}^ \ast }} \right]
\right\},\quad n,m \ge 0,
\end{align}
\begin{align}
S_{2nm}^{(2)}({r_2}) &= {\mathop{\rm Im}\nolimits} \left\{
{\left[ {(n + 1)a_n^{(2) \ast }{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{n + 1}} + n(n +
1)b_n^{(2) \ast }h_n^{(1) \ast }({k_v}{r_2})} \right.}
\right.\nonumber\\ & \quad \left. { - {{( - 1)}^n}n{{\left(
{\dfrac{{{r_2}}}{d}} \right)}^n}\sum\limits_{k = 0}^\infty
{a_k^{(1) \ast }{C_{kn}}\xi_1^{k + 1}} + n(n + 1){j_n}(k_v^
\ast {r_2})\sum\limits_{k = 1}^\infty
{{{(b_k^{(1)}{G_{nk}})}^ \ast }} } \right]\nonumber\\
& \quad \times \left[ {a_m^{(2)}{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{m + 1}} -
b_m^{(2)}[h_m^{(1)}({k_v}{r_2}) +
{k_v}{r_2}h_m^{(1)/}({k_v}{r_2})]} \right.\nonumber\\
& \quad \left. \left. + {{( - 1)}^m}{{\left(
{\dfrac{{{r_2}}}{d}} \right)}^m}\sum\limits_{k = 0}^\infty
{a_k^{(1)}{C_{km}}\xi_1^{k + 1}}\right.\right.\nonumber\\
&\left. \left.\quad - [{j_m}({k_v}{r_2}) +
{k_v}{r_2}j_m^/({k_v}{r_2})]\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{mk}}} \right] \right\},\quad n,m \ge 1,
\end{align}
\begin{align} Q_{1nm}^{(2)}({r_2}) &=
{\mathop{\rm Im}\nolimits} \left\{ {\left[ {(n +
1)a_n^{(2)}{{\left( {\dfrac{{{R_{20}}}}{{{r_2}}}}
\right)}^{n + 1}} + n(n + 1)b_n^{(2)}h_n^{(1)}({k_v}{r_2})}
\right.} \right.\nonumber\\ & \quad \left. { - {{( -
1)}^n}n{{\left( {\dfrac{{{r_2}}}{d}}
\right)}^n}\sum\limits_{k = 0}^\infty
{a_k^{(1)}{C_{kn}}\xi_1^{k + 1}} + n(n +
1){j_n}({k_v}{r_2})\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{nk}}} } \right]\nonumber\\ & \quad \times
\left[ { - (m + 3)a_m^{(2) \ast }{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{m + 1}} + b_m^{(2)
\ast }{{[2h_m^{(1)}({k_v}{r_2}) -
{{({k_v}{r_2})}^2}h_m^{(1)/{/}}({k_v}{r_2})]}^ \ast }}
\right.\nonumber\\ & \quad \left. \left. + {{( -
1)}^m}(m - 2){{\left( {\dfrac{{{r_2}}}{d}}
\right)}^m}\sum\limits_{k = 0}^\infty {a_k^{(1) \ast
}{C_{km}}\xi_1^{k + 1}}\right.\right.\nonumber\\
&\left.\left. \quad + {{[2{j_m}({k_v}{r_2}) -
{{({k_v}{r_2})}^2}j_m^{/{/}}({k_v}{r_2})]}^ \ast
}\sum\limits_{k = 1}^\infty {{{(b_k^{(1)}{G_{mk}})}^ \ast
}} \right] \right\},\quad n \ge 0,m \ge 1,
\end{align}
\begin{align} Q_{2nm}^{(2)}({r_2}) &=
{\mathop{\rm Im}\nolimits} \left\{ {\left\{
{a_n^{(2)}{{\left( {\dfrac{{{R_{20}}}}{{{r_2}}}}
\right)}^{n + 1}} - b_n^{(2)}[h_n^{(1)}({k_v}{r_2}) +
{k_v}{r_2}h_n^{(1)/}({k_v}{r_2})]} \right.}
\right.\nonumber\\ & \quad \left. { + {{( - 1)}^n}{{\left(
{\dfrac{{{r_2}}}{d}} \right)}^n}\sum\limits_{k = 0}^\infty
{a_k^{(1)}{C_{kn}}\xi_1^{k + 1}} - [{j_n}({k_v}{r_2}) +
{k_v}{r_2}j_n^/({k_v}{r_2})]\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{nk}}} } \right\}\nonumber\\ & \quad \times
\left\{ {a_m^{(2) \ast }{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{m + 1}} - b_m^{(2)
\ast }{{[h_m^{(1)}({k_v}{r_2}) +
{k_v}{r_2}h_m^{(1)/}({k_v}{r_2})]}^ \ast }}
\right.\nonumber\\
& \quad \left. \left. + {{( - 1)}^m}{{\left( {\dfrac{{{r_2}}}{d}}
\right)}^m}\sum\limits_{k = 0}^\infty {a_k^{(1) \ast}{C_{km}}\xi_1^{k + 1}} - [{j_m}({k_v}{r_2})\right. \right.\nonumber\\
&\left.\left. \quad + {k_v}{r_2}j_m^/({k_v}{r_2})]^{\ast} \sum\limits_{k = 1}^\infty {{{(b_k^{(1)}{G_{mk}})}^ \ast }} \right\}
\right\},\quad n,m \ge 1.\end{align}Equations (E3), (E4), (E9) and (E10) can be represented as
\begin{align}{v_{Sr}}({r_j},{\theta _j}) = \sum\limits_{l = 0}^\infty {V_{Srl}^{(j)}({r_j}){P_l}({\mu _j})} ,\end{align}
\begin{align}{v_{S\theta }}({r_j},{\theta _j}) = \sum\limits_{l = 1}^\infty {V_{S\theta l}^{(j)}({r_j})P_l^1({\mu _j})} ,\end{align}
where 
$V_{Srl}^{(j)}({r_j})$ is calculated by
\begin{align}V_{Srl}^{(j)}({r_j}) &=
\frac{{2l + 1}}{{4\omega r_j^3}}\int_{ - 1}^1 \left[
\sum\limits_{n,m = 0}^\infty
{S_{1nm}^{(j)}({r_j}){P_n}({\mu_j}){P_m}({\mu_j})}\right.\nonumber\\ &\left. \quad +
\sum\limits_{n,m = 1}^\infty
{S_{2nm}^{(j)}({r_j})P_n^1({\mu_j})P_m^1({\mu_j})}
\right]{P_l}({\mu _j})\,\textrm{d}{\mu _j}
,\end{align}
while 
$V_{S\theta l}^{(j)}({r_j})$ is calculated by
\begin{align} V_{S\theta l}^{(j)}({r_j}) &= \dfrac{{2l + 1}}{{4l(l +
1)\omega r_j^3}}\nonumber\\
& \quad \times \int_{ - 1}^1 \left\{
\displaystyle\sum_{\substack{n = 0\\ m = 1}}^\infty
{Q_{1nm}^{(j)}({r_j}){P_n}({\mu_j})P_m^1({\mu_j})}\right.\nonumber\\
&\left. \quad + \sqrt {1 - \mu_j^2} \sum\limits_{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})\,\textrm{d}{\mu _j} .
\end{align} 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\theta l}^{(j)}({r_j})$,
\begin{align}
V_{Srl}^{(j)}({r_j}) &= \dfrac{{2l + 1}}{{2\omega
r_j^3}}\left\{ {\sum\limits_{n = 0}^\infty {\sum\limits_{m
= |l - n|}^{l + n} {\dfrac{{{{(C_{l0n0}^{m0})}^2}}}{{2m +
1}}} S_{1nm}^{(j)}({r_j})} } \right.\nonumber\\
& \quad \left. + \sum\limits_{n = 1}^\infty \sum\limits_{m = 1}^{[(n +
1)/2]} \sum\limits_{k = |l - (n - 2m + 1)|}^{l + n - 2m +
1} \dfrac{{(2n - 4m + 3){{(C_{l0(n - 2m + 1)0}^{k0})}^2}}}{{2k + 1}}\right.\nonumber\\
&\left. \quad \left[ {\dfrac{{(k + 1)(k +
2)}}{{2k + 3}}S_{2n(k + 1)}^{(j)}({r_j}) - \dfrac{{k(k -
1)}}{{2k - 1}}S_{2n(k - 1)}^{(j)}({r_j})} \right] \right\},
\end{align}
\begin{align}
V_{S\theta l}^{(j)}({r_j}) &= \dfrac{{2l + 1}}{{2\sqrt {l(l
+ 1)} \omega r_j^3}}\left\{ {\sum\limits_{n = 0}^\infty
{\displaystyle\sum_{\substack{m = |n - l|\\m \ge 1}}^{n +
l} {\dfrac{{\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.\nonumber\\
& \quad \left. { + \sum\limits_{n = 2}^\infty {\sum\limits_{m = 1}^{[n/2]}
{\displaystyle\sum_{\substack{k = |n - 2m - l|\\ k \ge
1}}^{n - 2m + l} {\dfrac{{(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\}.
\end{align}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)/}({r_j})$ is used; see (C22) and (C29). It is calculated by
\begin{align}
V_{S\theta l}^{(j)/}({r_j}) & ={-} \dfrac{{3V_{S\theta
l}^{(j)}({r_j})}}{{{r_j}}}\nonumber\\
&\quad + \dfrac{{2l + 1}}{{2\sqrt {l(l
+ 1)} \omega r_j^3}}\left\{ {\sum\limits_{n = 0}^\infty
{\displaystyle\sum_{\substack{m = |n - l|\\ m \ge 1}}^{n +
l} {\dfrac{{\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.\nonumber\\
& \quad \left. { + \sum\limits_{n = 2}^\infty
{\sum\limits_{m = 1}^{[n/2]}
{\displaystyle\sum_{\substack{k = |n - 2m - l|\\ k \ge
1}}^{n - 2m + l} {\dfrac{{(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\}\!,
\end{align}where
\begin{align}
Q_{1nm}^{(1)/}({r_1}) &= {\mathop{\rm Im}\nolimits} \left\{
{\left[ { - \dfrac{{{{(n +
1)}^2}}}{{{R_{10}}}}a_n^{(1)}{{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{n + 2}} + n(n +
1)b_n^{(1)}{k_v}h_n^{(1)/}({k_v}{r_1})} \right.} \right.\nonumber\\
& \quad \left. { - \dfrac{{{n^2}}}{d}{{\left(
{\dfrac{{{r_1}}}{d}} \right)}^{n - 1}}\sum\limits_{k =
0}^\infty {a_k^{(2)}{{( - 1)}^k}{C_{kn}}\xi_2^{k + 1}} +
n(n + 1){k_v}j_n^/({k_v}{r_1})\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{nk}}} } \right]\nonumber\\
& \quad \times \left[ { - (m + 3)a_m^{(1) \ast }{{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{m + 1}} + b_m^{(1)
\ast }{{[2h_m^{(1)}({k_v}{r_1}) -
{{({k_v}{r_1})}^2}h_m^{(1)/{/}}({k_v}{r_1})]}^ \ast }}
\right.\nonumber\\
& \left. \quad + (m - 2){{\left(
{\dfrac{{{r_1}}}{d}} \right)}^m}\sum\limits_{k = 0}^\infty
{a_k^{(2) \ast }{{( - 1)}^k}{C_{km}}\xi_2^{k + 1}}\right. \nonumber\\
&\left. \quad + {{[2{j_m}({k_v}{r_1}) -
{{({k_v}{r_1})}^2}j_m^{/{/}}({k_v}{r_1})]}^ \ast
}\sum\limits_{k = 1}^\infty {{{(b_k^{(2)}{H_{mk}})}^ \ast
}} \right]\nonumber\\
& \quad + \left[ {(n + 1)a_n^{(1)}{{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{n + 1}} + n(n +
1)b_n^{(1)}h_n^{(1)}({k_v}{r_1})} \right. \nonumber\\
& \quad \left. { - n{{\left( {\dfrac{{{r_1}}}{d}}
\right)}^n}\sum\limits_{k = 0}^\infty {a_k^{(2)}{{( -
1)}^k}{C_{kn}}\xi_2^{k + 1}} + n(n +
1){j_n}({k_v}{r_1})\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{nk}}} } \right]\nonumber\\
& \quad \times \left[ \dfrac{{(m + 1)(m + 3)}}{{{R_{10}}}}a_m^{(1) \ast
}{{\left( {\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{m + 2}}\right. \nonumber\\
&\left.\quad + b_m^{(1) \ast }k_v^ \ast {{[2h_m^{(1)/}({k_v}{r_1}) -
2{k_v}{r_1}h_m^{(1)/{/}}({k_v}{r_1}) -
{{({k_v}{r_1})}^2}h_m^{(1)/{/}/}({k_v}{r_1})]}^ \ast }\right.\nonumber\\
& \quad + \dfrac{{m(m - 2)}}{d}{\left(
{\dfrac{{{r_1}}}{d}} \right)^{m - 1}}\sum\limits_{k =
0}^\infty {a_k^{(2) \ast }{{( - 1)}^k}{C_{km}}\xi _2^{k +
1}} \nonumber\\
&\quad \left. \left. +\, k_v^ \ast
[2j_m^/({k_v}{r_1}) - 2{k_v}{r_1}j_m^{/{/}}({k_v}{r_1})\right.\right.\nonumber\\
&\left. \left. - {{({k_v}{r_1})}^2}j_m^{/{/}/}({k_v}{r_1})]^\ast
\sum\limits_{k = 1}^\infty {{{(b_k^{(2)}{H_{mk}})}^ \ast
}}\right] \right\},\quad n \ge 0,\;m \ge 1, \end{align}
\begin{align}
Q_{2nm}^{(1)/}({r_1}) &= {\mathop{\rm Im}\nolimits} \left\{
{\left\{ { - \dfrac{{n + 1}}{{{R_{10}}}}a_n^{(1)}{{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{n + 2}} -
b_n^{(1)}{k_v}[2h_n^{(1)/}({k_v}{r_1}) +
{k_v}{r_1}h_n^{(1)/{/}}({k_v}{r_1})]} \right.} \right.\nonumber\\
&\quad + \dfrac{n}{d}\left. {{{\left( {\dfrac{{{r_1}}}{d}}
\right)}^{n - 1}}\sum\limits_{k = 0}^\infty {a_k^{(2)}{{( -
1)}^k}{C_{kn}}\xi_2^{k + 1}} - {k_v}[2j_n^/({k_v}{r_1}) +
{k_v}{r_1}j_n^{/{/}}({k_v}{r_1})]\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{nk}}} } \right\}\nonumber\\
& \quad \times \left\{
{a_m^{(1) \ast }{{\left( {\dfrac{{{R_{10}}}}{{{r_1}}}}
\right)}^{m + 1}} - b_m^{(1) \ast }{{[h_m^{(1)}({k_v}{r_1})
+ {k_v}{r_1}h_m^{(1)/}({k_v}{r_1})]}^ \ast }} \right.\nonumber\\
& \quad \left. { + {{\left( {\dfrac{{{r_1}}}{d}}
\right)}^m}\sum\limits_{k = 0}^\infty {a_k^{(2) \ast }{{( -
1)}^k}{C_{km}}\xi_2^{k + 1}} - {{[{j_m}({k_v}{r_1}) +
{k_v}{r_1}j_m^/({k_v}{r_1})]}^ \ast }\sum\limits_{k =
1}^\infty {{{(b_k^{(2)}{H_{mk}})}^ \ast }} } \right\}\nonumber\\
& \quad + \left\{ {a_n^{(1)}{{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{n + 1}} -
b_n^{(1)}[h_n^{(1)}({k_v}{r_1}) +
{k_v}{r_1}h_n^{(1)/}({k_v}{r_1})]} \right.\nonumber\\
& \quad + \left. {{{\left( {\dfrac{{{r_1}}}{d}}
\right)}^n}\sum\limits_{k = 0}^\infty {a_k^{(2)}{{( -
1)}^k}{C_{kn}}\xi_2^{k + 1}} - [{j_n}({k_v}{r_1}) +
{k_v}{r_1}j_n^/({k_v}{r_1})]\sum\limits_{k = 1}^\infty
{b_k^{(2)}{H_{nk}}} } \right\}\nonumber\\
& \quad \times \left\{ { -
\dfrac{{m + 1}}{{{R_{10}}}}a_m^{(1) \ast }{{\left(
{\dfrac{{{R_{10}}}}{{{r_1}}}} \right)}^{m + 2}} - b_m^{(1)
\ast }k_v^ \ast {{[2h_m^{(1)/}({k_v}{r_1}) +
{k_v}{r_1}h_m^{(1)/{/}}({k_v}{r_1})]}^ \ast }} \right.\nonumber\\
& \quad + \dfrac{m}{d}{\left( {\dfrac{{{r_1}}}{d}} \right)^{m
- 1}}\sum\limits_{k = 0}^\infty {a_k^{(2) \ast }{{( -
1)}^k}{C_{km}}\xi _2^{k + 1}} \nonumber\\
& \quad \left. {\left. { - k_v^ \ast {{[2j_m^/({k_v}{r_1}) +
{k_v}{r_1}j_m^{/{/}}({k_v}{r_1})]}^ \ast }\sum\limits_{k =
1}^\infty {{{(b_k^{(2)}{H_{mk}})}^ \ast }} } \right\}}
\right\},\quad n,m \ge 1, \end{align}
\begin{align}
Q_{1nm}^{(2)/}({r_2}) &= {\mathop{\rm Im}\nolimits} \left\{
{\left[ { - \dfrac{{{{(n +
1)}^2}}}{{{R_{20}}}}a_n^{(2)}{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{n + 2}} + n(n +
1)b_n^{(2)}{k_v}h_n^{(1)/}({k_v}{r_2})} \right.} \right.\nonumber\\
& \quad \left. { - \dfrac{{{{( - 1)}^n}{n^2}}}{d}{{\left(
{\dfrac{{{r_2}}}{d}} \right)}^{n - 1}}\sum\limits_{k =
0}^\infty {a_k^{(1)}{C_{kn}}\xi_1^{k + 1}} + n(n +
1){k_v}j_n^/({k_v}{r_2})\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{nk}}} } \right]\nonumber\\
& \quad \times \left[ { -
(m + 3)a_m^{(2) \ast }{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{m + 1}} + b_m^{(2)
\ast }{{[2h_m^{(1)}({k_v}{r_2}) -
{{({k_v}{r_2})}^2}h_m^{(1)/{/}}({k_v}{r_2})]}^ \ast }}
\right.\nonumber\\
& \quad \left. + {{( - 1)}^m}(m - 2){{\left(
{\dfrac{{{r_2}}}{d}} \right)}^m}\sum\limits_{k = 0}^\infty
{a_k^{(1) \ast }{C_{km}}\xi_1^{k + 1}} +
[2{j_m}({k_v}{r_2})\right.\nonumber\\
&\left.\quad -
{{({k_v}{r_2})}^2}j_m^{/{/}}({k_v}{r_2})]^ \ast
\sum\limits_{k = 1}^\infty {{{(b_k^{(1)}{G_{mk}})}^ \ast
}} \right]\nonumber\\
& \quad + \left[ {(n + 1)a_n^{(2)}{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{n + 1}} + n(n +
1)b_n^{(2)}h_n^{(1)}({k_v}{r_2})} \right.\nonumber\\
& \quad \left.
{ - {{( - 1)}^n}n{{\left( {\dfrac{{{r_2}}}{d}}
\right)}^n}\sum\limits_{k = 0}^\infty
{a_k^{(1)}{C_{kn}}\xi_1^{k + 1}} + n(n +
1){j_n}({k_v}{r_2})\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{nk}}} } \right]\nonumber\\
& \quad \times \left[
\dfrac{{(m + 1)(m + 3)}}{{{R_{20}}}}a_m^{(2) \ast
}{{\left( {\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{m + 2}} +
b_m^{(2) \ast }k_v^ \ast [2h_m^{(1)/}({k_v}{r_2})
\right.\nonumber\\
&\left. \quad - 2{k_v}{r_2}h_m^{(1)/{/}}({k_v}{r_2}) - {{({k_v}{r_2})}^2}h_m^{(1)/{/}/}({k_v}{r_2})]^ \ast
\right.\nonumber\\
& \quad + \dfrac{{{{( - 1)}^m}m(m -
2)}}{d}{\left( {\dfrac{{{r_2}}}{d}} \right)^{m -
1}}\sum\limits_{k = 0}^\infty {a_k^{(1) \ast }{C_{km}}\xi
_1^{k + 1}} \nonumber\\ & \quad \left. \left. + k_v^ \ast
[2j_m^/({k_v}{r_2}) - 2{k_v}{r_2}j_m^{/{/}}({k_v}{r_2})\right.\right. \nonumber\\
&\left.\left.\quad -
{{({k_v}{r_2})}^2}j_m^{/{/}/}({k_v}{r_2})]^ \ast
\sum\limits_{k = 1}^\infty {{{(b_k^{(1)}{G_{mk}})}^ \ast
}} \right] \right\},\quad n \ge 0,\;m \ge 1,
\end{align}
\begin{align}
Q_{2nm}^{(2)/}({r_2}) &= {\mathop{\rm Im}\nolimits} \left\{
{\left\{ { - \dfrac{{n + 1}}{{{R_{20}}}}a_n^{(2)}{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{n + 2}} -
b_n^{(2)}{k_v}[2h_n^{(1)/}({k_v}{r_2}) +
{k_v}{r_2}h_n^{(1)/{/}}({k_v}{r_2})]} \right.} \right.\nonumber\\
& \quad \left. { + \dfrac{{{{( - 1)}^n}n}}{d}{{\left(
{\dfrac{{{r_2}}}{d}} \right)}^{n - 1}}\sum\limits_{k =
0}^\infty {a_k^{(1)}{C_{kn}}\xi_1^{k + 1}} -
{k_v}[2j_n^/({k_v}{r_2}) + {k_v}{r_2}j_n^{/{/}}({k_v}{r_2})]\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{nk}}} } \right\}\nonumber\\
& \quad \times \left\{ {a_m^{(2) \ast }{{\left( {\dfrac{{{R_{20}}}}{{{r_2}}}}
\right)}^{m + 1}} - b_m^{(2) \ast }{{[h_m^{(1)}({k_v}{r_2})
+ {k_v}{r_2}h_m^{(1)/}({k_v}{r_2})]}^ \ast }} \right.\nonumber\\
& \quad \left. { + {{( - 1)}^m}{{\left( {\dfrac{{{r_2}}}{d}}
\right)}^m}\sum\limits_{k = 0}^\infty {a_k^{(1) \ast
}{C_{km}}\xi_1^{k + 1}} - {{[{j_m}({k_v}{r_2}) +
{k_v}{r_2}j_m^/({k_v}{r_2})]}^ \ast }\sum\limits_{k =
1}^\infty {{{(b_k^{(1)}{G_{mk}})}^ \ast }} } \right\}\nonumber\\
& \quad + \left\{ {a_n^{(2)}{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{n + 1}} -
b_n^{(2)}[h_n^{(1)}({k_v}{r_2}) +
{k_v}{r_2}h_n^{(1)/}({k_v}{r_2})]} \right.\nonumber\\
& \quad \left. { + {{( - 1)}^n}{{\left( {\dfrac{{{r_2}}}{d}}
\right)}^n}\sum\limits_{k = 0}^\infty
{a_k^{(1)}{C_{kn}}\xi_1^{k + 1}} - [{j_n}({k_v}{r_2}) +
{k_v}{r_2}j_n^/({k_v}{r_2})]\sum\limits_{k = 1}^\infty
{b_k^{(1)}{G_{nk}}} } \right\}\nonumber\\
& \quad \times \left\{ { -
\dfrac{{m + 1}}{{{R_{20}}}}a_m^{(2) \ast }{{\left(
{\dfrac{{{R_{20}}}}{{{r_2}}}} \right)}^{m + 2}} - b_m^{(2)
\ast }k_v^ \ast {{[2h_m^{(1)/}({k_v}{r_2}) +
{k_v}{r_2}h_m^{(1)/{/}}({k_v}{r_2})]}^ \ast }} \right.\nonumber\\
& \quad + \dfrac{{{{( - 1)}^m}m}}{d}{\left(
{\dfrac{{{r_2}}}{d}} \right)^{m - 1}}\sum\limits_{k =
0}^\infty {a_k^{(1) \ast }{C_{km}}\xi _1^{k + 1}} \nonumber\\
& \quad \left. {\left. { - k_v^ \ast {{[{2j_m^/({k_v}{r_2}) +
{k_v}{r_2}j_m^{/{/}}({k_v}{r_2})} ]}^ \ast }\sum\limits_{k
= 1}^\infty {{{(b_k^{(1)}{G_{mk}})}^ \ast }} } \right\}}
\right\},\quad n,m \ge 1.
\end{align}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:
\begin{equation}r_s^{(j)} = {R_{j0}} + s_0^{(j)}\cos (2{\rm \pi} {f_0}t) + s_2^{(j)}{P_2}(\cos \theta )\; \cos (2{\rm \pi} {f_0}t).\end{equation}
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.
With 
${R_{j0}} = 50\;\mathrm{\mu }\textrm{m}$, 
$s_0^{(j)} = 10\;\mathrm{\mu }\textrm{m}$, 
$s_2^{(j)} = 10\;\mathrm{\mu }\textrm{m}$ and 
${f_0} = 30\;\textrm{kHz}$. The initially deformed bubble interface is visible in figure 14(b).
