Self radiation force on a moving monopolar source

The radiation force exerted on an object by an acoustic wave is a widely studied phenomenon since the early work of Rayleigh, Langevin and Brillouin and has led in the last decade to tremendous developments for acoustic micromanipulation. Despite extensive work on this phenomenon, the expressions of the acoustic radiation force applied on a particle have so far been derived only for a steady particle, hence neglecting the effect of its displacement on the radiated wave. In this work we study the acoustic radiation force exerted on a monopolar source translating at a constant velocity small compared to the sound speed. We demonstrate that the asymmetry of the emitted field resulting from Doppler effect induces a radiation force on the source opposite to its motion.


Introduction
Since the seminal work by Rayleigh (1902Rayleigh ( , 1905 and Langevin (work reported later by Biquard (1932a,b)), much effort has been devoted to the derivation of theoretical expressions of the acoustic radiation force exerted by an acoustic wave on a particle. Brillouin (1925b,a) was the first to recognize the tensorial nature of the acoustic radiation force, which is not necessarily orthogonal to the insonified interface. Later on, King (1934) derived an expression of the axial acoustic radiation force exerted on a rigid spherical particle by a plane wave. This expression was extended to the case of a compressible fluid sphere and of an elastic particle by Yosika & Kawasima (1955) and Hasegawa & Yosika (1969) respectively. The case of an incident focused wave was treated by Embleton (1954) for a rigid particle and Chen & Apfel (1996) for an elastic particle. Later on, the more general case of the axial force exerted by a Bessel beam was addressed by Marston (2006Marston ( , 2009. In parallel, a general expression of the acoustic radiation force exerted by an arbitrary wavefield on a spherical particle in the Long Wavelength Regime (LWR i.e. when ka 1 with k the wavenumber and a the particle radius) was obtained by Gork'ov (1962). It was shown up to 3rd order (in ka) that the radiation force is proportional to the gradient of an acoustic potential, which is proportional to the difference between the time-averaged potential and kinetic acoustic energy weighed respectively by the monopole and dipole scattering coefficients. This expression was extended to the 6th order by Sapozhnikov & Bailey (2013), which is necessary when the average potential and kinetic energy are uniform in space (e.g. for plane propagating waves). This expression was also extended by Doinikov (1997a,b,c) to consider the effect of the viscous and thermal boundary layers in some asymptotic limits (of the boundary layer size compared to the particle size) and in the general case by Settnes & Bruus (2012) and Karlsen & Bruus (2015). The case of nonspherical particles such as disks and spheroids was treated by Keller (1957) and Silva & Drinkwater (2018) respectively, while the case of transient acoustic fields was recently addressed by Wang et al. (2021). We can also note that some expression for the secondary radiation force (inter-particle force) have been derived by Silva & Bruus (2011). The specific case of the radiation force exerted on a vibrating bubble known as the primary Bjerkness force was treated separately by Bjerknes (1906), Blake (1949), Eller (1968) and Crum (1975). Indeed, bubbles have some specificity: owing to their strong compressibility compared to the surrounding liquid, their monopolar resonance appears in the LWR. Hence bubbles can be attracted to the nodes or anti-nodes of a standing wave depending whether they are forced below or above their monopolar resonance frequency (see Eller (1968)).
Recently, there have been some renewed interest in the calculation of the acoustic radiation force with the development of selective acoustical tweezers (see  for a review of the subject). Acoustical tweezers rely on the acoustic radiation force to move objects. To reach selectivity, i.e. the ability to manipulate a single object independently of other neighboring objects, it is necessary to localize the acoustic energy close to the target particle to only affect it (see e.g. Baresch et al. (2016), Baudoin et al. (2019), ). Hence, such selectivity cannot be reached in the LWR. In addition, to calculate the restoring force (i.e. the force which brings back the particle toward the trap center), it is necessary to compute the radiation force when the particle is out-centered from the trap position. Yet, all the aforementioned expressions were either limited to the calculation of the radiation force for an axisymmetric configuration or to the LWR. To cope with this issue, general expressions of the radiation force exerted by an arbitrary incident field on a spherical particle, without restriction on the particle size compared to the wavelength, were obtained by Sapozhnikov & Bailey (2013) with an angular spectrum based method and by Silva (2011), and Baresch et al. (2013) with a multipole expansion method. The equivalence between the formula obtained with the different approaches was demonstrated by . Note also that some general expressions have been proposed to compute the acoustic radiation torque exerted by an arbitrary acoustic field on a particle of arbitrary size by Silva et al. (2012) and .
Yet, in all the theoretical developments mentioned so far, the calculation of the radiation force is made for a steady particle, hence neglecting the effect of its motion. Some account of the interaction between the oscillatory and translational motion of bubbles can be found in the literature. Following Saffman (1967), Benjamin & Ellis (1990) showed that even in an inviscid fluid, self-propulsion of a bubble can be achieved by nonlinear interactions between adjacent surface deformation modes. In their work however, the surface modes deformation are supposed to be known a priori. This work was extended later on by Mei & Zhou (1991) to account for the parametric excitation of surface modes by the isotropic volume mode and Feng & Leal (1995) who considered the direct coupling between translational motion, volume and shape modes. Finally, Doinikov (2004) obtained an expression of this coupling whatever the shape modes, their natural frequency, and the type of excitation (parametric forcing by the volume mode or direct excitation through externally induced pressure gradients at the surface of the bubble). In parallel it was also shown by Watanabe & Kukita (1993) and Doinikov (2002) that even if the bubble oscillations remain spherical, some complex coupling between volumetric Figure 1: Sketch illustrating the asymmetry of the acoustic field synthesized by a translating monopolar source. The normalized field is calculated with eq. (2.21) and for the sake of illustration, the asymmetry is magnified by choosing a Mach number M = 0.5.
oscillations and translational motion leading to erratic motion of the bubble can still occur when the bubble is excited by acoustic standing waves of high intensity. Finally, we can mention the work of Magnaudet & Legendre (1998) who computed how the viscous drag applied on a translating bubble is modified by its oscillation.
But none of these works considered the effect of the asymmetry of the acoustic wave radiated by a translating source on the acoustic radiation force. In this paper, we consider a monopolar source translating in a quiescent inviscid fluid at a constant velocity U along a fixed axis and demonstrate that the asymmetry of the acoustic field due to Doppler effect ( Figure 1) induces a self-induced radiation force on the source resisting its motion. This result is obtained by inserting the well known solution of the wavefield radiated by a moving monopolar source into a far-field integral expression of the radiation force exerted on a moving source, and finally computing this integral within the approximation of slow translating speed compared to the sound speed.

Wavefield radiated by a translating monopolar source
The first step to compute the self-induced radiation force exerted on a moving monopolar source is to compute the wavefield radiated by this source in a fixed reference frame. This classic calculation can be found in the acoustics textbook of Morse & Ingard (1968). In this first section, we recall the main steps of the derivation. Here we suppose the fluid to be inviscid.

Wave equation for a translating monopolar source
In acoustics, a monopolar source can be seen as a source of mass, whose strength is specified by the instantaneous mass flow rate q(t) created by this source. In the following, the source is supposed to be periodic of period T. For a punctual source translating at a velocity U = U x = M c o x along a fixed axis x, the mass and momentum conservation equations become: Galilean reference frame, M the mach number, c o the sound speed, ρ the density, v the fluid velocity,σ the stress tensor equal to −pĪ for an inviscid fluid,Ī the identity tensor and p the pressure. If we (i) make the classic asymptotic development of equations (2.1) and (2.2) up to first order: with 1, and obtain the linearized mass and momentum balance: (ii) introduce the sound speed: where s is the entropy and (iii) combine the time derivative of equation (2.4) with the divergence of equation (2.5), we obtain the wave equation: (2.7) Note that to ease the resolution of this problem, it is convenient to introduce the velocity potential ψ 1 defined by v 1 = − 1 ρ0 ∇ψ 1 such that p 1 = ∂ψ1 ∂t , which enables to suppress the time derivative in the rhs of equation (2.7): and will make the future change of variables easier.

Resolution of the wave equation and Lorentz transformation
The solution of the wave equation (2.8) is well-known for a fixed monopolar source (M = 0): with r = x 2 + y 2 + z 2 the radial distance. To solve the problem for the moving source, the idea is to rewrite equation (2.8) in a reference frame wherein the source is fixed and the wave equation remains unchanged. This can be achieved by using the invariance of the wave equation by the Lorentz transformation, which is at the core of special relativity: where γ defined by γ −1 = √ 1 − M 2 is the "Lorentz acoustic boost". With this transformation, the wave equation (2.8) becomes: Since the rhs of equation (2.11) is null when x = 0 for all t and using δ(x /γ) = γδ(x ), it comes: If we now introduce a second set of variables: the wave equation becomes: (2.14) which now resembles the static monopolar source problem and whose solution is: If we now perform the inverse transformations to obtain the potential as a function of (x, y, z, t), we obtain: with: If we introduce the distance R ± between the emission and the observation points: then the solution becomes: Figure 2: S represents the source surface, varying over time. The surface S ∞ is centered on the source and moves with it at the velocity U in R. R * is the frame of the source.
with R = R + in the subsonic case (M < 1).

Integral expression of the radiation stress in the far field
The next step is to derive a far-field integral expression of the radiation stress exerted on a moving source. Indeed in acoustics, a monopolar point source constitutes a far-field approximation of a real source of finite extent and hence the above expressions are only valid in the far field.

Far field expression of the radiation force for a moving source
The acoustic radiation stress exerted on an object of surface S(t) is by definition the time average of the surface integral of the stress exerted by the acoustic wave on its surface: where f = 1 T t+T t f (t)dt is the time average of the function f and T the period of the function f (t). In general, there are two difficulties when computing this integral: (i) the surface of the object is vibrating and hence depends on time (S = S(t)) and (ii) an expression of the wave scattered by the object in the near field must be known. Hence generally this integral is converted into an integral over a closed surface at rest surrounding the object in the far field by using the divergence theorem and Reynolds transport theorem (see e.g. the review by  for details of this process). Here an additional difficulty comes from the fact that the particle, in addition to its vibration, is translating at a constant velocity U . To solve this issue, we will transpose our integral of the stress on the surface of the object into an integral over a spherical surface S ∞ of radius r ∞ λ, centered on the source, and hence translating at the velocity U in R (see Figure 2), with λ = c 0 /f the wavelength and f the frequency. The volume between S and S ∞ is named V. The integral of equation (2.2) over V gives: Using the divergence theorem, this volume integral turns into: with n and n ∞ the outgoing normal vectors to the surface S and S ∞ respectively (see Figure 2). Another equation can be obtained by applying the Reynolds transport theorem to the momentum density ρv: since the surface S follows the object surface displacement (equal to the fluid displacement at the interface due to the continuity condition) and the surface S ∞ is translating at a constant velocity v. If we study the steady regime, then: Finally, if we substract the time average of equation (3.25) to the time average of equation (3.24) and take into account (3.26), we obtain: (3.27)

Expression as a function of the first order acoustic field
In order to compute the previous integrals and since the time average of first order terms are equal to 0 for a periodic signal, we need to express the terms appearing in equation (3.27) up to second order. Here we suppose again that the fluid is inviscid and henceσ = −pĪ. If we push the asymptotic development up to second order: the momentum balance (2.2) at second order in the volume V becomes: If we take the time-average of this equation, we obtain: since ∂v2 ∂t = 0 in the steady regime. From the first order momentum conservation (2.5), we have ∂v1 ∂t = − 1 ρ0 ∇p 1 and from the state equation (2.6) we have ρ 1 = p 1 /c 2 0 leading to: And since the acoustic field is by definition irrotationnal: Finally, by replacing (3.31) and (3.32) in (3.30), we obtain: From this equation, we obtain the second order averaged stress tensor: Then at second order, we have: so that if we replace equations (3.34), (3.35) and (3.36) in (3.27), we obtain the final expression of the radiation force as a function of the first order acoustic fields: (3.37)

Expression of the radiation force for a slowly translating source
The final step is to compute this integral from the first order acoustic field.

Pressure, density and velocity fields at first order
The first step to compute equation (3.37) is to determine the pressure, density and velocity fields at first order from the expression of the velocity potential determined previously (equation (2.21)). First, (4.38) Moreover, the coordinate of the velocity v 1i along the direction x i = (x, y, z) is: (4.39) Of course, we have: ρ 1 = p 1 /c 2 0 . Here R and R 1 are time and space dependent functions. To simplify the calculation of the derivatives and the force integral, we perform the change of variable corresponding to the Galilean transformation from R to R * : R 1 and R become: Figure 3: We make the change of variables corresponding to the Galilean transformation from R to R * and then use the local spherical coordinates (r * , θ * , ϕ * ).
The Jacobian matrix J * of this transformation is: (4.42) We use the spherical coordinates (r * , θ * , ϕ * ) to compute the integral (see Figure 3 for the notations): The key point in the following derivation of the self radiation force is that we consider M 1, which means that the source translates at low velocity compared to the sound speed. This enables us to get simpler expressions of the pressure and velocity fields. First, we derive the expressions of R 1 and R at the first order in M : (4.45) In order to get the velocity and the pressure fields at the first order in M , we compute at first the different time derivatives of R 1 and R: (4.46) with: Hence we obtain: (4.48) Then: Finally, the acoustic pressure is: (4.50) We choose a sphere S ∞ with a radius r ∞ huge compared to all the other lengths of the problem. We then write the acoustic pressure in the far field approximation. At first order in M we obtain: In order to derive the velocity field we compute the different space derivatives of R 1 and R: and (4.53) We obtain: We finally write the velocity field in the far field approximation: (4.55) (4.56) (4.57)

Calculation of the integrals
With the change of variables (4.40) and since the Jacobian of this change of variable is equal to 1, the integral expression of the radiation force (3.37) becomes: (4.58) with S ∞ * the surface defined by r * = r * ∞ a constant radius, the infinitesimal element of surface dS * = r * ∞ 2 sin θ * dθ * dϕ * , θ * ∈ [0, π] and ϕ * ∈ [0, 2π]. Finally, n * = n ∞ * = e * r , with (e * r , e * θ , e * ϕ ) the spherical coordinates unit vector of R * . Since we have already expressed the pressure, velocity and density fields as a function of the spherical coordinates (r * , θ * , ϕ * ), we can now compute integral (4.58). In the following paragraphs, we compute separately each term of this integral.

Potential energy term
The integration over ϕ * along the y * and z * axis is null. Hence only the term along x * remains: (4.60) If we swap the time and space integrals and since the function q(t) is periodic, we obtain: Since the last integral is equal to 4M/3, we obtain: (4.62)

Kinetic energy term
Using the same arguments than previously we have: Since the last integral is equal to 4M/3, we obtain:

Convective term
Due to the symmetry of the problem (invariance by rotation over ϕ * angle), no force can exist along y and z directions. Hence, we only need to compute the following components of v 1 ⊗ v 1 : Then due to the dependence of these terms over ϕ * given in equations (4.55) to (4.57) and since e * r = cos θ * x + sin θ * cos ϕ * y + sin θ * sin ϕ * y, we obtain: The integral term is equal to 8M/3, so that:

Source translation term
The term due the translation of the sphere S ∞ is: Since at leading order, the last integral is equal to 2/3, we finally obtain:

Final expression of the self-radiation force
If we now replace (4.62), (4.65), (4.68) and (4.69) in (4.58), we obtain the final expression of the self acoustic radiation force exerted on a monopolar source: There are many interesting things to notice in the above calculations and final expression. First, we see that the potential and kinetic energy terms cancel at leading order in M, so that the self radiation force is solely due to the convective and translation terms. Second, the radiation force is proportional to the radiated intensity, and inversely proportional to the sound speed square, which is classical in radiation force calculations. In addition, here the self-radiation force is proportional to the hydrodynamic Mach number M , which is expected since the force results from the asymmetry of the radiated field due to the translation of the source. Finally and most importantly F rad · U is always negative, which means that this force always slows down the movement of the bubble.

An example of a monopolar oscillator: a vibrating bubble
In the following section, we estimate this force for a translating and oscillating bubble in a liquid, which constitutes one example of acoustic monopolar source. Indeed, bubbles are exceptional resonators, which exhibit strong monopolar resonances in the long wavelength regime. Let's consider a spherical bubble of mean radius r b in a liquid of density ρ 0 and sound speed c 0 vibrating periodically at its resonance frequency, called the Minnaert frequency: with γ the heat capacity ratio of the gas in the bubble and p 0 the pressure of the surrounding fluid. At resonance, λ/r b = 2πc0 ω M r b = 2πc 0 ρ0 3γp0 1 (of the order of 5 × 10 2 for an air bubble in water). Since the bubble is very small compared to the wavelength in this case, it can be considered as a point source a few wavelengths away from the bubble surface. The oscillation of this bubble creates a periodic mass flux q(t) = Q cos(ω M t), whose magnitude Q is basically equal to the surface of the bubble 4πr 2 b , times the surrounding liquid mass density ρ 0 , times the amplitude of the oscillations αr b (where α designates a dimensionless parameter fixing the magnitude of the bubble oscillation), times the pulsation ω b : Consequently, we have: and: For small bubbles, it is interesting to compare this radiation force to the Stokes drag: with C = 4 for a bubble in a pure liquid moving at low Reynolds number, C = 6 if the surface is polluted so that the slip boundary condition is turned into a no slip boundary condition (see e.g. the book Kim & Karilla (2005)) and C = 12 for an undeformed bubble at large Reynolds numbers (see Moore (1963)). If we compare equation (5.74) to (5.75), we obtain: | F rad | |F Sto | ∼= 18 C γ 2 p 2 0 α 2 µρ 0 c 3 0 r b . (5.76) For a bubble rising in water at ambient temperature and pressure, we have ρ 0 ∼1 × 10 3 kg m −3 , µ ∼1 × 10 −3 Pa s, c 0 ∼1.5 × 10 3 m s −1 , γ ∼ 1.4, p 0 =1 × 10 5 Pa, so that with α ∼ 0.5, we obtain: | F rad | |F Sto | ∼ 7r b . (5.77) Hence, the self radiation force would be small compared to the Stokes drag for a millimetric bubble. However the self radiation force could become significant in cryogenic liquids such as liquid nitrogen or superfluid helium. In liquid nitrogen at T 77 K, ρ 0 =8 × 10 2 kg m −3 , µ =2 × 10 −4 Pa s, c 0 =8 × 10 2 m s −1 , γ ∼ 1, p 0 =1 × 10 5 Pa, so that for α ∼ 0.5 we obtain: (5.78) which means that for a bubble of a few millimeters in radius, the two phenomena would be of the same orders of magnitude. Note first that this calculation constitutes a rough comparison of the radiation force and Stokes drag since (i) for large bubble oscillations, the bubble dynamics becomes nonlinear, and (ii) the Stokes drag is also modified by the bubble oscillations as demonstrated by Magnaudet & Legendre (1998). Note also that the case of bubbles moving in a liquid constitutes just one possibility over the many configurations covered by equation (5.70) since the above theory applies for an arbitrary monopolar source moving in an arbitrary fluid as soon as (i) the monopolar source emits a signal in the long wavelength regime and (ii) that the source is moving at low speed compared to the sound speed.

Conclusion
In this paper, we calculated the theoretical expression of the radiation force exerted on a translating monopolar source by its own acoustic field. We showed that the asymmetry of the radiated field due to the motion of the source creates a self-induced radiation force opposite to its motion. This first calculation opens many perspectives to calculate this effect for different types of sources, beyond the long wavelength regime and potentially unveil situations wherein this force is strong and could be in the opposite direction leading to a propulsive force, and hence with a scatterer surfing on its own acoustic field.