## Abstract

The optical forces of a (2 + 1)D Airy beam on a microsphere are studied in the ray optics regime. The ray model of a (2 + 1)D Airy beam is derived from its Fourier angular spectrum using a stable aggregate of the flexible elements theory. Numerical results demonstrate that the microsphere can be trapped by the transverse optical force and pulled towards the beam major lobe. Longitudinal optical forces further push the microsphere towards the positive *z*-direction. The trend for the movement of a microsphere in an Airy beam is clearly demonstrated as the stream line of optical forces, which is consistent with the observed phenomena in optical trapping experiments. In the meantime, both the transverse and the longitudinal optical forces increase when the relative refractive index of the trapped microsphere increases. Calculation of optical forces on microspheres with larger size reveals that the optical forces contributed by rays on each hemisphere are actually different due to the asymmetry of the Airy beam. The force difference could cause natural torque on the trapped objects if they are ovals or other asymmetric shapes.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

By use of counter-propagating Gaussian beams or a highly focused Gaussian beam, a dielectric particle or a metallic nanoparticle illuminated by the incident light can be stably trapped in the optical potential [1]. In addition, with the development of wavefront shaping technique, light beams possessing fascinating optical features play increasingly important roles in optical trapping [2,3]. Apart from stable trapping, light beam with angular momentum is able to rotate a rod-like particle or guide trapped particles along a designed trajectory. Thereinto, self-accelerating light beam has been developed into a powerful tool for optical guiding. The most widely used self-accelerating light beam is the Airy beam which was first predicted as an exact solution for the free-particle nonlinear Schrödinger equation (NLSE) [4]. Experimental Airy beam was first produced in the optical domain via spectrum modulation exerted by high resolution spatial light modulator (SLM) [5]. Such Airy beam possesses nondiffracting, self-healing and self-accelerating features [6]. The beam lobes automatically bend along a parabolic trajectory while the beam propagates, this feature makes it more attractive in many optical applications including the particle clearing [7], optical path clearing [8] and optical redistribution of microparticles and cells [9]. In those experiments, particles or cells were firstly dragged into the area of beam lobes, and then moved along the beam’s bending trajectory under the action of optical force. Therefore, theoretical analysis of the optical force in Airy beam helps us understand the experimental phenomena, and design Airy beams with novel features including higher trapping quality.

Prior studies of the optical force rely on the integration of the Maxwell stress tensor (MST) for the scattering electromagnetic field in the frame of electrodynamics (EM). EM models include the dipole approximation [10], arbitrary beam theory [11] and numerical method such as T-matrix [12] and finite-difference time-domain method (FDTD) [13]. In general, they can be applied to both Rayleigh [14,15] and Mie particles [16] for Airy beam, however, complex mathematical procedures are involved when the particle size increases. When the size of trapped particle is much larger than the wavelength, optical force calculation can be done in the framework of geometrical optics also known as ray optics (RO). Due to its simplicity and immediacy, the RO model has been widely developed in the simulation of optical trapping by fundamental Gaussian beam [17-19], and provides reasonable results. Although, it is difficult to find appropriate ray models to represent the characteristics for some newly developed light beams, the ray model of Airy beam can be found under the framework of combination of geometrical optics and catastrophe optics [20,21], which discusses the relationship between optical field and light rays accompanied with the caustic curve formation. Moreover, the ray model of Airy beam can be also derived by the use of the stable aggregates of flexible elements (SAFE) theory [22,23].

Here, we use RO model to analyze the trapping force on microspheres in (2 + 1)D Airy beam. The ray model of (2 + 1)D Airy beam is derived from the Fourier angular spectrum of Airy beam according to the SAFE theory. Using this ray model, the optical force acted on a microparticle under Airy beam can be fully calculated in 3D. The stream line map of trapping force suggests that the moving tendency of trapped particle. The simulation results are consistent with the experimental phenomena. As the particle size increases, the force efficiency of Airy beam shows significant splitting due to the asymmetric energy flow. This phenomenon attributes to the asymmetric illumination of the Airy beam on the microsphere, and is well explained by the RO model. Finally, we also intuitively predict that the asymmetric illumination and optical force on the ellipsoid particles will induce optical torque and alignment.

