Helical tractor beam : analytical solution of Rayleigh particle dynamics

We analyze particle dynamics in an optical force field generated by helical tractor beams obtained by the interference of a cylindrical beam with a topological charge and a co-propagating temporally de-phased plane wave. We show that, for standard experimental conditions, it is possible to obtain analytical solutions for the trajectories of particles in such force field by using of some approximations. These solutions show that, in contrast to other tractor beams described before, the intensity becomes a key parameter for the control of particle trajectories. Therefore, by tuning the intensity value the particle can describe helical trajectories upstream and downstream, a circular trajectory in a fixed plane, or a linear displacement in the propagation direction. The approximated analytical solutions show good agreement to the corresponding numerical solutions of the exact dynamical differential equations. © 2015 Optical Society of America OCIS codes: (260.2110) Electromagnetic optics; (260.5430) Polarization; (140.7010) Laser trapping; (350.4855) Optical tweezers or optical manipulation. References and links 1. A. Ashkin, “Acceleration and trapping of particles by radiation presure,” Phys. Rev. Lett. 24, 156–159 (1970). 2. T. Cizmar, V. Garces-Chavez, K. Dholakia, and P. Zemanek, “Optical conveyor belt for delivery submicron objects,” Appl. Phys. Lett. 86, 174101 (2005). 3. T. Cizmar, M. Siler, and P. Zemanek, “An optical nanotrap array movable over a milimetre range,” Appl. Phys. B 84, 197–203 (2006). 4. S.-H. Lee, Y. Roichman, and D. G. Grier, “Optical solenoid beams,” Opt. Express 18(7), 6988–6993 (2010). 5. D. B. Ruffner and D. G. Grier, “Optical conveyors: A class of active tractor beams,” Phys. Rev. Lett. 109, 163903 (2012). 6. O. Brzobohaty, V. Karasek, M. Sailer, L. Chavatal, T. Cizmar, and P. Zemanek, “Experimental demonstration of optical transport sorting and sef-arragement using a tractor beam,” Nat. Photonics 7, 123–127 (2013). 7. V. Shvedov, A. R. Davoyan, C. Hnatovsky, N. Engheta, and W. Krolikowski, “A long-range polarizationcontrolled optical tractor beam,” Nat. Photonics 8(11), 846–850 (2014). 8. Huajin Chen, Shinyang Liu, Jiang Zi, and Zhifang Lin, “Fano resonance-induced negative optical scattering force on plasmonic nanoparticles,” ACS Nano 9(2), 1926–1935 (2015). 9. S. Sukhov and A. Dogariu, “Negative nonconservative forces: optical ‘tractor beams’ for arbitrary objects,” Phys. Rev. Lett. 107, 203602 (2011). URL http://link.aps.org/doi/10.1103/PhysRevLett.107.203602. 10. L. Carretero, P. Acebal, and S. Blaya, “Three-dimensional analysis of optical forces generated by an active tractor beam using radial polarization,” Opt. Express 22, 3284–3295 (2014). 11. J. Durnin, “Exact solutions for nondiffracting beams. I. The scalar theory,” J. Opt. Soc. Am. A 4(4), 651–654 (1987). 12. H. Chen, N. Wang, W. Lu, S. Liu, and Z. Lin, “Tailoring azimuthal optical force on lossy chiral particles in Bessel beams,” Phys. Rev. A 90, 043850 (2014). URL http://link.aps.org/doi/10.1103/PhysRevA.90.043850. #242346 Received 4 Jun 2015; revised 9 Jul 2015; accepted 10 Jul 2015; published 28 Jul 2015 © 2015 OSA 10 Aug 2015 | Vol. 23, No. 16 | DOI:10.1364/OE.23.020529 | OPTICS EXPRESS 20529 13. L. Carretero, P. Acebal, C. Garcia, and S. Blaya, “Periodic trajectories obtained with an active tractor beam using azimuthal polarization: design of particle exchanger,” IEEE Photon. J. 7(1), 1–12 (2015). 14. C. Vetter, T. Eichelkraut, M. Ornigotti, and A. Szameit, “Generalized radially self-accelerating helicon beams,” Phys. Rev. Lett. 113, 183901 (2014). 15. P. C. Chaumet and M. Nieto-Vesperinas, “Time-averaged total force on a dipolar sphere in an electromagentic fiel,” Opt. Lett. 25, 1065–1067 (2000). 16. M. Nieto-Vesperinas, J. S??enz, R. Gomez-Medina, and L. Chantada, “Optical forces on small magnetodielectric particles,” Opt. Express 18, 11428–11443 (2010). 17. N. Wang, J. Chen, S. Liu, and Z. Lin, “Dynamical and phase-diagram study on stable optical pulling force in Bessel beams,” Phys. Rev. A 87, 063812 (2013).


