# Self-Propulsion of Two Contacting Bubbles Due to the Radiation Interaction Force

Alexander A Doinikov, Thomas Micol, Cyril Mauger, Philippe Blanc-Benon,
Claude Inserra

► To cite this version:

Alexander A Doinikov, Thomas Micol, Cyril Mauger, Philippe Blanc-Benon, Claude Inserra. Self-Propulsion of Two Contacting Bubbles Due to the Radiation Interaction Force. Micromachines, 2023, 14 (8), pp.1615. 10.3390/mi14081615. [hal-04878659]

HAL Id: hal-04878659

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

Submitted on 10 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Self-Propulsion of Two Contacting Bubbles Due to the Radiation Interaction Force

Alexander A. Doinikov <fnref n="1" />, Thomas Micol <fnref n="1" />, Cyril Mauger <fnref n="1" />, Philippe Blanc-Benon <fnref n="1" /> and Claude Inserra <fnref n="2" />,*

**Abstract:** In this paper, we consider a new bubble-based microswimmer composed of two contacting bubbles. Under the action of an acoustic field, both bubbles are oscillating, and locomotion of the two-bubble system is observed. A theory is developed that allows one to calculate the acoustic radiation interaction forces between two gas bubbles in an incompressible viscous liquid for any small separation distance between the bubbles. This theory is used to demonstrate that two acoustically excited bubbles can create a self-propelled microswimmer due to a nonzero net force experienced by the bubbles when they come in contact. Experimental evidence of the creation of such a swimmer and of its motion is provided.

**Keywords:** gas bubble; radiation interaction force; viscous liquid; self-propulsion; acoustic manipulation

## 1. Introduction

Interest in artificial self-propelled microswimmers, which could be used in microfluidic and biomedical applications, has existed for decades. The range of possible applications of microswimmers is very wide: cargo transport, micromixing, sensing, targeted drug delivery, indirect manipulation of cells, and other microscopic objects, such as microsurgery, etc.

Different types of microswimmers have been proposed [1], which are actuated by various external energy sources such as light [2], electric [3], magnetic [4] or acoustic [5] fields. One of the promising ways is the development of acoustically controlled bubble-based microswimmers [6–8], which have attracted great attention due to their non-invasiveness and cheap implementation. Acoustically-driven swimmers have been used to capture and move micro-objects [9], to mix fluids [10], and to drive the assembly of microparticles, cells, or microstructures [11]. Bubble-based swimmers are usually composed of a rigid microcapsule that contains an air bubble inside [6,7]. Upon the application of an acoustic field, periodic oscillations of the bubble are induced, resulting in the motion of the capsule. Publications devoted to these investigations commonly attribute the source of locomotion to acoustic microstreaming generated by the bubble. With the aim of accurately manipulating and controlling acoustically-powered swimmers, microrobots composed of two bubbles have been recently designed [8,12]. In addition to the translational motion of the robot, a rotational motion is observed that results from the asymmetry in the bubble oscillations amplitudes and the induced microstreaming. Such asymmetry is generated by switching the frequency of the driving acoustic field between the resonant frequencies of each of the two bubbles. For a two-bubble robot composed of bubbles with equilibrium radii of about 45 µm and driven at a frequency of 40 kHz [8], it has been shown that the translational velocity can reach 8 mm/s, which corresponds to almost 50 body lengths per second. The predominant propulsion mechanism was established to be the bubble-induced microstreaming, while the radiation pressure exerted by the acoustic field on the swimmer was shown to be negligible.

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

### check for updates

**Citation:** Doinikov, A.A.; Micol, T.; Mauger, C.; Blanc-Benon, P.; Inserra, C. Self-Propulsion of Two Contacting Bubbles Due to the Radiation Interaction Force. *Micromachines* **2023**, *14*, 1615. https://doi.org/10.3390/mi14081615

Academic Editor: Hongjie An

Received: 12 July 2023

Revised: 3 August 2023

Accepted: 8 August 2023

Published: 16 August 2023

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

**Copyright:** © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ University Lyon, École Centrale de Lyon, INSA de Lyon, CNRS, LMFA UMR 5509, F-69002 Lyon, France</footnote>

<footnote n="2">² University Lyon, Université Claude Bernard Lyon 1, Centre Léon Bérard, INSERM, UMR 1032, LabTAU, F-69003 Lyon, France</footnote>

<footnote>* Correspondence: claude.inserra@inserm.fr</footnote>
In this paper, we propose a new acoustically powered system composed of two contacting bubbles driven by an acoustic field. We demonstrate that two acoustically excited contacting bubbles can create a self-propelled microswimmer due to the radiation interaction force acting between the bubbles. It is worth noting that a similar geometry has been considered by Pak et al. [13] in relation to two contacting rigid spheres that rotate about their axis of symmetry. This structure was called "snowman" by the authors. In our paper, we will refer to the robot based on two contacting bubbles as a "bubbleman".

When two bubbles are subjected to an acoustic field, the scattered wave generated by one bubble produces a time-averaged force on the other bubble, and vice versa. These forces make the bubbles attract or repel each other. This effect was first reported by C. A. Bjerknes and his son V. F. K. Bjerknes [14], and since then it has been well known in acoustics. In modern literature, this force is referred to by several names: radiation interaction force, secondary radiation force, Bjerknes force, or secondary Bjerknes force [15,16]. C. A. Bjerknes has derived an analytical expression for the time-averaged interaction force between two bubbles, assuming that the bubbles pulsate radially in an incompressible nonviscous liquid and that the separation distance between the bubbles is large compared to the bubble radii. Under such conditions, the system of two interacting bubbles is conservative, which means that the forces on the bubbles are equal and opposite. Doinikov [17,18] has shown that in a viscous liquid, the interaction forces experienced by two bubbles are no longer equal and opposite as the liquid viscosity breaks the conservatism of the system of two interacting bubbles. However, the assumption that the separation distance between the bubbles is large compared to the bubble radii, which was used in his derivation, does not allow one to apply his results to bubbles in contact.

In the present paper, a theory is developed that allows one to calculate the radiation interaction forces between two bubbles in an incompressible viscous liquid for any small separation distance between the bubbles. This theory reveals that if the bubbles come in contact, their agglomerate experiences a nonzero net force, which causes locomotion of the bubbleman. We also provide experimental evidence of this effect.

## 2. Theory

### 2.1. *Interaction Force between Two Bubbles in an Acoustic Field*

We begin with the calculation of the radiation interaction force, also referred to as the secondary radiation force, which acts between two spatially separated gas bubbles in an acoustically excited liquid. This calculation is based on the theory developed by Doinikov et al. [19] for acoustic microstreaming generated by two interacting bubbles. The bubbles are assumed to be immersed in a viscous, incompressible liquid and undergo axisymmetric oscillations. We used two spherical coordinate systems, $(r_1, \theta_1, \epsilon_1)$ and $(r_2, \theta_2, \epsilon_2)$, which originate at the equilibrium centers of bubbles 1 and 2, respectively, and in which the direction $\theta_1 = \theta_2 = 0$ corresponds to the z axis; see Figure 1. The distance between the equilibrium centers of the bubbles is denoted by $d$.

We assume that each bubble can undergo several oscillation modes, which can have different frequencies because some of the modes can be excited parametrically. By definition, the radiation interaction force is a time-averaged quantity. Since time averaging annihilates terms that arise from the interaction of modes with different frequencies, a contribution to the force can only come from the interaction of modes with the same frequency. Therefore, the modes can be divided into groups in accordance with their frequencies and the contributions of such groups can be considered independently.
<figure><img src="image_5.png" /><figcaption>Figure 1. Coordinate systems used in calculations.</figcaption></figure>

Let us assume that the jth bubble (j = 1, 2) undergoes $N_j$ modes with a frequency $\omega$. Let these modes have numbers $M_1^{(j)}$, $M_2^{(j)}$, ..., $M_{N_j}^{(j)}$. Then, the perturbation of the bubble surface produced by these modes can be represented by

$$r_s^{(j)} = R_{j0} + e^{-i\omega t} \sum_{n=M_1^{(j)}}^{M_{N_j}^{(j)}} s_n^{(j)} P_n(\mu_j), \quad (1)$$

