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

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.


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

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 , θ 1 , ε 1 ) and (r 2 , θ 2 , ε 2 ), which originate at the equilibrium centers of bubbles 1 and 2, respectively, and in which the direction θ 1 = θ 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. Let us assume that the jth bubble (j = 1, 2) undergoes N j modes with a frequency ω.

Let these modes have numbers
Then, the perturbation of the bubble surface produced by these modes can be represented by where r (j) s 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, µ j = cos θ j , and s (j) n 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] where S j0 is the surface of the jth bubble at rest, n j is the unit outward normal to S j0 , ρ is the liquid density, v is the first-order liquid velocity generated by the bubbles, η is the dynamic liquid viscosity, v E is the Eulerian velocity of acoustic streaming, p E is the time-averaged pressure corresponding to v E , and means the time average. The total first-order liquid velocity generated by both bubbles is given by where v (j) is the first-order liquid velocity generated by the jth bubble. The components of v (j) are shown in [19] to be calculated by v (j) where h (1) n is the spherical Hankel function of the first kind, k v = (1 + i)/δ is the viscous wavenumber, δ = √ 2ν/ω is the viscous penetration depth, ν = η/ρ is the kinematic liquid viscosity, the prime denotes the derivative with respect to the argument in brackets (such as the term h (1)/ n k v r j in Equation (5)), and P 1 n is the associated Legendre polynomial of the first order and degree n. Expressions for the constants a (j) n and b (j) n , 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 , θ 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 where V (j) Here, C nm = (n + m)!/(n!m!), ξ j = R j0 /d, j n is the spherical Bessel function, and 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 θn 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 where Ψ (j) n (r j ) and Ψ (j)/ n (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 The contribution of the second term is calculated as follows: Making use of the equation the calculation of the third term results in 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 θ 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.
where the constants C (j) 410 are calculated by Equations (C30)-(C33) of [19]. With the help of Equations (18)-(20), one obtains Combining Equations (14), (15), (17) and (21) brings the following expression for the radiation interaction force on the jth bubble: where C (j) 110 and C (j) 410 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.

Linear Scattering Coefficients When Parametric Excitation Is Absent
If the parametric excitation of shape modes is absent, a The linearized equations of an incompressible viscous liquid are given by [26] ∇ · v = 0, ∂v ∂t where p is the first-order liquid pressure.
In the case of two bubbles, a solution for v is sought as where with the scalar, ϕ (j) , and the vector, ψ (j) , velocity potentials defined by [27] where e εj is the unit azimuth vector of the jth bubble. Note that axial symmetry allows us to set ε 1 = ε 2 and e ε1 = e ε2 . Substituting Equation (27) into Equation (25) and taking into account that ψ (j) obeys the equation (∆ + k 2 v )ψ (j) = 0 [27], one finds that the first-order scattered pressure generated by the jth bubble is given by 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, 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)πR 3 j0 is the equilibrium volume of the jth bubble, V j (t) is the instantaneous volume of the jth bubble, γ is the specific hear ratio of the gas, P 0 is the hydrostatic pressure in the liquid, P ac (t) = P a e −iωt is the imposed acoustic pressure, and p (j) st is the pressure of surface tension on the surface of the jth bubble, given by [28] p (j) where σ is the surface tension coefficient. Accurate to first-order terms, V j (t) is calculated by To apply Equation (31) and thus calculate a n , we need to know ϕ (1) and v (1) r in the coordinates (r 2 , θ 2 ) and ϕ (2) and v (2) r in the coordinates (r 1 , θ 1 ). To this end, we use Equations (A5), (A6), (A11) and (A13) of [19], which give With the help of Equations (30) and (32)-(37), applying Equation (31) at j = 1 and j = 2, one obtains for n = 0, and for n ≥ 2, where x j = k v R j0 , is the natural frequency of the radial mode of the jth bubble, is the natural frequency of the nth shape mode of the jth bubble (n ≥ 2), It also follows from Equation (31) that 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, and for n ≥ 1, where the coefficients f n , α jnm , β jnm and γ jnm are calculated by Equations (A25)-(A34) of [19].

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.

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: ρ = 1000 kg/m 3 , η = 0.001 Pa s, σ = 0.0727 N/m, γ = 1.4, P 0 = 101.3 kPa, P a = 10 kPa, f = ω/2π = 30 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 (j) n for two contacting bubbles (bubbleman) with R 10 = 20 µm and R 20 = 50 µ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. n for a bubbleman with R 10 = 20 µm and R 20 = 50 µm in the case that parametric excitation is absent. The calculation was carried out by the equations derived in Section 2.2. 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 η. 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 → 0 .  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, η = 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, η = 0.001 Pa s and η = 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 (j) 1 given by Equation (14), and the long-dash curve shows the contribution of the terms T  (15) and (17). As one can see, at η = 0.001 Pa s, the dominant contribution comes from the terms T    (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) η = 0.001 Pa s and (b) η = 0.004 Pa s. R 10 = 20 µm, R 20 = 50 µm, f = 30 kHz.

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. 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. 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 used 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.

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.

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

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