Introduction
Optical trapping of particles [1] has become an essential tool for the manipulation of small objects in many disciplines such as physics or biology.In this way, new advances in laser beam control have led to the experimental realization of tractor beams that can draw particles toward light sources [2][3][4][5][6][7].Recently, it was pointed out that Fano resonance of a composite metallic nanoparticle can induce a negative pulling force [8].Unlike conventional optical tweezers, which rely on adjusting the positions of intensity gradients to move objects, tractor beams can exert non-conservative pulling forces on particles by continuously dragging them towards the beam sources [9].In these tractor beams, the dynamics of particles has also been analyzed; in particular, our team has theoretically studied the case of radially polarized zero-order Bessel tractors beam [10], analyzing the trajectories that, in general, are linear [2,3], [5].As found from many theoretical and experimental studies, non-diffracting beams are promising candidates for realizing tractor beams due to their unique properties of maintaining both intensity and spatial extent in the direction of propagation [11].Thus, a great variety of optical fields, ranging from the fundamental Gaussian beam, Laguerre-Gaussian beam, Airy beam, and Bessel beam (BB) to the plasmon-based optical field, has been employed to achieve optical micromanipulations [12] .In this sense, more complex trajectories can be obtained using, for example, azimuthally polarized beams [13].Moreover, recently, an interesting class of beam that could be used for sorting, mixing, or cell extraction applications could be obtained using generalized radially self-accelerating helicon beams [14] due to the three-dimensional spiraling trajectories that could be generated.In fact, a similar procedure was used in [4] for obtaining solenoidal waves through a particular superposition of many m-th order Bessel beams that do not differ in their relative temporal phase.In relation to these points, in this paper we propose to generate a helical tractor beam by means of the interference of a plane wave and a cylindrical beam (that differ in their relative temporal phase) and it is a paraxial solution to the Helmholtz equation.For this, we are going to analytically and numerically analyze the dynamics of Rayleigh particles inside this tractor beam.

Theoretical background
We consider a linearly polarized electric field E, obtained by a superposition of a plane wave (PW) and a paraxial solution of the Helmholtz equation in (r,ψ,z) cylindrical coordinates (i.e., a cylindrical beam (CB) like Bessel beams or Laguerre-Gaussian beams, etc.) propagating in the ẑ direction given by: where k 1 and k 2 are the wave numbers of the beams in a medium of refractive index n, G(r) is the radial solution to the Helmholtz equation in cylindrical coordinates, A is the electric field amplitude, 0 < η < 1 is the beam ratio between the plane wave and CB beam.The PW and the CB beams differ in their relative phases not only in the helical phase given by the topological charge m, but also in a temporal linear function ξ t that is critical for tractor working [5].
It should be noted, that the solenoidal beam described by Eq.( 1) operates due to the relative phase difference between the PW and CB beams (ξ t), and it is quite different to that showed in [4], because the second one uses periodic axial intensity gradients with discrete propagation invariance to achieve forward scattering from the interference between the incident field and the dipole radiation field of an illuminated object [5].
The optical forces acting on a particle in the Rayleigh regime (particle radius r p << λ ) as a consequence of an electric field given by Eq. ( 1) can be expressed according to [15], as: where ℜ denotes the real part of the expression, * is the complex conjugate value, and α is the complex polarizability of the particle: α R and α I are the real and imaginary parts of the polarizability, respectively, which for the dielectric particles can be obtained, for example, from [16].
The deterministic trajectories of particles can be obtained by solving the equation: where R(t) is the position vector of the particle at time t, m 0 is the particle mass, and −κ d R dt is the frictional force.It is important to remark that Eq. ( 4) does not take into account the thermal effects, that could be included by using the Langevin equation instead [10].

