Acoustically mediated long-range interaction among multiple spherical particles exposed to a plane standing wave

After the pioneering discovery of optical trapping [1, 2], Burns et al [3] observed a series of optically-induced bound states between two polystyrene spheres. Different from the optical trapping generated by the gradient of an external field directly, the optical binding effect occurs even in a plane wave field without any gradient. Since then, extensive efforts have been devoted to the light-mediated particle-particle interaction [4–12]. The optical binding force alternates between attraction and repulsion with the inter-particle separation, which leads to a light-mediated self-organization of particle clusters. The interaction between objects can also be induced by acoustic waves, as reported theoretically [13–24] and experimentally [25–28]. Comparing with its optical counterpart, less progress has been made on the sound-mediated interaction. To the best of our knowledge, most of the current theoretical studies focus on the acoustic interaction between two objects [14–20].

The acoustic radiation force (ARF), as a second order quantity of the acoustic field, relies on an accurate calculation of the self-consistent field distribution in the presence of the particles. The self-consistent sound field can be solved by various techniques, such as the finite-element method, the finite difference time domain simulation, and the multiple-scattering theory (MST). For spherical particles, the MST approach has been proved to be the most efficient since the sound field can be precisely captured by a finite number of spherical basis functions. It has been successfully developed to calculate the acoustic interaction between two spherical particles [14, 15]. Although a system involving multiple particles [21–23] can also be handled, the MST method becomes inefficient as the growth of the particle number, owing to the overloaded memory and time consumption. Recently, for the particles much smaller than the acoustic wavelength (i.e., Rayleigh particles), a simple formula of the pairwise interaction force has been derived based on a scalar potential theory [24]. Together with a mean-field approximation, this approach is further demonstrated to be powerful in treating a great number of the Rayleigh particles. It is worth pointing out that, the acoustic interaction among the Rayleigh particles is very weak, which could be hidden in the gradient force induced by a tiny defect of the external sound field.

Based on a single prior scattering approximation, here we present a theoretical study on the sound-mediated interactions among the multiple spherical particles distributed sparsely in the same node or antinode plane of a plane standing wave (PSW) field. Particularly, we focus on the particle size comparable with the acoustic wavelength, in which the acoustic interaction is anticipated to be much stronger than that in the Rayleigh situation. We derive a concise form for the total interaction force exerting on any given particle. It consists of a series of independent pairwise interactions between the particle and the others. Therefore, the multi-particle system is reduced into a two-body problem, which greatly simplifies the computation comparing with the rigorous MST method. Interestingly, we find that each pairwise interaction is a conservative force and thus the whole system could be described by a potential energy. This result is unusual for the particle beyond the Rayleigh limit. Besides, the pairwise force is oscillatory at a factor of cosine function and decays at a power exponent of n = −1 or −2. The oscillatory long-range interaction, controlled by the external sound field, could facilitate the self-organization of a two-dimensional (2D) cluster involving a large number of particles.

The remainder of this paper is organized as follows. In section 2, an analytical formula responsible for the long-distance acoustic interaction is derived, followed by a detailed discussion on the formula. In section 3 we check first the accuracy of our formula by some few-particle systems, and then study the system involving a large number of particles. Finally, a brief summary is made in the last section.

PSW is a good external field for exploring the acoustically mediated particle–particle interaction. The node or antinode plane forms a natural potential well to confine the particles tightly in the same plane without any transversal field gradient. Assume that the acoustic PSW field is described by ${{\rm{\Psi }}}_{{\rm{in}}}^{({e})}({\bf{r}})={{\rm{\Psi }}}_{0}({{\rm{e}}}^{{\rm{i}}kz}+{{\rm{e}}}^{-{\rm{i}}kz}),$ where ${{\rm{\Psi }}}_{0}$ characterizes the amplitude of the PSW field. In this case, the interaction force survives only in the xy plane owing to the symmetry of the external field.

As schematically depicted in figure 1, for an arbitrary spherical particle $i,$ the total velocity potential function can be written as a superposition of the incident wave ${{\rm{\Psi }}}_{{\rm{in}}}({\bf{r}})$ and the scattering wave ${{\rm{\Psi }}}_{{\rm{sc}}}({\bf{r}}),$ i.e.,

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\rm{in}}}({\bf{r}})=\displaystyle \sum _{lm}{a}_{l,m}^{i}{j}_{l}(k{r}_{i}){Y}_{l,m}({\hat{r}}_{i}),\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\rm{sc}}}({\bf{r}})=\displaystyle \sum _{lm}{b}_{l,m}^{i}{h}_{l}(k{r}_{i}){Y}_{l,m}({\hat{r}}_{i}),\end{eqnarray} \tag{ 1b }$$