## 2. Ray model of (2 + 1)D Airy beam

For a finite energy (2 + 1)D Airy beam at $z = 0$ plane, the Fourier angular spectrum $A(u,v,0)$ is proportional to [24]

*x*

_{0}(

*y*

_{0}) and

*a*(

*b*) are the length scale and non-negative decay factors along

*x*-axis (

*y*-axis),

*u*and

*v*are the directional cosine of the inclined angle of wave vectors along

*x*-axis and

*y*-axis. Since Airy beam is the exact solution of paraxial Helmholtz wave equation [4,24], its Fourier spectrum $A(u,v,z)$ at any given

*z*plane can be determined by Fresnel diffraction [25]. Then, $A(u,v,z)$ is derived as

*x*and

*y*of a point

**X**locating on that ray at any given

*z*coordinate can be determined by Fourier angular spectrum of Airy beam, and they are given by

In homogeneous medium, each ray can be regarded as an infinite line which launches from the reference point **X** and propagates along its direction **S**. Thus, the parametric equation of each ray is ${\textbf X}^{\prime} = {\textbf X} + t \cdot {\textbf S}$, where **X**’ is an arbitrary point locating on the ray and *t* is the parameter.

Usually, the coordinates of a reference point **X** for a ray at a given axial coordinate $z = \upsilon$ are pre-known, for example ${\textbf X} = (\zeta ,\eta ,\upsilon )$. In order to find the direction vectors of rays that pass through point **X**, we substitute the coordinates of **X** $(\zeta ,\eta ,\upsilon )$ into Eq. (3) to figure out parameters *u* and *v*. There are two solutions of Eq. (3a) for parameter *u*, while there are also two solutions for parameter *v* in Eq. (3b). The combination of *u* and *v* has four pairs which indicate that 4 rays are passing through the given point **X** simultaneously. Their direction vectors **S** are respectively given by

The light intensity of (2 + 1)D Airy beam at any given points $(x,y,z)$ can be found from Eq. (2). The closed-form expression of the intensity distribution is [24]

*C*is the power normalization constant which is given by

Figure 1 demonstrates the normalized intensity distribution of (2 + 1)D Airy beam with the corresponding ray model (positions and transverse directions). Here, ${\lambda _0} = 1.064$ μm, ${x_0} = {y_0} = 2$ μm, $a = b = 0.1$ and ${n_m} = 1.33$ for water. The intensity of (2 + 1)D Airy beam at $z = 0$ plane is mainly located in the third quadrant [Fig. 1(a)]. There are 4 rays (white arrow) launching from each given point **X** (gray point). The length of arrow is proportional to the sine of the angle between the ray and the *z* axis, which also indicates the magnitude of ray’s inclination. After beam propagating at $z = 90$ μm plane as shown in Fig. 1(b), the Airy beam is accelerating towards the first quadrant, rays’ directions are also inclined towards first quadrant. Figure 1(c) shows the ray model of Airy beam at $z = 180$ μm plane, the arrows of rays become longer than that of in Fig. 1(b) and Fig. 1(a) which implies the accelerating property of Airy beam.

## 3. Analyses of trapping forces

Assuming that a microsphere with refractive index *n _{b}* centered at ${\textbf c} = ({x_c},{y_c},{z_c})$, the surface of this sphere is described by parametric equations: {$x = {x_c} + \rho \sin \theta \cos \varphi $, $y = {y_c} + \rho \sin \theta \sin \varphi $, $z = {z_c} + \rho \cos \theta $}, where

*ρ*is radius of the microsphere, altitude $\theta \in [0,\pi ]$ and azimuth $\varphi \in [0,2\pi ]$. When a microsphere is illuminated by an Airy beam, there are light rays passing through points on the microsphere’s surface. The light intensities at those points are given by Eq. (5). In this situation, we choose the point on microsphere’s surface as the reference point

