Conservative and nonconservative forces for Mie particles in acoustic trapping

Optical tweezers (OTs), a well-established, non-invasive, contactless manipulation technique, were proposed in 1986 [1]. Subsequently, the technique was applied to manipulate microscopic, or even nanoscopic objects, such as colloids and living biological particles [2–4]. Nevertheless, OTs still have some shortcomings such as the heat they generate on absorbing particles and the minute force they can exert. As an alternative, acoustic tweezers (ATs) can also serve as versatile, non-invasive tools for trapping and manipulating particles [5–15]. And ATs have been implemented using various techniques. For instance, they have been achieved through the utilization of single beam [9], hologram [12], ultrasonic attenuation [16], beating sound waves [17], and other methods [10]. Due to the orders of magnitude slower wave speed, a sound wave is generally able to exert much larger forces than its optical counterpart at the same intensity. Also, a particle that absorbs light may not absorb acoustic waves, thus allowing ATs to serve as a complementary tool of OTs. While both techniques are operable for microparticles, nanoparticles are better trapped by OTs, whereas ATs are better suited for millimeter-sized particles [10, 18, 19]. Moreover, acoustic waves typically have better biocompatibility since it is more transparent to living issues than optical waves [19, 20].

The acoustic radiation force arises from momentum transfer during the scattering of sound waves. The system consisting of a particle immersed in an acoustic field is non-conservative and non-Hermitian since sound waves can also exchange energy with the particle owing to its open nature [21, 22]. A general force field ${\mathbf{F}}({\boldsymbol{r}})$, according to the Helmholtz's theorem, can always be decomposed into a curl-less conservative force field ${{\mathbf{F}}^{\text{G}}}(r) = - \nabla U(r)$ fulfilling $\nabla \times {{\mathbf{F}}^{\text{G}}} = {\mathbf{0}}$, and a divergence-less non-conservative force field ${{\mathbf{F}}^{\text{S}}}(r) = \nabla \times {\mathbf{g}}(r)$ fulfilling $\nabla \cdot {{\mathbf{F}}^{\text{S}}} = 0$, where $U({\boldsymbol{r}})$ and ${\mathbf{g}}({\boldsymbol{r}})$ are scalar and vector potential fields, respectively. We call ${{\mathbf{F}}^{\text{G}}}(r)$ the gradient force (GF) and ${{\mathbf{F}}^{\text{S}}}(r)$ the scattering force (SF). GF and SF have very different mathematical and physical properties and applications. GF is a conservative force that is responsible for particle trapping, while SF is a non-conservative force that can transport or revolve particles. Note that SF cannot trap particles by just itself as stated by Earnshaw's theorem [23]. However, writing down the analytical expression of GF and SF for an arbitrary particle is quite difficult. It is possible for Rayleigh particles [24, 25], but not for larger particles. And for deeply subwavelength Rayleigh particles where the SF is vanishing small, the acoustic force field can be considered conservative [17] and is associated with the Gor'kov potential [26, 27]. For larger particles, we also emphasize the crucial role of SF, which can greatly affect the trapping stability. Moreover, researches on GF and SF in acoustics [6, 9, 25–29] are significantly fewer than that in optics [23, 30–49].

There exists plenty of other methods for calculating the total acoustic radiation force, which include the finite element method [50, 51], the boundary element method [52–54], the finite volume method [55, 56], the lattice Boltzmann method [57, 58], and the finite-difference time-domain method [59]. Furthermore, various techniques for computing the acoustic radiation force applied to general structures [60–63] can be found in the literature. Here, we investigate the decomposed GF and SF of acoustic forces exerted on a particle with an arbitrary size and reveal the corresponding physics.

2.1. Computing GF and SF from total force using fast Fourier transform

One way to calculate the GF and SF of an arbitrarily sized spherical particle is the Helmholtz decomposition theorem [64]. However, there are several drawbacks that make this approach very challenging. For instance, the three-dimensional (3D) integration with differential operators is quite difficult to compute even numerically. Additionally, the physics behind is obscured by the math and thus not transparent. Another numerical approach [35, 45] adopts the fast Fourier transform (FFT) to calculate the GF and SF:

$$\begin{align} {{\mathbf{F}}^{\text{G}}}\left({\boldsymbol{r}}\right) & = \int {\frac{{{\boldsymbol{q}}\left( {{\boldsymbol{q}} \cdot {\tilde{\mathbf F}}\left({\boldsymbol{q}}\right)} \right)/{q^2}}}{{{{\left( {2\pi } \right)}^{3/2}}}}{{\text{e}}^{{\text{i}}{\boldsymbol{q}} \cdot {\boldsymbol{r}}}}{{\text{d}}^3}{\boldsymbol{q}}} , \nonumber \\ {{\mathbf{F}}^{\text{S}}}\left({\boldsymbol{r}}\right) & = \int {\frac{{\left( {{\boldsymbol{q}} \times {\tilde{\mathbf F}}\left({\boldsymbol{q}}\right)} \right) \times {\boldsymbol{q}}/{q^2}}}{{{{\left( {2\pi } \right)}^{3/2}}}}{{\text{e}}^{{\text{i}}{\boldsymbol{q}} \cdot {\boldsymbol{r}}}}{{\text{d}}^3}{\boldsymbol{q}}} ,\end{align} \tag{ 1 }$$