where $k$ is the wavenumber in the host fluid, the vector ${{\bf{r}}}_{i}={r}_{i}{\hat{r}}_{i}$ refers to any probe position ${\bf{r}}$ measured from the center of a given sphere $i,$ and $\{{j}_{l}\},$ $\{{h}_{l}\}$ and $\{{Y}_{l,m}\}$ are the spherical Bessel functions, the spherical Hankel functions of the first kind and the spherical harmonics, respectively.

Zoom In Zoom Out Reset image size

The incident coefficient ${a}_{lm}^{i}$ can be decomposed into two parts,

$$\begin{eqnarray}&&{a}_{l,m}^{i}={a}_{l,m}^{i(e)}+\displaystyle \sum _{j\ne i}{a}_{l,m}^{i(j)}.\end{eqnarray} \tag{ 2 }$$

The first part ${a}_{l,m}^{i(e)}$ stems from the external incidence (PSW field here),

$$\begin{eqnarray}&&{a}_{l,m}^{i(e)}={{\rm{\Psi }}}_{0}4\sqrt{(2l+1)\pi }{\delta }_{m0}\,{\rm{\cos }}\left(\alpha +\displaystyle \frac{l\pi }{2}\right),\end{eqnarray} \tag{ 3 }$$

with ${\delta }_{m0}$ being a Kronecker delta. The value of $\alpha$ depends on the location of the particle ensemble: $0\,$ for the antinode plane and $\pi /2$ for the node plane of the acoustic PSW field. The second part stems from the scattering of the other particles $j\ne i.$ Each component ${a}_{l,m}^{i(j)}$ can be derived according to the addition theorem of the spherical functions [29], i.e.,

$$\begin{eqnarray}&&{a}_{l,m}^{i(j)}=\displaystyle \sum _{l^{\prime} m^{\prime} }{b}_{l^{\prime} ,m^{\prime} }^{j}{G}_{l^{\prime} m^{\prime} lm}^{ij}.\,\end{eqnarray} \tag{ 4 }$$

Here ${G}_{l^{\prime} m^{\prime} lm}^{ij}=4\pi \displaystyle {\sum }_{l^{\prime\prime} }{{\rm{i}}}^{l+l^{\prime\prime} -l^{\prime} }{C}_{lml^{\prime\prime} (m^{\prime} -m)}^{l^{\prime} m^{\prime} }{h}_{l^{\prime\prime} }(k{R}_{ij}){Y}_{l^{\prime\prime} ,m^{\prime} -m}({\hat{R}}_{ij})$ with the coefficient ${C}_{lml^{\prime\prime} (m^{\prime} -m)}^{l^{\prime} m^{\prime} }=\displaystyle {\int }_{0}^{2\pi }\displaystyle {\int }_{0}^{\pi }{Y}_{l^{\prime} ,m^{\prime} }(\theta ,\phi ){Y}_{l,m}^{* }(\theta ,\phi ){Y}_{l^{\prime\prime} ,m^{\prime} -m}^{* }(\theta ,\phi ){\rm{\sin }}\,\theta {\rm{d}}\theta {\rm{d}}\phi .$ ${{\bf{R}}}_{ij}={R}_{ij}{\hat{R}}_{ij}$ represents the position vector of the spherical center $i$ measured from the spherical center $j$.

For a homogeneous spherical particle, the incident and scattering coefficients can be related by a diagonal scattering matrix $\{{t}_{l}^{i}{\delta }_{ll^{\prime} }\},$ i.e., ${b}_{l,m}^{i}={t}_{l}^{i}{a}_{l,m}^{i},$ according to the continuum boundary condition on the surface of the particle. Substituting this relation into equation (4) and using equation (2), we obtain a linear equation system

$$\begin{eqnarray}&&\displaystyle \sum _{jl^{\prime} m^{\prime} }({\delta }_{ij}{\delta }_{l{l}^{\prime }}{\delta }_{mm^{\prime} }-{t}_{l^{\prime} }^{j}{G}_{l^{\prime} m^{\prime} lm}^{ij}){a}_{l^{\prime} ,m^{\prime} }^{j}={a}_{l,m}^{i(e)}.\end{eqnarray} \tag{ 5 }$$

Solving this linear problem gives the incident coefficient ${a}_{l,m}^{i}$ and further gives the total self-consistent acoustic field according to equation (1). In this procedure, multiple scatterings among the particles are taken into account rigorously. However, the computational time and memory cost increase rapidly with the growth of the particle number. If the particles are distributed sparsely $k{R}_{ij}\gg 1,$ a single prior scattering approximation can be used to avoid solving the linear problem. Now the scattering contribution from the sphere $j$ to $i$ is dominated by a single prior scattering event, i.e., ${b}_{l,m}^{j}\approx {t}_{l}^{j}{a}_{l,m}^{j(e)},$ which finally gives rise to

