Calculation of optical forces for arbitrary light beams using the Fourier ray method

Meng Shao; Meng Shao; Shuhe Zhang; Shuhe Zhang; Jinhua Zhou; Yu-Xuan Ren

doi:10.1364/OE.27.027459

1. Introduction

Since photons carry momentum, the changes of their momentum induce forces on their interacting object. Utilizing such optical force, we can trap microscopic objects and guide them along specific trajectories. Depending on the transferred momentum, the optical force is categorized into scattering force and gradient force. The resultant forces of a highly focused beam can stably trap microparticles in three dimensions. This was named as optical tweezers which was invented by A. Ashkin [1], and is now widely used in many interdisciplinary fields. In life science, optical tweezers can be used to explore the motility of motor proteins [2], and to measure the elasticity of cells and biomacromolecule such as DNA and RNA [3], and protein [4]. Therefore, many investigations focused on how to establish an optical trapping system with high quality or advanced properties. With the development of wavefront shaping technique [5], light beams can possess fascinating optical features, and play increasingly important roles in optical manipulations. Apart from stably trapping a particle, light beams with self-reconstructing property, like Bessel beams, are able to trap many particles in multiple planes at the same time [6]. Self-accelerating light beam has been developed into a powerful tool for optical guiding [7,8]. In general, light beams with novel properties expand the applications of optical manipulation and have attracted extensive researches.

Simulation of optical force for a light beam acting on a particle gives us a way to understand its dynamical behaviors. The optical force can be calculated in both Electromagnetic (EM) model and Geometrical optics model, and the latter is also known as ray-optics model (RO model). The EM model relies on the integration of the Maxwell stress tensor for scattered electromagnetic field of incident light beam, including many different types of analytic and numerical methods [9–12]. Technically, they are able to calculate the optical force of trapped particles that are illuminated by any kind of beams. However, computational complexity scales up with particle size and the complexity of structured beam.

When the size of trapped particles is much larger than the wavelength of incident light beams, RO model is appropriate to calculate the optical force. Due to its simplicity and immediacy, the RO model has been widely developed in the simulation of optical trapping by fundamental Gaussian beam. Nevertheless, there are two fatal disadvantages for the RO model. First, the RO model ignores the field’s diffraction effect. So it fails to correctly calculate the optical force produced by diffractive optical field [13]. Second, the RO model fails to simulate the forces for the complex structured light beams. Since the use of RO model relies on ray tracing, finding an appropriate ray model for the incident beam, especially, for beams with complex structured optical fields, becomes the most challenging part in RO model, although there are different ways to model the incident beam by rays.

The traditional way to build a ray model for a light beam is based on the eiknoal equation in geometrical optics [14]. If the incident beam has a slow varying complex amplitude which is $U(x,y,z) = A(x,y,z)\exp [i\varphi (x,y,z)]$, the ray’s directional vector at point $(x,y,z)$ is given by $\nabla \varphi (x,y,z)/k$(k is the wave number) since the ray family are perpendicular to the wavefronts of the beam [14]. Accordingly, we can model some simple light beams. For example, for a paraxial Gaussian beam propagating along z-axis [Fig. 1(a)], it can be regarded as a family of spherical waves, and the centers of which locate on different positions in z-axis. The radii of the wavefront are given by $R(z) = z[1 + {({z_r}/z)^2}]$. Assuming that a ray launches from point C and strikes on a given point M. In order to find the direction of this incident ray, the coordinates of point C should be determined by the Pythagorean Theorem in the right triangle CMM’. Hence, line CM describes the direction of incident ray at point M. This model is suitable for paraxial Gaussian beam and gives reasonable simulation results in many previous works [15–17].

Fig. 1. Ray model of Gaussian beam. (a) Paraxial Gaussian beam. (b) Focused Gaussian beam through a high numerical objective lens.

Download Full Size | PDF

For a highly-focused Gaussian beam [13], since the size of the focal spot is very small compared with that of trapped particle, it can be regarded as a tiny dot. Hence, the focused Gaussian beam can be regarded as a part of a spherical wave converging to the focus, as shown in Fig. 1(b) [13]. In this way, the trajectories of light rays are given by converging lines from the wavefront to the focus. This model is known as ray-cone model or ray-pencil model, and is widely used and improved in simulation of force for highly-focused beams [18–21]. Notably, the ray-cone model is suitable for force simulation of speckle field as well, and gives reasonable results [22]. However, this model is not suitable for Gaussian beam in paraxial regime, and also breaks the fact that laser beam cannot be perfectly focused at a tiny point.

In addition, the ray model of the beam can be constructed according to its geometric shape. Since the envelope of the Gaussian beam is a hyperboloid of one sheet, the skew line ray model is proposed to model the nonparaxial Gaussian beam [23]. However, this model is not able to describe structured Gaussian beams. At present, there is no entity of RO model for arbitrary beams.

Apart from the spatial domain, the ray model for a light beam can be also constructed in the Fourier spectrum domain. If the incident beam has a slow varying Fourier angular spectrum which is $\tilde{U}(p,q,z) = B(p,q,z)\exp [i\Phi (p,q,z)]$, the ray with direction vector $\left( {p,q,\sqrt {1 - {p^2} - {q^2}} } \right)$ will pass through point $\nabla \Phi (p,q)/k$ [24]. Since the Fourier angular spectrum of self-accelerating Airy beam has a considerably simple mathematical expression [25], optical force on particles under this Airy beam can be analyzed according to RO model [26]. Unfortunately, there are large numbers of light beams that do not have abovementioned closed-form expressions neither in spatial nor in Fourier spectrum regimes, for example, the Hermite-Gaussian beam, Laguerre-Gaussian beam and Bessel beam [27]. It seems to be that the RO model is no longer able to calculate the optical force from those structured beams, since an appropriate ray scheme fails to model those intricate optical fields.

Notably, the RO model developed by Roosen [28] in 1977 is a very traditional method with unique simplicity, and is indeed a powerful technique for analyzing optical force on large particle. Here, we propose the Fourier ray (FR) method to calculate the optical force on a spherical particle for an arbitrary structured beam. Meanwhile, this method is suitable for calculation of the optical force in diffractive field as well. In the following texts, we start from the definition of FRs and demonstrate how these rays can be used to calculate the optical force. The FR method is validated by comparing the simulation results with previous method. We further take the ring beam, Hermite-Gaussian, and an arbitrary 3D structured beam as examples to demonstrate that the FR method is able to calculate optical force of particles in any incident beams.

2. Fourier rays

Let’s consider a scalar monochromatic light beam propagating in a homogeneous medium along positive z-direction. The optical field can be expressed by the superposition of plane waves (PW) through its Fourier angular spectrum $\tilde{U}(p,q,0)$ as [29]:

(1)$$U({x,y,z} )= \frac{{{C_0}}}{{{\lambda ^2}}}\iint\limits_{{p^2} + {q^2} \le 1} {\tilde{U}({p,q,0} )\exp \left[ {ik\left( {xp + yq + z\sqrt {1 - {p^2} - {q^2}} } \right)} \right]\textrm{d}p\textrm{d}q} ,$$