**X**for each group of 4 relevant rays. Note that not all ray equations determined by Eq. (4) would be considered as effective ray on this point. Rays with first intersection point close to the light source would be considered as an effective ray. For example, ray

**L**

_{1}strikes on point

**X**

_{1}and

**X**

_{4}in sequence [Fig. 2(a)], thus, it is an effective ray on the sphere at point

**X**

_{1}. In contrast, ray

**L**

_{2}strikes on point

**X**

_{3}first before it hits on point

**X**

_{2}, so that ray

**L**

_{2}is neglected for force calculation at point

**X**

_{2}.

We take advantage of the parametric equations of rays to distinguish which rays are effective for force calculation at arbitrary points on the microsphere. Recall that in homogeneous medium, ray trajectory is a straight line, and its parametric equation is ${\textbf X}^{\prime} = {\textbf X} + t \cdot {\textbf S}$. In this condition reference point of rays **X** is chosen as the point on the surface of the microsphere. Substituting the coordinate of **X** into Eq. (4), we then get four directions of ray. The intersection point between light ray and the microsphere can be found by solving the combination of ray equation and sphere equation, which gives a quadratic equation about parameter *t*. The roots to the quadratic equation correspond to two intersection points between the light ray and the sphere with parameters ${t_1} = 0$ (point **X**) and ${t_2}$ (second intersection point between ray and sphere). For example, in Fig. 2(a), ray **L**_{1} (**L**_{2}) has two intersection points on spherical surface: point **X**_{1} (**X**_{2}) where ${t_1} = 0$, and point **X**_{4} (**X**_{3}) where ${t_2} > 0$ (${t_2} < 0$). Obviously, the effective ray **S** at point **X** satisfies

**X**must be the first intersection point between ray and sphere when we trace the ray along its propagating direction. It is an important criterion to determine which rays are valid for the force calculation.

As shown in Fig. 2(b), when a valid incident ray **L**_{I} strikes on an infinitesimal region $\mbox{d}\sigma $ (red curve) that close to point **X** in the sphere surface, due to the mismatch of the refractive indexes between medium and microbead, **L**_{I} will produce reflected (not shown) and refractive ray **L**_{R}. $\alpha $ is the incident angle of ray **L**_{I} at point **X**, while $\beta $ is the refraction angle of **L**_{R}. In the meantime, the momentum of photon changes, and generates optical force on that bead. The optical force can be decomposed into a scattering component $\mbox{d}{{\textbf F}_{\textrm s}} = {n_m}/c{Q_s}\mbox{d}P{{\hat{\textbf F}}_{\textrm s}}$ and a gradient component $\mbox{d}{{\textbf F}_{\textrm g}} = {n_m}/c{Q_g}\mbox{d}P{{\hat{\textbf F}}_{\textrm g}}$, where *c* is the speed of light, d*P* is the light power

**X**[18], and

*P*

_{0}is the measured total power of the Airy beam. Note that there are four rays intersecting at point

**X**with the same weight, so that the intensity should be uniformly assign to each of those four rays, which leads to a constant ratio of 1/4 in the right-hand side of Eq. (8).

The direction of $\mbox{d}{{\textbf F}_{\textrm s}}$ is along unit vector ${{\hat{\textbf F}}_{\textrm s}}$ which is parallel to the incident ray, and the direction of $\mbox{d}{{\textbf F}_{\textrm g}}$ is along unit vector ${{\hat{\textbf F}}_{\textrm g}}$ which is perpendicular to the incident ray and toward outside of the sphere. Let the surface normal **N**_{1} points outward at incident point **X**, ${{\hat{\textbf F}}_{\textrm s}}$ and ${{\hat{\textbf F}}_{\textrm g}}$ can be found by the using spatial analytic geometry method according to the directions of **N**_{1} and **L**_{I} [19]. The scattering and gradient trapping efficiencies *Q*_{s} and *Q*_{g} are respectively given by [17],