$$\begin{eqnarray}&&{a}_{l,m}^{i(j)}\approx \displaystyle \sum _{l^{\prime} m^{\prime} }{t}_{l^{\prime} }^{j}{a}_{l^{\prime} ,m^{\prime} }^{j(e)}{G}_{l^{\prime} m^{\prime} lm}^{ij}.\end{eqnarray} \tag{ 6 }$$

Combining the equations (1), (2) and (6) with the relation ${b}_{l,m}^{i}={t}_{l}^{i}{a}_{l,m}^{i},$ we can obtain the field distribution in the presence of particles. The ARF exerted on any particle $i$ can be calculated by integrating the time-averaged radiation stress tensor $\langle \mathop{{\bf{S}}}\limits^{\leftrightarrow}\rangle$ over an arbitrary surface $S$ enclosing the particle, i.e.,

$$\begin{eqnarray}&&{{\bf{F}}}^{i}=\displaystyle {\oiint }_{s}\langle \mathop{{\boldsymbol{S}}}\limits^{\leftrightarrow}\rangle \cdot {\rm{d}}{\bf{A}}.\end{eqnarray} \tag{ 7 }$$

Specifically, the time-averaged radiation stress tensor $\langle \mathop{{\bf{S}}}\limits^{\leftrightarrow}\rangle$ can be written as

$$\begin{eqnarray}&&\langle \mathop{{\bf{S}}}\limits^{\leftrightarrow}\rangle =\displaystyle \frac{1}{2}{\rho }_{0}{\rm{Re}}({\bf{v}}* {\bf{v}})-\left(\displaystyle \frac{{\rho }_{0}{| {\bf{v}}| }^{2}}{4}-\displaystyle \frac{{| p| }^{2}}{4{\rho }_{0}{c}_{0}^{2}}\right)\mathop{{\bf{I}}}\limits^{\leftrightarrow},\end{eqnarray} \tag{ 8 }$$

with $\mathop{{\bf{I}}}\limits^{\leftrightarrow}$ being a unit tensor, ${\rho }_{0}$ and ${c}_{0}$ being the static mass density and sound velocity of the fluid background, and ${\bf{v}}$ and $p$ denoting the first-order velocity and pressure fields, respectively. If an infinitely large contour is selected, the integral form of the ARF ${{\bf{F}}}^{i}=({F}_{x}^{i},\,{F}_{y}^{i})$ can be rewritten as a series [30, 31]

$$\begin{eqnarray}{F}_{x}^{i}+{\rm{i}}{F}_{y}^{i} & = & \displaystyle \frac{{\rm{i}}{\rho }_{0}}{4}\displaystyle \sum _{lm}[{\mu }_{l+1,m-1}(2{b}_{l+1,m-1}^{i}{b}_{l,m}^{i* }+{b}_{l+1,m-1}^{i}{a}_{l,m}^{i* }+{a}_{l+1,m-1}^{i}{b}_{l,m}^{i* })\\ & & +{\mu }_{l+1,-m-1}(2{b}_{l,m}^{i}{b}_{l+1,m+1}^{i* }+{b}_{l,m}^{i}{a}_{l+1,m+1}^{i* }+{a}_{l,m}^{i}{b}_{l+1,m+1}^{i* })],\end{eqnarray} \tag{ 9 }$$

where the coefficient ${\mu }_{l,m}=\sqrt{[(l-m)(l-m-1)]/[(2l-1)(2l+1)]}.$ Using the relation ${b}_{l,m}^{i}={t}_{l}^{i}{a}_{l,m}^{i},$ equation (9) can be reorganized into

$$\begin{eqnarray}&&{F}_{x}^{i}+{\rm{i}}{F}_{y}^{i}=\displaystyle \frac{{\rm{i}}{\rho }_{0}}{4}\displaystyle \sum _{lm}[{T}_{l}^{i}{\mu }_{l+1,m-1}{a}_{l,m}^{i* }{a}_{l+1,m-1}^{i}+{T}_{l}^{i* }{\mu }_{l+1,-m-1}{a}_{l,m}^{i}{a}_{l+1,m+1}^{i* }],\end{eqnarray} \tag{ 10 }$$

where the factor ${T}_{l}^{i}=2{t}_{l+1}^{i}{t}_{l}^{i* }+{t}_{l+1}^{i}+{t}_{l}^{i* }$ characterizes the scattering property of the particle $i.$ Substituting equation (2) into equation (10), the acoustic interaction can be divided into three parts

