Conservative and nonconservative forces for Mie particles in acoustic trapping

A general acoustic force field can be decomposed into a conservative gradient force (GF) and a non-conservative scattering force (SF), which have very different physical and mathematical properties. However, the profiles of such forces for Mie particles are unknown, let alone their underlying physics. Here, by using a fast Fourier transform approach, we calculated the GF and SF for spherical particle of various sizes and various incident waves. For the same focused incident waves, the normalized GF and SF are similar for different particle sizes, while the total force can be quite different owing to the varying relative strength between the GF and SF. GF and SF possess symmetries that are not found in the incident waves, indicating that these physically and mathematically distinct forces have symmetries that are hidden from the beam profile. For a vortex beam carrying a well-defined topological charge, acoustic forces alone cannot trap particles.


Introduction
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][3][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][6][7][8][9][10][11][12][13][14][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 F(r), according to the Helmholtz's theorem, can always be decomposed into a curl-less conservative force field F G (r) = −∇U(r) fulfilling ∇ × F G = 0, and a divergence-less non-conservative force field F S (r) = ∇ × g(r) fulfilling ∇ • F S = 0, where U(r) and g(r) are scalar and vector potential fields, respectively.We call F G (r) the gradient force (GF) and F 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][26][27][28][29] are significantly fewer than that in optics [23,[30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][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][53][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][61][62][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.

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: where ´F(r)e −iq•r d 3 r is the Fourier transform of the total force field F(r), and q is the coordinate in the Fourier space.In fact, F(r) can be any force, including optical and acoustic force.The Fourier transform in equation (1) requires F(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.

Calculating GF and SF for focused acoustic beams
Throughout this paper, a spherical expanded polystyrene particle with radius r p , density ρ p = 29 kg m −3 and longitudinal sound speed is c p = 900 m 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, θ, φ ) be the spherical coordinate centered at the focus (indicated by the red dot), then the OSS only spans from θ = −π + α 1 to θ = −π + α 0 .Its outer radius is R 0 = 825 mm and inner radius is R 1 = 375 mm.The focal length of the OSS is r 0 = 1050 mm indicated by the green dashed lines.The existence of the opening will reduce the forward radiation pressure, thus enhancing the trapping.
We assume a velocity distribution u(φ ) 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 φ , an additional phase of e im ′ φ is imposed, where m ′ is an integer that determines the topological charge, and thus the azimuthal dependence of the velocity distribution is 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 ′ , standing and non-standing fields consisting of PWs are also taken into account.

Comparing with monopole and dipole approximation
For a 1 mm-radius particle manipulated by a focused beam with m ′ = 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 with Here, p is the pressure field, v is the velocity field, is the compressibility, ρ b , c b are the density and sound speed of the background medium, respectively.k b = ω/c b is the wavenumber in the background medium, and ω is the angular frequency.
3a 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 p of 1 mm at a wavelength of 8.7 mm, where k b r p = 0.72.
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][67][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.

GF and SF for focused acoustic beams
We begin by examining the GF and SF exerted on a polystyrene particle by the focused m ′ = 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 F(r) is based on the highly efficient and accurate Mie theory and translational addition theorem outlined in the appendix.After applying equation (1) to F(r), we can obtain F G (r) and F S (r).
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 ′ = 0 focused acoustic beam, the total force, GF, and SF are respectively given by F(ρ, φ , z) = (F tot ρ , 0, Since m ′ = 0, no force depends on φ and F φ = 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 G ρ and F G z ) 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 S z is not only the dominating component of SF, but also greater than GF near the beam axis, making stable trapping with m ′ = 0 impossible.
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 ′ = 1, where the particle radius is 1 mm and 3 mm.The intensities of m ′ = 1 beam in the xy-plane and xz-plane are shown in figure 5.With a nonzero m ′ , the azimuthal force SF emerges.The particle may be pushed near the focal plane due to the dominated GF, while the increased F 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 ρ , 0, F tot z ) (a), (b) for 1 mm-radius particle and (g), (h) for 3 mm-radius particle, GF (F G ρ , 0, F G z ) (c), (d) for 1 mm-radius particle and (i), (j) for 3 mm-radius particle, and SF (F S ρ , 0, F S z ) (e), (f) for 1 mm-radius particle and (k), (l) for 3 mm-radius particle in the ρz-plane.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 ′ , 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].
When the topological charge increases from m ′ = 1 to m ′ = 3, the low-intensity region near the focus expands, and F φ 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.
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 ′ ) 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].

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

Conclusion
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.and p in = −ik b ρ b c b ϕ in .R is the distance between an arbitrary point P(r 0 , θ 0 , φ 0 ) on the transducer's surface and the target point M(r, θ, φ ) with R 2 = r 2 + r 2 0 − 2rr 0 cos γ.The velocity distribution u on the transducer surface in the configuration of the focused beam is defined by a pupil function: Here, α 0 and α 1 are denoted in figure 1(a).The constant U 0 is the real amplitude of the vibration velocity.The topological charge m ′ denotes the phase variation, while the variation of polar angle is not taken into account by letting f(θ 0 ) = 1.
Meanwhile, the potential expression can be expanded using with cos γ = sin θ sin θ 0 cos(φ − φ 0 ) + cos θ cos θ 0 .The incident pressure field can then be expressed as Therefore, the BSC in pressure expansion is