where $k = 2\pi {n_0}/\lambda$, $\lambda$ is the wavelength of the beam, ${n_0}$ is the refractive index of the medium, p and q are the directional cosines of the plane wave in x- and y-directions respectively, ${C_0}$ is a normalization coefficient of light power satisfying $\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{{|{U(x,y,0)} |}^2}\textrm{d}x\textrm{d}y} } = 1$. The complex amplitude of the plane wave that propagates in directional vector ${\textbf v}({p,q} )= (p,q,\sqrt[{}]{{1 - {p^2} - {q^2}}})$ is $A(x,y,z;p,q) = \tilde{U}(p,q,0)\exp [ik(xp + yq + z\sqrt {1 - {p^2} - {q^2}} )]$.

In geometrical optics, each plane wave can be regarded as a family of countless parallel rays spreading all over the real space. Figure 2 shows two sets of the plane waves (PW₁ and PW₂), both of which intersect across a given point C. Each plane wave is regarded as a family of parallel rays (red dashed line for PW₁ and blue dashed line for PW₂). Since these rays inlay on the plane wave weighted by the Fourier angular spectrum, we name these rays as the Fourier rays (FR). According to Eq. (1), the amplitudes of the FRs are given by $|{\tilde{U}(p,q,0)} |$, and remain unchanged with those Fourier rays transmitting.

Fig. 2. Sketch of Fourier rays.

Download Full Size | PDF

3. Optical force of Fourier ray family

An incident beam is decomposed into superposition of plane waves propagating in different directions. Since there are infinite numbers of FRs belonging to each specific plane wave, how can we utilize these FRs to calculate the optical force? Here, we stipulate that there is only one FR belonging to the corresponding plane wave and passing through a given point. Once the incident beam is decomposed by FRs, the optical force can be calculated according to following four steps.

First, effective FRs should be determined. Assuming that a microsphere with refractive index of ${n_b}$ is centered at point $c({x_c},{y_c},{z_c})$ [Fig. 3(a)]. Point $M({x_m},{y_m},{z_m})$ locates on the surface of the microsphere, and its coordinates can be described by the parametric equations of the sphere, $\{ {x_m} = {x_c} + r\sin \eta \cos \xi ,{y_m} = {y_c} + r\sin \eta \sin \xi ,{z_m} = {z_c} + r\cos \eta \} $, where r is the radius of the microsphere, parameters $\eta \in [0,\pi ]$ and $\xi \in [0,2\pi ]$. At point M, there are many plane waves passing this point. Figure 3(a) shows two sets of those plane waves, PW₁ and PW₂, which are perpendicular to xoz plane. For PW₁, only one FR (${L_1}$) travels across point M. Similarly, there is only one FR (${L_2}$) traveling across point M and belongs to PW₂ as well. Note that not all FRs can effectively strike on point M. For example, ${L_2}$ in Fig. 3(a) has struck on point S before it arrived at point M. In this case, ray ${L_2}$ cannot be taken into account for force calculation. To distinguish which rays can effectively strike on microsphere, we can use the following criterion [26]:

(2)$$rp\sin \eta \cos \xi + rq\sin \eta \sin \xi + r\sqrt {1 - {p^2} - {q^2}} \cos \eta < 0.$$

Only if the parameters η, ξ, p and q satisfy with Eq. (2), this FR of direction of ${\textbf v}({p,q} )$ is effective at point M. In physical insight, Eq. (2) implies that point M must be the first intersection point between ray and sphere when we trace the ray along its traveling direction.

Fig. 3. Ray tracing for Fourier rays. (a) Sketch of two Fourier rays, L₁ (blue line) and L₂ (red line) strike on the surface of the microsphere at point M. The blue dashed line and red dashed line are the wavefronts of two plane waves passing through point M. (b) Correlative vectors for calculating trapping forces. Incident point is M and the corresponding outward normal is unit vector N₁. Unit vectors e and N₂ indicate directions of ray polarization, and normal of the incident plane, respectively.

Download Full Size | PDF

Second, the optical force of an incident FR is then given according to the change of momentum. When a valid FR ${{\textbf L}_\textrm{I}}$ strikes on a tiny region (red curve) of point M with incident angle $\alpha$, as shown in Fig. 3(b), it will produce reflected (not shown) and refractive ray ${{\textbf L}_\textrm{R}}$ with refraction angle $\beta$. Due to the mismatch of the refractive indices between medium and microbead, the momentum of photon changes. Then the change of momentum generates optical force on the microsphere simultaneously. The optical force can be decomposed into a scattering component $\textrm{d}{{\textbf F}_s}$ which is parallel to the incident ray, and a gradient component $\textrm{d}{{\textbf F}_g}$ which is perpendicular to the incident ray. $\textrm{d}{{\textbf F}_s}$ and $\textrm{d}{{\textbf F}_g}$ are given by [13,26]

(3)$$\textrm{d}{{\textbf F}_s} = \frac{{{n_m}\textrm{d}P}}{c}{Q_s}{\hat{{\textbf F}}_s},$$

(4)$$\textrm{d}{{\textbf F}_g} = \frac{{{n_m}\textrm{d}P}}{c}{Q_g}{\hat{{\textbf F}}_g},$$

respectively, where c is the speed of light. dP is the light power of the FR. ${\hat{{\textbf F}}_s}$ and ${\hat{{\textbf F}}_g}$ are directional unit vector, which can be obtained from spatial analytic geometry [30] or quaternions [31] methods. ${Q_s}$ and ${Q_g}$ are the trapping efficiencies of scattering and gradient components, they are [13,26]

(5)$${Q_s} = 1 + R\cos ({2\alpha } )- {T^2}\left[ {\frac{{\cos ({2\alpha - 2\beta } )+ R\cos ({2\alpha } )}}{{1 + {R^2} + 2R\cos ({2\beta } )}}} \right],$$

(6)$${Q_g} = - R\sin ({2\alpha } )+ {T^2}\left[ {\frac{{\sin ({2\alpha - 2\beta } )+ R\sin ({2\alpha } )}}{{1 + {R^2} + 2R\cos ({2\beta } )}}} \right],$$

where R and T are the Fresnel reflection and transmission coefficients of energy flow. For a linear polarized ray they are separately given by [14]

(7)$$R = \frac{{{{\tan }^2}({\alpha - \beta } )}}{{{{\tan }^2}({\alpha + \beta } )}}{\sin ^2}\gamma + \frac{{{{\sin }^2}({\alpha - \beta } )}}{{{{\sin }^2}({\alpha + \beta } )}}{\cos ^2}\gamma ,$$

(8)$$T = \frac{{\sin ({2\alpha } )\sin ({2\beta } )}}{{{{\sin }^2}({\alpha + \beta } ){{\cos }^2}({\alpha - \beta } )}}{\sin ^2}\gamma + \frac{{\sin ({2\alpha } )\sin ({2\beta } )}}{{{{\sin }^2}({\alpha + \beta } )}}{\cos ^2}\gamma ,$$