$$\begin{eqnarray}&&{F}_{x}^{i}+{\rm{i}}{F}_{y}^{i}=\displaystyle \frac{{\rm{i}}{\rho }_{0}}{4}({Q}_{{\rm{e}}}+{Q}_{{\rm{s}}}+{Q}_{{\rm{c}}}),\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&{Q}_{{\rm{e}}}=\displaystyle \sum _{lm}[{T}_{l}^{i}{\mu }_{l+1,m-1}{a}_{l,m}^{i({e})}{a}_{l+1,m-1}^{i({e})}+{T}_{l}^{i* }{\mu }_{l+1,-m-1}{a}_{l,m}^{i({e})}{a}_{l+1,m+1}^{i({e})}],\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&{Q}_{{\rm{s}}}=\displaystyle \sum _{lm}\displaystyle \sum _{j\ne i,h\ne i}[{T}_{l}^{i}{\mu }_{l+1,m-1}{a}_{l,m}^{i(j)* }{a}_{l+1,m-1}^{i(h)}+{T}_{l}^{i* }{\mu }_{l+1,-m-1}{a}_{l,m}^{i(j)}{a}_{l+1,m+1}^{i(h)* }],\end{eqnarray} \tag{ 12b }$$

$$\begin{eqnarray}{Q}_{{\rm{c}}} & = & \displaystyle \sum _{lm}\displaystyle \sum _{j\ne i}[{T}_{l}^{i}{\mu }_{l+1,m-1}({a}_{l,m}^{i(e)}{a}_{l+1,m-1}^{i(j)}+{a}_{l+1,m-1}^{i(e)}{a}_{l,m}^{i(j)* })\\ & & +{T}_{l}^{i* }{\mu }_{l+1,-m-1}({a}_{l,m}^{i(e)}{a}_{l+1,m+1}^{i(j)* }+{a}_{l+1,m+1}^{i(e)}{a}_{l,m}^{i(j)})].\end{eqnarray} \tag{ 12c }$$

The first term ${Q}_{{\rm{e}}}$ involves only the external field ${a}_{l,m}^{i(e)},$ which is always zero because ${a}_{l,m}^{i(e)}{a}_{l+1,m\pm 1}^{i(e)}=0$ (see equation (3)). This stems inherently from the translational invariance of the external PSW field in the xy plane. The second term ${Q}_{{\rm{s}}}$ involving only ${a}_{l,m}^{i(j)}$ is contributed by the scattering from the other particles $j\ne i,$ and the third term ${Q}_{{\rm{c}}}$ is a cross term that involves the external incidence ${a}_{l,m}^{i(e)}$ and the scattering ${a}_{l,m}^{i(j)}$ simultaneously. It is straight to prove that ${a}_{l,m}^{i(j)}$ becomes negligibly small with respect to ${a}_{l,m}^{i(e)}$ as the inter-particle separation increases. Therefore, the cross term ${Q}_{{\rm{c}}}$ is dominant for the dilute particle suspension.

Using the approximation ${h}_{l}(x)=\tfrac{{(-i)}^{l}{{\rm{e}}}^{{\rm{i}}x}}{{\rm{i}}x}+\tfrac{(l+1)l}{2}\tfrac{{(-i)}^{l}{{\rm{e}}}^{{\rm{i}}x}}{{x}^{2}}+O\left(\tfrac{1}{{x}^{2}}\right)$ for $x\gg 1,$ ${a}_{l,m}^{i(j)}$ in equation (6) can be reshaped as

$$\begin{eqnarray}&&{a}_{l,m}^{i(j)}={S}_{l,m}^{ij(1)}\displaystyle \frac{{{\rm{e}}}^{{\rm{i}}k{R}_{ij}}}{k{R}_{ij}}+{S}_{l,m}^{ij(2)}\displaystyle \frac{{{\rm{e}}}^{{\rm{i}}k{R}_{ij}}}{{(k{R}_{ij})}^{2}}+O\left(\displaystyle \frac{1}{{(k{R}_{ij})}^{2}}\right),\end{eqnarray} \tag{ 13 }$$

where the coefficients

$$\begin{eqnarray}&&{S}_{l,m}^{ij(1)}=\displaystyle \sum _{l^{\prime} }\left[{t}_{l^{\prime} }^{j}{a}_{l^{\prime} ,0}^{j(e)}4\pi {{\rm{i}}}^{l-l^{\prime} -1}\displaystyle \sum _{l^{\prime\prime} }{C}_{lm{l}^{\prime\prime }(-m)}^{l^{\prime} 0}{Y}_{l^{\prime\prime} ,-m}\left(\displaystyle \frac{\pi }{2},{\phi }_{ij}\right)\right],\end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray}&&{S}_{l,m}^{ij(2)}=\displaystyle \sum _{l^{\prime} }\left[{t}_{l^{\prime} }^{j}{a}_{l^{\prime} ,0}^{j(e)}4\pi {{\rm{i}}}^{l-l^{\prime} }\displaystyle \sum _{l^{\prime\prime} }\displaystyle \frac{(l^{\prime\prime} +1)l^{\prime\prime} }{2}{C}_{lml^{\prime\prime} (-m)}^{l^{\prime} 0}{Y}_{l^{\prime\prime} ,-m}\left(\displaystyle \frac{\pi }{2},{\phi }_{ij}\right)\right].\end{eqnarray} \tag{ 14b }$$

