Feasibility of Multiple Micro-Particle Trapping—A Simulation Study

Both optical tweezers and acoustic tweezers have been demonstrated for trapping small particles in diverse biomedical applications. Compared to the optical tweezers, acoustic tweezers have deeper penetration, lower intensity, and are more useful in light opaque media. These advantages enable the potential utility of acoustic tweezers in biological science. Since the first demonstration of acoustic tweezers, various applications have required the trapping of not only one, but more particles simultaneously in both the axial and lateral direction. In this research, a method is proposed to create multiple trapping patterns, to prove the feasibility of trapping micro-particles. It has potential ability to electronically control the location and movement of the particles in real-time. A multiple-focus acoustic field can be generated by controlling the excitation of the transducer elements. The pressure and intensity of the field are obtained by modeling phased array transducer. Moreover, scattering force and gradient force at various positions are also evaluated to analyze their relative components to the effect of the acoustic tweezers. Besides, the axial and lateral radiation force and the trapping trajectory are computed based on ray acoustic approach. The results obtained demonstrate that the acoustic tweezers are capable of multiple trapping in both the axial and lateral directions.


Introduction
Optical tweezers have been applied to biology and physical research after they were first introduced in 1986 [1], and have been applied in manipulating micro/nano particles such as bacteria, cells, molecular motors, and single molecules [2]. Although optical tweezers have excellent precision as demonstrated in many applications [3], they have several limitations. Optical tweezers are largely limited to transparent media. The trapping ability is significantly reduced in the presence of opaque media, which scatters and attenuates the optical ray on its path to the particles. In addition, optical tweezers are associated with poor penetration in the biotic environment, such as skin and other human tissues. Moreover, previous studies [4][5][6] have demonstrated that optical tweezers can cause physiological damage to biological objects due to laser-induced heating.
Previous studies have also demonstrated the feasibility of manipulating particles by acoustic waves [7][8][9][10][11]. Acoustic tweezers can propagate deeper than optical tweezers in opaque biomedical environments, giving them many advantages in biological applications, such as drug delivery. In contrast with highly focused laser or electron beams, ultrasound has been proven to be safe to biomedical tissues and has been used in diagnostic imaging for over 60 years [12,13]. These two major advantages have driven the development of acoustic tweezers.
Many studies have demonstrated the feasibility of acoustic waves for the non-contact and non-invasive manipulation of microscopic particles. The first report of acoustic tweezers was given by Wu [7] in 1991. A stable force potential well generated by two identical focusing transducers (3.5 MHz) was demonstrated to trap latex particles and clusters of frog eggs. In 1995, Hertz [8] described a 3D trapping mechanism for microscopic particles in a standing wave field produced by two collimated ultrasound transducers. Courtney et al. [9] demonstrated that 5 μm polystyrene particles could be trapped and moved in one dimension using a pair of acoustically matched piezoelectric transducers (5.25 MHz) that generate the standing wave field. However, acoustic trappings described above were achieved using acoustic waves produced by multiple transducers; the implementation of a multiple-transducer system is complicated, which involves small trapping spacing between two collimated transducers that restricts its application in one side open configuration.
To address this problem, single-transducer acoustic trapping mechanisms have been proposed. Lee et al. [14] applied a highly focused ultrasonic transducer (100 MHz) to generate radiation forces to trap target particle along axial direction in Mie regime (in which particle size is six times larger than wavelength). Kang and Yeh [11] designed an acoustic tweezer based on a four-element planar transducer with low driving frequency (1 MHz). In his simulation, Laguerre-Gaussian model is generated by the four-element planar transducer for confining small particles in Rayleigh regime (in which particles are much smaller than the wavelength).
Since the first demonstration of acoustic tweezers, various potential applications have required the simultaneous trapping of many particles rather than one particle only, including sorting small particles [15], controlling multi-molecular complexes [16], particle-patterning technology [17] etc.
While many optical tweezers have now been designed with multiple trapping capability [18], to the best of our knowledge, we were the first to describe single-transducer acoustic tweezers that have multiple trapping ability in a three-dimensional setting [19]. In this paper, we further improve the previous system in the area of scattering and gradient force analysis, the effect of acoustic streaming and resistance due to viscosity of medium in trapping performance.
In this paper, we describe a multiple trapping acoustic tweezers with four focal points, each of which can control a single particle. To demonstrate the feasibility of this system, the radiation forces of a particle along and across the beam propagation direction are analyzed in one of the four trapping points. Further, the trapping trajectory along and across the beam propagation directions was determined with the consideration of acoustic streaming effect and the resistance due to the viscosity of water. Scattering force and gradient force components at various positions were also evaluated to analyze their relative contribution to the effect of the multiple trapping acoustical tweezers. We demonstrate that the proposed acoustic tweezer can generate a set of trapping points in different positions.