where $R + T \equiv 1$, $\gamma$ is the angle between the polarization direction e of FR and the normal ${{\textbf N}_2}$ of the incident plane Ω, i.e., ${{\textbf N}_2}\ =\ {{\textbf L}_\textrm{I}} \times {\hat{{\textbf F}}_g}$. Particularly, for normal incident, there is no need to determine vector N₂ since the distinction between the parallel and perpendicular components disappears. In this situation, $R = {(n - 1)^2}/{(n + 1)^2}$ and $T = 1 - R$, where $n = {n_b}/{n_m}$ is the relative refractive index. In addition, for circularly polarized or unpolarized (scalar) ray, $\gamma \equiv \pi /4$. For simplicity, the FR are unpolarized in the following analysis.

The third step is to determine the quantity of dP in Eq. (3) and Eq. (4), which denotes the light power of each FR arriving at point M. Since FR is weighted by the Fourier angular spectrum, dP is the light power of corresponding plane wave $A(x,y,z;p,q)$ at point M when the FR propagates in direction ${\textbf v}({p,q} )$.

Assuming the total light power of the incident beam be ${P_0}$, and the light power at point M is $\textrm{d}{P_M}$, $\textrm{d}{P_M}$ is given as

(9)$$\textrm{d}{P_M} = {P_0} \cdot I({{x_M},{y_M},{z_M}} )\cdot |{\cos \eta } |\cdot \textrm{d}\sigma ,$$

where I is the power-normalized intensity at point M, and can be determined by Eq. (1). At point M the area of the tiny region has $\textrm{d}\sigma = {r^2}\sin \eta \textrm{d}\eta \textrm{d}\xi$. Let the light power of plane wave $A(x,y,z;p,q)$ be $\textrm{d}{P_{\textrm{PW}}}$. If the proportion of $\textrm{d}{P_{\textrm{PW}}}$ in ${P_0}$ is $\tau $, $\tau $ is given as

(10)$$\tau = \frac{{\textrm{d}{P_{\textrm{PW}}}}}{{{P_0}}} = C_0^2{|{\tilde{U}({p,q,0} )} |^2}/{\lambda ^2}$$

according to Eq. (1) and the Parseval's theorem. Notice that this proportion remains invariant for each plane wave at all points in the incident beam. Since the $\textrm{d}P$ is the light power of plane wave $A(x,y,z;p,q)$ at point M, it can be evaluated by replacing ${P_0}$ by $\tau {P_0}$ in Eq. (9). Hence, combining Eq. (9) and Eq. (10), $\textrm{d}P$ is given as

(11)$$\textrm{d}P = {P_M} \times \frac{{\textrm{d}{P_{\textrm{PW}}}}}{{{P_0}}} = \frac{{{P_0}C_0^2{r^2}}}{{{\lambda ^2}}}I({{x_M},{y_M},{z_M}} ){|{\tilde{U}({p,q,0} )} |^2}|{\cos \eta } |\sin \eta \textrm{ d}p\textrm{d}q\textrm{d}\eta \textrm{d}\xi .$$

Finally, the total optical forces on the trapped microsphere are the sum of scattering force $\textrm{d}{{\textbf F}_s}$ and gradient force $\textrm{d}{{\textbf F}_g}$ at all points on microsphere for effective FRs. In analyses of trapping force, we are more accustomed to use trapping efficiencies, which are represented by $\textrm{d}{{\textbf Q}_\chi } = c/({{n_m}{P_0}} )\textrm{d}{{\textbf F}_\chi },\ \chi = s,g$. According to Eq. (3) to Eq. (11), the total trapping efficiencies at any given position $c({x_c},{y_c},{z_c})$ are given by

(12)$$\begin{aligned}{{\textbf Q}_\chi }({{x_c},{y_c},{z_c}} ) &= \frac{{C_0^2{r^2}}}{{{\lambda ^2}}}\iint\limits_{{p^2} + {q^2} \le 1} \left[ {\iint\limits_\Sigma {{Q_\chi }I({{x_M},{y_M},{z_M}} )|{\cos \eta } |\sin \eta } } \right. \\ &\left. \times {{{|{\tilde{U}({p,q,0} )} |}^2}{{\hat{{\textbf F}}}_\chi }\textrm{d}\xi \textrm{d}\eta \ } \vphantom{\iint\limits_{{p^2}}}\right]\textrm{d}p\textrm{d}q,\ \chi = s,g. \end{aligned}$$

Notation Σ in Eq. (12) denotes the set of all points on the microsphere surface that satisfy the criterion in Eq. (2). The integral of Eq. (12) can be expressed by two steps. First, we integrate the force produced by the a given plane wave at all points (ξ, η) satisfying condition (2) on the surface of the sphere. Then we integrate on (p, q) and calculate the resultant force produced by all plane waves. In the following simulations, ${Q_x}$, ${Q_y}$, and ${Q_z}$ are used to denote the components of total trapping efficiency ${{\textbf Q}_s} + {{\textbf Q}_g}$ along x, y, and z directions, respectively. The calculation of Eq. (12) is a numerical integration. Hence, we should divide many grids for parameters ξ, η, as well as p and q. The number of the grids should be chosen appropriately in order to achieve the correct calculate and less time consuming. The standard for the choosing of number of grids is given in the Appendix.

4. Verification of FR method

In order to validate the FR method, simulation results from FR method are compared with the results from traditional RO model under the same parameters. Here, we take fundamental Gaussian beam and Airy beam as examples, since people are more familiar with these two kinds of beams.

4.1 Force of the fundamental Gaussian beam

Fundamental Gaussian beam is a typical beam for optical manipulation system. Lots of works have been done to systematically study its optical force in both experiment [15,17,32] and theory [10,13,16,17]. Hence, fundamental Gaussian beam is an appropriate example for us to validate the FR method. The traditional ray model of fundamental Gaussian beam is simply demonstrated in Fig. 1(a), and also detailed in Ref. [16]. The intensity distribution of fundamental Gaussian beam is [33]

(13)$$I(x,y,z) = C_0^2\frac{{\omega _0^2}}{{{\omega ^2}(z)}}\exp \left[ { - 2\frac{{({x^2} + {y^2})}}{{{\omega^2}(z)}}} \right],$$

where ${C_0} = {[2/(\pi \omega _0^2)]^{1/2}}$, ${\omega _0}$ is the radius of beam waist and satisfies the relation of $\omega (z) = {\omega _0}{(1 + {z^2}/z_r^2)^{1/2}}$, and ${z_r} = k\omega _0^2/2$ is the Rayleigh distance. According to the Fourier transform, the Fourier angular spectrum of fundamental Gaussian beam in $z = 0$ plane is

(14)$$\tilde{U}({p,q,0} )= \pi \omega _0^2\exp \left[ { - \frac{1}{4}\omega_0^2{k^2}({{p^2} + {q^2}} )} \right].$$

Substituting Eq. (13) and Eq. (14) into Eq. (12), we can calculate the optical forces for fundamental Gaussian beam by FR method.