*R* and *T* in Eq. (9) are the Fresnel reflection and transmission coefficients of energy flow for a single ray. According to conservation of energy, $R + T \equiv 1$, however, they are related to many factors, including incident position on microsphere, incident angle and beam polarization. For a linearly polarized beam, the electric unit vector is ${\textbf E} = ({\cos \xi ,\sin \xi ,0} )$ which is perpendicular to the ray, and the angle *ξ* describes the polarization direction relative to the *x*-axis. When ray **L**_{I} hits the sphere surface at point **X** [Fig. 2(b)], the direction normal **N**_{2} of the incident plane Ω is determined by ${{\textbf N}_2} = {{\textbf L}_{\textrm I}} \times {{\hat{\textbf F}}_{\textrm g}}$. Defining that *γ* is the angle between vector **E** and **N**_{2}, the reflection and transmission coefficients *R* and *T* of linear polarized beam are given by [26]

For circularly polarized or unpolarized beam, $\gamma \equiv \pi /4$. In the simulation presented in this paper, the rays are considered to be with circular polarization.

From Eqs. (8) to (10), the trapping efficiencies of a single ray at a surface element are

**Q**

_{s}and

**Q**

_{g}on the microsphere at any position ${\textbf c} = ({x_c},{y_c},{z_c})$ are the sum of scattering efficiency d

**Q**

_{s}and gradient efficiency d

**Q**

_{g}. They are given by

**X**on the microsphere surface that satisfy the validation in Eq. (7). We use

*Q*,

_{x}*Q*, and

_{y}*Q*to denote the components of total trapping efficiency ${{\textbf Q}_{\textrm s}} + {{\textbf Q}_{\textrm g}}$ along

_{z}*x*,

*y*, and

*z*directions.

## 4. Results and discussion

The transverse pattern of Airy beam shifts along the diagonal line during propagation (Fig. 1). We calculate the optical forces of a micron-sized particle in a circularly polarized Airy beam. The parameters for the Airy beam are chosen as ${\lambda _0} = 1.064$μm, ${x_0} = {y_0} = 2$ μm, and $a = b = 0.1$. The refractive indexes of the medium and the microsphere are ${n_m} = 1.33$ and ${n_b} = 1.46$, respectively. In order to elucidate the trapping efficiency distribution in transverse planes at different propagation distances where $z = 0$, $z = 90$ μm and $z = 130$ μm.

Figures 3(a)–(c) show the transverse beam profiles at corresponding locations. We are specifically interested in trapping performances along two lines in the transverse beam patterns [Figs. 3(a)–(c)]. The first line along *l*_{1} [where $y = {z^2}/(4{k^2}y_0^3)$] is the caustics of Airy beam parallel to *x*-axis, while the second line *l*_{2} [where $y = x$] is the symmetry axis of transverse beam profile. Analysis of the trapping efficiencies along these two lines helps us understand how the microspheres are manipulated and guided in 3D in the Airy beam field. We thus evaluate the trapping performances along these two lines at three axial locations for a microsphere with radius $\rho = 2$μm.

When the microsphere is placed along *l*_{1}, the transverse trapping efficiencies ${Q_x}$ oscillate between positive and negative *x-*directions, and the sphere restores at stable trapping points with ${Q_x} = 0$. The oscillation has been verified along all the transverse planes with different propagation distances, while representative curves at $z = 0$, $z = 90$ μm and $z = 130$ μm are shown in Fig. 3(d). These restoring forces compel the microsphere to move toward each beam lobes. As the *z* increases, the maximum restoring force decreases and the stable trapping point shifts together with the central position of the main lobe. The transverse trapping efficiencies *Q _{y}* along

*y*-direction are, however, non-positive [Fig. 3(e)], which attracts the microsphere toward the beam region and prevents the microsphere from escaping. The axial trapping efficiencies