Methods
The proposed multiple trapping system is based on a multiple-focus field generated by a phased array transducer. Each focal point in this acoustic tweezers can be used to manipulate an individual particle. The schematic of the system is shown as Figure 1. In Section 2.1, we describe the approach to calculate the excitation vector for each phased array element. In Section 2.2, the multiple-focus field patterns are synthesized with continuous excitation from every element in the phased array transducer instead of assumed Gaussian intensity distribution [10]. shows the Gaussian beams generated by a phased array transducer (gray plane). The orange line is the intensity profile of Gaussian beam. The red ball is the objective particle to be trapped by the focal points in the center of beams. Green arrows denote wave propagation direction in Gaussian beam.
The theory of acoustic tweezers based on ray acoustic approach is introduced to calculate radiation force on the axial and lateral direction in Section 2.3. Resistance force and acoustic streaming effect are calculated in Section 2.4. The overall mechanism of the multiple trapping acoustic tweezers is outlined in the flow chart shown in Figure 2.

Determine Excitation for Each Phased Array Element
A 50 MHz two-dimensional (2D) phased array transducer with 25 × 25 elements, which is shown as Figure 3a was used to generate an acoustic field with multiple focal points for multiple trapping. In this section, a method is introduced to calculate the excitation vector of the every array element to generate desired multiple-focus field patterns. A matrix form propagation operator [20] generated by the discretized Rayleigh-Sommerfeld was applied in this work. The acoustic pressure from source point ′ to a target point q in Figure 2a can be expressed as: where ( ) is acoustic pressure at point q, = √−1, and are the density and the speed of sound in the medium, respectively, is sound wavelength, is time dependence term, is acoustic wave frequency, is acoustic wavenumber, q and q′ are, respectively, the target and source points. denotes the surface of source. is excitation of the source point. For the phased array we used to generate multiple trapping energy field, Equation (1) can be rewritten as: where n is the th element surface, is the excitation of the th element. ′ represents source points on the th element. Suppose the pressure at the M focal spots that used to trap target particle is known, Equation (2) can be rewritten as： or in matrix form: is a vector representing pressure at M focal spots in the multiple trapping energy field. Vector = [ , , … , ] denotes a vector including excitation of N elements. H is a matrix denoting acoustic wave forward propagation, its element can be expressed mathematically as: in which ( , ) represents acoustic wave propagation from source space on the th element to the th focal spots. In Equation (4), can be obtained by inverting the propagation matrix using the generalized inverse approach, expressed as: = * ( * ) (6) * denotes the pseudoinverse of , which can be obtained by calculating decomposition value of matrix H using SVD. To increase the array excitation efficiency, an iterative weighting algorithm is applied based on Equation (6): where is an N × N real, positive definite weighting matrix. It could be initialized as an identity matrix first. Then an iterative weighting algorithm [20] is taken to yield . Finally, given the pressure at the M focal spots, vector , which includes excitation vector of every phased array element that makes specified power at the focal points, can be resolved.

Acoustic Field Modeling
The pressure at a point in the field can be obtained by summarizing all the acoustic pressure from each infinitesimal element. Figure 3b illustrates the definition of the coordinate system and geometry. Δw and Δh indicate the width and height of the phased array element. On, the center of the nth array element, is signified by (xn, yn) in coordinate system (x, y, z). A secondary coordinate system denoted by (x0, y0) has the center of the th array element as its origin. Every phased array element is divided into infinitesimal elements, each of which acts as a small source. dx0dy0 is the surface of the infinitesimal element. For a uniformly moving radiation source with time dependence, the acoustic pressure [21] is given by Equations (8): where represents the distance from the center of infinitesimal element to the observation point in the field, and are the density and the speed of sound in the medium, respectively. The vector p(x, y, z, t) indicates the pressure at observation point (x, y, z) at time t. The excitation of the nth array element can be obtained by Equation (7). Based on Equation (8), the intensity value I(x, y, z, t) for each point at time t in the field is expressed as [22]: where and * are, acoustic pressure in Equation (8) and its conjugate. and are the density and the speed of sound in the medium, respectively.