Figure 4 shows the optical forces of particles under fundamental Gaussian beam. In the simulation, $\lambda = 1.064$ µm, ${n_0} = 1.33$ and ${n_b} = 1.59$. Figure 4(a) shows the axial trapping efficiency ${Q_z}$ for different sizes of microspheres. The radius of beam waist is ${\omega _0} = 2$ µm. Since the Gaussian beam is not tightly focused, the axial force pushes the microsphere toward positive z-axis. For radius $r = 2,\ 3,\ 6,$ and $8$ µm, the results from FR method are plotted in black full lines, while the results from the ray optics [16] are plotted by red dashed lines. The tendency and the value of ${Q_z}$ from FR method are accordance with the results using ray optics [16]. As for microsphere with $r = 3$ µm when the beam waist varies with ${\omega _0} = 2,\ 4,\ 6$ µm, ${Q_z}$ from FR method are still consistent with the results calculated by ray optics [16], as shown in Fig. 4(b).

Fig. 4. Optical forces of fundamental Gaussian beam. For different radius of microsphere r, (a) is the comparison of two groups of axial trapping efficiency Q_z that are calculated by the Fourier ray method and the method in Ref. [16], respectively. While (b) is the comparison of the results with different beam waist ω₀. Red dashed lines in both (a) and (b) are calculated by the method in Ref. [16], black dashed lines in (a) and (b) are calculated by the Fourier ray method. (c) and (d) are Q_z and Q_x calculated by Fourier ray method with different refractive indexes of microspheres, respectively.

Download Full Size | PDF

Setting ${\omega _0} = 2$ µm and $r = 3$ µm for microspheres with different refractive indices, Figs. 4(c) and 4(d) show calculation results of ${Q_z}$ and ${Q_x}$ by the FR method, respectively. Fourier ray method still works well with different refractive indices. Therefore, Fourier ray method is available for the simulation of optical force for fundamental Gaussian beam.

4.2 Force of the Airy beam

Airy beam is a type of self-accelerating beam which is used to optically guide particle to move along a parabolic trajectory. Study of optical force in Airy beam helps us understand the experimental phenomena, and design Airy beams with better features including higher trapping quality. Prior studies of the optical force were mainly based on EM model [12,34]. Recently, ray analysis for Airy beam has been reported [26]. In this subsection, we apply the FR method to Airy beam, and demonstrate that this method is also suitable for calculating the optical force of Airy beam.

The intensity distribution of a finite energy (2 + 1)D Airy beam is [25]

(15)$$\begin{aligned}I({x,y,z} ) &= C_0^2 \cdot {\left|{\textrm{Airy}\left( {\frac{x}{{{x_0}}} - \frac{{{z^2}}}{{4{k^2}x_0^4}} + i\frac{{az}}{{kx_0^2}}} \right)\textrm{Airy}\left( {\frac{y}{{{y_0}}} - \frac{{{z^2}}}{{4{k^2}y_0^4}} + i\frac{{bz}}{{ky_0^2}}} \right)} \right|^2}\\ \ &\times \exp \left( {\frac{{2ax}}{{{x_0}}} + \frac{{2by}}{{{y_0}}} - \frac{{a{z^2}}}{{{k^2}x_0^4}} - \frac{{b{z^2}}}{{{k^2}y_0^4}}} \right), \end{aligned}$$

where ${x_0}$ and ${y_0}$ are the length scale along x-axis and y-axis, respectively. a is non-negative decay factor along x-axis, while b is the decay factor along y-axis. The normalization coefficient of light power is ${C_0} = {[8\pi /({x_0}{y_0})\sqrt {ab} ]^{1/2}}\exp [ - 1/3({a^3} + {b^3})]$.The Fourier angular spectrum of the finite energy (2 + 1)D Airy beam at $z = 0$ plane is [25]

(16)$$\tilde{U}({p,q,0} )= {x_0}{y_0}\exp \left[ {\frac{1}{3}{{({a - ik{x_0}p} )}^3}} \right]\exp \left[ {\frac{1}{3}{{({b - ik{y_0}q} )}^3}} \right].$$

Substituting Eq. (15) and Eq. (16) into Eq. (12), we can calculate the optical force of Airy beam using FR method.

Figure 5(a) shows the transverse intensity pattern of Airy beam in $z = 0$ plane. Since line l₀ of $y = x$ is the symmetry axis of the transverse intensity pattern, we calculate the optical force of Airy beam in different z along line $y = x$, and compare the calculation results from FR method with the ray optics calculation [26]. In the simulations, ${n_b} = 1.46$ and $r = 2$ µm, $x = y = 2$ µm and $a = b = 0.1$. Figure 5(b) shows the tendency of ${Q_x}$ in $z = 0$, $90$ µm and $130$ µm. The results from FR method are plotted with black full lines, while the results from the method in Ref. [26] are plotted with red dashed lines. Although the magnitude of trapping efficiency from FR method is weaker than that from the method in Ref. [26], the whole tendencies of the ${Q_x}$ at different z are consistent with each other for both two methods. For ${Q_z}$ in Fig. 5(c), the simulations from two methods have similar tendency along line l₀ in the plane with different z-value.

Fig. 5. Optical forces of Airy beam. (a) is the normalized intensity pattern of Airy beam in z = 0 plane. (b) and (c) are comparison of Q_x and Q_z for the FR method and the ray optics [26], where r = 2 µm. (d) is the normalized intensity pattern in y = x plane. Green line is the caustic of Airy beam, and gray lines are light ray of the ray model [26]. (e) and (f) are comparison of Q_x for the FR method and ray optics [26] in the vicinity of the caustic. r = 5 µm in (e), while r = 1 µm in (f).

Download Full Size | PDF

In addition, the curves of trapping efficiency from FR method smoothly decrease to 0 at the neighbor of the beam’s caustic (for example, see curves of Q_x in Fig. 5(b) and Q_z in Fig. 5(c) at the neighbor of x = 0 when z = 0), while the curves from method in Ref. [26] vary dramatically. The main reason is that the traditional RO method fails at the vicinity of the caustic. Figure 5(d) shows the traditional ray model of the (2 + 1)D Airy beam in y = x plane. In this model, there is no ray (gray lines) at the right-hand side of the caustic (green curve). However, there are optical fields at that region due to the diffraction phenomenon. So that the traditional RO method cannot be used to calculate the optical force at that region. When the size of the particle is much larger than wavelength, this failure causes small calculation error which can be ignored as shown in Fig. 5(e) in the case of $r = 5$ µm. While the size of the particle is approach to the wavelength, this failure will lead considerable calculation error in the case of $r = 1$ µm, as shown in Fig. 5(f). Fortunately, the FR method avoids this failure, and can give valid way to calculate optical force produced by diffractive field.

The FR method gives reasonable simulation results of optical force under fundamental Gaussian beam and Airy beam, and is able to calculate the optical force that is produced by diffractive wave fields. In general, the FR method can be used to simulate the optical forces for arbitrary incident light beams in ray optics regime, since it relies on the Fourier angular spectrum of incident beam and any optical field can be expressed by Fourier angular spectrum superposition.

5. Validation of Fourier Ray method