where ${\tilde{\mathbf F}}({{\boldsymbol{q}}}) = {\left( {2\pi } \right)^{ - 3/2}}\int {{\mathbf{F}}({{\boldsymbol{r}}}){{\text{e}}^{ - {\text{i}}{{\boldsymbol{q}}} \cdot {{\boldsymbol{r}}}}}{{\text{d}}^3}{{\boldsymbol{r}}}}$ is the Fourier transform of the total force field ${\mathbf{F}}({\boldsymbol{r}})$, and ${\boldsymbol{q}}$ is the coordinate in the Fourier space. In fact, ${\mathbf{F}}({\boldsymbol{r}})$ can be any force, including optical and acoustic force. The Fourier transform in equation (1) requires ${\mathbf{F}}({\boldsymbol{r}})$ to be either a rapidly decaying field or a periodic field. This (total) force field is inputted to our FFT approach, where the conservative and nonconservative parts are separately computed. Nevertheless, FFT automatically assumes that the input force field is periodic. While the FFT-algorithm has no problem with periodic sound field, non-periodic field will suffer from truncation error, and this error is tolerable only when the sound field is decaying sufficiently rapidly.

Here, we present rigorous calculations of acoustic GF and SF for Mie particles using equation (1). We examined both focused acoustic beams, which have a rapidly decaying force field away from the focus, and multiple plane wave (PW) fields, which are periodic. We considered the variation of the forces with particle sizes and topological charges of focused beams. The underlying physics and symmetries are also discussed.

2.2. Calculating GF and SF for focused acoustic beams

Throughout this paper, a spherical expanded polystyrene particle with radius ${r_{\text{p}}}$, density ${\rho _\text{p}}{\text{ = 29 kg}}\;{{\text{m}}^{ - {\text{3}}}}$ and longitudinal sound speed is ${c_\text{p}} = {\text{900 m}}\;{{\text{s}}^{ - 1}}$. In OTs, a high numerical aperture objective lens can be used to focus a paraxial beam. Here, we consider ATs in air, using the 40 kHz acoustic wave, corresponding to a wavelength of 8.7 mm in air at 25 °C. The sound wave can be focused using a concave spherical transducer. We employ the open spherical sector (OSS) transducer similar to the experiment in [65], which is illustrated in figure 1(a). The lake blue in figure 1(b) denotes the continuous sound source with a particular phase such that the acoustic beam can be focused on the red dot. The OSS is a part of a hemisphere but with an empty opening at the bottom. Let $(r,\theta ,\varphi )$ be the spherical coordinate centered at the focus (indicated by the red dot), then the OSS only spans from $\theta = - \pi + {\alpha _1}$ to $\theta = - \pi + {\alpha _0}$. Its outer radius is ${R_0} = 825{\text{ mm}}$ and inner radius is ${R_1} = 375{\text{ mm}}$. The focal length of the OSS is ${r_0} = 1050{\text{ mm}}$ indicated by the green dashed lines. The existence of the opening will reduce the forward radiation pressure, thus enhancing the trapping.

Zoom In Zoom Out Reset image size

We assume a velocity distribution $u(\varphi )$ on the transducer surface with the amplitude ${U_0}$ normally incident on the transducer shown in equation (2). To generate a vortex beam, at each point represented by the coordinate $\varphi$, an additional phase of ${{\textrm{e}}^{{\textrm{i}}m^{\prime}\varphi }}$ is imposed, where $m^{\prime}$ is an integer that determines the topological charge, and thus the azimuthal dependence of the velocity distribution is