*Q*are always positive in the beam region, and push the microsphere towards positive

_{z}*z*-direction [Fig. 3(f)]. To point out that in Figs. 3(d), (e) and (f), the peak [or valley in Fig. 3(e)] values decrease for large propagation distance and shift towards positive

*x*-direction due to self-acceleration of the Airy beam.

Next, we examine the trapping performance along the diagonal line *l*_{2}. The trapping efficiencies *Q _{x}* also oscillate between positive and negative

*x-*direction as in Fig. 3(g). Since

*l*

_{2}is the symmetry axis of intensity pattern, the magnitude of

*Q*is naturally the same as the magnitude of

_{y}*Q*. The direction of

_{x}*Q*shows the same oscillation behavior between positive and negative

_{y}*y-*direction as the

*Q*[Fig. 3(h)]. The maximum trapping efficiency is found near the major lobe of the beam, which is on the diagonal line

_{x}*l*

_{2}. The axial efficiency

*Q*is always non-negative [Fig. 3(i)] and shows the same trend as that of in Fig. 3(f). The peak of

_{z}*Q*is found near the main lobe of beam. The magnitude of axial trapping efficiency decreases along z- direction. Owing to the self-acceleration, the maximum position shifts towards the corresponding self-accelerating direction of the Airy beam.

_{z}To analyze the tendency of motion of a trapped microsphere, we mapped the force distribution in two typical planes. Figure 4(a) shows the force stream lines in the transverse plane where $z = 0$ (the *xoy* plane). The background profile indicates the normalized intensity distribution of the Airy beam in that plane. The stream lines (red lines) of the transverse forces predicts the probable motion trajectories of the microsphere exposed by the Airy beam. Under the transverse optical forces, the microsphere is pulled toward the beam region in the transverse force field. To examine how the microsphere moves in the axial plane, we mapped the force distribution in the vertical plane determined by *z*-axis and diagonal line $y = x$ (*l*_{2}). Stream lines of optical force show that the microsphere is pulled to local lobes at low power region, and further pulled to the major lobes due to the interplay of Brownian motion and the transverse trapping forces when the laser power (thus optical force) increases [Fig. 4(b)]. The axial motion of the particle follows the parabolic trajectory determined by the bending direction of Airy beam. To provide a clear 3D picture, 3D stream lines of optical force are mapped to show the possible motion trajectory of microsphere that immerged in this (2 + 1)D Airy beam [Fig. 4(c)].

The 3D motion of trapped microsphere can be classified into two categories. First, if the microsphere is initially placed at the major lobe, it will be pushed along the self-bending trajectory of main lobe. Second, if the microsphere is placed at the neighboring side lobes, it will experience two stages, initial lifting along the nearby side lobe, and the subsequent attraction into major lobe before it is further pushed by the beam in the major lobe. This phenomenon is intuitively explained by the fast intensity decay on the side lobes. As a result, the transverse optical forces at major lobe are much larger on magnitude than that experiences at neighboring side lobes. This is confirmed by the fact that the branch stream lines in Fig. 4(c) converge into one line on the main lobe after a certain propagation distance. The 3D stream lines explanation in Fig. 4(c) is consistent with the movement of microparticle in optical guiding experiment by Airy beam and expounds the experimental phenomena [7].

In real trapping experiment, the particles have variant physical properties, including the refractive index. For example, the polystyrene microsphere has a refractive index of 1.57, while the refractive index for a silica microsphere is 1.46 at 1064 nm wavelength. The trapping ability will be greatly affected by these physical properties. We thus evaluate the trapping performance for particles with different physical properties. Figure 5 shows the variations of trapping efficiencies along the caustic line *l*_{1} at initial position at $z = 0$ in three cases of refractive index. In our simulations, the refractive indexes are chosen to be 1.46, 1.57 and 1.65 respectively. The trapping efficiency for all three orthogonal directions becomes stronger for particles with higher refractive indexes since the change of momentum of incident rays is enhanced. Meanwhile, the positions of the peak of trapping efficiency distribution keep the same while refractive index varies, since the size of microsphere, beam parameters *x*_{0} and *y*_{0} remain unchanged. For particles with higher relative refractive index, ray optics approach is also applicable [17,27].