In this section, we will take three specific structured light beams as examples to demonstrate that the FR method is appropriate to calculate optical forces, and gives reasonable prediction of the moving tendency of the microsphere emerged in these beams.

5.1 Force in the diffraction pattern of light ring

As mentioned in Subsection 4.2, FR method can be used to calculate the optical force produced by diffractive fields. To further demonstrate the applicability, we calculate the optical force of the diffraction field of a light ring. In the cylindrical coordinate system $(\rho ,\varphi ,z)$, where, $x = \rho \cos \varphi $ and $y = \rho \sin \varphi $, the field of the light ring is

(17)$$U({\rho ,\varphi ,0} )= {C_0}\exp \left[ { - \frac{{{{({\rho - {\rho_0}} )}^2}}}{{w_0^2}}} \right].$$

${w_0}$ in Eq. (17) is a positive real number that controls the wide of the light ring, ${r_0}$ is an arbitrary non-negative number. The normalization coefficient of light power is ${C_0} = {2^{ - 1/4}}{\pi ^{ - 3/4}}/\sqrt {{\rho _0}{w_0}}$. Figure 6(a) shows the intensity pattern of this light ring in $z = 0$ plane with ${\rho _0} = 20$ µm and ${w_0} = 0.5$ µm.

Fig. 6. Trapping forces of the light ring with ${\rho _0} = 20$ µm and ${w_0} = 0.5$ µm. The intensity profile at (a) $z = 0$ and (b) $z = 20$ µm, respectively. Inset in (b) shows the zoom-in of the central lobes. (c) Q_x for a trapped microsphere placed along x-direction at different propagation distances. (d) Q_z for a trapped microsphere placed along z-axis. Other parameters are r = 3 µm, n_m = 1.33, n_m = 1.59, and λ = 1.064 µm.

Download Full Size | PDF

It is known that the far-field Fraunhofer diffraction of the light ring forms Bessel-like beam [35], thus, we demonstrate that the Fourier ray method is able to calculate the force in diffractive fields. When the light ring propagates at $z = 20$ µm, the intensity pattern in the center of the beam is similar to Bessel beam with finite energy [Fig. 6(b)]. The Fourier spectrum of this light ring can be calculated numerically by applying Hankel transform to Eq. (17). In this simulation, we choose r = 3 µm, n_m = 1.33, n_m = 1.59, and λ = 1.064 µm. When the microsphere is placed along x-direction at $z = 20$, $25$ and $30$ µm planes, ${Q_x}$ is calculated using the Fourier ray method [Fig. 6(c)]. One stable trapping point appears at $x = 0$ where the microsphere can be stably trapped along x-direction. Due to the rotational symmetry of the diffraction pattern, the microsphere can be stably trapped in the transverse plane at different axial distance. The longitudinal trapping efficiency Q_z increases with distance [Fig. 6(d)], which implies that this beam always pushes the microsphere along positive z-axis.

In addition, due to the self-healing property, the major lobe of this beam will automatically reconstruct itself after partial blocking of the central part. We block (green line) the region of $\rho \le 5$ µm at $z = 40$ µm plane [Fig. 7(a)], and then calculated the Q_x along four lines which are l₁, l₂, l₃ and l₄ locating at $z = 45$ µm, $z = 50$ µm, $z = 60$ µm and $z = 70$ µm, respectively. When the microsphere is placed along l₁, since it locates in the shadow area, ${Q_x}$ is relatively weak, and so does ${Q_x}$ along l₂ [Fig. 7(b)]. When the microsphere is placed along l₃ or along l₄, where the beam recovers, the microsphere can be stably trapped at the center of the main lobe. In general, due to the self-healing property, this beam is able to stably trap many microspheres in multiple transverse planes with different axial position at the same time. Our simulation results are consistent with previous experiment [6].

Fig. 7. Effect of self-healing property on Q_x. (a) Intensity pattern of light ring in xoz plane. The region of $\rho \le 5$ µm in $z = 40$ µm plane is blocked. Due to the self-healing, the intensity lobe automatically reconstruct itself. (b) Q_x for a trapped microsphere placed along l₁, l₂, l₃ and l₄, respectively. Other parameters are r = 3 µm, n_m = 1.33, n_m = 1.59, and λ = 1.064 µm.

Download Full Size | PDF

5.2 Force in the Hermite-Gaussian beam

Hermit-Gaussian beam is a typical high-order transverse mode. Since it has an array-like intensity pattern with $m + 1$ column in x-direction and $n + 1$ row in y-direction, it becomes a good candidate for the optical tweezers array, which is suitable for trapping multiple microspheres simultaneously. In this subsection, we demonstrate the transverse optical force for a high-order Hermite-Gaussian beam using the Fourier ray method.

The intensity pattern of scalar Hermite-Gaussian beam at a point $(x,y,z)$ is [36]

(18)$${I_{mn}}({x,y,z} )= C_0^2\frac{{\omega _0^2}}{{{\omega ^2}(z )}}{{\mathop{\textrm {H}}\nolimits} _m}{\left[ {\frac{{\sqrt 2 x}}{{{\omega^2}(z )}}} \right]^2}{{\mathop{\textrm{H}}\nolimits} _n}{\left[ {\frac{{\sqrt 2 y}}{{{\omega^2}(z )}}} \right]^2}\exp \left[ { - \frac{{2({{x^2} + {y^2}} )}}{{{\omega^2}(z )}}} \right],$$

where m and n is the horizontal and vertical orders, ${\textrm{H}_m}$ is the mth-order Hermite polynomial, and ${C_0} = {(\pi {2^{m + n - 1}}m!n!\omega _0^2)^{ - 1/2}}$ [36]. Due to the self-similarity property of Hermite-Gaussian beam, the Fourier angular spectrum is also a Hermite-Gaussian like pattern which is [36]

(19)$${\tilde{U}_{mn}}({p,q,0} )= \pi \omega _0^2{{\mathop{\textrm {H}}\nolimits} _m}\left[ {\frac{{k{\omega_0}p}}{{\sqrt 2 }}} \right]{{\mathop{\textrm{H}}\nolimits} _n}\left[ {\frac{{k{\omega_0}q}}{{\sqrt 2 }}} \right]\exp \left[ { - \frac{1}{4}\omega_0^2{k^2}({{p^2} + {q^2}} )} \right].$$

Substituting Eq. (18) and Eq. (19) into Eq. (12), we can calculate the optical forces of Hermite-Gaussian beam using the FR method.

In following simulations, we choose ${n_b} = 1.59$, $r = 3$ µm, $m = 4$ and $n = 0$. The intensity pattern of this Hermite-Gaussian beam has five lobes on the x-axis [Fig. 8(a)]. The optical forces acting on the microsphere are calculated with different beam waist ${\omega _0}$. For ${\omega _0} = 2$ µm, there are 2 stable trapping positions on the x-axis [Fig. 8(b)]. When ${\omega _0} = 4.5$ µm, there are 4 positions where the beam can stably trap the microsphere in x-direction [Fig. 8(c)]. Since there are 5 lobes on x-axis, the maximum number of stable positions in x-direction is also 5 for ${\omega _0} = 6$ µm [Fig. 8(d)] and for ${\omega _0} = 9$ µm [Fig. 8(e)]. With waist ${\omega _0}$ increasing, the distance between arbitrary two adjacent intensity lobes also increase, and so does the number of the stable points.