$$\begin{equation}u\left(\varphi \right) = \left\{ {\begin{array}{*{20}{c}} {{U_0}{{\text{e}}^{{\text{i}}m^{\prime}\varphi }}{\text{,}}}&{{\text{if }}\pi - {\alpha _0} \unicode{x2A7D} \theta \unicode{x2A7D} \pi - {\alpha _1},} \\ 0&{{\text{otherwise}}.} \end{array}} \right.\end{equation} \tag{ 2 }$$

The field represented by equation (2) is focused by the OSS. In our formalism, the focused field is mathematically expanded by the acoustic Mie theory outlined in appendix A and B. The expansion series in the Mie theory is sufficient for accurately calculating the acoustic force, as discussed in appendix F. The incident pressure field of the focused beam and the PW field can be described by beam-shape coefficients in equations (B7) and (B9), respectively. The maximum of the pressure amplitudes of the focused or multiple PW fields are normalized to 3000 Pa. We calculate the acoustic force of the target particle placed in focused fields with different $m^{\prime}$, standing and non-standing fields consisting of PWs are also taken into account.

3.1. Comparing with monopole and dipole approximation

For a 1 mm-radius particle manipulated by a focused beam with $m^{\prime} = 0$, we compare the FFT results (dashed lines) with the forces calculated by the monopole and dipole approximation (solid lines) [9, 25], where the GF and SF are

$$\begin{align} {\mathbf{F}}_{{\text{Anal}}{\text{.}}}^{\text{G}} & = \text{Re} \left( {{\alpha _{\text{m}}}} \right)\nabla {W^{\left(p\right)}} + \text{Re} \left( {{\alpha _{\text{d}}}} \right)\nabla {W^{\left({\mathbf{v}}\right)}}, \nonumber\\ {\mathbf{F}}_{{\text{Anal}}{\text{.}}}^{\text{S}} & = 2\omega \left[ {\text{Im} \left( {{\alpha _{\text{m}}}} \right){{\mathbf{P}}^{\left(p\right)}} + \text{Im} \left( {{\alpha _{\text{d}}}} \right){{\mathbf{P}}^{\left({\mathbf{v}}\right)}}} \right] - \frac{{k_{\text{b}}^3\omega }}{{6\pi }}\text{Re} \left( {{\alpha _{\text{m}}}} \right)\text{Re} \left( {{\alpha _{\text{d}}}} \right){{\mathbf{P}}^{\left(p\right)}}, \end{align} \tag{ 3 }$$

with

$$\begin{align} W & = \frac{1}{4}\left( {\beta {{\left| p \right|}^2} + {\rho _{\text{b}}}{{\left| {\mathbf{v}} \right|}^2}} \right) \equiv {W^{\left(p\right)}} + {W^{\left({\mathbf{v}}\right)}}, \nonumber \\ {\mathbf{P}} &= \frac{1}{{4\omega }}\text{Im} \left( {\beta {p^*}\nabla p + {\rho _{\text{b}}}{{\mathbf{v}}^*} \cdot \nabla {\mathbf{v}}} \right) \equiv {{\mathbf{P}}^{\left(p\right)}} + {{\mathbf{P}}^{\left({\mathbf{v}}\right)}} .\end{align} \tag{ 4 }$$

Here, $p$ is the pressure field, v is the velocity field, $\beta = \frac{1}{{{\rho _{\text{b}}}c_{\text{b}}^2}}$ is the compressibility, ${\rho _{\text{b}}}$, ${c_{\text{b}}}$ are the density and sound speed of the background medium, respectively. ${k_b} = \omega /{c_\text{b}}$ is the wavenumber in the background medium, and $\omega$ is the angular frequency. ${\alpha _{\text{m}}} = - \frac{{4\pi i}}{{k_\text{b}^3}}{a_0}$ and ${\alpha _{\text{d}}} = - \frac{{4\pi i}}{{k_\text{b}^3}}3{a_1}$ are the polarizabilities of monopole and dipole, respectively, with ${a_0}$ and ${a_1}$ being the first two Mie scattering coefficients. Deviations between the analytical multipole expansion method and the more accurate numerical FFT method are observed in figure 2, highlighting the limitations of using equation (3) to decompose GF and SF, even for particles with a radius ${r_\text{p}}$ of 1 mm at a wavelength of 8.7 mm, where ${k_\text{b}}{r_\text{p}} = 0.72$.

Zoom In Zoom Out Reset image size

Here, we compared the decomposed force fields obtained through our FFT approach with those calculated using the analytical multipole expansion method. The multipole expansion method fails to accurately calculate the force fields in this case, because the truncated higher order multipoles are not negligible, and one must resort to the full Mie theory. To provide evidence for the accuracy of our force calculation method, we present the force fields of a spherical particle with a radius of 1 mm (λ = 8.7 mm) in appendix D (figure A10), where the force values are determined using Mie theory (our method), multipole expansion (taking into account only monopole and dipole contributions), the finite element method (COMSOL Multiphysics), and Multiple Scattering package in Julia language (an open-source program based on T-matrix) [66–68]. The results obtained by the Mie theory, the finite element method, and the Julia language are identical, while they deviate from the multiple expansion, which only considers the monopole and dipole contributions.

3.2. GF and SF for focused acoustic beams

We begin by examining the GF and SF exerted on a polystyrene particle by the focused $m^{\prime} = 0$ acoustic beam. Due to the axial symmetry of the beam, the intensity distribution on the xy-plane (figure 3(a)) and xz-plane (figure 3(b)) can fully represent the focused beam in all space. The field of the beam center marked with green crosses is strong, meeting the rapidly decaying presumption specified by equation (1). As long as a sufficiently large 3D supercell is selected, such that the force acting on the particle decays to zero at the boundary, we can numerically compute GF and SF. The calculation of the acoustic force ${\mathbf{F}}({\boldsymbol{r}})$ is based on the highly efficient and accurate Mie theory and translational addition theorem outlined in the appendix. After applying equation (1) to ${\mathbf{F}}({\boldsymbol{r}})$, we can obtain ${{\mathbf{F}}^{\text{G}}}(r)$ and ${{\mathbf{F}}^{\text{S}}}(r)$.

Zoom In Zoom Out Reset image size

For a spherical particle with a radius of 1 mm (shown in figures 4(a)–(f)) and 3 mm (shown in figures 4(g)–(l)) manipulated by $m^{\prime} = 0$ focused acoustic beam, the total force, GF, and SF are respectively given by ${\mathbf{F}}(\rho ,\varphi ,z) = (F_\rho ^{{\text{tot}}},0,F_z^{{\text{tot}}})$, ${{\mathbf{F}}^{\text{G}}}(\rho ,\varphi ,z) = (F_\rho ^{\text{G}},0,F_z^{\text{G}})$, ${{\mathbf{F}}^{\text{S}}}(\rho ,\varphi ,z) = (F_\rho ^{\text{S}},0,F_z^{\text{S}})$. Since $m^{\prime} = 0$, no force depends on $\varphi$ and ${F_\varphi } = 0$. Because of the axial symmetry, we only display the force in the ρz-plane, while the complete force can be obtained by revolving this plane. The force (both $F_\rho ^{\text{G}}$ and $F_z^{\text{G}}$) repel the particle away from the center, as expected for a solid particle located near the intensity maximum of the focused beam in air. On the other hand, the $F_z^{\text{S}}$ is not only the dominating component of SF, but also greater than GF near the beam axis, making stable trapping with $m^{\prime} = 0$ impossible.

Zoom In Zoom Out Reset image size

We find that when the particle size increases from 1 mm-radius to 3 mm-radius, the SF will increase much more than GF. The profiles of both GF and SF are quite similar for particles with different sizes, indicating that they are mostly determined by the beam profile (intensity for GF and both intensity and phase for SF). However, their total forces are very different, because the relative strength of GF and SF varies with particle size. The symmetries in GF and SF are higher than that in the total force, indicating that these physically and mathematically distinct forces have hidden symmetries that are not revealed in the beam profile [42].

Figure 6 presented the GF and SF for a focused acoustic beam with $m^{\prime} = 1$, where the particle radius is 1 mm and 3 mm. The intensities of $m^{\prime} = 1$ beam in the xy-plane and xz-plane are shown in figure 5. With a nonzero $m^{\prime}$, the azimuthal force SF emerges. The particle may be pushed near the focal plane due to the dominated GF, while the increased $F_\varphi^{{\text{tot}}}$ will rotate the particle around the propagation axis as long as it moves a little away from the beam axis. Meanwhile, along the beam axis, SF surpasses the GF in the longitudinal direction. Nevertheless, trapping along the axial direction may still be possible at an appropriate level of incident power, given that the mass of the particle can balance out the axial SF. We note that an angular momentum-carrying beam, such as those with nonzero $m^{\prime}$, will drive the particle to rotate thus continuously transferring angular momentum and energy to the particle, and eventually allowing the particle to escape. For stable trapping, such angular momentum and energy have to be damped out by ambient damping, and it is not always possible, especially for large particles [65, 69].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When the topological charge increases from $m^{\prime} = 1$ to $m^{\prime} = 3$, the low-intensity region near the focus expands, and $F_\varphi ^{\text{S}}$ also increases accordingly, as shown in figure 7. When the particle is placed on the focal spot, it is almost free of the acoustic force. Nevertheless, it is not a stable equilibrium point, i.e. once the particle leaves a little distance away from the focus, it will escape from that area due to the off-centered transverse force.

Zoom In Zoom Out Reset image size

In acoustic trapping, the stability of the trapped particle is not solely determined by the conservative GF, which is typically designed to trap the particle at the center of the beam. On the other hand, SF also affects the stability of the trapped particle, but it cannot trap any particle by itself. For this reason, SF is commonly regarded as a force that undermines the stability of the trapped particle. In our study, we employ the FFT approach to decompose the total force field into GF and SF, allowing for a comparison between them. By examining the decomposed profiles of the GF and SF, one can gain a clearer understanding of how the trapped particle is manipulated in acoustic trapping. In the context of acoustic trapping using a focused beam, we have examined various topological charges ($m^{\prime}$) and observed their significant impact on both the GF and SF profiles. As the topological charge goes non-zero, the GF turns the focus from a repulsive one to a trapping one. Nonetheless, the vortex-like SF aligned along the z-direction can still lead to the instability in trapping. Such observation aligns with the experiments in [65].

3.3. GF and SF for standing and non-standing wave fields produced by interfering PWs

Consider the interaction between a 5 mm-radius sphere and the acoustic PW fields in the xy-plane for a three-PW and a four-PW configuration. Figure 8(a) (figure 9(a)) shows three (four) equally angular PWs propagating on the xy-plane. The three-PW constitutes a non-standing field while the four-PW forms a standing field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In both kinds of fields, we plotted the acoustic intensity and particle forces in the xy-plane. The magnitude of the GF is a little larger than that of the SF in the non-standing field, and the SF completely vanished in the standing field. This indicates that SF is closely tied to the energy flux of the light, thus it takes a larger value in the non-standing field than in the standing wave case.

The total force for a millimeter-radius particle immersed in a 40 kHz acoustic field has been rigorously calculated and presented (within a square pillar domain). This total force field is then inputted to our FFT-based decomposition algorithm, where the conservative and nonconservative parts are separately computed. Nevertheless, FFT automatically assumes that the input force field is periodic. Therefore, while the FFT-algorithm has no problem with periodic sound field, non-periodic field will suffer from truncation error, and this error is tolerable only when the sound field is decaying sufficiently rapidly. There is a decreasing truncation error associated with increasing sample size. In appendix E (figure A11), we showed the calculated SF at different L, and they deviate only slightly from each other. The errors are relatively strong near the boundaries of the box, owing to the artificial discontinuity of the total force introduced by the FFT-algorithm. The force fields near the center display a higher level of consistency.

For the focused beam configuration we considered, no zero-force position is found along the beam axis owing to the large pushing SF. Nevertheless, stable axial trapping is still possible since the particle mass can balance out the SF under appropriate conditions. GF and SF were found to be insensitive to particle size but depend strongly on the incident fields with different topological charges, as well as standing or non-standing waves.

The GF and SF also show symmetries that are not found in the total force or incident wave field. We conclude that this work may serve as a map for acoustic manipulation. One can select appropriate configurations to optimize or attenuate GF and SF to harmonize with specific applications.

This work was supported by Stable Support Plan Program of Shenzhen Natural Science Fund (Grant No. 20200925152152003), National Natural Science Foundation of China (Grant No. 12074169), and Guangdong Province Talent Recruitment Program (Grant No. 2021QN02C103). Huimin Cheng and Xixi Zhang contributed equally to this work.

The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.

To obtain the acoustic force acting on a spherical particle, the radiation stress tensor is often applied (see appendix C), which requires the incident and scattered fields. The Mie theory offers a semi-analytic series solution to calculate the scattered fields for spherical particles [70–75]. We present a brief introduction to Mie theory.

The incident pressure field can be expanded as:

$$\begin{align}{p_{\text{in}}}\left(r,\theta ,\phi \right) = {p_0}\sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {a_n^m} } {j_n}\left({k_\text{b}}r\right)Y_n^m\left(\theta ,\phi \right),\end{align} \tag{ A1 }$$