We further study how the trapping efficiency varies with the microsphere radius $\rho $. Figure 6 shows the trapping efficiency for a microsphere with ${n_b} = 1.46$ along line $y = z = 0$ [*l*_{1} in Fig. 3(a)]. The transverse trapping efficiency *Q _{x}* gets stronger when radius increases, and the peak position of

*Q*shifts away from main lobe, due to the fact that large microsphere is easier to be irradiated by the major lobe as compared with the small ones [Fig. 6(a)]. When particle size increases, we observed a splitting of the peak on the axial trapping efficiency

_{x}*Q*curve [Fig. 6(b)]. This splitting is attributed to the asymmetry of Airy beam. To reveal this phenomenon, the total

_{z}*Q*is decomposed into

_{x}*Q*and

_{x,l}*Q*, where

_{x,r}*Q*denotes the summation of

_{x,l}*x*-direction trapping efficiencies among those rays that strike on the surface of the left hemisphere of the microsphere, and

*Q*is

_{x,r}*x*-direction trapping efficiencies yielded by rays that strike on the right hemisphere. Meanwhile,

*Q*is also decomposed into

_{z}*Q*, the axial trapping efficiency originating from rays that hit the surface of left hemisphere, and

_{z,l}*Q*, on the surface of right hemisphere. Figure 6(c) shows the definition of left and right hemisphere in our simulations.

_{z,r}Figures 6(d)–(f) show the decomposed *Q _{x}* from contributing rays hitting on the left and right hemispheres. We found that there is a separation between maximum of modulus of

*Q*and

_{x,r}*Q*along

_{x,l}*x*-axis. Assume a microsphere moves from left to right and experiences radiation from the main lobe of Airy beam, the right hemisphere first enters the main lobe and is followed by its left hemisphere. The

*x*-direction trapping force contributed by rays on right hemisphere

*Q*then drags the sphere toward positive

_{x,r}*x*-direction. Meanwhile, the trapping force on left hemisphere

*Q*pulls the sphere toward negative

_{x,l}*x*-direction with smaller magnitude. When the microsphere enters the main lobe, the modulus of

*Q*and

_{x,r}*Q*starts to increase owing to larger intensity of main lobe than nearby side lobes. Note that the right hemisphere enters the main lobe first, the modulus of

_{x,l}*Q*will reach its maximum value compared with modulus of

_{x,r}*Q*, which creates the highest peak of curve

_{x,l}*Q*. After that, the right hemisphere is about to leave the main lobe, so that

_{x,r}*Q*begins to decrease and the modulus of

_{x,r}*Q*will subsequently increase to its maximum. When the left hemisphere is leaving the main lobe, the modulus of

_{x,l}*Q*begins to decrease. In general, there will be a lag of

_{x,l}*x*-position between two poles of

*Q*and

_{x,r}*Q*, and the distance of the lag is proportional to the size of radius of microsphere [Figs. 6(d)–(f)].

_{x,l}The same analyses would also apply to *Q _{z}* which is the sum of

*Q*and

_{z,l}*Q*. Figures 6(g)–(i) show the decomposed

_{z,r}*Q*. There is also a separation of

_{z}*x*-coordinates between two peaks of

*Q*and

_{z,l}*Q*[see Figs. 6(g)–(i)]. Since the

_{z,r}*Q*and

_{z,l}*Q*are both toward positive

_{z,r}*z*-direction, there will be only one peak in the line

*Q*when the

_{z}*x*-distance between the two peaks is relatively small. For example, when $\rho = 2$ μm, there is only one peak in line

*Q*as shown in Fig. 6(b). When $\rho = 4$ μm or 6 μm, the positional separation along