The Theory of Acoustic Tweezers
Since acoustic waves possess similar physical properties as optical waves, the theory of calculation of the radiation force for acoustic tweezers is analogous to optical tweezers in the Mie regime, where particle diameters are at least six times larger than the wavelength. Figure 4a shows an overall view of a single Gaussian beam acoustic tweezers and a target particle (red ball) with radius suspending in a liquid. The wave propagation direction (green arrows) convergent in front of the trapping point and divergent beyond it. Pink arrow is an individual incident ray hit on the target particle. Red ball is the target particle; (b) 2D zooming geometry for the analysis of the radiation force when a sphere located at the arbitrary location. Blue line represents the acoustic beam. The incident ray assumed from point with angle α to the beam axis z hits the spherical particle (red line) at point Q. The absolute distance between and is , which is also the radius of the spherical wave front. L is the axial distance between and .
and z axis has angle to each other. and indicate the incident angle and refractive angle at the interacting point respectively. ̂ and are the unit vectors representing the directions of the scattering and the gradient forces. Beam waist size is . Figure 4b shows the 2D geometry for the analysis of the radiation force. In a Gaussian distributed beam, an individual incident ray with power strikes a small region on the surface of the particle at point Q. According to the law of conservation, the radiation force caused by an individual ray can be decomposed into two orthogonal components known as scattering force component dFs , which mainly arise from reflection momentum transfer, and gradient force component dFg, which mainly arise from refraction momentum transfer. Commonly they can be written as [10,23]: where and c represent the average power of the incident ray and the acoustic speed in the medium, respectively. and q are dimensionless fractions of the peak scattering force and gradient force [24] transferred to the sphere by the emergent ray, respectively. and are the incident and transmitted angles. and are the Fresnel reflection coefficient and transmission coefficient at the surface of the sphere, they are given by: where and are the acoustic impedances in the surrounding medium and the particles.Then the radiation force along z direction and along y direction on Figure 4a can be estimated by decomposing dFs and dFg, they can be written as Equations (14) and (15): = + (15) in which, and indicate the z component of the scattering force and gradient force , likewise, and represent the y component of the scattering force and gradient force , respectively. As integrating and from the differential area over which the incident ray hits the particle, it leads to the total force and produced by the entire beam. Finally, in the Mie regime where particle size is larger than the wavelength of incident ray, the complete formulations of radiation force in both axial and lateral directions can be described as follows: where is the speed of sound in the medium; is the angle between the beam axis and line .
indicates the incident angle at the interacting point Q. is the radius of the sphere. Term represents an individual ray with power entering the particle across a differential area. ̂ and are the unit vectors representing the directions of the scattering and the gradient forces shown as Figure 4b.
Here, they are decomposed into the y and z directions and denoted as , = − (20) where is the wavelength of the incident ray; and are coordinates of point. represents the vector from to . | | is the absolute distance from to . Other signs in Equations (18)

Resistance Force and Acoustic Streaming Effect
In the numerical analysis of the object motion, the resistance due to viscosity of medium is also taken into consideration in the form of Stokes' formula as Equation (21) [25]: where μ is the viscosity coefficient of medium; and v are the radius and velocity of the spherical particle, respectively. Drag force caused by acoustic streaming is also taken into account, since experiments [7] has found that the main disruption in acoustic trapping is the competition force due to acoustic streaming. To investigate the acoustic streaming effect on multiple-trapping acoustic tweezers, the drag force due to acoustic streaming is taken into account by [26]: where is the dynamic viscosity of medium, is the beam width. is the axial streaming velocity, which can be estimated by the approach proposed by Nowicki [27].

Results and Discussion
The theoretical method to create acoustic field and radiation force for the multiple trapping has been introduced in Section 2. In this section, the results will be shown to demonstrate the feasibility of this multiple trapping acoustic tweezers model. A 2D phased array transducer consisted of 625 elements with kerf width 5 μm, and pitch size 27 μm is used to generate and control the multiple-focus acoustic field as shown in Figure 3a. The simulation parameters including the detailed characteristics of the phased array and material parameters are summarized in Table 1. Here a lipid or fat sphere was chosen as the object, supposing its homogeneous refractive index is higher than that of the surrounding medium.