Fig. 8. Q_x for Hermite-Gaussian beam with m = 4 and n = 0. (a) The intensity pattern in the xoy plane. The radii of beam waist are (b) ω₀ = 2 µm, (c) ω₀ = 4.5 µm, (d) ω₀ = 6 µm, and (e) ω₀ = 9 µm. Other parameters are r = 3 µm, n_m = 1.33, n_m = 1.59, and λ = 1.064 µm.

Download Full Size | PDF

When ${\omega _0} = 9$ µm, the stream lines (red lines) of the transverse optical forces $({Q_x},{Q_y})$ in xoy plane are shown in Fig. 9(a). The gray background indicated the normalized intensity distribution. The stream lines (red lines) of the transverse optical forces can be regarded as the moving tendency of trapped microsphere. Under the transverse optical forces, microsphere is dragged toward the beam region and then moved toward the center of 5 lobes. For high-order Hermite-Gaussian beam in two cases, the stream lines of transverse optical forces are separately plotted in Fig. 9(b) ($m = 5$ and $n = 3$) and Fig. 9(c) ($m = 2$ and $n = 6$). As shown in Fig. 9(b), the beam has $6 \times 4$ stable points in the center of each intensity lobe. This means that the beam can trap at least 24 microspheres transversely at the same time. In Fig. 9(c), this beam has $3 \times 7$ stable points, then it can trap at least 21 microspheres transversely at the same time.

Fig. 9. Stream lines of transverse optical forces for Hermite-Gaussian beams with (a) m = 4 and n = 0, (b) m = 5 and n = 3, and (c) m = 2 and n = 6. The beam waist is ω₀ = 9 µm. n_b= 1.59 and r = 3 µm.

Download Full Size | PDF

5.3 Force in a dual-arcs light beam

Complex structured light provides complex manipulation to the microscopic objects, and is potential colloidal assembly. In section 4 and subsections 5.1 and 5.2, we calculated optical forces for beams with pre-known initial field in spatial domain. However, FR method, in general, is suitable to special structured light beam with only close-form expression in its Fourier domain.

Here, we demonstrate the performance of optical force for a special light beam which only has the close-form expression of its Fourier angular spectrum

(20)$$\tilde{U}({p,q,0} )= \exp \left( { - \frac{{{p^2} + {q^2}}}{{{w_0}}}} \right)\exp \left\{ {i{k^3}w_0^3\cos \left[ {\arctan \left( {\frac{q}{p}} \right)} \right]} \right\},$$

where ${w_0}$ is an arbitrary positive real number. The distribution of optical field can be obtained by Eq. (1). The power normalization coefficient ${C_0} = \sqrt {2{\lambda ^2}/(\pi w_0^2)}$.

Figure 10(a) shows the intensity distribution of this beam with ${w_0} = 0.4$ µm. Due to the geometric shape of intensity profile, we call this structured beam as a dual-arcs beam. Figure 10(b) shows the transverse intensity distribution at $z = 70$ µm. Two main lobes locate along y-direction. Combining Fig. 10(a) and Fig. 10(b), the intensity of this beam for $z > 0$ mainly distributes in yoz plane. While $z < 0$, the intensity of this beam mainly distributes in xoz plane according to Fig. 10(a) and Fig. 10(c). The intensity distribution of this beam at $z = 0$ µm is shown in Fig. 10(d). It has a cross-like intensity lobe. In general, this beam is so-called quasi-nondiffracting optical field, and is self-accelerating.

Fig. 10. Intensity distributions of the dual-arc beam according to Eq. (20). (a) Spatial intensity distribution. (b), (c) and (d) are transverse intensity patterns in z = 70 µm, z = -70 µm, and z = 0 µm planes, respectively.

Download Full Size | PDF

To evaluate the optical force under the dual-arc beam, we calculated the optical forces in the part of $z < 0$. According to FR method, Eq. (20) is substituted into Eq. (1). The parameter of the microsphere is ${n_b} = 1.59$ and $r = 3$ µm. Figure 11 shows ${Q_x}$ and ${Q_z}$ at different height. For ${Q_x}$ in xoz plane, there are two main stable trapping positions along two main intensity lobes at three heights. Self-accelerating property makes the two main stable points gradually approach to each other with z increasing. Figure 11(b) shows that ${Q_z}$ has two main peaks and is always along positive z-axis, which can push the microsphere. Due to the self-acceleration, the distance between the two main peaks of ${Q_z}$ becomes closer with z increasing.

Fig. 11. Simulations of optical forces for dual-arcs beam. (a) and (b) are Q_x and Q_z, respectively.

Download Full Size | PDF

Finally, the stream lines (red lines) of the transverse optical forces ${\textbf Q} = ({Q_x},{Q_y})$ in xoy plane and ${\textbf Q} = ({Q_x},{Q_z})$ in xoz plane are shown in Fig. 12(a) and (b), respectively. The background gray patterns indicate the normalized intensity distribution in those planes. Optical forces with $|{\textbf Q} |< \textrm{|}{{\textbf Q}_{\max }}\textrm{|/20}$ are not shown since they are negligible. The stream lines (red lines) of the transverse optical forces denote the moving tendency of trapped microsphere. Figure 12(a) shows the stream lines in the neighbor of the main lobe in the right-hand side of the yoz plane. It is obvious that under the transverse optical forces, the microsphere is dragged toward the beam region and then moved toward the one of the main lobes of this beam. Meanwhile, the microsphere under the forces of ${{\textbf Q}_x} + {{\textbf Q}_z}$ is dragged toward the beam’s main lobe, and then moves along a trajectory of the distribution of the main lobe [Fig. 12(b)]. The same analyses can also be applied to optical forces for the left-hand side main lobe due to the symmetry of the intensity profile.

Fig. 12. Stream lines of optical forces. (a) in transverse plane at z = -70 µm. (b) in xoz plane. The background gray patterns denote normalized intensity distribution in corresponding planes. Optical forces with $\left| {\textbf Q} \right|<\left| {{\textbf Q}_{max}} \right|/20$ are not shown since they are negligible. The ring of attraction in (a) is caused by optical forces although the light intensity at those regions are relatively weak.

Download Full Size | PDF

6. Fourier rays for vectorial light beam

Although the previous simulations only involve the scalar optical field, it is believed that the FR method can be also applied to vectorial optical field. The crucial step is to obtain the Fourier angular spectrum. For instance, let ${E_x}(x,y,0)$ and ${E_y}(x,y,0)$ of the electric component E for an optical wave be initially known. Then the Fourier angular spectrum is given by Fourier transform [29]:

(21)$${\tilde{E}_{x,y}}({p,q,0} )= \int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {{E_{x,y}}({x,y,0} )\exp [{ - ik({px + qy} )} ]\textrm{d}x\textrm{d}y} } .$$