B2. Plane wave field
A plane wave can be expanded as where p = p 0 e ik•r represents the plane wave pressure field p 0 is the characteristic amplitude of pressure, both the plane wave vector k = (k b , θ k, ϕ k ) and the position vector r = (r, θ, ϕ ) are expressed in a spherical coordinate system.Thus, BSCs for plane waves are b m n = (−1) We remark that to move the coordinate center by d, it is sufficient to add a constant e ik•d to every BSC:

Appendix C: Acoustic force calculation
We calculate the acoustic force by integrating the radiation stress tensor over the scatterer's surface.The time-averaged radiation stress tensor ⟨ ↔ T⟩ is given by where ⟨• • • ⟩ stands for the time-average of the corresponding argument, v = |v| and ↔ I is the identity matrix in three-dimension.The acoustic radiation force can be expressed by where is the characteristic energy density, e x , e y and e z are the Cartesian unit vectors [74], and where

Appendix D: Comparing the acoustic force obtained from Mie theory with other methods
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][67][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.

Appendix E: Verification of the numerical accuracy
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.

Appendix F: Numerical accuracy of the acoustic obtained from Mie theory
We note that the expansion series in the Mie theory are truncated with the empirical formula [76,77] N max = k b r p + 4(k b r p ) 1/3 + 2, where r p and k 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.

Figure 1 .
Figure 1.Schematic of the concave spherical transducer in the form of an open spherical sector (OSS).(a) With an outer radius of R0 = 825 mm, an inner radius of R1 = 375 mm, and a focal length of r0 = 1050 mm (as illustrated by the green dashed lines).(b) Illustration of the focusing effect.The red dots denote the coordinate origin and the focus.

Figure 2 .
Figure 2. The disagreement of the forces calculated by analytical multipole expansion method (Anal., solid lines) and the more accurate numerical FFT method (FFT, dashed lines).(a) Force distribution on the radial axis and z = 0 mm.(b) Force distribution on the beam axis and ρ = 0 mm.

Figure 3 .
Figure 3. Intensity distribution of the focused acoustic beam with m ′ = 0 in the xy-plane (a) and the xz-plane (b).The green crosses mark the coordinate center.a.u.denotes arbitrary unit.

Figure 5 .
Figure 5. Intensity distribution of the focused acoustic beam with m ′ = 1 in the xy-plane (a) and the xz-plane (b).The green crosses mark the coordinate center.a.u.denotes arbitrary unit.

Figure 6 .
Figure 6.GF and SF for a spherical particle with a radius of 1 mm or 3 mm manipulated by a focused acoustic beam with m ′ = 1.Total force field (F tot ρ , F tot φ , F tot z ) (a)-(c) for 1 mm-radius particle and (j)-(l) for 3 mm-radius particle, GF (F G ρ , Fφ G , F G z ) (d)-(f) for 1 mm-radius particle and (m)-(o) for 3 mm-radius particle, and SF (F S ρ , Fφ S , F S z ) (g)-(i) for 1 mm-radius particle and (p)-(r) for 3 mm-radius particle in the ρz-plane.

Figure 7 .
Figure 7. GF and SF for a spherical particle with a radius of 1 mm manipulated by a focused acoustic beam with m ′ = 3. (a)-(b) Intensity distribution of the focused acoustic beam in the xy-plane (a) and the xz-plane (b).Total force field (F tot ρ , F tot φ , F tot z ) (c)-(e), GF (F G ρ , Fφ G , F G z ) (f)-(h), and SF (F S ρ , Fφ S , F S z ) (i)-(k) in the ρz-plane.

Figure 8 .
Figure 8.The acoustic radiation force field exerted on a 5 mm-radius spherical particle positioned in a three-beam acoustic plane wave is investigated.(a) The schematic of the three plane waves is shown.(b) Intensity distribution of the interfered plane waves in the xy-plane.Total force field (F tot x (x, y), F tot y (x, y)) (c), (d), GF (F G x (x, y), F G y (x, y)) (e), (f), and SF (F S x (x, y), F S y (x, y)) (g), (h) in the xy-plane, where (x, y) denote the coordinate of the particle in the xy-plane.

Figure 9 .
Figure 9.The acoustic radiation force field exerted on a 5 mm-radius spherical particle positioned in a four-beam acoustic plane wave is investigated.(a) The schematic of the four plane waves is shown.(b) Intensity distribution of the interfered plane waves in the xy-plane.Total force field (F tot x (x, y), F tot y (x, y)) (c), (d), GF (F G x (x, y), F G y (x, y)) (e), (f), and SF (F S x (x, y), F S y (x, y)) (g), (h) in the xy-plane, where (x, y) denote the coordinate of the particle in the xy-plane.

Figure A10 .
Figure A10.Acoustic force calculated by different methods.A 1 mm-radius particle is manipulated by three coherent plane waves (λ = 8.7 mm), as illustrated in (a).We show the y-direction force obtained by Mie theory (our method, red), dipole/monopole expansion (dipole/monopole, green), the finite element method (COMSOL Multiphysics, blue), and Multiple Scattering package in Julia language (an open-source program, orange).

Figure A12 .
Figure A12.Acoustic force obtained from Mie theory versus the truncation order N.A focused beam with a topological charge m = 0 is considered.Different particle positions (a) (0, 0, 0) mm, (b) (1, 0, 0) mm, and (c) (1, 1, 0) mm are considered.The dashed line denotes the empirical truncation order we adopted, which is given by Nmax = k b rp + 4(k b rp) 1/3 + 2, where rp and k b denote the particle radius and the wavenumber, respectively.