Open Access
Add to BibSonomyAbstract
The computation details related to computing the optical radiation pressure force on various objects using a 2-D grid FDTD algorithm are presented. The technique is based on propagating the electric and magnetic fields through the grid and determining the changes in the optical energy flow with and without the trap object(s) in the system. Traces displayed indicate that the optical forces and FDTD predicted object behavior are in agreement with published experiments and also determined through other computation techniques. We show computation results for a high and low dielectric disc and thin walled shell. The FDTD technique for computing the light-particle force interaction may be employed in all regimes relating particle dimensions to source wavelength. The algorithm presented here can be easily extended to 3-D and include torque computation algorithms, thus providing a highly flexible and universally useable computation engine.
©2005 Optical Society of America
1. Introduction
In the early part of 1970 A. Ashkin experimentally demonstrated that optical radiation pressure could be used to suspend and manipulate a small dielectric object against a gravitational force [1]. He labeled this type of light-particle interaction optical levitation. This lead to the study of optical levitation in the research laboratory and in the theoretical modeling of the interaction. In 1986, this same researcher introduced optical trapping of a dielectric object in a highly focused region of a laser beam [2]. The gradient force pulling the object towards the focus having been made stronger than the scattering force due to the high focus, resulted in the particle being trapped using light only and no longer required the gravitational counterpart to balance forces. This type of light-particle interaction was labeled optical trapping. In the basic optical trapping experimental configuration the trapped object will follow the beam focus and as such the object may be translated in 3-D using the non-contact optical force. In several adaptations of the basic trapping system an optical torque can be generated through either an alteration of the Gaussian beam or through the shape and composition of the object [3,4]. The optical generation of force and torque on objects has lead to the trapping system being used to trap and activate micro-optic gears [5,6], align and sort cells[7,8], study of micro-system forces [9], and many other applications in all branches of the sciences.
Considerable effort has been directed towards the computation and simulation of the optical levitation and trapping experimental configurations. In general the objects under study may range from a small fraction of the wavelength to several tens of wavelengths, passing through the range where the object is comparable to the wavelength. In the two extremes of this size range quite adequate and continuously evolving optical models exist to compute the light-particle interaction and predict the dynamic behavior of the object [10–17]. In this paper we develop and present an optical force computation technique based on a 2-D Finite-Difference-Time-Domain (FDTD) algorithm. The technique is applied to computing the force on spherical 2-D objects in one or more optical beams.
2. FDTD technique
The present FDTD technique for solving Maxwell’s time dependent equations can be traced back to Yee in 1966[18]. The technique involves expanding Maxwell’s vector curl equations and isolating the field component mixed scalar equations that result. Considerations as to the nature of the media involved, the source excitation and dimensionality of the space are usually made at this stage in order to simplify as much as possible the resultant expressions. The computation domain is then discretized into a fine mesh of points and finite difference expression developed for each of the scalar equations. A stability criterion is applied which sets an upper bound between the grid spacing of the mesh and the time step interval. In order to minimize the effects of reflections at the boundary of the mesh grid a perfectly matched layer (PML) is usually defined and encloses the entire computation domain [19,20]. Through a leapfrog arrangement the time updated values for the electric and magnetic fields are computed throughout the grid. Given the nature of the source (plane wave, Gaussian beam, …) and objects in the grid domain, the light-particle interaction can be computed as a function of time. The resulting E and H field grid information can be analyzed in order to extract the radiation pressure generated force on the grid domain objects.
Fig. 1. (X, Y) computation domain showing a single sphere axially offset and displaced in the divergent region of the source plane. The computation domain is discretized into an (I, J) grid. The boundary is composed of a perfectly matched layer 15 grid points wide.
Figure 1 shows the 2-D computation domain that is enclosed by a PML of 15 grid points thick. The X-axis runs from the top downwards and the Y-axis runs from left to right, Z-direction is out of page. A spherical object of relative dielectric εs, index of refraction ns=(εs)½, and radius rs is shown centered on the point (dx, dy). The straight line on the left side indicates the optical source that is characterized by its wavelength λo, waist Wo, amplitude profile Axy and propagation direction. The computation domain is discretized into a fine (X, Y) mesh. In order to assure sufficient resolution of the optical system the grid should posses at least 10 grid points per wavelength in any direction and at the same time provide at least 40 grid points with respect to the smallest object dimension. The fine grid ensures that the stability criterion is met and that the objects have a low level of “stair casing”. The efficient computerization of the FDTD equations permits large grid areas to be discretized and analyzed in a reasonable computational time and memory requirement. In 2-D, the polarization of the source can be selected as either TM with (Ez, Hx, Hy) components, or, TE with (Hz, Ex, Ey) components. The computation results presented are for a TM polarization source since extending the computations to include TE introduces no additional computation complexities.
Computation of the radiation pressure force using the FDTD technique requires that the time average of the Poynting vector be computed for the object under analysis [21,22]. Several approaches exist for obtaining this result but the most efficient and easiest to implement is to first propagate the beam without the desired object in place until a computation state is reached where all grid points are (E, H) field excited. Then the time average of the Poynting vector (Sxo, Syo, Szo) is obtained using all grid points that are directly adjacent to the PML but not part of the PML. This averaging can be taken over several source cycles. The Poynting vector obtained represents the energy flow of the source without the object of interest in the system. The second computation is performed with the object reintroduced in the grid and the grid fields re-initialized. The beam is again propagated through the grid until all grid points are excited and a speudo-steady state in the fields is obtained. The time average of the Poynting vector (Sx, Sy, Sz) is again obtained using the same grid points as for the reference. The averaging is taken over the same time interval as the reference (object free) time interval. The reference and object Poynting vector values calculated can be used to obtain the optical forces on the object. The advantages of using this two step computation technique lies in the fact that; the single reference reading can be used for many objects and positions; the system can contain many objects of which only the desired object is removed for the referencing and replaced for force computation; the energy flow is computed over the entire grid boundary and as a result includes all light-particle interactions; there is no need to extend computation results to the far field; there is no need to track the interactions taking place at or near the objects surface. In a 2-D computation domain with the Δx and Δy grid space equal, the X and Y force components can be computer through:
with L equal to the number of points over which the Poynting value is computed multiplied by the axial grid point spacing. As a final step in the determination of the optical force present on the object the force values are normalized to a source providing one mW of optical power. This is accomplished by computing the Poynting vector over the source points using Axy and the grid point spacing. In some instances an additional computation step is required as most published FDTD computation algorithms utilize normalized E and H field values [18]. The FDTD force computation technique does not lend itself easily to the computation of the dynamic trap [23] since the time scales involved are incompatible. The FDTD computation time scale is in the order of pico-seconds while that of the dynamic behavior of a particle in the beam is in the order of milli-seconds to seconds. The difference in time scales could be overcome by introducing a hybrid time scale leapfrog model in which force computations occur in one time scale followed by micro-displacements in the other time scale.
3. Computation results
The FDTD algorithm used to propagate the beam through the optical system is one available by Sullivan [24] and rewritten in a graphic display friendly language of Visual Basic. The computer code hierarchy implemented for our computations requires that the objects, source and PML be defined first, then the system is discretized into a fine mesh. The optical source is propagated through the mesh and after a sufficient number of computation iterations the field information is processed and the optical force information is extracted. In the computations, a soft source is employed with each grid point of the source contributing to the amplitude of the incident E and H fields. A diverging Gaussian beam propagating left to right across the grid can be defined at the source location in Fig. 1 by specifying that the amplitude factor Axy be Gaussian in nature:
and that the Ez and Hx components on the source grid points be evaluated as:
for normalized E and H field components and source optical frequency ω. (xc, yc) refer to the center of the minimum waist location of the optical source. Defining Ez and Hx ensures that the cross product gives a left to right directed beam. A plane wave propagating from left to right would have unit amplitude factor, Axy=1, and be defined over a line of grid points extending from the lower to upper range of X values. If a focusing beam is required at the trap object, it can be obtained by imaging either the plane wave or Gaussian beam through a positive lens or by defining the beam over a convergent wavefront.
The computation results are presented using the trapping efficiency defined by the dimensionless function Qx and Qy expressed as:
where c is the speed of light in vacuum, Fx,y are the radial and axial forces present on a object, and P is the optical power of the source. Figures 2 and 3 show results for a dielectric disc subjected to a divergent Gaussian intensity profile beam (Wo=1 µm, λo=0.6328 µm). In these figures the (radial, axial) offsets correspond to the objects coordinates (dx, dy). In Fig. 2 the disc has a higher index, ns=1.58, than the ambient medium, na=1.33 and a radius rs=λo.
Fig. 2. (a)Efficiency factor Qx for radial force on a high dielectric sphere in water. (square-1 µm, circle-4 µm from minimum waist). (b)Efficiency factor Qy for the axial force on the same sphere. Sphere and beam parameters defined in text. The FDTD analysis predicts that axial trapping is expected and that the sphere is pushed in the direction of beam propagation.
The radial force acts to pull the disc in alignment with the beam center while the axial force tends to push the disc along the beam propagation direction, a behavior in agreement with experiments and published results [1]. Interpretation of Fig. 3 indicates that the radial force computed using the same Gaussian beam pushes the low dielectric disc, ns=1.00, out of axial alignment with the beam. A behavior consistent with the properties of low dielectric objects subjected to a focused laser beam in trapping experiments [1].
Fig. 3. (a) Efficiency factor Qx for radial force on a low dielectric sphere in water. (square-1 µm, circle-4 µm from minimum waist). (b) Efficiency factor Qy for the axial force on the same sphere. Sphere and beam parameters defined in text. The FDTD analysis correctly predicts that the sphere is pushed out of the beam and along the beam propagation direction.
Figure 4 shows the FDTD grid domain when a positive lens is used to generate a focusing beam. The index of refraction and diameter of the lens are chosen as nℓ=1.6 and rℓ=5 µm giving a computed focal length of 14.81 µm when placed in water and a beam minimum waist of Wo=0.59 µm using the ABCD method. If desired the focal length and minimum waist can be determined from the propagated beam characteristics of the FDTD algorithm.
Fig. 4. FDTD computation domain traces. Center trace is a gray scale representation of the electric field plotted for each (I, J) grid point. White represents high electric field values. In this trace a plane wave is incident from the left onto a positive lens (large circle) and is focused before diverging. The top trace shows the electric field profile down the Y-axis centerline. The envelope maximum of the E field profile corresponds to the focal point of the lens and can be quantified knowing the grid dimensions or source wavelength in the various mediums. The lower trace is a top plot of the E field profile through the focal region. From this trace the minimum waist can be determined as well the level of aberrations present in the optical trap system.
The top trace of Fig. 4 shows the electric field amplitude as a function of the J grid point along the Y-axis centerline of the computation domain. Knowing the dimensions of the grid in microns, or equivalently the source wavelength in the various dielectric mediums allows accurate measurement of various system parameters. For instance the focal point of the lens corresponds to the minimum waist location of the focused plane incident wave. This focal point corresponds to the maximum of the electric field value as indicated in the top trace. The focal point can then be related back to the exit surface of the lens and through elementary matrix optics manipulations, related to the principal or nodal points. The lower trace shows the electric field profile for a top down cut through the minimum waist location of the focused beam. Measurements on this profile enable the determination of the minimum waist and also provide information on the beam’s aberrations at the focus. Computing the focusing beam parameters through the FDTD technique has the advantage that the lens system can be rendered optically equivalent to the lens system used in the laser trap experiment. Then the aberrations introduced by the lens system are naturally included in the FDTD field computations. Adding the lens system to the trap environment does not increase the computation complexity or computation time. The system referencing of the Poynting vector is accomplished using the entire optical system less the trap object. When the trap object is reintroduced for the second beam propagation, the difference in the Poynting vectors is related to the optical force on the re-introduced object. Figure 5 shows the efficiency factor, Qy, computation result when the high dielectric disc (rs=2 µm, ns=1.45) is translated through the focused beam. In the region of the focus, the negative nature of the trapping efficiency factor indicates optical trapping is possible in this focused beam along the centerline. This is the feature of the single beam trap that has been exploited in numerous experimental laser systems where light is used to manipulate objects. The other structure on the curve results from the incident plane wave being poorly focused through a single element lens. In the focus region the beam is not Gaussian and contains aberrations. The structure in the curve has been observed by others in their theoretical models and is interpreted based on a ray-particle interaction versus object position in the focal region [25–27].
Fig. 5. Efficiency factor Qy for the high dielectric disc axially aligned with the beam center and propagated through the focal region of the beam. Focusing system and E field profile shown in figure 4. In the minimum waist region Qy is less than zero indicating that laser trapping can be accomplished with these beam-object-focusing-system parameters.
The FDTD technique presented lends itself to the study of the optical trapping forces on all types of objects. The requirement is that the object’s parameters be suitably defined such that the FDTD grid points, which make up the discretized object, represent the objects electromagnetic properties. Any optical trapping system configuration and light source beam profile can be explored and the system optimized through FDTD analysis. Figure 6 shows the efficiency factor, Qx, for a dielectric thick walled shell of 5 µm outer radius, 3 µm inner radius. The shell has an index of refraction 1.45, the inner region has an index of refraction of 1.00, and the ambient medium has an index of refraction of 1.33. The beam is chosen to have a minimum waist of 1 µm and a wavelength of 0.785 µm. The figure shows that an unstable equilibrium exists when the shell’s center coincides with the beam’s central propagation axis. Axial trapping is possible when the shell is offset from the beam’s central axis and occurs at the position where the zero crossing of Qx is observed. This behavior is consistent with the dielectric object minimizing its energy by placing the high dielectric region(s) in the strong E field region. The system parameters used here for the thick shell and beam were also used in a previously published paper by Gauthier[28] where the optical force properties were modeled based on a hybrid ray-wave-EM model. The FDTD computation and the hybrid model agree on the nature of the optical forces present for this system.
Fig. 6. Radial trapping efficiency (Qx) for the thick walled dielectric shell. The central region has 1.00 index, the shell has 1.45, and the ambient index is 1.33. The shell is pushed out of axial alignment with the beam but may be axially trapped, offset from the beam axis as indicated in the figure.
Of interest in several system configurations is the dual beam optical trap. It is simpler to implement than the single beam trap since the object is held in 3-D using a push-push effect of the diverging counter propagating beams [17,29–31]. The FDTD technique has been applied to the modeling of the dual beam system by introducing a pair of counter propagating Gaussian beams. The left to right beam has its source location on the left and defined as above through Eqs. (1) and (2), the other beam placed on the right is defined as above except that the Hx source contribution sign is made negative by placing negative sign in front of the amplitude Axy. The cross product of E and H gives a negative Y directed beam. The two beams have a wavelength of λ0=0.6328 µm and minimum waist of 1.5 µm. The ambient medium has an index of refraction of 1.33. Figure 7(a), shows the optical configuration when the beam minimum waist separation corresponds to a 4F system with F=3.99 µm determined as the focal length of the polystyrene sphere introduced between the beams (ns=1.58, rs=1.2656 µm). Both beams have been propagated a short distance in the optical system. As shown, the propagated beam on the left is being focused by the sphere and is converging, while the beam on the right has not yet encountered the sphere and is diverging. Figure 7(b) shows the axial efficiency factor Qy on the sphere in the 4F system. The central mid point between the two sources corresponds to an unstable trap position as indicated by the zero crossing of the Qy factor and positive slope of the Qy line in the zero region. Figure 7(c) shows the axial efficiency factor Qy on the sphere in the 6F system. The central mid point between the two sources corresponds to an stable trap position as indicated by the zero crossing of the Qy factor and negative slope of the Qy line in this zero region. In the dual beam optical trap configuration, the presence of the ball acts as a thick lens and if properly positioned between the sources and with correct source-to-source separation, the minimum waist of one beam may be imaged onto the other beam’s minimum waist. This imaging suggests a means of enhancing fiber-to-fiber coupling should the beams be generated by single mode optical fibers. The FDTD technique would facilitate the calculation of the coupling in such a system as the algorithm could use the beam profiles defined by the optical fibers and include a full fiber-to-fiber propagation system. Figure 8 shows this optical configuration where the sphere is placed midway between the end faces of two single mode optical fibers. The ratio of the input and output (E, H, S) values can be related to optical coupling and loss in the coupling system. Equivalently, FDTD based field profiles could be used to compute the usual modal field overlap integrals in order to obtain the coupling ratio. Experimental confirmation of the power coupling enhancement in a dual beam optical fiber trap configuration has been published by Gauthier [17] et al.
Fig. 7. (a) FDTD field profile for the dual beam counter-propagating trap configuration. Left and right Gaussian beams have been propagated a short distance in the system. Beam on the left is focusing after passing through the lens while the beam on the right is still diverging. (b) Axial force profile obtained when the beam separation corresponds to 4F. The mid point axial location corresponds to an unstable equilibrium for this sphere between these beams. (c) Axial force profile obtained when the beam separation corresponds to 6F. The mid point axial location corresponds to a stable equilibrium for this sphere between these beams.
Fig. 8. FDTD E-field trace of the fiber-to-fiber coupling through a dual beam trapped ball lens. The FDTD trapping computation engine developed permits the optical system to be modeled, the field and beams propagated and the resultant E fields permit the determination of the coupling efficiency.
4. Conclusion
In this paper we have shown that the 2-D FDTD computation technique for propagating electric and magnetic fields can be used to determine the radiation pressure force on various dielectric objects. The technique is based on computing the energy flow at the boundary of the FDTD grid domain with and without the trap object present. The computation technique does not require information on the fields in the far field or at the object surface. The computation technique is shown to reproduce the standard observations of high and low dielectric spheres and the thick wall shell. Due to the grid discretization process of the FDTD technique, experimental optical systems can be accurately represented and the propagating fields will include effects such as optical aberrations, mechanical misalignments, losses, … Steps are being taken to refine the computation engine such that a 3-D FDTD grid domain can be defined for the analysis of 3-D optical systems and the introduction of computation routines to include optical torque.
References
1. A. Ashkin, “Acceleration and trapping of particles by radiation pressure,” Phys. Rev. Lett. 24, 156–159 (1970). [CrossRef]
2. A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, “Observation of a single-beam gradient force optical trap for dielectric particles,” Opt. Lett. 11, 288–290 (1986). [CrossRef] [PubMed]
3. P. Galajda and Pal Ormos, “Orientation of flat particles in optical tweezers by linearly polarized light,” Opt. Express 11, 446–451 (2003). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-5-446 [CrossRef] [PubMed]
4. S. J. Parkin, T. A. Nieminen, N. R. Heckenberg, and H. Rubinsztein-Dunlop, “Optical measurement of torque exerted on an elongated object by a noncircular laser beam,” Phys. Rev A 70, 023816 (2004). [CrossRef]
5. R. C. Gauthier, R. N. Tait, and M. Ubriaco, “Activation of microcomponents with light for micro-electro-mechanical systems and micro-optical-electro-mechanical systems applications,” Opt. Lett. 41, 2361–2367 (2002).
6. E. Higurashi, H. Ukita, H. Tanaka, and O. Ohguchi, “Optically indiced rotation of anisotropic micro-objects fabricated by surface micromachining,” Appl. Phys. Lett. 64, 2209–2210 (1994). [CrossRef]
7. R. W. Applegate Jr., J. Squier, T. Vestad, J. Oakey, and D. W. M. Marr, “Optical trapping, manipulation, and sorting of cells and colloids in microfluidic systems with diode laser bars,” Opt. Express 12, 4390–4398 (2004). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-19-4390 [CrossRef] [PubMed]
8. M. Goksor, J. Enger, and D. Hanstrop, “Optical manipulation in combination with multiphoton microscopy for single-cell studies,” Appl. Opt. 43, 4831–4837 (2004). [CrossRef] [PubMed]
9. F. Qian, S. Ermilov, D. Murdock, W. E. Brownell, and B. Anvari, “Combining optical tweezers and patch clamp for studies of cell membrane electromechanics,” Rev. Sci. Inst. 75, 2937–2942 (2004). [CrossRef]
10. K. F. Ren, G. Grehan, and G. Gouesbet, “Prediction of reverse radiation pressure by generalized Lorentz-Mie theory,” Appl. Opt. 35, 2702–2710 (1996). [CrossRef] [PubMed]
11. J. Lock, “Calculation of the radiation trapping force for laser tweezers by use of generalized Lorentz-Mie theory. I. Localized model description of an on-axis tightly focused laser beam with spherical aberrations,” Appl. Opt. 43, 2532–2544 (2004). [CrossRef] [PubMed]
12. D. Ganic, X. Gan, and M. Gu, “Exact radiation trapping force calculation based on vectorial diffraction theory,” Opt. Express 12, 2670–2675 (2004). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-12-2670 [CrossRef] [PubMed]
13. R. Goussgard and T. Lindmo, “Calculation of the trapping force in a strongly focused laser beam,” J. Opt. Soc. Am B 9, 1922–1930 (1992). [CrossRef]
14. A. Ashkin, “Forces of a single-beam gradient laser trap on a dielectric sphere in the ray optics regime,” J. Biophys. 61, 569–582 (1992). [CrossRef]
15. R. C. Gauthier and A. Frangioudakis, “Theoretical investigation of the optical trapping properties of a micro-optic cube glass structure,” Appl. Opt. 39, 3060–3070 (2000). [CrossRef]
16. R. C. Gauthier, “Optical levitation and trapping of a micro-optic inclined end-surface cylindrical spinner,” Appl. Opt. 40, 1961–1973 (2001). [CrossRef]
17. R. C. Gauthier, M. Friesen, T. Gerrard, W. Hassouneh, P. Koziorowski, D. Moore, K. Oprea, and S. Uttamalingam, “Self-centering of a ball lens by laser trapping: fiber-to-fiber coupling analysis,” Appl. Opt. 42, 1610–1619 (2002). [CrossRef]
18. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas and Propagat. 14, 302–307 (1966). [CrossRef]
19. J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1994). [CrossRef]
20. A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain, Boston: Artech House, 1995.
21. W. L. Collett, C. A. Ventrice, and S. M. Mahajan, “Electromagnetic wave technique to determine radiation torque on micromachines driven by light,” Appl. Phys. Lett. 82, 2730–2732 (2003). [CrossRef]
22. D. Zhang, X.-C. Yuan, S. C. Tjin, and S. Krishnan, “Rigorous time domain simulation of momentum transfer between light and microscopic particles in optical trapping,” Opt. Express 12, 2220–2230 (2004). [CrossRef] [PubMed]
23. R. C. Gauthier and M. Ashman, “Simulated dynamic behavior of single and multiple spheres in a trap region of focused laser beams,” Appl. Opt. 37, 6421–6431 (1998). [CrossRef]
24. D. Sullivan, Electromagnetic simulation using the FDTD method, IEEE Press Series on RF and Microwave Technology, New York, 2000.
25. K. Visscher and G. J. Brakenhoff, “Theoretical study of optically induced forces on spherical particles in a single beam trap I: Rayleigh scatterers,” Optik 89, 174–180 (1992).
26. K. Visscher and G. J. Brakenhoff, “Theoretical study of optically induced forces on spherical particles in a single beam trap II: Mie scatterers,” Optik 90, 57–60 (1992).
27. S. Nemoto and H. Togo, “Axial force acting on a dielectric sphere in a focus laser beam,” Appl. Opt. 37, 6386–6394 (1998). [CrossRef]
28. R. C. Gauthier, “Laser-trapping properties of dual component spheres,” Appl. Opt. 41, 7135–7144 (2002). [CrossRef] [PubMed]
29. G. Roosen, “A theoretical and experimental study of the stable equilibrium positions of spheres levitated by two horizontal laser beams,” Opt. Commun. 21, 189–194 (1977). [CrossRef]
30. A. Constable, J. K. Kim, J. Mervis, F. Zarinetchi, and M. Prentiss, “Demonstration of a fiber-optic light-force trap,” Opt. Lett. 18, 1867–1869 (1993). [CrossRef] [PubMed]
31. E. Sidick, S. D. Collins, and A. Knoesen, “Trapping forces in a multiple-beam fiber- optic trap,” Opt. Lett. 36, 6423–6433 (1997).