Approximated solution
Assuming that the inertia force is smaller than the drag force (over-damped case), the left-hand side of Eq.(4) vanishes, so introducing Eqs.(1,2) and 3 in Eq. (4), motion particle equations can be written as: where ṗ, (p = r, ψ, z) represents the temporal derivative of the p-coordinate, the function describes the temporal coupling of the axial and azimuthal variables, and the constant matrix M is defined as: where we have used the parameter Ω = A 2 /(2κ) = I 0 /(κ n c ε 0 ) (I 0 being the intensity of the field I 0 = n c ε 0 A 2 /2 ).The motion equation for the radial coordinate is a typical equation for an optical trap, where the particles are quickly trapped in the radial positions r q that satisfy, independently of their axial and azimuthal positions, the equation Ġ(r)| r=r q = 0. Therefore, at these points r q , matrix elements M 11 , M 12 and M 13 are null, so Eq. ( 5) can be rewritten as: where we have introduced the parameters: Parameters v z0 and v ψ0 represent the axial and angular velocity that the particles would obtain in absence of interference, namely, if the electric fields in Eq.( 1) were incoherent, while v z and v ψ correspond to the axial and angular velocity due to the interference, respectively.In order to solve the general case given by Eqs.(8,9), it is necessary to uncouple the ψ(t) and z(t) variables.To do this, removing function H(t) from Eqs. (8,9) obtains that axial velocity ż and angular velocity ψ fulfill the following ellipse equation: where δ = α z − α ψ .Using the definitions of α z and α ψ given in Eq. ( 10), the key parameters Cos(δ ) and Sin(δ ) = [1 −Cos(δ ) 2 ] 1/2 can be written as: We center our study on the case of Rayleigh dielectric particles so we can assume that α R >> α I [3] and then, using Eq. ( 12), we can approximate: where we have assumed that the wave number of PW is higher than the one of CB (k 1 > k 2 ).For Rayleigh particles, |Sin(δ )| << 1.Therefore, Eq.( 11) can be simplified to: Thus, by solving differential equation 14, the temporal relation between axial and azimuthal coordinates is given by: where: with ψ 0 and z 0 the initial position coordinates of the particle.Therefore, by using the temporal relation between the axial and azimuthal coordinates (Eq.( 15)), it is possible to express H(t) only as a function of the time and the z(t) or ψ(t) coordinate, so differential equations 8-9 can be written as: Previous equations show that trapped particles are axially and azimuthally translated with tractor velocities: Analytical solutions to Eqs.(19) are given by: where f z and f ψ are the functions: ) c 1 and c 2 are constants that take boundary conditions into account.Moreover, it must be noted that f z and f ψ are bounded functions that verify f z < v zt t and f ψ < v ψt t, so Eqs.(22) can be approximated to: As a consequence, the position's particle vector R(t) can be written as: where p = 2πv zt v ψt .Equation 26shows that particles describe a helical motion where ẑ is the axis of the helix, r q is the radius, p is the pitch and v ψt is the angular frequency of the particle vector.Then, R(t) traces a helix whose end point moves around a cylinder once every 2π v ψt units of time and moves a distance p in the direction ẑ for every revolution around the cylinder.

