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]

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]

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.

 

Fig. 1. Ray model of (2 + 1)D Airy beam. Each point X is crossed by four rays at the same time. The transverse directions of rays are represented by arrows, the length of arrow is proportional to the sine of the angle between the ray and the z axis. (a) $z = 0$; (b) $z = 90$ μm and (c) $z = 180$ μm. Along the diagonal $y = x$, the magnitude of arrows becomes stronger when the beam propagates, which implies the self-acceleration of Airy beam (see Visualization 1).

Download Full Size | PDF

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₁ strikes on point X₁ and X₄ in sequence [Fig. 2(a)], thus, it is an effective ray on the sphere at point X₁. In contrast, ray L₂ strikes on point X₃ first before it hits on point X₂, so that ray L₂ is neglected for force calculation at point X₂.

 

Fig. 2. Sketch of single ray tracing. (a) Ray validation. L₁ is valid (green tick) and L₂ is invalid (red cross). (b) Physical vectors associated with force calculation. Incident point, X, the outward normal unit vector, N₁. Unit vectors E and N₂ indicate directions of beam polarization, and the incident plane normal direction, respectively.

Download Full Size | PDF

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₁ (L₂) has two intersection points on spherical surface: point X₁ (X₂) where ${t_1} = 0$, and point X₄ (X₃) where ${t_2} > 0$ (${t_2} < 0$). Obviously, the effective ray S at point X satisfies

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, dP is the light power

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₁ 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₁ 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₂ 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₂, 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

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₁ [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₂ [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.

 

Fig. 3. Trapping efficiency of microbead along specific paths. Transverse intensity profiles of the Airy beam at propagation distance of (a) $z = 0$, (b) $z = 90$ μm and (c) $z = 130$ μm. White straight lines l₁ and l₂ indicate the positions under study. Trapping efficiencies along respective lines [(d), (e) and (f) along l₁, and (g), (h) and (i) along l₂] for Q_x [(d) and (g)], Q_y [(e) and (h)], and Q_z [(f) and (i)].

Download Full Size | PDF

When the microsphere is placed along l₁, 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_z are always positive in the beam region, and push the microsphere towards positive 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₂. The trapping efficiencies Q_x also oscillate between positive and negative x-direction as in Fig. 3(g). Since l₂ is the symmetry axis of intensity pattern, the magnitude of Q_y is naturally the same as the magnitude of Q_x. The direction of Q_y shows the same oscillation behavior between positive and negative y-direction as the Q_x [Fig. 3(h)]. The maximum trapping efficiency is found near the major lobe of the beam, which is on the diagonal line l₂. The axial efficiency Q_z is always non-negative [Fig. 3(i)] and shows the same trend as that of in Fig. 3(f). The peak of Q_z 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.

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₂). 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)].

 

Fig. 4. Stream lines of optical force in (a) xoy plane, (b) vertical plane along line y = x (l₂), and (c) (2 + 1)D. The caustic surface is demonstrated by a translucent surface in the 3D stream line plot (c). Optical forces with $|{\textbf Q} |< |{{{\textbf Q}_{{\textrm {max}}}}} |/50$ are not shown due to negligible magnitude.

Download Full Size | PDF

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₁ 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₀ and y₀ remain unchanged. For particles with higher relative refractive index, ray optics approach is also applicable [17,27].

 

Fig. 5. Trapping efficiency on microsphere with $\rho = 2\; \mu \mbox{m}$ in three cases of refractive index. The microsphere is placed at $z = 0$ plane along line l₁ [Fig. 3(a)]. The (a) Q_x, (b) Q_y and (c) Q_z components on microsphere with refractive indexes of 1.46 (black), 1.57 (red) and 1.65 (blue), respectively.

Download Full Size | PDF

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₁ in Fig. 3(a)]. The transverse trapping efficiency Q_x gets stronger when radius increases, and the peak position of Q_x 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 Q_z curve [Fig. 6(b)]. This splitting is attributed to the asymmetry of Airy beam. To reveal this phenomenon, the total Q_x is decomposed into Q_x,l and Q_x,r, where Q_x,l denotes the summation of x-direction trapping efficiencies among those rays that strike on the surface of the left hemisphere of the microsphere, and Q_x,r is x-direction trapping efficiencies yielded by rays that strike on the right hemisphere. Meanwhile, Q_z is also decomposed into Q_z,l, the axial trapping efficiency originating from rays that hit the surface of left hemisphere, and Q_z,r, on the surface of right hemisphere. Figure 6(c) shows the definition of left and right hemisphere in our simulations.

 

Fig. 6. Splitting trapping efficiency along line l₁. The microsphere with ${n_b} = 1.46$ is placed at $z = 0$ plane. The Q_x (a) and Q_z (b) components of the trapping efficiency are summation of trapping efficiency on left and right hemispheres (c) in three cases of radius ρ with 2 μm, 4 μm and 6 μm. (d), (e) and (f) are decomposed trapping efficiency Q_x on left and right sphere for ρ=2 μm, 4 μm and 6 μm. (g), (h) and (i) are decomposed trapping efficiency Q_z on left and right sphere for ρ=2 μm, 4 μm and 6 μm.

Download Full Size | PDF

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_x,r and Q_x,l along 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_x,r then drags the sphere toward positive x-direction. Meanwhile, the trapping force on left hemisphere Q_x,l pulls the sphere toward negative x-direction with smaller magnitude. When the microsphere enters the main lobe, the modulus of Q_x,r and Q_x,l 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 Q_x,r will reach its maximum value compared with modulus of Q_x,l, which creates the highest peak of curve Q_x,r. After that, the right hemisphere is about to leave the main lobe, so that Q_x,r begins to decrease and the modulus of Q_x,l will subsequently increase to its maximum. When the left hemisphere is leaving the main lobe, the modulus of Q_x,l begins to decrease. In general, there will be a lag of x-position between two poles of Q_x,r and Q_x,l, and the distance of the lag is proportional to the size of radius of microsphere [Figs. 6(d)–(f)].

The same analyses would also apply to Q_z which is the sum of Q_z,l and Q_z,r. Figures 6(g)–(i) show the decomposed Q_z. There is also a separation of x-coordinates between two peaks of Q_z,l and Q_z,r [see Figs. 6(g)–(i)]. Since the Q_z,l and Q_z,r are both toward positive z-direction, there will be only one peak in the line Q_z when the x-distance between the two peaks is relatively small. For example, when $\rho = 2$ μm, there is only one peak in line Q_z as shown in Fig. 6(b). When $\rho = 4$ μm or 6 μm, the positional separation along x axis between two peaks becomes remarkable [Fig. 6(h) and Fig. 6(i)], and the Q_z curve will then display two peaks which are relatively higher than other peaks.

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

 

Fig. 7. Optical torque on ellipsoid object mediated by asymmetric force in Airy beam. The axial optical force on right part of the object is stronger than on the left part, which will induce an optical torque on the object. (a) initially, the object lays horizontal in the transverse plane. (b) the object gradually rotates under the asymmetric axial optical forces. (c) The object reaches the final vertical configuration.

Download Full Size | PDF

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.