with ${\phi }_{ij}$ denoting the azimuthal angle of ${\hat{R}}_{ij}.$ Especially, one may prove that ${S}_{l,m}^{ij(1)}=0$ if the particles are positioned on the node plane of the PSW field, since ${a}_{l^{\prime} ,m^{\prime} }^{j(e)}=0$ for even $l^{\prime}$ (see equation (3)) and $\displaystyle {\sum }_{l^{\prime\prime} }{C}_{lml^{\prime\prime} (-m)}^{l^{\prime} 0}{Y}_{l^{\prime\prime} ,-m}\left(\tfrac{\pi }{2},\phi \right)=0$ for odd $l^{\prime} .$ Therefore, the total interaction force exerting on the sphere $i$ can be approximately written as

$$\begin{eqnarray}\begin{array}{rcl}{F}_{x}^{{i}}+{\rm{i}}{F}_{y}^{{i}} & \approx & \displaystyle \sum _{j\ne {i}}\displaystyle \sum _{l}\displaystyle \frac{{\rm{i}}{\rho }_{0}}{4}\displaystyle \frac{1}{{(k{R}_{{i}j})}^{n}}[{\mu }_{l+1,-1}{a}_{l,0}^{{i}({e})}({T}_{l}^{{i}}{{\rm{e}}}^{{\rm{i}}k{R}_{{i}j}}{S}_{l+1,-1}^{{i}j(n)}+{T}_{l}^{{i}* }{{\rm{e}}}^{-{\rm{i}}k{R}_{{i}j}}{S}_{l+\mathrm{1,1}}^{{i}j(n)* })\\ & & +{\mu }_{l+\mathrm{1,0}}{a}_{l+\mathrm{1,0}}^{{i}({e})}({T}_{l}^{{i}}{{\rm{e}}}^{-{\rm{i}}k{R}_{{i}j}}{S}_{l,1}^{{i}j(n)* }+{T}_{l}^{{\rm{i}}* }{{\rm{e}}}^{{\rm{i}}k{R}_{{i}j}}{S}_{l,-1}^{{i}j(n)})],\end{array}\end{eqnarray} \tag{ 15 }$$

where n = 1 and n = 2 correspond to the antinode and node planes, respectively.

Note that here only the terms ${S}_{l,1}^{ij(n)}$ and $\,\,{S}_{l,-1}^{ij(n)}$ are involved, because ${a}_{l,m}^{i(e)}$ is nonzero only for m = 0 (see equation (3)).

To further simplify equation (15), we substitute the definition of the spherical harmonics ${Y}_{l,m}(\theta ,\phi )={\zeta }_{l}{P}_{l}^{m}({\rm{\cos }}\,\theta ){{\rm{e}}}^{{\rm{i}}m\phi }$ into equation (14), with $\{{P}_{l}^{m}\}$ being the associated Legendre polynomials and ${\zeta }_{l}=\sqrt{[(2l+1)(l+1)!]/[4\pi (l-1)!]}.$ This gives rise to

$$\begin{eqnarray}&&{S}_{l,1}^{ij(n)}={U}_{l}^{j(n)}{{\rm{e}}}^{-{\rm{i}}{\phi }_{ij}},\,{S}_{l,-1}^{ij(n)}=-{U}_{l}^{j(n)}{{\rm{e}}}^{{\rm{i}}{\phi }_{ij}}.\end{eqnarray} \tag{ 16 }$$

with

$$\begin{eqnarray}&&{U}_{l}^{j(n)}=\displaystyle \sum _{l^{\prime} }\left[{t}_{l^{\prime} }^{j}{a}_{l^{\prime} ,0}^{j(e)}4\pi {{\rm{i}}}^{l-l^{\prime} -{\delta }_{n1}}\displaystyle \sum _{l^{\prime\prime} }\left.{\left[\displaystyle \frac{(l^{\prime\prime} +1)l^{\prime\prime} }{2}\right]}^{{\delta }_{n2}}{C}_{l1l^{\prime\prime} (-1)}^{l^{\prime} 0}{\zeta }_{l^{\prime\prime} }{P}_{l^{\prime\prime} }^{-1}(0)\right]\right..\end{eqnarray} \tag{ 17 }$$

Combining equation (15) with equation (16), we derive a relatively simple and compact form

$$\begin{eqnarray}&&{F}_{x}^{{i}}+{\rm{i}}{F}_{y}^{{i}}=\displaystyle \sum _{j\ne {i}}\displaystyle \frac{{A}_{{i}j}\,{\rm{\cos }}(k{R}_{{i}j}+{\varphi }_{{i}j})}{{(k{R}_{{i}j})}^{n}}{{\rm{e}}}^{{\rm{i}}{\phi }_{{i}j}}.\end{eqnarray} \tag{ 18 }$$