where $r_s^{(j)}$ is the radial coordinate of the surface of the jth bubble, $R_{j0}$ is the equilibrium radius of the jth bubble, $P_n$ is the Legendre polynomial of degree $n$, $\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, measured experimentally [19], or evaluated theoretically [20,21].

The time-averaged acoustic radiation force on the jth bubble can be represented as [22]

$$F_j = \int_{S_{j0}} \{-\rho \langle (\mathbf{n}_j \cdot \mathbf{v})\mathbf{v} \rangle + 2\eta (\mathbf{n}_j \cdot \nabla)\mathbf{v}_E + \eta \mathbf{n}_j \times (\nabla \times \mathbf{v}_E) - p_E \mathbf{n}_j\} dS_{j0}, \quad (2)$$

where $S_{j0}$ is the surface of the jth bubble at rest, $\mathbf{n}_j$ is the unit outward normal to $S_{j0}$, $\rho$ is the liquid density, $\mathbf{v}$ is the first-order liquid velocity generated by the bubbles, $\eta$ is the dynamic liquid viscosity, $\mathbf{v}_E$ is the Eulerian velocity of acoustic streaming, $p_E$ is the time-averaged pressure corresponding to $\mathbf{v}_E$, and $\langle \rangle$ means the time average.

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

$$\mathbf{v} = \mathbf{v}^{(1)} + \mathbf{v}^{(2)}, \quad (3)$$

where $\mathbf{v}^{(j)}$ is the first-order liquid velocity generated by the jth bubble. The components of $\mathbf{v}^{(j)}$ are shown in [19] to be calculated by

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

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

where $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, the prime denotes the derivative with respect to the argument in brackets (such as the term $h_n^{(1)}/(k_v r_j)$ in Equation (5)), and $P_n^1$ is the associated Legendre polynomial of the first order and degree $n$. Expressions for the constants $a_n^{(j)}$ and $b_n^{(j)}$, called linear scattering coefficients, are provided in Appendix A of [19].

Equations (4) and (5) give the components of $v^{(j)}$ in the coordinates of the jth bubble. However, to make use of Equation (2), we also need to know how these components are expressed in the coordinates of the other bubble, i.e., how the components of $v^{(1)}$ are expressed in the coordinates ($r_2, \theta_2$) of bubble 2 and vice versa. Necessary equations are provided in Appendix A of [19]; see Equations (A11), (A13), (A15) and (A16) therein. By using these equations, the components of the total first-order liquid velocity are represented by

$$v_r(r_j, \theta_j) = e^{-i\omega t} \sum_{n=0}^{\infty} V_{rn}^{(j)}(r_j) P_n(\mu_j), \quad (6)$$

$$v_{\theta}(r_j, \theta_j) = e^{-i\omega t} \sum_{n=1}^{\infty} V_{\theta n}^{(j)}(r_j) P_n^1(\mu_j), \quad (7)$$

where $V_{rn}^{(j)}(r_j)$ and $V_{\theta n}^{(j)}(r_j)$ are calculated by

$$V_{rn}^{(1)}(r_1) = -\frac{n+1}{r_1} \left[ a_n^{(1)} \left(\frac{R_{10}}{r_1}\right)^{n+1} + b_n^{(1)} n h_n^{(1)}(k_v r_1) \right] + \frac{n}{r_1} \left(\frac{r_1}{d}\right)^n \sum_{m=0}^{\infty} (-1)^m C_{nm} \xi_2^{m+1} a_m^{(2)} \\ -(2n+1) i^n \sqrt{n(n+1)} \frac{j_n(k_v r_1)}{r_1} \sum_{m=1}^{\infty} \frac{\sqrt{m(m+1)}}{(2m+1) i^m} b_m^{(2)} \sum_{l=0}^{\infty} i^{-l} (2l+1) C_{n0l0}^{m0} C_{n1l0}^{m1} h_l^{(1)}(k_v d), \quad (8)$$

$$V_{rn}^{(2)}(r_2) = -\frac{n+1}{r_2} \left[ a_n^{(2)} \left(\frac{R_{20}}{r_2}\right)^{n+1} + b_n^{(2)} n h_n^{(1)}(k_v r_2) \right] + \frac{(-1)^n}{r_2} \left(\frac{r_2}{d}\right)^n \sum_{m=0}^{\infty} C_{nm} \xi_1^{m+1} a_m^{(1)} \\ -(2n+1) i^n \sqrt{n(n+1)} \frac{j_n(k_v r_2)}{r_2} \sum_{m=1}^{\infty} \frac{\sqrt{m(m+1)}}{(2m+1) i^m} b_m^{(1)} \sum_{l=0}^{\infty} (-1)^l i^{-l} (2l+1) C_{n0l0}^{m0} C_{n1l0}^{m1} h_l^{(1)}(k_v d), \quad (9)$$