_{z}*x*axis between two peaks becomes remarkable [Fig. 6(h) and Fig. 6(i)], and the

*Q*curve will then display two peaks which are relatively higher than other peaks.

_{z}Since the optical force contributed by rays from left and right hemispheres are quite different, they would cause a natural torque on the trapped objects as long as the shape is non-spherical. The object will rotate under the asymmetry-mediated torque. For instance, if the shape of the particle is an ellipsoid, and initial configuration of the particle suspended in the sample chamber is with long axis horizontal to the transverse plane [Fig. 7(a)]. Under the action of Airy beam, the right-hand side will experience stronger scattering force than the left part [Fig. 7(b)]. Therefore, the right side will be pushed more forward, creating the rotation of the particle until it reaches the balanced rotational configuration [Fig. 7(c)].

## 5. Concluding remarks

The ray model of (2 + 1)D paraxial Airy beam is derived using the Fourier spectrum of Airy beam and SAFE theory. This construction applies for any optical field with slowly varying spectral amplitude. The ray model was confirmed by calculation of the trapping efficiencies along *x*, *y* and *z* directions. The calculation shows that the Airy beam is able to trap particles towards its intensity lobes transversely and to push particles axially. In the action of axial pushing force, the trapped particles shift along the bending trajectory of Airy beam. The simulation results are consistent with that of in experiment. We also demonstrate that the trapping efficiency increases with the refractive index of microsphere. Due to the asymmetry of Airy beam, the axial trapping efficiency would have two relatively higher peaks when the radius of microsphere is greater than the size of beam main lobe. Although the ray optics method neglects the diffraction effect of optical field, it can still predict the trapping performance on microspheres with radius $\rho < 10\lambda $ [17], both axially and transversely. In conclusion, the ray optics method applies for the simulation of optical force in self-accelerating Airy beam, and gives a rather straight-forward approach to study the trapping force on particles with size larger than wavelength.

## Funding

Scientific Research Foundation of the Institute for Translational Medicine of Anhui Province (2017zhyx25); Excellent Young Talents Fund Program of Higher Education Institutions of Anhui Province (KJ2016A361); Anhui Medical University (XJ201518).

## References

**1. **K. Svoboda and S. M. Block, “Optical trapping of metallic Rayleigh particles,” Opt. Lett. **19**(13), 930–932 (1994). [CrossRef]

**2. **K. Dholakia and T. Čižmár, “Shaping the future of manipulation,” Nat. Photonics **5**(6), 335–342 (2011). [CrossRef]

**3. **M. Woerdemann, C. Alpmann, M. Esseling, and C. Denz, “Advanced optical trapping by complex beam shaping,” Laser Photon. Rev. **7**(6), 839–854 (2013). [CrossRef]

**4. **M. V. Berry and N. L. Balazs, “Nonspreading wave packets,” Am. J. Phys. **47**(3), 264–267 (1979). [CrossRef]

**5. **G. A. Siviloglou, J. Broky, A. Dogariu, and D. N. Christodoulides, “Observation of Accelerating Airy Beams,” Phys. Rev. Lett. **99**(21), 213901 (2007). [CrossRef]

**6. **J. Broky, G. A. Siviloglou, A. Dogariu, and D. N. Christodoulides, “Self-healing properties of optical Airy beams,” Opt. Express **16**(17), 12880–12891 (2008). [CrossRef]

**7. **J. Baumgartl, M. Mazilu, and K. Dholakia, “Optically mediated particle clearing using Airy wavepackets,” Nat. Photonics **2**(11), 675–678 (2008). [CrossRef]

**8. **J. Baumgartl, T. Čižmár, M. Mazilu, V. C. Chan, A. E. Carruthers, B. A. Capron, W. McNeely, E. M. Wright, and K. Dholakia, “Optical path clearing and enhanced transmission through colloidal suspensions,” Opt. Express **18**(16), 17130–17140 (2010). [CrossRef]