Multiple-Focus Acoustic Field
First, four focal points are created in the ultrasonic field using the pseudoinverse method to achieve the desired field patterns and intensity levels at the control points. All four control points are set as the same weight = [1, 1, 1, 1] t in Equation (7) at the same depth 670 μm in the plane perpendicular with the Z direction in the XYZ coordinate system, as shown in Figure 1. By detecting the four focal points in the 3D energy field, we can get the coordinates of the four control points, which are (0.18, −0.18, 0.67) mm, (−0.18, 0.18, 0.67) mm, (−0.18, −0.18, 0.67) mm, and (0.18, 0.18, 0.67) mm, respectively. Then the pseudoinverse approach with weighting matrix described by Equation (7) is used to determine the excitation vector to synthesize the desired multiple-focus acoustic field with the maximum array excitation efficiency. After that, the fields with four focal points are obtained by combining the excitation vector with field modeling method expressed by Equations (8) and (9). The results are shown in Figure 5. Figure 5a is the diagonal section of the instantaneous pressure field, which is through points (−0.18, −0.18, 0.67) mm and (0.18, 0.18, 0.67) mm along z axis. Within the instantaneous pressure field, the wave front converges to two focal points and diverges afterwards. So it satisfies the geometry and behavior of a Gaussian beam. According to Equation (9), the diagonal plane of ultrasonic intensity field with desired intensity level at control points is obtained, which is shown in Figure 5b, and used as a parameter to calculate the radiation force. Figure 5c shows the section of acoustic intensity field across the z axis. It can be seen that the field intensity around the four focal points is symmetrical, also the intensity at the four focal points are equal to each other. Figure 5d presents the intensity profile across the Gaussian wave propagation direction, showing the peak at the waist. It can be seen that the spatial peak intensity in the beam is 150 W/cm 2 , the beam width is defined as the spatial extent between the half power points, which is 22 μm here. Hence, the four focal points can be regarded as four beams with 22 μm beam waist.

Radiation Force Analysis for the Particles in Different Location
The acoustic intensity field with four multiple-focus points ( Figure 5) has been created by a phased array transducer with parameters as mentioned above. It can be seen that the four focal points and the field around them are symmetry. So the radiation force in one of the beam is the same with others. Here, the directions along and across the beam propagation direction are defined as axial direction (z-axis), and lateral direction (y-axis) respectively. And the focal point o is set as origin of coordinates yoz, shown in Figure 5b. The axial and lateral radiation force are calculated in one of the four focal points following the steps indicated by the flowchart in Figure 2.
In coordinate system zoy that shown in Figure 5b, the axial and lateral radiation force for a sphere with 240 μm diameter at various locations are shown in Figure 6. The axial radiation force, when the sphere moves along z-axis, is shown in Figure 6a. It can be seen that positive axial radiation force was first exerted to push the target particle forward, afterwards, negative axial radiation force is applied to draw particle back. Therefore, positive and negative force exists alternatively, it makes the particle be trapped. Assuming a sphere moves from (−0.2, 0.02) mm to (1.2, 0.02) mm, parallels to z axis, but 20 μm distance off it, the axial radiation force is plotted in Figure 6b. It is interesting to find that even its motion is off the z-axis, it still produces negative radiation force for acoustic trapping, though the magnitude is smaller than that on the beam axis. Figure 6c shows the axial radiation force when the particle moves along the parallel line with farther distance 50 μm off the z-axis. It can be seen that, although the gradient force component has similar performance with Figure 6a, the scattering force with positive value counteracts the negative force produced by intensity gradient, only the positive axial radiation force can be found in Figure 6c, which results in the failure of acoustic trapping. From Figure 6a-c, it can be deduced that the negative radiation forces along the z axis get smaller and smaller as increased distance off the z-axis. Also, the scattering force mainly contributes to the positive axial force which pushes forward, whereas gradient force gives majority negative force which attracts the objective particle back.   Figure 6d shows the lateral radiation force as the particle moves along the line that parallels to y-axis and passes through point (0.13, 0) mm in the zoy coordinate system. The alternative positive and negative lateral forces indicate that once the particle locates off the beam axis, the lateral force will tend to pull it back to the z-axis. Also, it can be seen that the dash-dot line nearly overlaps with the solid line, while the dashed line is close to the zero. Since the scattering force has the same direction with the incident ray, the major part of the scattering force is in axial direction, while the part of the scattering force for lateral force is very small. Figure 6e shows lateral radiation force as the particle moves along the line which parallels to y-axis and passes through point (0.23, 0) mm. Compared with Figure 6d, the results in Figure 6e show that the lateral radiation force try to keep the particle at z-axis similarly, however, as farther from the y-axis, the magnitude of the lateral force gets smaller. From the gradient and scattering forces component shown in Figure 6a-e, it can be deduced that the gradient force caused by refraction plays a much more critical role than scattering force produced by reflection in producing lateral trapping force