$$V_{\theta n}^{(1)}(r_1) = \frac{1}{r_1} \left\{ a_n^{(1)} \left(\frac{R_{10}}{r_1}\right)^{n+1} - b_n^{(1)} \left[ h_n^{(1)}(k_v r_1) + k_v r_1 h_n^{(1)}/(k_v r_1) \right] \right\} + \frac{1}{r_1} \left(\frac{r_1}{d}\right)^n \sum_{m=0}^{\infty} (-1)^m C_{nm} \xi_2^{m+1} a_m^{(2)} \\ - \frac{(2n+1) i^n}{\sqrt{n(n+1)}} \frac{j_n(k_v r_1) + k_v r_1 j_n'(k_v r_1)}{r_1} \sum_{m=1}^{\infty} \frac{\sqrt{m(m+1)}}{(2m+1) i^m} b_m^{(2)} \sum_{l=0}^{\infty} i^{-l} (2l+1) C_{n0l0}^{m0} C_{n1l0}^{m1} h_l^{(1)}(k_v d), \quad (10)$$

$$V_{\theta n}^{(2)}(r_2) = \frac{1}{r_2} \left\{ a_n^{(2)} \left(\frac{R_{20}}{r_2}\right)^{n+1} - b_n^{(2)} \left[ h_n^{(1)}(k_v r_2) + k_v r_2 h_n^{(1)}/(k_v r_2) \right] \right\} + \frac{(-1)^n}{r_2} \left(\frac{r_2}{d}\right)^n \sum_{m=0}^{\infty} C_{nm} \xi_1^{m+1} a_m^{(1)} \\ - \frac{(2n+1) i^n}{\sqrt{n(n+1)}} \frac{j_n(k_v r_2) + k_v r_2 j_n'(k_v r_2)}{r_2} \sum_{m=1}^{\infty} \frac{\sqrt{m(m+1)}}{(2m+1) i^m} b_m^{(1)} \sum_{l=0}^{\infty} (-1)^l i^{-l} (2l+1) C_{n0l0}^{m0} C_{n1l0}^{m1} h_l^{(1)}(k_v d). \quad (11)$$

Here, $C_{nm} = (n+m)!/(n!m!)$, $\xi_j = R_{j0}/d$, $j_n$ is the spherical Bessel function, and $C_{l_1 m_1 l_2 m_2}^{LM}$ are the Clebsch-Gordan coefficients [23–25]. The superscript (j) in $V_{rn}^{(j)}$ and $V_{\theta n}^{(j)}$ emphasizes that these quantities are taken in the coordinates of the jth bubble. According to the theory developed in [19], Equations (8)–(11) are valid for $r_j < d$. This is sufficient for our purpose because we will see below that the expression for the interaction force only contains the values of $V_{rn}^{(j)}$ and $V_{\theta n}^{(j)}$ at $r_j = R_{j0}$.

Equations (6) and (7) allow us to calculate the first term in the integrand of Equation (2). To calculate the other terms, we again use the results of [19].

According to Equations (2.17) and (2.18) of [19], the components of the Eulerian streaming velocity $v_E$ are given by

$$v_{Er}(r_j, \theta_j) = -\frac{1}{r_j} \sum_{n=1}^{\infty} n(n+1) \Psi_n^{(j)}(r_j) P_n(\mu_j), \quad (12)$$

$$v_{E\theta}(r_j, \theta_j) = -\frac{1}{r_j} \sum_{n=1}^{\infty} \left[ \Psi_n^{(j)}(r_j) + r_j \Psi_n^{(j)}/(r_j) \right] P_n^1(\mu_j), \quad (13)$$
where $\Psi_n^{(j)}(r_j)$ and $\Psi_n^{(j)\prime}(r_j)$ are calculated by Equations (C8) and (C17) of [19]. Equations (12) and (13) allow us to calculate the terms of Equation (2) that depend on $v_E$ and $p_E$.

With the help of Equations (6), (7) and (A1)–(A4) from Appendix A of the present paper, the first term in the integrand of Equation (2) brings

$$ 
\begin{aligned}
T_1^{(j)} &= -\int_{S_{j0}} \rho \langle (\mathbf{n}_j \cdot \mathbf{v}) \mathbf{v} \rangle dS_{j0} = -\pi \rho R_{j0}^2 \mathbf{e}_z \mathrm{Re} \int_0^\pi v_r^*(v_r \cos \theta_j - v_\theta \sin \theta_j) \sin \theta_j d\theta_j \\
&= -\pi \rho R_{j0}^2 \mathbf{e}_z \mathrm{Re} \left\{ \sum_{n,m=0}^\infty V_{rn}^{(j)}(R_{j0}) V_{rm}^{(j)*}(R_{j0}) \int_{-1}^1 P_n(\mu_j) P_m(\mu_j) \mu_j d\mu_j \right. \\
&\qquad \left. - \sum_{\substack{n=1 \\ m=0}}^\infty V_{\theta n}^{(j)}(R_{j0}) V_{r m}^{(j)*}(R_{j0}) \int_{-1}^1 P_n^1(\mu_j) P_m(\mu_j) \sqrt{1 - \mu_j^2} d\mu_j \right\} \\
&= -2\pi \rho R_{j0}^2 \mathbf{e}_z \mathrm{Re} \left\{ V_{r0}^{(j)}(R_{j0}) V_{r1}^{(j)*}(R_{j0}) + \sum_{n=1}^\infty \frac{1}{(2n+1)^2} \left\{ V_{rn}^{(j)}(R_{j0}) [nV_{r(n-1)}^{(j)*}(R_{j0}) + (n+1)V_{r(n+1)}^{(j)*}(R_{j0})] \right. \right. \\
&\qquad \left. - n(n+1)V_{\theta n}^{(j)}(R_{j0}) [V_{r(n+1)}^{(j)*}(R_{j0}) - V_{r(n-1)}^{(j)*}(R_{j0})] \right\}.
\end{aligned} 
\quad (14) $$

The contribution of the second term is calculated as follows:

$$ 
\begin{aligned}
T_2^{(j)} &= 2\eta \int_{S_{j0}} (\mathbf{n}_j \cdot \nabla) \mathbf{v}_E dS_{j0} = 4\pi\eta R_{j0}^2 \mathbf{e}_z \int_{-1}^1 \left( \frac{\partial v_{Er}}{\partial r_j} \mu_j - \frac{\partial v_{E\theta}}{\partial r_j} \sqrt{1 - \mu_j^2} \right)_{r_j=R_{j0}} d\mu_j \\
&= 4\pi\eta R_{j0}^2 \mathbf{e}_z \left\{ -\sum_{n=1}^\infty n(n+1) \left[ \frac{\Psi_n^{(j)\prime}(R_{j0})}{R_{j0}} - \frac{\Psi_n^{(j)}(R_{j0})}{R_{j0}^2} \right] \int_{-1}^1 P_n(\mu_j) \mu_j d\mu_j \right. \\
&\qquad \left. - \sum_{n=1}^\infty \left[ \frac{\Psi_n^{(j)}(R_{j0})}{R_{j0}^2} - \frac{\Psi_n^{(j)\prime}(R_{j0})}{R_{j0}} - \Psi_n^{(j)\prime\prime}(R_{j0}) \right] \int_{-1}^1 P_n^1(\mu_j) \sqrt{1 - \mu_j^2} d\mu_j \right\} \\
&= \frac{16}{3} \pi\eta R_{j0}^2 \mathbf{e}_z \left[ \frac{2\Psi_1^{(j)}(R_{j0})}{R_{j0}^2} - \frac{2\Psi_1^{(j)\prime}(R_{j0})}{R_{j0}} - \Psi_1^{(j)\prime\prime}(R_{j0}) \right].
\end{aligned} 
\quad (15) $$

Making use of the equation

$$ \mathbf{n}_j \times (\nabla \times \mathbf{v}_E) = \frac{e_{\theta j}}{r_j} \left[ \frac{\partial v_{Er}}{\partial \theta_j} - \frac{\partial (r_j v_{E\theta})}{\partial r_j} \right], \quad (16) $$

the calculation of the third term results in

$$ 
\begin{aligned}
T_3^{(j)} &= \eta \int_{S_{j0}} \mathbf{n}_j \times (\nabla \times \mathbf{v}_E) dS_{j0} \\
&= 2\pi\eta R_{j0} \mathbf{e}_z \sum_{n=1}^{\infty} \left[ \frac{n(n+1)}{R_{j0}} \Psi_n^{(j)}(R_{j0}) - 2\Psi_n^{(j)\prime}(R_{j0}) - R_{j0} \Psi_n^{(j)\prime\prime}(R_{j0}) \right] \int_{-1}^{1} P_n^1(\mu_j) \sqrt{1 - \mu_j^2} d\mu_j \\
&= -\frac{8}{3}\pi\eta R_{j0}^2 \mathbf{e}_z \left[ \frac{2\Psi_1^{(j)}(R_{j0})}{R_{j0}^2} - \frac{2\Psi_1^{(j)\prime}(R_{j0})}{R_{j0}} - \Psi_1^{(j)\prime\prime}(R_{j0}) \right].
\end{aligned} 
\quad (17) $$
To calculate the fourth term of Equation (2), we first need to calculate $p_E$. In order to perform this calculation, we only need the expression of $p_E$ at $r_j = R_{j0}$. In addition, it is reasonable to assume that the time-averaged gas pressure inside the jth bubble has the same value at all points of the gas medium, which means that $p_E$ should be constant, independent of $\theta_j$, at $r_j = R_{j0}$. It is easy to check that, in this case, the fourth term of Equation (2) does not contribute to the force.

To calculate $\Psi_1^{(j)}(R_{j0})$, $\Psi_1^{(j)\prime}(R_{j0})$ and $\Psi_1^{(j)\prime\prime}(R_{j0})$, we use Equations (C8), (C17) and (C18) of [19],

$$ \Psi_1^{(j)}(R_{j0}) = C_{110}^{(j)} + \frac{C_{210}^{(j)}}{R_{j0}^2} + R_{j0}C_{310}^{(j)} + R_{j0}^3 C_{410}^{(j)}, \quad (18) $$

$$ \Psi_1^{(j)\prime}(R_{j0}) = -\frac{2C_{210}^{(j)}}{R_{j0}^3} + C_{310}^{(j)} + 3R_{j0}^2 C_{410}^{(j)}, \quad (19) $$

$$ \Psi_1(j) \parallel (Rj0) = 6C210(j)Rj04 + 6Rj0C410(j), \quad (20) $$

where the constants $C_{110}^{(j)}-C_{410}^{(j)}$ are calculated by Equations (C30)–(C33) of [19].

With the help of Equations (18)–(20), one obtains

$$ \frac{\Psi_1^{(j)}(R_{j0})}{R_{j0}^2} - \frac{\Psi_1^{(j)\prime}(R_{j0})}{R_{j0}} - \frac{\Psi_1^{(j)\prime\prime}(R_{j0})}{2} = \frac{C_{110}^{(j)}}{R_{j0}^2} - 5R_{j0}C_{410}^{(j)}. \quad (21) $$

Combining Equations (14), (15), (17) and (21) brings the following expression for the radiation interaction force on the jth bubble:

$$ F_j = \frac{16\pi}{3}\eta(C_{110}^{(j)} - 5R_{j0}^3 C_{410}^{(j)}) - 2\pi\rho R_{j0}^2 \text{Re} \left\{ V_{r0}^{(j)}(R_{j0}) V_{r1}^{(j)*}(R_{j0}) + \sum_{n=1}^{\infty} \frac{1}{(2n+1)^2} \left\{ V_{rn}^{(j)}(R_{j0}) \left[ nV_{r(n-1)}^{(j)*}(R_{j0}) + (n+1)V_{r(n+1)}^{(j)*}(R_{j0}) \right] - n(n+1)V_{\theta n}^{(j)}(R_{j0}) \left[ V_{r(n+1)}^{(j)*}(R_{j0}) - V_{r(n-1)}^{(j)*}(R_{j0}) \right] \right\} \right\}, \quad (22) $$

where $C_{110}^{(j)}$ and $C_{410}^{(j)}$ are constants that are calculated by Equations (C30) and (C33) of [19], Re means “the real part of”, and the asterisk denotes complex conjugate. Note that the force is directed along the line joining the equilibrium centers of the bubbles, i.e., along the z axis in Figure 1. Therefore, $F_j > 0$ means that the force acts in the positive direction of the z axis, while $F_j < 0$ means that the force acts in the negative direction of the z axis.

In this subsection we have derived an analytical formula that allows one to calculate the interaction forces on two bubbles in the case that the bubbles undergo strong shape modes, which are excited parametrically, and for which we know the magnitudes and the phases of these modes.

## 2.2. Linear Scattering Coefficients When Parametric Excitation Is Absent

If the parametric excitation of shape modes is absent, $a_n^{(j)}$ and $b_n^{(j)}$ can be expressed in terms of the amplitude of the imposed acoustic pressure. In this case, Equation (1) is transformed to

$$ r_s^{(j)} = R_{j0} + e^{-i\omega t} \sum_{n=0}^{\infty} s_n^{(j)} P_n(\mu_j), \quad (23) $$

The linearized equations of an incompressible viscous liquid are given by [26]

$$ \nabla \cdot \mathbf{v} = 0, \quad (24) $$

$$ \frac{\partial \mathbf{v}}{\partial t} = -\frac{1}{\rho} \nabla p + \nu \Delta \mathbf{v}, \quad (25) $$
where *p* is the first-order liquid pressure.

In the case of two bubbles, a solution for *v* is sought as

$$v = v^{(1)} + v^{(2)}, \qquad (26)$$

where $v^{(j)}$ is represented by

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

with the scalar, $\varphi^{(j)}$, and the vector, $\varphi^{(j)}$, velocity potentials defined by [27]

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

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

where $e_{\varepsilon_j}$ is the unit azimuth vector of the jth bubble. Note that axial symmetry allows us to set $\varepsilon_1 = \varepsilon_2$ and $e_{\varepsilon_1} = e_{\varepsilon_2}$.

Substituting Equation (27) into Equation (25) and taking into account that $\varphi^{(j)}$ obeys the equation $(\Delta + k_v^2)\varphi^{(j)} = 0$ [27], one finds that the first-order scattered pressure generated by the jth bubble is given by

$$p^{(j)} = i\omega\rho\varphi^{(j)}. \qquad (30)$$

Note that the total first-order scattered pressure is equal to $p = p^{(1)} + p^{(2)}$.

We use the boundary condition for normal stress at the surface of the jth bubble,

$$P_{gj} \left( \frac{V_{j0}}{V_j(t)} \right)^\gamma = p|_{r_j=R_{j0}} - 2\eta \left. \frac{\partial v_r(r_j, \theta_j)}{\partial r_j} \right|_{r_j=R_{j0}} + p_{st}^{(j)} + P_{ac}(t) + P_0, \qquad (31)$$

where the gas pressure within the bubbles is assumed to be spatially homogeneous and to obey the adiabatic law, $P_{gj}$ is the equilibrium gas pressure inside the jth bubble, $V_{j0} = (4/3)\pi R_{j0}^3$ is the equilibrium volume of the jth bubble, $V_j(t)$ is the instantaneous volume of the jth bubble, $\gamma$ is the specific heat ratio of the gas, $P_0$ is the hydrostatic pressure in the liquid, $P_{ac}(t) = P_a e^{-i\omega t}$ is the imposed acoustic pressure, and $p_{st}^{(j)}$ is the pressure of surface tension on the surface of the jth bubble, given by [28]

$$p_{st}^{(j)} = \frac{2\sigma}{R_{j0}} + e^{-i\omega t} \frac{\sigma}{R_{j0}^2} \sum_{n=0}^{\infty} (n-1)(n+2)s_n^{(j)} P_n(\mu_j), \qquad (32)$$

where $\sigma$ is the surface tension coefficient.

Accurate to first-order terms, $V_j(t)$ is calculated by

$$V_j(t) = \int_0^{r_s^{(j)}} r^2 dr \int_0^{\pi} \sin \theta_j d\theta_j \int_0^{2\pi} d\varepsilon_j = \frac{2\pi}{3} \int_{-1}^{1} (r_s^{(j)})^3 d\mu_j \approx V_{j0} \left( 1 + e^{-i\omega t} \frac{3s_0^{(j)}}{R_{j0}} \right). \qquad (33)$$

To apply Equation (31) and thus calculate $a_n^{(j)}$ and $b_n^{(j)}$, we need to know $\varphi^{(1)}$ and $v_r^{(1)}$ in the coordinates $(r_2, \theta_2)$ and $\varphi^{(2)}$ and $v_r^{(2)}$ in the coordinates $(r_1, \theta_1)$. To this end, we use Equations (A5), (A6), (A11) and (A13) of [19], which give

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

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

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

With the help of Equations (30) and (32)–(37), applying Equation (31) at $j=1$ and $j=2$, one obtains for $n=0$,

$$s_0^{(1)} = \frac{\omega}{\omega_{10}^2 R_{10}} \left[ (2\tau_1 - i)a_0^{(1)} - i \sum_{m=0}^{\infty} (-1)^m \xi_2^{m+1} a_m^{(2)} - \frac{P_a}{\rho\omega} \right], \quad (38)$$

$$s_0^{(2)} = \frac{\omega}{\omega_{20}^2 R_{20}} \left[ (2\tau_2 - i)a_0^{(2)} - i \sum_{m=0}^{\infty} \xi_1^{m+1} a_m^{(1)} - \frac{P_a}{\rho\omega} \right], \quad (39)$$

for $n=1$,

$$(6\tau_1 - i)a_1^{(1)} + 2\tau_1[h_1^{(1)}(x_1) - x_1h_1^{(1)\prime}(x_1)]b_1^{(1)}$$

$$-i\xi_1 \sum_{m=0}^{\infty} (-1)^m (m+1) \xi_2^{m+1} a_m^{(2)} + 3\sqrt{2i\tau_1} [j_1(x_1) - x_1j_1'(x_1)] \sum_{m=1}^{\infty} \kappa_{1m}^{(2)} b_m^{(2)} = 0, \quad (40)$$

$$(6\tau_2 - i)a_1^{(2)} + 2\tau_2[h_1^{(1)}(x_2) - x_2h_1^{(1)\prime}(x_2)]b_1^{(2)}$$

$$+i\xi_2 \sum_{m=0}^{\infty} (m+1)\xi_1^{m+1} a_m^{(1)} + 3\sqrt{2i\tau_2} [j_1(x_2) - x_2j_1'(x_2)] \sum_{m=1}^{\infty} \kappa_{1m}^{(1)} b_m^{(1)} = 0, \quad (41)$$

and for $n \ge 2$,

$$\frac{\omega_{1n}^2 R_{10}}{(n+1)\omega} s_n^{(1)} = [(n+1)(n+2)\tau_1 - i]a_n^{(1)} + n(n+1)\tau_1[h_n^{(1)}(x_1) - x_1h_n^{(1)\prime}(x_1)]b_n^{(1)} \\ + [n(n-1)\tau_1 - i]\xi_1^n \sum_{m=0}^{\infty} (-1)^m C_{nm}\xi_2^{m+1}a_m^{(2)} + \tau_1 i^n (2n+1)\sqrt{n(n+1)}[j_n(x_1) - x_1j_n'(x_1)]\sum_{m=1}^{\infty} \kappa_{nm}^{(2)}b_m^{(2)}, \quad (42)$$

$$\frac{\omega_{2n}^2 R_{20}}{(n+1)\omega} s_n^{(2)} = [(n+1)(n+2)\tau_2 - i]a_n^{(2)} + n(n+1)\tau_2[h_n^{(1)}(x_2) - x_2h_n^{(1)\prime}(x_2)]b_n^{(2)} \\ + (-1)^n [n(n-1)\tau_2 - i]\xi_2^n \sum_{m=0}^{\infty} C_{nm}\xi_1^{m+1}a_m^{(1)} + \tau_2 i^n (2n+1)\sqrt{n(n+1)}[j_n(x_2) - x_2j_n'(x_2)]\sum_{m=1}^{\infty} \kappa_{nm}^{(1)}b_m^{(1)}, \quad (43)$$

where $x_j = k_v R_{j0}$,

$$\omega_{j0} = \frac{1}{R_{j0}} \sqrt{\frac{1}{\rho} \left( 3\gamma P_{gj} - \frac{2\sigma}{R_{j0}} \right)} \quad (44)$$
is the natural frequency of the radial mode of the jth bubble,

$$ \omega_{jn} = \sqrt{(n^2 - 1)(n + 2)\frac{\sigma}{\rho R_{j0}^3}} \qquad (45) $$

is the natural frequency of the nth shape mode of the jth bubble (n ≥ 2),

$$ \tau_j = \frac{2\nu}{\omega R_{j0}^2}, \qquad (46) $$

$$ \kappa_{nm}^{(j)} = \frac{\sqrt{m(m+1)}}{(2m+1)i^m} \sum_{l=0}^{\infty} (-1)^{jl} i^{-l} (2l+1) C_{n0l0}^{m0} C_{n1l0}^{m1} h_l^{(1)}(k_v d). \qquad (47) $$

It also follows from Equation (31) that

$$ P_{gj} = P_0 + \frac{2\sigma}{R_{j0}}. \qquad (48) $$

To make the system of Equations (38)–(43) closed, we use Equations (A19)–(A38) of [19], which were obtained there from the boundary conditions for the normal component of the liquid velocity and the liquid tangential stress. The above equations give for n = 0,

$$ a_0^{(1)} = i\omega R_{10} s_0^{(1)}, \qquad (49) $$

$$ a_0^{(2)} = i\omega R_{20} s_0^{(2)}, \qquad (50) $$

and for n ≥ 1,

$$ a_n^{(1)} + f_n^{(1)} b_n^{(1)} - n \sum_{m=0}^{\infty} \alpha_{2nm} a_m^{(2)} + \sum_{m=1}^{\infty} \beta_{2nm} b_m^{(2)} = \frac{i\omega R_{10}}{n+1} s_n^{(1)}, \qquad (51) $$

$$ a_n^{(2)} + f_n^{(2)} b_n^{(2)} - n \sum_{m=0}^{\infty} \alpha_{1nm} a_m^{(1)} + \sum_{m=1}^{\infty} \beta_{1nm} b_m^{(1)} = \frac{i\omega R_{20}}{n+1} s_n^{(2)}, \qquad (52) $$

$$ a_n^{(1)} + g_n^{(1)} b_n^{(1)} - \frac{n^2 - 1}{n + 2} \sum_{m=0}^{\infty} \alpha_{2nm} a_m^{(2)} + \sum_{m=1}^{\infty} \gamma_{2nm} b_m^{(2)} = 0, \qquad (53) $$

$$ a_n^{(2)} + g_n^{(2)} b_n^{(2)} - \frac{n^2 - 1}{n + 2} \sum_{m=0}^{\infty} \alpha_{1nm} a_m^{(1)} + \sum_{m=1}^{\infty} \gamma_{1nm} b_m^{(1)} = 0, \qquad (54) $$

where the coefficients $f_n^{(j)}$, $g_n^{(j)}$, $\alpha_{jnm}$, $\beta_{jnm}$ and $\gamma_{jnm}$ are calculated by Equations (A25)–(A34) of [19].

Combining Equations (38), (39), (49) and (50), one obtains

$$ \left(1 - \frac{\omega_{10}^2}{\omega^2} + 2i\tau_1\right) a_0^{(1)} + \sum_{m=0}^{\infty} (-1)^m \xi_2^{m+1} a_m^{(2)} = \frac{iP_a}{\rho\omega}, \qquad (55) $$

$$ \left(1 - \frac{\omega_{20}^2}{\omega^2} + 2i\tau_2\right) a_0^{(2)} + \sum_{m=0}^{\infty} \xi_1^{m+1} a_m^{(1)} = \frac{iP_a}{\rho\omega}. \qquad (56) $$

Combining Equations (42), (43), (51) and (52), one obtains for n ≥ 2,

$$ \left[1 - \frac{\omega_{1l}^2}{\omega^2} + i(n+1)(n+2)\tau_1\right] a_n^{(1)} + \left\{in(n+1)\tau_1\left[h_n^{(1)}(x_1) - x_1h_n^{(1)}/(x_1)\right] - \frac{\omega_{1l}^2}{\omega^2}f_n^{(1)}\right\}b_n^{(1)} \\ + \sum_{m=0}^{\infty} \left\{(-1)^m C_{nm}[n(n-1)i\tau_1 + 1]\xi_1^n \xi_2^{m+1} + n\frac{\omega_{1l}^2}{\omega^2}\alpha_{2nm}\right\}a_m^{(2)} $$
$$ + \sum_{m=1}^{\infty} \left\{ \tau_1 i^{n+1} (2n+1) \sqrt{n(n+1)} \left[ j_n(x_1) - x_1 j'_n(x_1) \right] \kappa_{nm}^{(2)} - \frac{\omega_{1n}^2}{\omega^2} \beta_{2nm} \right\} b_m^{(2)} = 0, \quad (57) $$

$$ \begin{aligned} & \left[ 1 - \frac{\omega_{2n}^2}{\omega^2} + i(n+1)(n+2)\tau_2 \right] a_n^{(2)} + \left\{ i n(n+1)\tau_2 \left[ h_n^{(1)}(x_2) - x_2 h_n^{(1)\prime}(x_2) \right] - \frac{\omega_{2n}^2}{\omega^2} f_n^{(2)} \right\} b_n^{(2)} \\ & + \sum_{m=0}^{\infty} \left\{ (-1)^n C_{nm} [n(n-1)i\tau_2 + 1] \xi_2^n \xi_1^{m+1} + n \frac{\omega_{2n}^2}{\omega^2} \alpha_{1nm} \right\} a_m^{(1)} \\ & + \sum_{m=1}^{\infty} \left\{ \tau_2 i^{n+1} (2n+1) \sqrt{n(n+1)} \left[ j_n(x_2) - x_2 j'_n(x_2) \right] \kappa_{nm}^{(1)} - \frac{\omega_{2n}^2}{\omega^2} \beta_{1nm} \right\} b_m^{(1)} = 0. \end{aligned} \quad (58) $$

Equations (40), (41), (53), and (54) do not contain the mode amplitudes $s_n^{(j)}$ and therefore remain unchanged.

Equations (40), (41), (53), (54) and (55)–(58) form a system of equations in the unknowns $a_n^{(j)}$ and $b_n^{(j)}$. The number of the equations, just as the number of the coefficients $a_n^{(j)}$ and $b_n^{(j)}$, is infinite. However, since $a_n^{(j)}$ and $b_n^{(j)}$ decrease for $n \to \infty$, their number can be truncated at some value of $n = N$. Doing so, we obtain a finite system of $4N + 2$ equations, which can be solved numerically,

$$ \left( 1 - \frac{\omega_{10}^2}{\omega^2} + 2i\tau_1 \right) a_0^{(1)} + \sum_{n=0}^{N} (-1)^n \xi_2^{n+1} a_n^{(2)} = \frac{iP_a}{\rho\omega}, \quad (59) $$

$$ \sum_{n=0}^{N} \xi_1^{n+1} a_n^{(1)} + \left( 1 - \frac{\omega_{20}^2}{\omega^2} + 2i\tau_2 \right) a_0^{(2)} = \frac{iP_a}{\rho\omega}, \quad (60) $$

$$ \begin{aligned} & (6\tau_1 - i)a_1^{(1)} - i\xi_1 \sum_{m=0}^{N} (-1)^m(m+1)\xi_2^{m+1}a_m^{(2)} \\ & + 2\tau_1[h_1^{(1)}(x_1) - x_1h_1^{(1)\prime}(x_1)]b_1^{(1)} + 3\sqrt{2}i\tau_1[j_1(x_1) - x_1j_1^{\prime\prime}(x_1)]\sum_{m=1}^{N} \kappa_{1m}^{(2)}b_m^{(2)} = 0, \end{aligned} \quad (61) $$

$$ \begin{aligned} & i\xi_2 \sum_{m=0}^{N} (m+1)\xi_1^{m+1}a_m^{(1)} + (6\tau_2 - i)a_1^{(2)} \\ & + 3\sqrt{2}i\tau_2[j_1(x_2) - x_2j_1^{\prime\prime}(x_2)]\sum_{m=1}^{N} \kappa_{1m}^{(1)}b_m^{(1)} + 2\tau_2[h_1^{(1)}(x_2) - x_2h_1^{(1)\prime}(x_2)]b_1^{(2)} = 0, \end{aligned} \quad (62) $$

$$ \begin{aligned} & \left[ 1 - \frac{\omega_{1n}^2}{\omega^2} + i(n+1)(n+2)\tau_1 \right] a_n^{(1)} + \sum_{m=0}^{N} \left\{ (-1)^m C_{nm} [n(n-1)i\tau_1 + 1] \xi_1^n \xi_2^{m+1} + n \frac{\omega_{1n}^2}{\omega^2} \alpha_{2nm} \right\} a_m^{(2)} \\ & + \left\{ i n(n+1)\tau_1 \left[ h_n^{(1)}(x_1) - x_1 h_n^{(1)\prime}(x_1) \right] - \frac{\omega_{1n}^2}{\omega^2} f_n^{(1)} \right\} b_n^{(1)} \\ & + \sum_{m=1}^{N} \left\{ \tau_1 i^{n+1} (2n+1) \sqrt{n(n+1)} \left[ j_n(x_1) - x_1 j'_n(x_1) \right] \kappa_{nm}^{(2)} - \frac{\omega_{1n}^2}{\omega^2} \beta_{2nm} \right\} b_m^{(2)} = 0, \quad 2 \le n \le N, \end{aligned} \quad (63) $$

$$ \begin{aligned} & \sum_{m=0}^{N} \left\{ (-1)^n C_{nm} [n(n-1)i\tau_2 + 1] \xi_2^n \xi_1^{m+1} + n \frac{\omega_{2n}^2}{\omega^2} \alpha_{1nm} \right\} a_m^{(1)} + \left[ 1 - \frac{\omega_{2n}^2}{\omega^2} + i(n+1)(n+2)\tau_2 \right] a_n^{(2)} \\ & + \sum_{m=1}^{N} \left\{ \tau_2 i^{n+1} (2n+1) \sqrt{n(n+1)} \left[ j_n(x_2) - x_2 j'_n(x_2) \right] \kappa_{nm}^{(1)} - \frac{\omega_{2n}^2}{\omega^2} \beta_{1nm} \right\} b_m^{(1)} \\ & + \left\{ i n(n+1)\tau_2 \left[ h_n^{(1)}(x_2) - x_2 h_n^{(1)\prime}(x_2) \right] - \frac{\omega_{2n}^2}{\omega^2} f_n^{(2)} \right\} b_n^{(2)} = 0, \quad 2 \le n \le N, \end{aligned} \quad (64) $$
$$a_n^{(1)} + \frac{1-n^2}{n+2} \sum_{m=0}^{N} \alpha_{2nm} a_m^{(2)} + g_n^{(1)} b_n^{(1)} + \sum_{m=1}^{N} \gamma_{2nm} b_m^{(2)} = 0, 1 \le n \le N, \qquad (65)$$

$$\frac{1-n^2}{n+2} \sum_{m=0}^{N} \alpha_{1nm} a_m^{(1)} + a_n^{(2)} + \sum_{m=1}^{N} \gamma_{1nm} b_m^{(1)} + g_n^{(2)} b_n^{(2)} = 0, 1 \le n \le N. \qquad (66)$$

Equations (51) and (52) at $n = 1$ are not used when $a_n^{(j)}$ and $b_n^{(j)}$ are calculated. They define the amplitudes of the translational oscillations of the bubbles, $s_1^{(1)}$ and $s_1^{(2)}$.

Changing $N$ in Equations (59)–(66), one can calculate the linear scattering coefficients $a_n^{(j)}$ and $b_n^{(j)}$ with any desired accuracy.

## 2.3. Net Force Experienced by Two Contacting Bubbles

The theory developed above allows one to calculate the forces on the bubbles, $F_1$ and $F_2$, at any distance $d$ between the bubbles. For the bubbles in contact, we should set $d = R_{10} + R_{20}$. The problem is that at this value of $d$ division by zero appears in equations obtained in [19] for acoustic microstreaming, which are used in the present paper. To overcome this problem, the forces on the bubbles can be calculated at a very small but nonzero distance between the bubbles' surfaces and then the sum $F_1 + F_2$ can be used as approximation of the net force acting on the two bubbles in contact. It is reasonable to assume that this approach should provide the net force on two contacting bubbles with acceptable accuracy.

## 3. Numerical Simulations

We consider the case without parametric excitation and calculate the linear scattering coefficients by means of the equations derived in Section 2.2. The following material parameters are used: $\rho = 1000 \text{ kg/m}^3$, $\eta = 0.001 \text{ Pa s}$, $\sigma = 0.0727 \text{ N/m}$, $\gamma = 1.4$, $P_0 = 101.3 \text{ kPa}$, $P_a = 10 \text{ kPa}$, $f = \omega/2\pi = 30 \text{ kHz}$. These parameters correspond to air bubbles in water.

When parametric excitation is absent, dominant oscillation modes are the radial mode (mode 0) and the translational mode (mode 1). As an example, Figure 2 shows the magnitudes of the mode amplitudes $s_n^{(j)}$ for two contacting bubbles (bubbleman) with $R_{10} = 20 \text{ µm}$ and $R_{20} = 50 \text{ µm}$. As one can see, for bubble 1 (smaller bubble), the translational mode is dominant, while for bubble 2 (bigger bubble), the radial mode is dominant.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** The magnitudes of the mode amplitudes $s_n^{(j)}$ for a bubbleman with $R_{10} = 20 \text{ µm}$ and $R_{20} = 50 \text{ µm}$ in the case that parametric excitation is absent. The calculation was carried out by the equations derived in Section 2.2.</figcaption></figure>
Figure 3 shows the net force experienced by the bubbleman. The equilibrium radius of bubble 2 is kept fixed, $R_{20} = 50$ µm, while the equilibrium radius of bubble 1 is varied. The calculation was made for three values of the liquid viscosity $\eta$. As one can see, the force magnitude decreases with increasing viscosity. For lower viscosity, the net force is always directed from the smaller bubble to the bigger bubble, while for higher viscosity, the force can act in the opposite direction. As expected, the net force vanishes for $R_{10} = R_{20}$ and for $R_{10} \to 0$.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Net force experienced by a bubbleman in the case that parametric excitation is absent. The equilibrium radius of bubble 2 is kept fixed, $R_{20} = 50$ µm, while the equilibrium radius of bubble 1 is varied.</figcaption></figure>

Figure 4 shows the magnitude of the shape modes developing on bubble 1 as a function of $R_{10}$. The equilibrium radius of bubble 2 is kept fixed, $R_{20} = 50$ µm, $\eta = 0.001$ Pa s, and the other parameters are as in Figure 3. The comparison of Figures 3 and 4 suggests that the trough in the upper curve of Figure 3 is caused by the resonance of mode 2 of bubble 1, as at $R_{10} = 29.065$ µm, the driving frequency of 30 kHz is equal to the natural frequency of the second shape mode; see Equation (45).

Figure 5 shows the contributions of the force components, given by Equations (14), (15) and (17), to the net force on the bubbleman as a function of the acoustic pressure amplitude $P_a$ for two values of the liquid viscosity, $\eta = 0.001$ Pa s and $\eta = 0.004$ Pa s. It is assumed that $R_{10} = 20$ µm, $R_{20} = 50$ µm and $f = 30$ kHz. The solid curve shows the net force $F_1 + F_2$, the short-dash curve shows the contribution of the terms $T_1^{(j)}$ given by Equation (14), and the long-dash curve shows the contribution of the terms $T_2^{(j)}$ and $T_3^{(j)}$ given by Equations (15) and (17). As one can see, at $\eta = 0.001$ Pa s, the dominant contribution comes from the terms $T_1^{(j)}$, while the terms $T_2^{(j)}$ and $T_3^{(j)}$, which arise due to acoustic microstreaming, reduce the net force slightly. At $\eta = 0.004$ Pa s, the contribution of the terms $T_2^{(j)}$ and $T_3^{(j)}$ becomes comparable to the contribution of the terms $T_1^{(j)}$, which results in a considerable reduction in the net force.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** The magnitude of the shape modes developing on bubble 1 as a function of $R_{10}$. The equilibrium radius of bubble 2 is kept fixed, $R_{20} = 50$ µm, $\eta = 0.001$ Pa s, and the other parameters are as in Figure 3.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** The contributions of the force components, given by Equations (14), (15) and (17), to the net force $F_1 + F_2$ on the bubbleman as a function of the acoustic pressure amplitude $P_a$ for (a) $\eta = 0.001$ Pa s and (b) $\eta = 0.004$ Pa s. $R_{10} = 20$ µm, $R_{20} = 50$ µm, $f = 30$ kHz.</figcaption></figure>

## 4. Experiments

The experimental design of the proposed microswimmer requires controlling the approach and the contact of two gas bubbles in an infinite medium. To this end, we take advantage of an experimental setup that allows one to control the coalescence of two bubbles in an acoustic levitation chamber. The setup and the experimental technique have already been described in detail in one of our previous studies [29], so here we only describe shortly the main experimental activities. Figure 6 depicts the experimental setup that is used for the creation of a two-bubble microswimmer (bubbleman) as well as for the capture of its oscillation dynamics and its translational trajectory in two orthogonal planes.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Schematic of the experimental setup used for the creation of a two-bubble microswimmer (bubbleman) and for the capture of its oscillation and translation dynamics.</figcaption></figure>

An 8 cm-edge cubic tank is filled with a viscous fluid medium constituted of sodium alginate power mixed with water (at a concentration of 2 g/L), whose viscosity has been measured by using a Anton Paar rheometer (MCR302 equipped with parallel plate) at 4.3 mPa s. Single bubbles are nucleated by short laser pulses using a Nd: YAG pulsed laser (λ = 532 nm New Wave Solo PIVIII), focused by a lens set. This laser-nucleation system allows generating microbubbles with a radius ranging from 20 µm to 50 µm. A 30.8-kHz Langevin transducer, which is in contact (by means of an echographic gel) with the underside of the tank, generates an acoustic standing wave field inside the tank, the driving frequency of which corresponds to a certain acoustic resonant mode of the cubic cavity. The standing wave has at least one pressure antinode that is located approximately at the middle-top part of the tank. This pressure antinode is a stable equilibrium location for bubbles with radii smaller than the resonance radius at the driving frequency, which is estimated by Equation (44) to be about 110 µm.

The process of the creation of a bubbleman is illustrated by Figure 7. A first bubble is nucleated by the laser at a distance of a few millimeters from the pressure antinode. Due to the driving acoustic field, this bubble oscillates spherically and is subject to a primary radiation force that makes it move to the nearest pressure antinode and settle there; see the trapped bubble on the right side of Figure 7a.

A second bubble is then nucleated; see the bubble on the left side of Figure 7a. The trajectory of this bubble is first determined by the primary radiation force and then by the secondary radiation (interaction) force that acts between the two bubbles when they are in close proximity. The interaction force exerted on the first (right) trapped bubble results in a small translational motion towards the approaching (left) bubble, as one can see in Figure 7a. At the final stage of the bubble approach, they can either coalesce to form a bigger bubble or come into contact and remain in this state. The exact mechanism underlying the coalescence and the contact of two bubbles is not yet well understood. The transition from coalescence to contacting bubbles surely depends on the velocities of approach of the bubbles, the amplitude of their spherical oscillations their equilibrium radii and their surface contamination. Anyway, although coalescence can happen in the
levitation chamber, its absence is regularly observed and contacting bubble pairs often occur. Figure 7b–e show the last steps of the approach and the contact of two bubbles, which result in the creation of a bubbleman. Once the bubbleman is created (see Figure 7d), it moves slowly to a new trapping location, shown in Figure 7e, which results from a balance between the primary radiation force on the bubbleman and the buoyancy force.

<figure><img src="image_11.png" /><figcaption>**Figure 7.** (a) Bubble trajectories prior to the contact. A first bubble (right side) is trapped at a stable location in the acoustic standing wave. Once a second bubble is nucleated (left side), it moves due to the primary radiation force towards the same stable location. When the bubbles come to proximity, the radiation interaction force makes them come into contact. (b,c) The last moments prior to the contact. (d) The bubbles come into contact and remain in this state. (e) As soon as the bubbles come into contact, the bubbleman moves to a new stable location.</figcaption></figure>

Since the bubbleman is located at a pressure antinode, both bubbles continuously undergo spherical oscillations and experience an interaction force, which is the source of a net propulsion force. However, if the interaction force does not exceed the trapping force caused by the primary radiation force, the bubbleman remains trapped. When the applied acoustic pressure is increased, the interaction force also increases and the net propulsion force can exceed an energy barrier that provides the trap of the contacting bubbles. As a result, the bubbleman begins to move and leaves the pressure antinode. Figure 8 shows the trajectory of the bubbleman captured from two orthogonal views. One camera (Phantom v12.1 equipped with a 12× Navitar objective lens, Vision Research, USA) captures the motion of the bubbleman from the side. A frame size of 1280 × 800 pixels (magnification of 4.6 µm/px) is used at a 60-Hz frame rate. Backlight illumination is provided by a continuous light-emitted diode (LED) light source. This allows capturing with high accuracy the shape of the bubbleman and ensures that it is spherically oscillating; see Figure 8b. Another camera (Basler ac640–750 µm) captures the top view of the bubbleman motion with a 640 × 480 frame size (magnification of 22 µm/px) at a 60-Hz frame rate. Due to experimental constraints, it is not possible to use backlight for the top view observation. The bubbleman is visualized thanks to the scattered light from the LED light source. The exposure time of the camera is therefore increased significantly in order to capture a bright spot corresponding to the propeller; see Figure 8a. The bubbleman under study is constituted of two bubbles with equilibrium radii of $R_{10} = 60$ µm and $R_{20} = 25$ µm. The propeller exhibits quasi-circular orbits around the pressure antinode with a predominant
motion in the plane orthogonal to the direction of gravity. As one can see from the side view, the largest travelled distance is ~1 mm. From the top view, the diameter of the largest orbit is around 4 mm. The propeller velocity has been estimated to reach 0.8 cm/s, which corresponds to 47 body lengths per second. The motion of the bubbleman can last as long as 10 s. For an acoustic frequency of 30.8 kHz, this corresponds to a bubbleman performing spherical oscillations for more than 300,000 acoustic periods. It is worth noting that such velocities correspond to the fastest acoustically-powered microswimmers mentioned in the literature. At the end of the motion, the intensity of the propulsion force decreases and the bubbleman becomes trapped again, returning to its equilibrium location. This loss of propulsion efficiency probably comes from a change in the bubble equilibrium radii due to rectified diffusion, a phenomenon that occurs on long timescales. Such a size modification can change drastically the magnitude of the net propulsion force, resulting then in a lack of motion. This effect can be prevented by coating bubble surfaces, which should enhance the size stability of the bubbles. At this stage, we cannot easily compare the experimental trajectories to the theoretical modelling. The main reason is that the theoretical modelling allows the calculation of the net propulsion force but not of the motion of the system within the levitation chamber. In order to describe the trajectories of the bubbleman, it is necessary to derive the equation of motion of the two-bubble system in the three-dimensional acoustic field by including the forces acting on the system: the buoyancy, the drag force, the radiation interaction force, and the primary radiation force (induced by the ultrasound field onto the bubbleman itself). This requires knowledge of the 3D pressure field around the trapping location, and relevant experiments are currently in progress in order to quantitatively determine the pressure field in the levitation chamber and hence to predict and control the trajectory of the bubbleman.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** (a) Trajectory of the bubbleman captured from above. The diameter of the quasi-circular trajectory can reach 4 mm. (b) A zoom of the trajectory captured from the side. The bubbleman is clearly visible, showing that both bubbles remain spherical during the motion.</figcaption></figure>

**5. Conclusions**

In this paper, a theory has been developed that allows one to analytically calculate the acoustic radiation interaction forces between two gas bubbles in an incompressible viscous liquid for any small separation distance between the bubbles. This theory has been used to demonstrate that two acoustically excited bubbles experience a nonzero net force when they come into contact. This result suggests a new mechanism that can be used for the development of artificial self-propelled microswimmers actuated and controlled by an acoustic field.
**Author Contributions:** Conceptualization, A.A.D. and C.I.; methodology, A.A.D. and C.I.; software, A.A.D.; validation, T.M. and C.I.; formal analysis, C.I.; investigation, T.M. and C.I.; resources, C.M. and C.I.; writing–original draft preparation, A.A.D. and C.I.; writing–review and editing, C.M., P.B.-B. and C.I.; visualization, A.A.D., T.M. and C.I.; supervision, C.M., P.B.-B. and C.I.; project administration, C.M. and C.I. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by l'Agence Nationale de la Recherche (ANR) under the project ANR-22-CE92-0062. For the purpose of open access, the authors have applied a CC-BY public copyright license to any Author Accepted Manuscript (AAM) version arising from this submission.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** This work was supported by the LabEx CeLyA of the University of Lyon (No. ANR-10-LABX-0060/ANR-11-IDEX-0007).

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A. Mathematical Formulas Used in Calculations

Here, we provide auxiliary mathematical formulas that are used in our derivation [23],

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

$$ \mu P_{n}(\mu) = \frac{n}{2n+1} P_{n-1}(\mu) + \frac{n+1}{2n+1} P_{n+1}(\mu), \qquad (A2) $$

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

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

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

## References

1. Bunea, A.-I.; Taboryski, R. Recent advances in microswimmers for biomedical applications. *Micromachines* **2020**, *11*, 1048. [CrossRef]

2. Sridhar, V.; Park, B.-W.; Sitti, M. Light-driven Janus hollow mesoporous TiO2-Au microswimmers. *Adv. Funct. Mater.* **2018**, *28*, 1704902. [CrossRef]

3. Alapan, Y.; Yigit, B.; Beker, O.; Demirörs, A.F.; Sitti, M. Shape-encoded dynamic assembly of mobile micromachines. *Nat. Mater.* **2019**, *18*, 1244–1251. [CrossRef]

4. Ceylan, H.; Yasa, I.C.; Yasa, O.; Tabak, A.F.; Giltinan, J.; Sitti, M. 3D-printed biodegradable microswimmer for theranostic cargo delivery and release. *ACS Nano* **2019**, *13*, 3353–3362. [CrossRef]

5. Louf, J.F.; Bertin, N.; Dollet, B.; Stephan, O.; Marmottant, P. Hovering microswimmers exhibit ultrafast motion to navigate under acoustic forces. *Adv. Mater. Interfaces* **2018**, *5*, 1800425. [CrossRef]

6. Bertin, N.; Spelman, T.A.; Stephan, O.; Gredy, L.; Bouriau, M.; Lauga, E.; Marmottant, P. Propulsion of bubble-based acoustic microswimmers. *Phys. Rev. Appl.* **2015**, *4*, 064012. [CrossRef]

7. Feng, J.; Yuan, J.; Cho, S.K. Micropropulsion by an acoustic bubble for navigating microfluidic spaces. *Lab Chip* **2015**, *15*, 1554–1562. [CrossRef]

8. Ahmed, D.; Lu, M.; Nourhani, A.; Lammert, P.E.; Stratton, Z.; Muddana, H.S.; Crespi, V.H.; Huang, T.J. Selectively manipulable acoustic-powered microswimmers. *Sci. Rep.* **2015**, *5*, 9744. [CrossRef]

9. Jeong, J.; Jang, D.; Kim, D.; Lee, D.; Chung, S.K. Acoustic bubble-based drug manipulation: Carrying, releasing and penetrating for targeted drug delivery using an electromagnetically actuated microrobot. *Sens. Actuators A* **2020**, *306*, 111973. [CrossRef]

10. Ahmed, D.; Mao, X.; Shi, J.; Juluri, B.K.; Huang, T.J. A millisecond micromixer via single-bubble-based acoustic streaming. *Lab Chip* **2009**, *9*, 2738–2741. [CrossRef]

11. Gao, Y.; Wu, M.; Luan, Q.; Papautsky, I.; Xu, J. Acoustic bubble for spheroid trapping, rotation and culture: A tumor-on-a-chip platform (ABSTRACT platform). *Lab Chip* **2022**, *22*, 805–813. [CrossRef]
12. Luo, T.; Wu, M. Biologically inspired micro-robotic swimmers remotely controlled by ultrasound waves. *Lab Chip* **2021**, *21*, 4095–4103. [CrossRef]

13. Pak, O.S.; Zhu, L.; Brandt, L.; Lauga, E. Micropropulsion and microrheology in complex fluids via symmetry breaking. *Phys. Fluids* **2012**, *24*, 103102. [CrossRef]

14. Bjerknes, V.F.K. *Fields of Force*; Columbia U. P.: New York, NY, USA, 1906.

15. Crum, L.A. Bjerknes forces on bubbles in a stationary sound field. *J. Acoust. Soc. Am.* **1975**, *57*, 1363–1370. [CrossRef]

16. Doinikov, A.A. Acoustic radiation forces: Classical theory and recent advances. In *Recent Research Developments in Acoustics*; Transworld Research Network: Trivandrum, India, 2003; Volume 1, pp. 39–67.

17. Doinikov, A.A. Dissipative effects on Bjerknes forces between two bubbles. *J. Acoust. Soc. Am.* **1997**, *102*, 747–751. [CrossRef]

18. Doinikov, A.A. Bjerknes forces between two bubbles in a viscous fluid. *J. Acoust. Soc. Am.* **1999**, *106*, 3305–3312. [CrossRef]

19. Doinikov, A.A.; Regnault, G.; Mauger, C.; Blanc-Benon, P.; Inserra, C. Acoustic microstreaming produced by two interacting gas bubbles undergoing axisymmetric shape oscillations. *J. Fluid Mech.* **2022**, *931*, A19. [CrossRef]

20. Shaw, S.J. Translation and oscillation of a bubble under axisymmetric deformation. *Phys. Fluids* **2006**, *18*, 072104. [CrossRef]

21. Guédra, M.; Inserra, C. Bubble shape oscillations of finite amplitude. *J. Fluid Mech.* **2018**, *857*, 681–703. [CrossRef]

22. Doinikov, A.A. Acoustic radiation pressure on a compressible sphere in a viscous fluid. *J. Fluid Mech.* **1994**, *267*, 1–21. [CrossRef]

23. Abramowitz, M.; Stegun, I.A. *Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables*, 10th ed.; Dover: New York, NY, USA, 1972.

24. Varshalovich, D.A.; Moskalev, A.N.; Khersonskii, V.K. *Quantum Theory of Angular Momentum*; World Scientific: Teaneck, NJ, USA, 1988.

25. Zwillinger, D. *Standard Mathematical Tables and Formulae*, 31st ed.; CRC: Boca Raton, FL, USA, 2003.

26. Landau, L.D.; Lifshitz, E.M. *Fluid Mechanics*, 2nd ed.; Pergamon Press: Oxford, UK, 1987.

27. Doinikov, A.A.; Cleve, S.; Regnault, G.; Mauger, C.; Inserra, C. Acoustic microstreaming produced by nonspherical oscillations of a gas bubble. I. Case of modes 0 and *m*. *Phys. Rev. E* **2019**, *100*, 033104. [CrossRef]

28. Prosperetti, A. Viscous effects on perturbed spherical flows. *Quart. Appl. Math.* **1977**, *34*, 339–352. [CrossRef]

29. Cleve, S.; Guédra, M.; Inserra, C.; Mauger, C.; Blanc-Benon, P. Surface modes with controlled axisymmetry triggered by bubble coalescence in a high-amplitude acoustic field. *Phys. Rev. E* **2018**, *98*, 033115. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.