The amplitude factor ${A}_{ij}=| {Z}^{ij(n)}|$ and the phase shift ${\varphi }_{ij}={\rm{\arg }}({Z}^{ij(n)}),$ where the complex quantity ${Z}^{ij(n)}$ is defined by

$$\begin{eqnarray}&&{Z}^{ij(n)}=-{\rm{i}}\displaystyle \frac{{\rho }_{0}}{2}\displaystyle \sum _{l}({\mu }_{l+1,-1}{a}_{l,0}^{i(e)}{U}_{l+1}^{j(n)}{T}_{l}^{i}+{\mu }_{l+\mathrm{1,0}}{a}_{l+\mathrm{1,0}}^{i(e)}{U}_{l}^{j(n)}{T}_{l}^{i* }).\end{eqnarray} \tag{ 19 }$$

Interestingly, for Rayleigh particles ${Z}^{ij(n)}$ can be simplified further by considering the lowest two scattering channels only. This gives rise to ${A}_{ij}=\tfrac{\pi }{72}{{\rm{\Psi }}}_{0}^{2}{\rho }_{0}{(k{D}_{i})}^{3}{(k{D}_{j})}^{3}{f}_{i,0}{f}_{j,0}$ and ${\varphi }_{ij}=\pi /2$ for the antinode plane, and ${A}_{ij}=\tfrac{\pi }{32}{{\rm{\Psi }}}_{0}^{2}{\rho }_{0}{(k{D}_{i})}^{3}{(k{D}_{j})}^{3}{f}_{i,1}{f}_{j,1}$ and ${\varphi }_{ij}=0$ for the node plane, where ${f}_{i,0}=({\kappa }_{i}-{\kappa }_{0})/{\kappa }_{i}$ and ${f}_{i,1}=2({\rho }_{i}-{\rho }_{0})/(2{\rho }_{i}+{\rho }_{0})$ are the monopole and dipole scattering factors of the Rayleigh particle $i,$ respectively. Here ${\rho }_{i},$ ${\kappa }_{i}$ and ${D}_{i}$ are the mass density, the bulk modulus and the diameter of the particle $i,$ respectively, and ${\kappa }_{0}$ is the bulk modulus of the background fluid. The result responsible for the Rayleigh limit is consistent with that derived by Silva and Bruus (see equations (22b) and (23b) in reference [24]).

Below we give a summary on the property of the acoustic interaction among the sparsely distributed multiple particles. (i) The total force (see equation (18)) exerting on the sphere $i$ can be viewed as a vector addition of all pairwise interactions labeled by $ij,$ since the quantity ${Z}^{ij(n)}$ in equation (19) depends on the scattering property of the particles $i$ and $j$ only. The additivity of the interaction does not hold naturally since the ARF is a second order quantity of the sound field. In this case, the conclusion stems from the negligible contribution of the interference among the scattering fields (see equation (12b)). As a consequence, a multi-body system can be simply decomposed into a series of independent two-body problems, which will greatly simplify the computation of a dilute particle ensemble involving a large number of particles. (ii) For each pairwise interaction, the orientation of the force component depends on the azimuthal angle of ${\hat{R}}_{ij},$ i.e., ${\phi }_{ij}.$ This means that the pairwise interaction is a 2D central force and thus is conservative, as long as the particle pair is tightly bounded on the node or antinode plane. Therefore, we can use a potential energy function to characterize the multi-particle system, considering the additivity of the pairwise force. (iii) As a whole, the pairwise interaction decays at a power exponential factor ${(k{R}_{ij})}^{-n},$ where $n=1\,$ or $2$ corresponds to the case of the antinode plane or the node plane. Similar to the Coulomb force or the gravitational force, the acoustically mediated interaction decays moderately and is a long-range force. Besides, the cosine-like oscillation factor ${\rm{\cos }}(k{R}_{ij}+{\varphi }_{ij})$ suggests a binding effect between the two particles: there exists a series of stable configurations occurring at $k{R}_{ij}+{\varphi }_{ij}=(4N+1)\pi /2$ with $N$ denoting an integer. The periodicity exhibited in the inter-particle separations of the bound states originate inherently from the distance dependent factor ${a}_{l,m}^{i(j)}$ (see equation (13)) involved in the dominant force component ${Q}_{{\rm{c}}}$ (see equation (12c)). Although a similar formula has been reported in optic systems [3, 5, 9], which involves the dipole–dipole interaction only, the present form of the pairwise acoustic interaction is not obvious since all scattering channels are considered here.