For a plane wave propagating along unit vector ${\textbf v} = (p,q,\sqrt[{}]{{1 - {p^2} - {q^2}}})$, the x and y components of the electric field $\tilde{{\textbf E}}$ are weighted by ${\tilde{E}_x}({p,q,0} )$ and ${\tilde{E}_y}({p,q,0} )$, respectively. According to the transversality condition, the z component of $\tilde{{\textbf E}}$ is [37]

(22)$${\tilde{E}_z} = - \frac{{{{\tilde{E}}_x}p + {{\tilde{E}}_y}q}}{{\sqrt {1 - {p^2} - {q^2}} }}.$$

Then the Fourier angular spectrums of the electric component are given by $\tilde{{\textbf E}} = ({\tilde{E}_x},{\tilde{E}_y},{\tilde{E}_z})$. Similarly, the Fourier angular spectrums of the magnetic component is $\tilde{{\textbf H}} = \sqrt {\varepsilon /\mu } {\textbf v} \times \tilde{{\textbf E}}$, where $\mu$ is the permeability and $\varepsilon$ is the permittivity of the medium. For one Fourier ray propagating along directional vector ${\textbf v}$, its light intensity can be derived from Eq. (21) and Eq. (22), and its polarization direction is ${\textbf e} = \tilde{{\textbf E}}/|\tilde{{\textbf E}}|$. Its coefficients R and T in Eq. (7) and Eq. (8) are then determined by the angle between ${\textbf e}$ and the normal of the incident plane. According to Eq. (11), light power of FRs at the point on the microsphere can be determined by the electromagnetic field of incident beam. Following this method, we believe that the FR method also can be used to calculate the optical force for vectorial filed, however, those are beyond the scope of this manuscript.

7. Concluding remarks

Combining the Fourier optics and ray optics, we proposed the Fourier ray (FR) method to calculate the optical force on particles in arbitrary light beam. The core idea of FR method is that arbitrary optical fields can be decomposed by the superposition of plane waves, and each plane wave can be further regarded as a family of parallel light rays that propagate in the same direction. For those parallel light rays, their intensities depend on the Fourier angular spectrum. By the virtue of Fourier analysis in optics, the FR method provides a different view for force calculation in the RO model, and can be used to evaluate the optical force generated by diffractive optical field, and even more complex structured light beams. Since any optical field has its Fourier angular spectrum, the FR method can be applied to calculate the optical force for arbitrary light beams.

Although the simulations demonstrated in this manuscript are only for spherical particle, we believe that the FR method can also be suitable for non-spherical objects like spheroid particles, since the only difference between FR method and traditional RO method is the calculation of light power for incident rays. The multiple reflecting of rays inside the particle can be calculated by using rotating matrix [38] or quaternions method [31].

Appendix: The appropriate number of grids for numerical integration of Eq. (12)

The calculation of Eq. (12) is a numerical integration. Hence, we should divide many grids for parameters ξ, η, p and q. The number of the grids should be chosen appropriately in order to achieve the correct calculate and less time consuming.

The number of grids for parameters ξ and η

Firstly, we determine the number of grids for parameters ξ and η. For simplicity, let N (a positive integer) be the number of grids for both ξ and η. We use the area of the surface of the microsphere as a criterion to choose the value of N. Let

(23)$${S_1} = {R^2}\sum\limits_{n = 1}^N {\sum\limits_{m = 1}^N {\sin {\eta _n}\varDelta \xi } } \varDelta \eta .$$

be the surface area of the microsphere calculated by the numerical integration, where $\varDelta \xi = 2\pi /N$ and $\varDelta \eta = \pi /N$. The total area is ${S_0} = 4\pi {R^2}$. If ${S_1}/{S_0} > 0.95$, it is considered that the value of N is sufficient for this numerical integration. The ratio of ${S_1}/{S_0}$ with respect to N is plotted in Fig. 13.

Fig. 13. Ratio of ${S_1}/{S_0}$ with respect to number of grids (N).

Download Full Size | PDF

As shown in Fig. 13, when $N \ge 50$, ${S_1}/{S_0} > 0.95$. In our simulation, there are 101 girds for parameter ξ in range of $(0 \le \xi \le 2\pi )$, while there are also 101 girds for parameter η in range of $(0 \le \eta \le \pi )$.

It is worth to note that this criterion cancels out the influence of the size of the sphere, however we are more interested in large microsphere since ray optics method only fits for large size of microsphere.

The number of grids for parameters p and q

Then, we determine the number of grids for parameters p and q. Similarly, let N (a positive integer) be the number of grids for both p and q. The total power of the incident beam is used as a criterion to choose the value of N.

According to the Parseval's theorem: the total light power in spatial domain equals to the power in Fourier spectrum domain, that is $\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{{|{U({x,y} )} |}^2}\textrm{d}x\textrm{d}y = } } \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{{|{\tilde{U}({p,q} )} |}^2}\textrm{/}{\lambda ^2}\textrm{d}p\textrm{d}q} }$. Let P₀ be the total power of the incident beam, and P₁ be the power calculated by

(24)$${P_1} = \frac{1}{{{\lambda ^2}}}\sum\limits_{n = 1}^N {\sum\limits_{m = 1}^N {{{|{\tilde{U}({p,q} )} |}^2}\varDelta p\varDelta q} } .$$

We define that the number of grids meeting ${P_1}/{P_0} > 0.95$ is sufficient for the numerical integration on p and q.

Since $\varDelta p$ and $\varDelta q$ denote the width of each grid for parameter p and q, respectively. If there are N girds for parameter p in ranges of $( - {p_l} \le p \le {p_r})$, $\varDelta p$ should be $\varDelta p = ({{p_r} + {p_l}} )/N$. Similarly, $\varDelta q$ should be $\varDelta q = ({q_r} + {q_l})/N$ if there are N girds for parameter q in ranges of $( - {q_l} \le q \le {q_r})$. Hence, our aim is to determine the appropriate integral range and number of grids for parameter p and q. Generally, we can simply set ${q_l} = {q_r} = {p_l} = {p_r} = 1$, since the evanescent wave (${p^2} + {q^2} > 1$) cannot propagate to far field. However, the values of parameters ${p_l}$, ${p_r}$, ${q_l}$, and ${q_r}$ should be chosen carefully for different types of incident beams in order to make the integral convergence precisely and quickly. Meanwhile, if the profile of ${|{\tilde{U}(p,q)} |^2}$ is rotational symmetry, the calculation of P₁ can be simplified as $2\pi \sum\nolimits_{n = 1}^N {|\tilde{U}(\rho ){|^2}} \textrm{/}{\lambda ^2}\rho \varDelta \rho$, where $p = \rho \cos \theta $, $q = \rho \sin \theta $, and $\theta = \arctan (q,p)$.

Taking the fundamental Gaussian beam as an example, the total light power of Gaussian beam can be calculated analytically as ${P_0} = \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\exp [{ - 2({{x^2} + {y^2}} )/\omega_0^2} ]\textrm{d}x\textrm{d}y = \pi \omega _0^2} } /2$. The Fourier angular spectrum of Gaussian beam is $\tilde{U}({p,q,0} )= \pi \omega _0^2\exp [ - \omega _0^2{k^2}({p^2} + {q^2})/4]$. We chose ${\omega _0} = 2$µm, $\lambda = 1.064$ µm, then the normalized distribution of ${|{\tilde{U}({p,q,0} )} |^2}$ is shown in Fig. 14(a).