Tractor velocity intensity dependence
In this section, we are going to assume that m ≥ 0. Parameters A 1 and A 2 play an important role in the movement of particles, so it is important to analyze their relative values.Taking into account the definitions given by Eqs.(10), it can be deduced that velocities v ψ0 , v ψ , v z0 , v z are greater than or equal to zero, and as a result, A 1 ≥ 0 and A 2 ≤ 0. It is important to remark that the A 1 parameter depends on the position of the maximum r q where the particle is located, and due to the dependence of velocities v ψ0 and v ψ on Ω, A 1 is proportional to the intensity.However, A 2 depends only on r q , while it does not on I 0 .According to this, we can write A 1 = a 1 (r q ) I 0 and A 2 = a 2 (r q ), with (a 1 (r q ) ≥ 0 and a 2 (r q ) ≤ 0).Taking into account this fact, the tractor velocities and pitch can be rewritten as: Previous equations show that pitch, axial, and azimuthal tractor velocities (and as a consequence trajectories given in Eq.( 26)) depend on the intensity, the topological charge, the local maximum where the particle is located, the relative phase difference between the PW and CB beam, and of course, the amplitude ratio η between the beams.It can be noted that in the limit case in which η = 0, there is no action of the tractor beam (ξ = 0) and the particle dynamic is controlled only by the Bessel beam.The stability of the negative pulling force for a particle in a Bessel beam has been discussed in the reference [17].In the other cases, the influence of amplitude ratio η in the particle movement is included in the a 1 parameter in a complex way.In this sense, a 1 shows a parabolic dependence on η, but the parabola's coefficients change basically with radial position, topological charge, and intensity.
If there is no topological charge (m = 0), then v ψ and v ψ0 are null, as a consequence a 1 (r q ) = a 2 (r q ) = 0, and then from Eq.( 27), it can be deduced that azimuthally tractor velocity v ψt is zero and the axial tractor velocity is equal to that obtained in reference [5] ).In this case, all particles are translated upstream or downstream (depending on the sign of the ξ parameter), and the helix trajectory given by Eq. ( 26) is reduced to a linear movement in the ẑ direction [5].Therefore, the lack of topological charge implies that axial tractor velocity does not depend on intensity I 0 .In the case of null relative phase difference ξ = 0 and non-null topological charge m = 0, it can be also deduced from Eqs.(27,26) that a helical trajectory is obtained with the tractor velocities and pitch: where, under these conditions, v zt > 0; then, particles move upstream describing a helix in the positive ẑ direction with velocities that are proportional to the intensity I 0 and vary for every r q position.It is interesting to point out that particles can not be moved downstream for any m parameter value.Furthermore, in this case, the pitch is constant and independent of the intensity and the radial position of particles r q .
The most interesting cases are obtained when the relative phase difference ξ and the topological charge m are both non-null.Eqs. ( 27) show that v zt can be positive, negative, or null, and v ψt is ever positive or null depending on the intensity I 0 and the position of the intensity maximum r q .Thus, the behavior of particle movement could be very different by adequately choosing the intensity and ξ ratio, so we can obtain: 1.All trapped particles describe a helix upstream.
In this case, it is necessary that v zt > 0, so, according to Eq.( 20): If ξ is positive, independently of the position trap r q , this inequality is verified for all particles.For negative values of the ξ parameter, this inequality is fulfilled for particles located at position traps r q that satisfy the condition |ξ | < m a 1 (r q )I 0 , so negative values of ξ can generate upstream movement in the helical tractor beams.
2. All trapped particles describe a helix downstream.
In this case, it is necessary that v zt < 0, so, according to Eq.( 20): This inequality is verified only if ξ is negative and |ξ | > m a 1 (r q )I 0 .
3. All trapped particles describe a helix except the corresponding ones located at the intensity maxima position r q = r qψ that describe a circular movement at the initial z 0 plane.
In this case, it is necessary that v zt (r qψ ) = 0, so, according to Eq.( 20): Then, intensity I qψ selects particles located at the radial trap r qψ and moves them in a circular trajectory at the initial z 0 plane, while all the other ones describe helical trajectories.Then, it is possible to control the type of particle' movement by using the adequate intensity.Note that the intensity I qψ is positive only if ξ is negative, so, in general, in this case, the helical movement will be downstream.
4. All trapped particles describe a helix upstream except the ones located at intensity maxima position r q = r qz that show a linear movement in the ẑ direction.
In this case, it is necessary that v zt > 0 for all particles (so according to the discussion of point 1, we choose ξ > 0) and v ψt (r qz ) = 0, so, from Eq.(20): Then, intensity I qz selects particles located at the radial trap r qz and moves them linearly in the ẑ direction, while all the other ones, describe upstream helical trajectories.Then, by using the right intensity, linear movement along the ẑ direction and an upstream helix can be chosen.Note that the intensity I qz is positive only when ξ is positive, so it is not possible to obtain a downstream helix movement.