Radiation Force and Trajectory Analysis for Particles with Various Sizes
Particles with various diameters are also used to evaluate the trapping performance. Meantime, according to Newton's second law of motion, the displacement trajectories of particles are determined after motion analysis in consideration of radiation force and the resistance of water viscosity. To simplify the analysis, the radiation force and displacement trajectory by axial and lateral direction are calculated separately. Figure 7a shows the axial radiation force of particles with varied diameter from 180 μm to 360 μm traveling along the z-axis in the coordinate system zoy. It is intriguing to note that bigger trapping force and longer attractive radiation force region are more likely to occur when using bigger particles.
To investigate the acoustic streaming effect on multiple-trapping acoustic tweezers, we introduce the approach studied by Nowicki [25] to estimate the axial streaming velocities. The streaming velocity along the axial direction is shown in Figure 7b. The result of drag force due to acoustic streaming velocity combined with radiation force along the axial direction is shown in Figure 7c. It can be noticed that the drag force is positive when the radiation force is negative in the trapping region, it means drag force intends to prevent drawing back target particle towards trapping point. The integrated force consisted with drag force due to acoustic streaming and radiation force along the axial direction is shown in Figure 7d. We can see even consider the acoustic streaming effect the integrated force along the axial direction is still negative in the trapping region, which means the particles can be trapped in this region. Figure 7e shows the plot of displacement trajectories along z-axis, which is calculated with the consideration of axial radiation force and the resistance of water viscosity. The initial speed is zero and the initial position is in the end of negative force points in Figure 7d. It is noted that although the initial positions are different, the particles are generally pulled back and finally captured within four seconds at their balanced positions. Assuming the particles are located on the lines passing through balanced positions in z-axis and paralleling to y-axis, and their locations varying along this line, the lateral force along this line could be calculated and is plotted in Figure 8a.
It can be noted that although the asymmetric lateral radiation force at the two sides of z-axis, the lateral radiation forces have opposite directions at two sides of z-axis, which means that once the particle is off z-axis, lateral radiation force will tend to trap it back.
Since there is almost zero radial acoustic streaming velocity at the ultrasonic transducer center, the drag force due to radial acoustic streaming effect is omitted when calculate the trajectory as particle move across the axial direction. Figure 8b shows the particles displacement trajectory along this line, which are calculated with the consideration of lateral radiation force and the resistance of water viscosity. The initial lateral positions are 150 μm off the beam axis and the initial speed are zero. It reveals that all the three particles can be trapped near the z-axis.  Given the specific particle size, the maximum axial trapping force, its axial trapping position, axial trapping regions and maximum lateral trapping force are listed in Table 2. It can be seen that the axial and lateral trapping force gets bigger when the size of particle increases. The magnitude of lateral force is much bigger than axial force, which indicates that the lateral trapping is easier to be achieved than the axial trapping. It also shows that the axial trapping range increases as the increased size of the particles.

Conclusions
A multiple trapping acoustic tweezers model is proposed in this paper. The pseudoinverse approach is used to produce desired intensity levels at the control points. Then, the pressure and intensity at any field point can be obtained by modeling the surface of the phased array. Besides, scattering force and gradient force at various positions are also evaluated to analysis their relative components to the effect of the acoustic tweezers. The radiation forces along and across the beam propagation direction for particles with various diameter and positions are calculated based on the acoustic ray approach. In addition, this paper reveals the trapping trajectory along and across the beam propagation directions with the consideration of acoustic streaming effect and the resistance due to viscosity of water. The results obtained demonstrate that the acoustic tweezers have feasibility for multiple trapping in both the axial and lateral directions. They could be used for non-contact and non-invasive manipulation of small particles such as cell sorting, trapping, and manipulation. Cancerous cells could be separated from healthy cells. The measurement of mechanical properties of the cell could also be achieved by acoustic tweezers by non-contact squeezing the cell. In addition, they could be used for drug delivery to drive particles along and across the wave propagation direction.