Fig. 14. (a) Normalized distribution of ${|{\tilde{U}(p,q,0)} |^2}$. (b) Ratio of ${P_1}/{P_0}$ with respect to N.

Download Full Size | PDF

We find that ${|{\tilde{U}(p,q,0)} |^2} \to 0$ for ${p^2} + {q^2} > 0.4$, hence we can choose ${q_l} = {q_r} = {p_l} = {p_r} = 0.4$. The ratio of ${P_1}/{P_0}$ with respect to N is plotted in Fig. 14(b). As shown in Fig. 14(b), the red curve converges faster than the black curve. When $N \ge 50$, ${P_1}/{P_0} > 0.95$. In all of our simulations, there are 201 girds for each parameter p or q, however, the ranges of integral are different.

Funding

Natural Science Foundation of Anhui Province (1908085MA14); Scientific Research Foundation of the Institute for Translational Medicine of Anhui Province (2017zhyx25).

References

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

2. N. J. Carter and R. A. Cross, “Mechanics of the kinesin step,” Nature 435(7040), 308–312 (2005). [CrossRef]

3. M. C. Williams and I. Rouzina, “Force spectroscopy of single DNA and RNA molecules,” Curr. Opin. Struct. Biol. 12(3), 330–336 (2002). [CrossRef]

4. D. J. Stevenson, F. Gunn-Moore, and K. Dholakia, “Light forces the pace: optical manipulation for biophotonics,” J Biomed Opt 15(4), 041503 (2010). [CrossRef]

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

6. V. Garcés-Chávez, D. McGloin, H. Melville, W. Sibbett, and K. Dholakia, “Simultaneous micromanipulation in multiple planes using a self-reconstructing light beam,” Nature 419(6903), 145–147 (2002). [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, 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]

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

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

11. S.-Y. Sung and Y.-G. Lee, “Trapping of a micro-bubble by non-paraxial Gaussian beam: computation using the FDTD method,” Opt. Express 16(5), 3463–3473 (2008). [CrossRef]

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

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

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

15. T. C. B. Schut, G. Hesselink, B. G. Degrooth, and J. Greve, “Experimental and theoretical investigations on the validity of the geometrical-optics model for calculating the stability of optical traps,” Cytometry 12(6), 479–485 (1991). [CrossRef]

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

17. S. Nemoto and H. Togo, “Axial force acting on a dielectric sphere in a focused laser beam,” Appl. Opt. 37(27), 6386–6394 (1998). [CrossRef]

18. Y.-R. Chang, L. Hsu, and S. Chi, “Optical trapping of a spherically symmetric sphere in the ray-optics regime: a model for optical tweezers upon cells,” Appl. Opt. 45(16), 3885–3892 (2006). [CrossRef]

19. A. Callegari, M. Mijalkov, A. B. Gököz, and G. Volpe, “Computational toolbox for optical tweezers in geometrical optics,” J. Opt. Soc. Am. B 32(5), B11–B19 (2015). [CrossRef]

20. J. Liu, C. Zhang, Y. W. Zong, H. L. Guo, and Z.-Y. Li, “Ray-optics model for optical force and torque on a spherical metal-coated Janus microparticle,” Photonics Res. 3(5), 265–274 (2015). [CrossRef]

21. A. Devi and A. K. De, “Theoretical investigation on optical Kerr effect in femtosecond laser trapping of dielectric microspheres,” J. Opt. 19(6), 065504 (2017). [CrossRef]

22. G. Volpe, L. Kurz, A. Callegari, G. Volpe, and S. Gigan, “Speckle optical tweezers: micromanipulation with random light fields,” Opt. Express 22(15), 18159–18167 (2014). [CrossRef]

23. S. H. Zhang, J. H. Zhou, and L. Gong, “Skew line ray model of nonparaxial Gaussian beam,” Opt. Express 26(3), 3381–3393 (2018). [CrossRef]

24. G. W. Forbes and M. A. Alonso, “What on earth is a ray and how can we use them best?” Proc. SPIE 3482, 22–31 (1998). [CrossRef]

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

26. S. H. Zhang, J. H. Zhou, and Y.-X. Ren, “Ray optics analysis of optical forces on a microsphere in a (2 + 1)D Airy beam,” OSA Continuum 2(2), 378–388 (2019). [CrossRef]

27. M. A. Alonso and M. R. Dennis, “Ray-optical Poincaré sphere for structured Gaussian beams,” Optica 4(4), 476–486 (2017). [CrossRef]

28. G. Roosen, B. F. de Saint Louvent, and S. Slansky, “Etude de la pression de radiation exercee sur une sphere creuse transparente par un faisceau cylindrique,” Opt. Commun. 24(1), 116–120 (1978). [CrossRef]

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

30. J. H. Zhou, H. L. Ren, C. Jin, 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]

31. S. H. Zhang, Z. Liang, and J. H. Zhou, “Using quaternions to analyze the trapping force of an ellipsoidal bead,” Acta. Phys. Sin. 66(4), 048701 (2017). [CrossRef]

32. A. Ashkin, “Acceleration and trapping of particles by radiation pressure,” Phys. Rev. Lett. 24(4), 156–159 (1970). [CrossRef]

33. M. Dietrich, Light transmission optics (Van Nostrand Reinhold, 1972).

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

35. J. C. Gutiérrez-Vega and M. A. Bandres, “Helmholtz–Gauss waves,” J. Opt. Soc. Am. A 22(2), 289–298 (2005). [CrossRef]

36. S. Feng and H. G. Winful, “Physical origin of the Gouy phase shift,” Opt. Lett. 26(8), 485–487 (2001). [CrossRef]

37. B. Andreas, G. Mana, and C. Palmisano, “Vectorial ray-based diffraction integral,” J. Opt. Soc. Am. A 32(8), 1403 (2015). [CrossRef]

38. J. H. Zhou, M.-C. Zhong, Z.-Q. Wang, and Y. M. Li, “Calculation of optical forces on an ellipsoid using vectorial ray tracing method,” Opt. Express 20(14), 14928–14937 (2012). [CrossRef]

Calculation of optical forces for arbitrary light beams using the Fourier ray method

Abstract

1. Introduction

2. Fourier rays

3. Optical force of Fourier ray family

4. Verification of FR method

4.1 Force of the fundamental Gaussian beam

4.2 Force of the Airy beam

5. Validation of Fourier Ray method

5.1 Force in the diffraction pattern of light ring

5.2 Force in the Hermite-Gaussian beam

5.3 Force in a dual-arcs light beam

6. Fourier rays for vectorial light beam

7. Concluding remarks

Appendix: The appropriate number of grids for numerical integration of Eq. (12)

The number of grids for parameters ξ and η

The number of grids for parameters p and q

Funding

References

Cited By

Figures (14)

Equations (24)

Optics Express