The sound-mediated interaction is closely related to the scattering property of the particles. In the following paper, the particles are assumed to be identical despite the fact that equation (18) can deal with an ensemble of particles with different geometry and material parameters. Intentionally, the water-immersed polystyrene particle is under consideration. The system carries rich scattering resonances at a wavelength comparable with the diameter $D,$ as revealed by figure 2, the size-dependent magnitudes of the scattering matrix elements $| {t}_{l}| .$ The material parameters involved are: the mass density $\rho =1050\,{\rm{kg}}\,{{\rm{m}}}^{-3},$ the longitudinal velocity ${v}_{{\rm{l}}}=2400\,{\rm{m}}\,{{\rm{s}}}^{-1},$ and the transverse velocity ${v}_{{\rm{t}}}=1150\,{\rm{m}}\,{\rm{s}}{}^{-1}$ for polystyrene; the mass density ${\rho }_{0}=1000\,{\rm{kg}}\,{{\rm{m}}}^{-3}$ and the sound speed ${c}_{0}=1490\,{\rm{m}}\,{{\rm{s}}}^{-1}$ for water. For the convenience of presentation, below all lengths are scaled by the operation wavelength $\lambda ,$ and the interaction force is scaled by ${F}_{0}={E}_{0}{S}_{0},$ where ${E}_{0}={\rho }_{0}{k}^{2}{{\rm{\Psi }}}_{0}^{2}/2$ is the energy density of a single plane wave, and $\,{S}_{0}=\pi {D}^{2}/4$ is the cross-section area of the spherical particle.

Zoom In Zoom Out Reset image size

3.1. The system involving few particles

We first provide a numerical validation for the formula of pairwise interaction. According to equation (18), the normalized interaction force ${F}_{{\rm{N}}}$ between a pair of particles $i$ and $j$ with a large separation $d$ can be expressed as

$$\begin{eqnarray}&&{F}_{{\rm{N}}}(d/\lambda )={A}_{{\rm{N}}}{(d/\lambda )}^{-n}\,{\rm{\cos }}(2\pi d/\lambda +\varphi ).\end{eqnarray} \tag{ 20 }$$

where the amplitude factor ${A}_{{\rm{N}}}=| {Z}^{ij(n)}| /[{(2\pi )}^{n}{F}_{0}],$ the phase shift $\varphi ={\rm{\arg }}({Z}^{ij(n)}),$ and the power attenuation factor $n=1$ or $2$ corresponds to the antinode or node plane. The positive (negative) sign of ${F}_{{\rm{N}}}$ corresponds to a repulsive (attractive) force. Figure 3 shows the pairwise interaction force plotted as a function of the dimensionless inter-particle separation $d/\lambda ,$ where the red dashed and black solid lines provide the comparative results obtained through the theoretical prediction and the rigorous MST approach, respectively. Three different particle sizes are considered, one close to the Rayleigh limit and the other two comparable with the acoustic wavelength. Depending on the stability of a single particle in the $z$ direction, the pair of particles are placed on the same node or antinode plane. It is observed that, the accuracy of equation (20) depends on the particle size and the inter-particle separation simultaneously. A remarkable deviation occurs at the first several wavelengths, especially for the big particles. As the inter-particle separation grows, the analytical results always capture well the MST data for all particle sizes. Without data presented here, we have also studied the pairwise interactions for the non-Rayleigh particles made of, e.g., steel, olive, and air, which carry strikingly different acoustic properties. Again, excellent accuracy of our analytical approach has been revealed by comparing with the MST results, as long as the inter-particle separation is larger than several wavelengths.

Zoom In Zoom Out Reset image size

Without resorting to strict MST calculations, the equation (20) enables a direct understanding of the long-range pairwise interaction. In addition to the well-defined power decay factor, the acoustic interaction can be described by two independent quantities, i.e., the dimensionless amplitude factor ${A}_{{\rm{N}}}$ and the phase shift $\varphi .$ Figure 4 shows ${A}_{{\rm{N}}}$ and $\varphi$ for a pair of polystyrene particles located on the same node plane (blue solid lines) or the antinode plane (red solid lines) of the PSW field, plotted as a function of the particle size $D/\lambda .$ For small particles $(D/\lambda \ll 1),$ as predicted in the Rayleigh limit (dashed lines), ${A}_{{\rm{N}}}$ grows fast because of the rapidly enhanced scattering, accompanying with a nearly constant $\varphi .$ As $D/\lambda$ increases, ${A}_{{\rm{N}}}$ becomes nonmonotonic and $\varphi$ varies sharply (where the latter suggests that the stable configuration is very sensitive to the particle size near the resonance). These features can be understood directly from the characteristics of the scattering matrix elements (see figure 2) and the incident coefficient ${a}_{lm}^{i(e)}$ (which is zero for either odd or even $l,$ depending on the position of the particle, see equation (3)). For comparison, we have also calculated the pairwise interaction by using the rigorous MST method, and fit the parameters ${A}_{{\rm{N}}}$ and $\varphi$ according to equation (20). The fitted data agree well with the analytical ones again, as shown in figure 4 by the blue and red circles.