**9. **J. Baumgartl, G. M. Hannappel, D. J. Stevenson, D. Day, M. Gu, and K. Dholakia, “Optical redistribution of microparticles and cells between microwells,” Lab Chip **9**(10), 1334–1336 (2009). [CrossRef]

**10. **A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, “Observation of a single-beam gradient force optical trap for dielectric particles,” Opt. Lett. **11**(5), 288–290 (1986). [CrossRef]

**11. **J. P. Barton, D. R. Alexander, and S. A. Schaub, “Theoretical determination of net radiation force and torque for a spherical particle illuminated by a focused laser beam,” J. Appl. Phys. **66**(10), 4594–4602 (1989). [CrossRef]

**12. **T. A. Nieminen, H. Rubinsztein-Dunlop, N. R. Heckenberg, and A. I. Bishop, “Numerical Modelling of Optical Trapping,” Comput. Phys. Commun. **142**(1-3), 468–471 (2001). [CrossRef]

**13. **Z. Fei, G. Xiaosong, X. Wendong, and G. Fuxi, “Comment on: Computation of the optical trapping force using an FDTD based technique,” Opt. Express **14**(25), 12494–12496 (2006). [CrossRef]

**14. **H. Cheng, W. Zang, W. Zhou, and J. Tian, “Analysis of optical trapping and propulsion of Rayleigh particles using Airy beam,” Opt. Express **18**(19), 20384–20394 (2010). [CrossRef]

**15. **W. Lu, H. Chen, S. Liu, S. Liu, and Z. Lin, “Rigorous full-wave calculation of optical forces on dielectric and metallic microparticles immersed in a vector Airy beam,” Opt. Express **25**(19), 23238–23253 (2017). [CrossRef]

**16. **Z. Zhao, W. Zang, and J. Tian, “Optical trapping and manipulation of Mie particles with Airy beam,” J. Opt. **18**(2), 025607 (2016). [CrossRef]

**17. **A. Ashkin, “Forces of a Single-Beam Gradient Laser Trap on a Dielectric Sphere in the Ray Optics Regime,” Biophys. J. **61**(2), 569–582 (1992). [CrossRef]

**18. **E. Sidick, S. D. Collins, and A. Knoesen, “Trapping forces in a multiple-beam fiber-optic trap,” Appl. Opt. **36**(25), 6423–6433 (1997). [CrossRef]

**19. **J. H. Zhou, H. L. Ren, J. Cai, and Y. M. Li, “Ray-tracing methodology: application of spatial analytic geometry in the ray-optic model of optical tweezers,” Appl. Opt. **47**(33), 6307–6314 (2008). [CrossRef]

**20. **D. Ludwig, “Wave propagation near a smooth caustic,” Bull. Amer. Math. Soc. **71**(5), 776–780 (1965). [CrossRef]

**21. **M. V. Berry and C. Upstill, “Catastrophe optics: morphologies of caustics and their diffraction patterns,” Prog. Opt. **18**, 257–346 (1980). [CrossRef]

**22. **M. A. Alonso and G. W. Forbes, “Using rays better. II. Ray families to match prescribed wave fields,” J. Opt. Soc. Am. A **18**(5), 1146–1159 (2001). [CrossRef]

**23. **M. A. Alonso and G. W. Forbes, “Stable aggregates of flexible elements give a stronger link between rays and waves,” Opt. Express **10**(16), 728–739 (2002). [CrossRef]

**24. **G. A. Siviloglou and D. N. Christodoulides, “Accelerating finite energy Airy beams,” Opt. Lett. **32**(8), 979–981 (2007). [CrossRef]

**25. **J. W. Goodman, * Introduction to Fourier Optics* (McGraw-Hill, New York, 1968).

**26. **M. Born and E. Wolf, * Principles of Optics* (Cambridge University Press, 1999).

**27. **L. A. Ambrosio and H. E. Hernández-Figueroa, “Inversion of gradient forces for high refractive index particles in optical trapping,” Opt. Express **18**(6), 5802–5808 (2010). [CrossRef]