where ${p_0}$ is the characteristic pressure amplitude, $a_n^m$ is the beam-shape coefficient (BSC) to be determined, ${j_n}$ is the nth order spherical Bessel function, $Y_n^m$ is the spherical harmonic function of nth order and mth degree. The time dependence ${{\text{e}}^{ - {\text{i}}\omega t}}$ equation (A1) is omitted for the sake of simplicity, as we are calculating the time-averaged acoustic force exerted by a single-frequency wave.

In equation (A1), the incident field is expanded on the transducer's focus. To solve the acoustic scattering field, it is necessary to expand the incident field at center of the target sphere where the coordinates are $O^{\prime}(r^{\prime},\theta ^{\prime},\varphi ^{\prime})$:

$$\begin{align}{p_{{\text{in}}}}\left(r^{\prime},\theta ^{\prime},\phi ^{\prime}\right) = {p_0}\sum\limits_{\nu = 0}^\infty {\sum\limits_{\mu = - \nu }^\nu {a_\nu ^\mu } } {j_\nu }\left({k_\text{b}}r^{\prime}\right)Y_\nu ^\mu \left(\theta ^{\prime},\phi ^{\prime}\right),\end{align} \tag{ A2 }$$

where $a_\nu ^\mu$ is the BSC about the coordinate system of $(r^{\prime},\theta ^{\prime},\varphi ^{\prime})$ [74]. Using the translational addition theorem, the translational BSCs are related to the original BSCs by a separation matrix $S_{\nu n}^{\mu m}({k_\text{b}}{\boldsymbol{d}})$:

$$\begin{align}a_\nu ^\mu = \sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {a_n^m} } S_{\nu n}^{\mu m}\left({k_\text{b}}{\boldsymbol{d}}\right).\end{align} \tag{ A3 }$$

Here ${\boldsymbol{d}}$ is a vector for the displacement of the coordinate origin, ${\boldsymbol{r}} = {\boldsymbol{r}}^{\prime} + {\boldsymbol{d}}$.

The scattered pressure field can be expanded similarly:

$$\begin{align}{p_{\text{sc}}}\left(r^{\prime},\theta ^{\prime},\phi ^{\prime}\right) = {p_0}\sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {s_n^m} } h_n^{\left(1\right)}\left({k_\text{b}}r^{\prime}\right)Y_n^m\left(\theta ^{\prime},\phi ^{\prime}\right),\end{align} \tag{ A4 }$$

where $s_n^m$ is the scattering coefficient, and $h_n^{(1)}$ is the nth order spherical Hankel function of first type. The scattering coefficient is determined by the boundary conditions [25], which states

$$\begin{align}{p_{{\text{in}}}} + {p_{{\text{sc}}}} = {p_{\text{p}}},{\text{ }}\frac{1}{{{\rho _{\text{b}}}}}\left( {\frac{{\partial {p_{{\text{in}}}}}}{{\partial r}} + \frac{{\partial {p_{{\text{sc}}}}}}{{\partial r}}} \right) = \frac{1}{{{\rho _{\text{p}}}}}\frac{{\partial {p_{\text{p}}}}}{{\partial r}},\end{align} \tag{ A5 }$$