Zoom In Zoom Out Reset image size

Now we check the vector additivity of the long-range acoustic interaction among multiple spherical particles. Without losing generality, we consider three identical polystyrene particles arranged into an equidistant linear array or a regular triangle array. As shown by the insets in figure 5, both configurations can be defined by the separation between adjacent particles. We consider the total force exerted on the third particle, which is along the x direction due to the symmetry of the system. It is easy to derived that for the linear and triangular configurations, the dimensionless total forces labeled by ${F}_{{\rm{N}}{\rm{L}}}$ and ${F}_{{\rm{N}}{\rm{T}}}$ can be expressed as

$$\begin{eqnarray}&&{F}_{{\rm{N}}{\rm{L}}}(d/\lambda )={F}_{{\rm{N}}}(d/\lambda )+{F}_{{\rm{N}}}(2d/\lambda ),\end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray}&&{F}_{{\rm{N}}{\rm{T}}}(d/\lambda )=\sqrt{3}{F}_{{\rm{N}}}(d/\lambda ),\end{eqnarray} \tag{ 21b }$$

where ${F}_{{\rm{N}}}(d/\lambda )$ is the pairwise interaction in equation (20). In figure 5 we present the analytical results (red dashed lines) for both configurations, together with the MST (black solid lines) for comparison. The particles with sizes ${D}/\lambda =0.5$ and ${D}/\lambda =0.8$ are considered, respectively. Excellent agreements between the two approaches confirm well the vector additivity of the sound-mediated long-range interaction.

Zoom In Zoom Out Reset image size

3.2. The system involving a large number of particles

Below we consider the systems involving many particles. Figure 6(a) shows 91 particles (diameter $D/\lambda =0.5)$ located on the node plane of the PSW field and arranged in a hexagonal lattice. We have calculated the total ARFs exerting on two of the particles (labeled by 1 and 2), one close to the center and the other near the boundary of the particle array, respectively. For these two particles, only the x-component is nonzero because of the mirror symmetry of the system. Figure 6(b) shows the comparative results for one hundred of configurations with different lattice constants $(a/\lambda ),$ evaluated by our analytical formula (solid lines) and the rigorous MST method (circles). Excellent agreement can be observed. (Note that in both calculations, the cutoff of the angular quantum number ${l}_{{\rm{\max }}}=8\,$ is used.) In the MST approach, ∼28 h is paid to finish the calculation, in which solving the linear problem is time-consuming and requires large memory storage. In contrast, only ∼1.0 s is cost in our analytical method, associated with a negligible memory requirement. It is worth pointing out that, the particle cluster will be eventually driven into a mechanically stable configuration by the nonzero force, as shown by the green circles in figure 6(a), which slightly deviates from the initial configuration of hexagonal lattice.

Zoom In Zoom Out Reset image size

The significant optimization in computation time and memory cost also allows the analytic method to handle a system involving a huge number of particles, which is unattainable by the conventional MST approach. For instance, we consider a system consisting of 101 × 101 particles $(D/\lambda =0.5),$ arranged by a square lattice in the node plane of the PSW field. The lattice constant is taken as $a/\lambda =10.36,$ one of the stable distances in the corresponding two-particle system. Obviously, the particle arrangement does not guarantee the stability of the whole system. For every particle, we have calculated the total interaction force exerted by the other particles. Figure 7(a) shows the force amplitudes for all particles. Interestingly, one can observe a much weaker force for the inner particle, comparing with the case near the boundary. Physically, the screening effect of the acoustic interaction can be understood by the cancelation of many pairwise interactions of relatively short distance. For example, as indicated in figure 7(a), the pairwise interaction $\langle j,i\rangle$ can be balanced by $\langle k,i\rangle$ if the position vector satisfies ${{\bf{R}}}_{ij}+{{\bf{R}}}_{ik}=0.$ A similar screening effect can be seen in figure 7(b) for the case of $D/\lambda =0.8.$ Comparing with figure 7(a), the force attenuation from the boundary is much slower, since the particles are located on the antinode plane now (n = 1).

Zoom In Zoom Out Reset image size

Starting from a single prior scattering approximation, we have derived a compact analytical formula for the sound-mediated interaction among sparsely distributed particles. The formula, confirmed by the rigorous MST method, reveals that a multi-particle system can always be reduced to a two-particle problem, irrelevant to the size of particles. A combination of our analytical method with the MST approach could be helpful to further handle a big particle cluster with arbitrary inter-particle separations. The possibility of such new type of oscillatory long-range interactions in creating thermodynamically stable 2D colloidal crystals is also of great interest⁴ .

This work is supported by the National Basic Research Program of China (Grant No. 2015CB755500); National Natural Science Foundation of China (Grant Nos. 11374233, 11534013, 11674250, and 11547310).