Comparison between exact numerical results and approximations
In this section, we are going to compare our approximated solutions described in previous sections and the corresponding exact numerical solution to Eq.( 4).solving coupled differential equations we have used a Runge-Kutta fourth-order method using the same boundary conditions in numerical calculations as in approximated solutions (with initial velocity particles equal to 0).We have also used a Bessel beam as the CB beam in Eq.( 1), then we take k 2 = k 1 Cos(θ ) and the radial function G(r) = J m (k 1 Sin(θ )r), with J m the Bessel function of the first kind of order m.Moreover, the constants κ = 6πη 1 r p and m 0 = 4 3 ρπr 3 p (with r p the particle radius) have been obtained using the water viscosity coefficient η 1 = 8.9 × 10 −4 (Pa s) and the material density (PMMA) ρ = 1.19 × 10 3 (Kg m −3 ).The other numerical values used in our study are given in table 1. Due to the large number of possible cases, we will focus on the results described by Eqs.(31,32).Figure 1 shows the axial length, axial velocity, angular velocity and pitch obtained at a fixed time t = 5 s for ξ = −10 Hz as a function of the intensity position maxima r q using three different topological charges (m = 1, 2, 3) by accomplishing v zt = 0 (Eq.(31)) at the first maximum position r 1 .Under this condition, the obtained intensities are I 1ψ = 0.71, 0.6, 0.56W /cm 2 for topological charges m = 1, 2, 3, respectively.As can be observed, good agreement is obtained between approximated analytical equations and exact numerical results.As was expected, particles located at the r 1 maxima position do not move axially, but describe a circular movement of angular velocity ψ = v ψt ≈ −ξ m .It can also be observed that, for a fixed topological charge, axial displacement increases as r q rises, although when the topological charge increases, axial length is lower for equivalent maxima r q .3D helical trajectories can be observed in Fig. 2 for the particular cases of a topological charge m = 2 and condition v zt = 0 obtained with an intensity I 0 = I 1ψ = 0.6W /cm 2 .It can be seen that particles located at maxima r 2 , r 3 and r 4 describe a downstream helix, with the numerical result very close to the obtained one by means of analytical Eqs.(26).As was previously mentioned, the particles located at the maximum r 1 describe a circular trajectory.Finally, we are going to realize the corresponding study when v ψt = 0 (Eq.(32)).Figure 3 shows the axial length, axial velocity, angular velocity, and pitch obtained at a fixed time t = 5 s for ξ = 10 Hz as a function of the intensity position maxima r q using three different topological charges (m = 1, 2, 3) with the condition given by Eq.(32) v ψt = 0 for the first maximum position r 1 .Under this condition, the intensities are I 1z = 4.9, 6.1, 6.8W /cm 2 for the topological charges m = 1, 2, 3 respectively.As can be observed, good agreement between approximated analytical equations and exact numerical results is shown, but small differences can be noted for the topological charge m = 1 case.From Fig. 3, as was mentioned, particles located at the r 1 maxima describe a linear movement on the z-axis because the angular velocity ψ = v ψt = 0. 3D helical trajectories can be observed in Fig. 4 for the particular cases of a topological charge m = 2, and the condition v ψt = 0 obtained by an intensity I 0 = I 1ψ = 6.1W /cm 2 .It can be seen that particles located at the maxima r 2 , r 3 , and r 4 describe an upstream helix, with the numerical result very close to that obtained by means of analytical Eqs.(26).

Table 1 .
Values of the parameters used in the numerical simulations