where ${p_{\text{p}}}$ is the pressure field inside the particle. After applying the boundary condition, one obtains

$$\begin{equation}s_n^m = {s_n}a_n^m,\end{equation} \tag{ A6 }$$

here ${s_n}$ is the Mie coefficient given by

$$\begin{equation}{s_n} = - \det \left[ {\begin{array}{*{20}{c}} {\gamma {j_n}\left({k_{\text{b}}}{r_{\text{p}}}\right)}&{{j_n}\left({k_{\text{p}}}{r_{\text{p}}}\right)} \\ {{{j_n^{^{\prime}}}}\left({k_{\text{b}}}{r_{\text{p}}}\right)}&{{{j_n^{^{\prime}}}}\left({k_{\text{p}}}{r_{\text{p}}}\right)} \end{array}} \right]\det {\left[ {\begin{array}{*{20}{c}} {\gamma {h_n}^{\left(1\right)}\left({k_{\text{b}}}{r_{\text{p}}}\right)}&{{j_n}\left({k_{\text{p}}}{r_{\text{p}}}\right)} \\[4pt] {{h_n}{{^{\left(1\right)}}^{\prime} }\left({k_{\text{b}}}{r_{\text{p}}}\right)}&{{j_n}^{\prime} \left({k_{\text{p}}}{r_{\text{p}}}\right)} \end{array}} \right]^{ - 1}},\end{equation} \tag{ A7 }$$

where $\gamma = \frac{{{\rho _{\text{b}}}{k_{\text{p}}}}}{{{\rho _{\text{p}}}{k_{\text{b}}}}}$, and ${\rho _{\text{b}}}$ (${\rho _{\text{p}}}$) and ${k_{\text{b}}}$ (${k_{\text{p}}}$) are the density and wavenumber of the background medium (sphere), respectively. The apostrophe indicates the derivative with respect to the argument.

B1. Focused beam field

The pressure amplitude in the focused beam field can be deduced from Rayleigh's integral, which performs the Huygens' principle and reveals that the pressure at a position ${\boldsymbol{r}}$ is the sum of wavelets generated on the surface of the transducer ${S_{\text{T}}}$ [73]:

$$\begin{align}{\phi _{\text{in}}}\left({\boldsymbol{r}}\right) = \frac{1}{{2\pi }}\int_{{S_{\text{T}}}} {u\frac{{{{\text{e}}^{{\text{i}}{k_\text{b}}R}}}}{R}} {\text{d}}{S_0},\end{align} \tag{ B1 }$$

and ${p_{\text{in}}} = - i{k_\text{b}}{\rho _\text{b}}{c_\text{b}}{\phi _{\text{in}}}$. $R$ is the distance between an arbitrary point $P({r_0},{\theta _0},{\varphi _0})$ on the transducer's surface and the target point $M(r,\theta ,\varphi )$ with ${R^2} = {r^2} + r_0^2 - 2r{r_0}\cos \gamma$. The velocity distribution $u$ on the transducer surface in the configuration of the focused beam is defined by a pupil function:

$$\begin{align}u\left({\theta _0},{\varphi _0}\right) = \left\{ {\begin{array}{*{20}{c}} {{U_0}{{\text{e}}^{{\text{i}}m^\prime{\varphi _0}f\left({\theta _0}\right)}},}&{{\text{if }}\pi - {\alpha _0} \unicode{x2A7D} {\theta _0} \unicode{x2A7D} \pi - {\alpha _1},} \\ 0&{{\text{otherwise}}.} \end{array}} \right.\end{align} \tag{ B2 }$$

Here, ${\alpha _0}$ and ${\alpha _1}$ are denoted in figure 1(a). The constant ${U_0}$ is the real amplitude of the vibration velocity. The topological charge $m^{\prime}$ denotes the phase variation, while the variation of polar angle is not taken into account by letting $f({\theta _0}) = 1$.

Meanwhile, the potential expression can be expanded using

$$\begin{align}\frac{{{{\text{e}}^{{\text{i}}{k_\text{b}}R}}}}{{i{k_\text{b}}R}} = \sum\limits_{n = 0}^\infty {\left(2n + 1\right)h_n^{\left(1\right)}} \left({k_\text{b}}{r_0}\right){j_n}\left({k_\text{b}}r\right){P_n}\left(\cos \gamma \right),\end{align} \tag{ B3 }$$

with $\cos \gamma = \sin \theta \sin {\theta _0}\cos (\varphi - {\varphi _0}) + \cos \theta \cos {\theta _0}$. The incident pressure field can then be expressed as

$$\begin{align}{p_{\text{in}}}\left(r,\theta ,\phi \right) & = {p_0}{\left({k_\text{b}}{r_0}\right)^2}\sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {{\delta _{m,m^{\prime}}}{{\left( - 1\right)}^{n + m}}h_n^{\left(1\right)}\left({k_\text{b}}{r_0}\right)N_n^mI_n^m{j_n}\left({k_\text{b}}r\right)Y_n^m\left(\theta ,\phi \right)} } ,\end{align} \tag{ B4 }$$

$$\begin{align}N_n^m &= 4\pi \sqrt {\frac{{\left(2n + 1\right)}}{{4\pi }}\frac{{\left(n - m\right)!}}{{\left(n + m\right)!}}} ,\end{align} \tag{ B5 }$$

$$\begin{align}I_n^m &= \int_{\cos {\alpha _0}}^{\cos {\alpha _1}} {P_n^m\left(x\right){\text{d}}x} .\end{align} \tag{ B6 }$$

Therefore, the BSC in pressure expansion is

$$\begin{align}a_n^m = {\delta _{m,m^{\prime}}}{\left( - 1\right)^{n + m}}{\left({k_\text{b}}{r_0}\right)^2}{h_n}^{\left(1\right)}\left({k_\text{b}}{r_0}\right)N_n^mI_n^m.\end{align} \tag{ B7 }$$

B2. Plane wave field

A plane wave can be expanded as

$$\begin{align}{{\text{e}}^{{\text{i}}k \cdot r}} = \sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {{{\left( - 1\right)}^m}} } 4\pi {i^n}{j_n}\left({k_\text{b}}r\right){Y_{n,m}}\left(\theta ,\phi \right){Y_{n, - m}}\left({\theta _k},{\phi _k}\right),\end{align} \tag{ B8 }$$

where $p = {p_0}{{\text{e}}^{{\text{i}}k \cdot r}}$ represents the plane wave pressure field ${p_0}$ is the characteristic amplitude of pressure, both the plane wave vector ${\boldsymbol{k}} = ({k_\text{b}},{\theta _{k,}}{\phi _k})$ and the position vector ${\boldsymbol{r}} = (r,\theta ,\phi )$ are expressed in a spherical coordinate system. Thus, BSCs for plane waves are

$$\begin{equation}b_n^m = {\left( - 1\right)^m}4\pi {i^n}{Y_{n, - m}}\left({\theta _k},{\phi _k}\right).\end{equation} \tag{ B9 }$$

We remark that to move the coordinate center by ${\boldsymbol{d}}$, it is sufficient to add a constant ${{\text{e}}^{{\text{i}}k \cdot d}}$ to every BSC:

$$\begin{equation}b_\nu ^\mu = b_n^m \cdot {{\text{e}}^{{\text{i}}k \cdot d}}.\end{equation} \tag{ B10 }$$

We calculate the acoustic force by integrating the radiation stress tensor over the scatterer's surface. The time-averaged radiation stress tensor $\langle {{\mathbf{\overset{{\smash{\scriptscriptstyle\leftrightarrow}}} {{\mathbf{T}}} }}} \rangle$ is given by

$$\begin{align}\left\langle {{\mathbf{\overset{{\smash{\scriptscriptstyle\leftrightarrow}}} {{\mathbf{T}}} }}} \right\rangle = {\rho _{b}}\left\langle {{\boldsymbol{vv}}} \right\rangle - \left( {\frac{{{\rho _{b}}\left\langle {{v^2}} \right\rangle }}{2} - \frac{{\left\langle {{p^2}} \right\rangle }}{{2{\rho _{b}}c_{b}^2}}} \right){\mathbf{\overset{{\smash{\scriptscriptstyle\leftrightarrow}}} {{\mathbf{I}}} }},\end{align} \tag{ C1 }$$

where $\left\langle \cdots \right\rangle$ stands for the time-average of the corresponding argument, $v = \left| {\boldsymbol{v}} \right|$ and ${\mathbf{\overset{{\smash{\scriptscriptstyle\leftrightarrow}}} {{\mathbf{I}}} }}$ is the identity matrix in three-dimension. The acoustic radiation force can be expressed by

$$\begin{align}{\mathbf{F}} = \pi r_\text{p}^2{E_0}\left({Y_x}{{\boldsymbol{e}}_x} + {Y_y}{{\boldsymbol{e}}_y} + {Y_z}{{\boldsymbol{e}}_z}\right),\end{align} \tag{ C2 }$$

where ${E_0} = \frac{{p_0^2}}{{2{\rho _b}c_b^2}}$ is the characteristic energy density, ${{\boldsymbol{e}}_x}$, ${{\boldsymbol{e}}_y}$ and ${{\boldsymbol{e}}_z}$ are the Cartesian unit vectors [74], and

$$\begin{align}{Y_x} + i{Y_y} = \frac{i}{{2\pi {{\left({k_{\text{b}}}{r_\text{p}}\right)}^2}}}\sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {\sqrt {\frac{{\left(n + m + 1\right)\left(n + m + 2\right)}}{{\left(2n + 1\right)\left(2n + 3\right)}}} } } \left({S_n}a_n^ma_{n + 1}^{m + {1^*}} + S_n^*a_n^{ - {m^*}}a_{n + 1}^{ - m - 1}\right),\end{align} \tag{ C3 }$$

$$\begin{align}{Y_z} = \frac{1}{{\pi {{\left({k_\text{b}}{r_\text{p}}\right)}^2}}}\text{Im} \sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {\sqrt {\frac{{\left(n - m + 1\right)\left(n + m + 1\right)}}{{\left(2n + 1\right)\left(2n + 3\right)}}} } } {S_n}a_n^ma_{n + 1}^{{m^*}},\end{align} \tag{ C4 }$$

where

$$\begin{align}{S_n} = {s_n} + s_{n + 1}^* + 2{s_n}s_{n + 1}^*.\end{align} \tag{ C5 }$$

Figure A10 depicts the y-direction force fields of a spherical particle with a radius of 1 mm (λ = 8.7 mm) as the particle's position varies along the x-axis. In this case, the particle is being manipulated by three coherent plane waves. The force values were determined using different methods: Mie theory (our method, indicated in red), multipole expansion (considering only monopole and dipole contributions, shown in green), the finite element method (COMSOL Multiphysics, shown in blue), and Multiple Scattering package in Julia language (an open-source program, shown in orange) [66–68]. Notably, the results obtained using Mie theory are consistent with those obtained using the finite element method and Julia language, but they differ from the results obtained using multipole expansion, which only accounts for monopole and dipole contributions. Evidently, the multipole expansion method, even with monopole and dipole contributions, has limitations in accurately calculating the 'total force', even for relatively small particles. To achieve accurate calculations of acoustic forces, it is necessary to include higher-order multipoles.

Zoom In Zoom Out Reset image size

In our approach, we utilize the fast Fourier transform (FFT) to decompose the total force field into two components: the conservative gradient force (GF) and the nonconservative scattering force (SF). It is required that the force field must exhibit either a 'rapidly decaying' or a 'periodic' behavior in order to be compatible with the FFT property.

In the case of a periodic force field, as long as the sampling points are sufficient, the resulting errors can be considered negligible. However, when dealing with a rapidly decaying force field, despite its rapid decay, the force field theoretically extends throughout the entire space. Therefore, in order to apply the FFT, it becomes necessary to 'truncate' the rapidly decaying force field using a 'box' domain, where the force field is sufficiently small outside its bounds. In our study, we used a cubic box domain with a volume ${L^3}$ to truncate the force field.

Intuitively, the errors that arise are highly dependent on the size of the truncation box domain. The larger the box domain size (L), the smaller the resulting errors. In figure A11, we calculated and illustrated the decomposed SF using different L. It has been observed that a slight inconsistency exists between each case. However, it is important to note that the errors primarily occur near the boundaries of the box, where the discontinuity happens in the FFT technique. On the other hand, the force fields in the central region display a higher level of consistency across each case.

Zoom In Zoom Out Reset image size

We note that the expansion series in the Mie theory are truncated with the empirical formula [76, 77] ${N_{\max }} = {k_\text{b}}{r_\text{p}} + 4{({k_\text{b}}{r_\text{p}})^{1/3}} + 2$, where ${r_\text{p}}$ and ${k_\text{b}}$ represents the particle radius and wavenumber, respectively. In figure A12, we present the results of our simulations, demonstrating the acoustic forces exerted on a particle with a radius of 3 mm as the number of expansion series, N, increases. The dashed red line indicates the value of ${N_{\max }}$. Notably, it is observed that the force values converge before reaching ${N_{\max }}$, which shows that ${N_{\max }}$ given by the empirical formula is sufficient to calculate the acoustic forces with a very high accuracy.

Zoom In Zoom Out Reset image size