1. Introduction

The first demonstration showing the ability of light to trap a tiny particle in three dimensions (3D) utilized two mildly focused Gaussian beams in a counterpropating geometry [1]. An interesting physical observation noted in that work was the need for a ‘positive’ separation between the two beam waists to achieve a stable optical potential- well at the midpoint of separation. Positive beam waist separation s means that the beam propagating to the right has its waist positioned to the left of the other beam’s waist. A variant of such a counterpropagating-beam (CB) trap that uses two coaxially aligned optical fibers [2, 3] also qualitatively suggests the same requirement on s.

In an attempt to rigorously model CB traps [4, 5], a more quantitative study on how trap stability depends on s has been carried out. It was found that to produce a stable CB trap (for a microsphere with diameter comparable to the beam waist) the two waists must have sufficiently large separation and that decreasing it to below a critical value s _c results in an unstable potential-well with possible emergence of additional multiple points of stable equilibrium. Furthermore, one can deduce in Ref. 5 that the value of s _c increases with microsphere diameter. Trap efficiencies were also investigated as functions of parameters such as the microsphere’s refractive index and radius. However, only the case of equally powerful opposing beams was analyzed.

In this paper, we report the extended calculations of the optical forces due to counterpropagating-beam traps of variable power ratios. This analysis is made specific to the case of CB traps generated via the generalized phase contrast (GPC) method [6–9], thereby providing a theoretical model for GPC-based CB traps for the first time. We calculate for the value of s _c corresponding to a symmetric GPC-based CB trap to find the operational beam waist separation and compare with our previous experimental implementations. GPC-based CB traps are real-time and interactively reconfigurable arrays of optical traps that rely on the energy-efficient intensity-visualization of phase patterns imparted on an incident plane wave by a spatial light modulator [10, 11]. The versatility of GPC-based CB traps has been exemplified by its ability to interactively manipulate a large number of microspheres or cells, and actuate microfabricated structures in 3D [6–9]. Interactive micromanipulation is truly achieved in real-time as the formation of arbitrary trapping patterns is direct and does not require iterative holographic computation.

The succeeding section describes the analytical derivation of expressions for the x, y, and z -components of the trapping force acting on a sphere due to a GPC-based CB trap. In section 3, we obtain the intensity distribution describing the propagation of one of the constituent beams through the host medium. In section 4, the analysis in sections 2 and 3 are used to numerically calculate the axial and transverse force curves as functions of beam waist separation, particle refractive index, and particle radius. We also propose a force constant (trap stiffness) versus s plot as a quantitative means of finding a critical separtation s _c . Finally, the axial force curves for counterpropagating beams with unequal strengths are used to estimate the CB trap’s dynamic range of axial position control.

2. Three-dimensional trapping force components

We consider a purely dielectric microsphere with diameter=2a illuminated by two coaxial (i.e. along z-axis) counterpropagating orthogonally polarized optical fields (Fig. 1(a)) corresponding to intensity distributions described by I ⁺(x, y, z) and I ^-(x, y, z) where +(-) denotes a beam propagation directed towards the positive (negative) z-axis. First, we examine the forces induced by I ⁺(x, y, z) on the particle. Guided by Fig. 1(b), two force components, namely the parallel force dF ⁺ _‖ and the perpendicular force dF ⁺ _⊥, may be resolved when a ray is incident to an infinitesimal region dS on the sphere’s surface:

where n ₁ is the refractive index of the host medium, dP is the differential power of the ray, c is the speed of light in vacuum, and ê_‖ and ê_⊥ are the unit vectors parallel and perpendicular to the incident ray, respectively. The relations of ê_‖ and ê_⊥ with surface normal n̂, and the angles of incidence α _i and refraction α _r , which are related by Snell’s law, are shown in Fig. 1(b). Factors q _‖ and q _⊥ denote the fractions of momentum, transferred to the particle by the ray, in the direction ê_‖ and ê_⊥, respectively. They are related to α _i and α _r by [4, 5]

where the particle surface reflectance R and transmittance T are carefully determined by noting the polarization of the incident field and the plane of incidence. The differential power dP through dS is given by

where the coordinates (x _p , y _p , z _p ) refer to the point of incidence P with respect to the unprimed rectangular coordinate system in Fig. 2. The effective forces acting along three Cartesian axes can be obtained from the integration of the scalar products of dF ⁺ _‖ dF ⁺ _⊥ with unit vectors x̂, ŷ, and ẑ over the illuminated sphere surface.

 

Fig 1. (a) Graphical illustration of a microsphere centered at (x _o ,y _o ,z _o ) illuminated by two counterpropagating beams I ⁺ and I ^-. (b) The microsphere’s cross-section in the plane of incidence associated with a ray hitting the surface at an angle α _i measured from the normal.

Download Full Size | PDF

 

Fig. 2. 3D illustration of an incident ray (red arrow) impinging at point P on the surface of a microsphere centered at point O, which is located at (x _o ,y _o ,z _o ) of the unprimed Cartesian coordinate system.

Download Full Size | PDF

Finding from our geometry that ê_‖=ẑ and ê_⊥=x̂cosφ+ŷsinφ and that α _i =θ(of the spherical coordinates of point P with respect to sphere center O positioned at (x _o, y _o , z _o )), we obtain

where, the relations x _p =x _o +asinθcosφ, y _p =y _o +asinθsinφ, and z _p =z _o -acosθ are used to determine the intensity at different points on the surface of integration. The Fresnel integral [12] is used to numerically obtain the intensity I ⁺(x _p , y _p , z _p ), which is a result of the free-medium propagation of the function e ⁺(x, y) describing the field distribution at the objective focal plane z=0.

The integrations in Eq. (4b) through Eq. (4d) are performed numerically to determine the force components. We also consider the incident beam to be linearly polarized along the x-axis. Since different planes of incidence exist for various points on the illuminated sphere surface, each of the force components is calculated for the TE and the TM components separately, and afterward the contributions are combined. Therefore, in addition to the proper calculation of R and T in Eq. (2a) and Eq. (2b), one must also multiply the integrands of Eqs. (4) a factor of sin²φ or cos² φ that corresponds to the fraction of the input power carried by the TE or TM component.

3. Tophat field distribution and propagation

A commonly applied transverse field distribution at z=0 when forming a GPC-based counterpropagating-beam trap may be approximated by a tophat profile

where beam power P ⁺ is uniformly distributed over the circle of radius R. The 3D intensity distribution I ⁺(x, y, z) for z>0 is obtained by numerically calculating the Fresnel integral [12]. This corresponds to the diffraction pattern that the field distribution described by Eq. (5) undergoes upon propagation within the medium of refractive index n ₁.

Figure 3 shows a plot of the calculated diffraction pattern in the xz-plane. In this sample calculation we used R=1.5 µm and P ⁺=5 mW, which are typical values realized in our past experiments. An observation one can make from Fig. 3 is that a global intensity maximum (about four times the intensity P ⁺/πR ² at z=0) occurs at z~3.8 µm which is consistent with the predicted position z _I,max=n ₁ R ²/λ, when considering that the beam (λ=830 nm) propagates in water - commonly used particle host medium with index n ₁=1.33. It is also worth to note that a plane defined by z=z _I,max separates the Fresnel and the Fraunhofer regions of the diffraction pattern. In the Fraunhofer region, one finds a central lobe that decreases in peak intensity but increases in width for increasing z position. In the section that follows, we use these observations to explain the obtained dependence of axial and transverse forces on the position of the particle within the CB trap.

 

Fig. 3. Pseudo-color plot of the xz-plane intensity distribution appearing at z>0 for a typical GPC-generated beam directed towards the positive z-axis. At z=0, the beam has a tophat field profile extending from -1.5 µm to 1.5 µm of the x-axis. Red color represents maximum intensity while blue represents minimum intensity.

Download Full Size | PDF

4. Numerical calculation of force curves

With the aid of the mathematical analysis outlined above, we performed simulations of typical trapping conditions we have implemented in previous GPC-based trapping experiments. But in order to describe the net axial and net transverse forces independently of the total CB power P _t =P ⁺+P ^-, we instead calculate for the trapping parameters given by

where $F_x^{-}$ and $F_z^{-}$ are the force components due to the other beam with power P ^- directed towards negative z-axis (see Fig. 1(a)). Numerical calculations of $F_x^{-}$ and $F_z^{-}$ similarly use Eqs. (4) and (5) except that P ⁺ in Eq. (5) is replaced by P ^- and z _o, which denotes the axial position of the microsphere, is replaced by s-z _o . As depicted in Fig. 1(a), the beam waist separation s is defined here as the distance separating the two microscope objective focal planes where the two GPC-generated tophat beams are imaged.

For a polystyrene microsphere (refractive index n ₂=1.59 ; radius a=1.5 µm) in water, we first look at the case where the opposing tophat beams are identical in power, P ⁺=P ^-, and in size R=1.5 µm. Figures 4(a) and 4(b) illustrate, for different waist separation s, the axial and transverse force profiles, respectively - axial (transverse) force as a function of axial z _o -s/2 (transverse x _o ) offset. Note that the force curves feature an odd-function property indicating a trapping potential with both axial and transverse symmetry about the separation midpoint. The plots in Fig. 4(a) also clearly show that, except for s=20 µm, a GPC-based CB trap provides a stable harmonic axial potential-well as characterized by the linear restoring force (i.e. negative slope) about the equilibrium position. For s=20 µm however, one obtains at the separation midpoint a positive slope in the axial force curve making it an axially unstable point of equilibrium. Although at the midpoint (in this case given by (0,0, s/2=10 µm) of the xyz-coordinate system) the trap potential appears to be transversely stable and harmonic (Fig. 4(b)), a small perturbation in the particle’s axial position will cause it to be ejected into either direction of the z-axis. These results are reminiscent of those for a fiber-based CB trap with a Gaussian mode profile [5]. Despite of finding that the region where we expect to trap the particle belongs to the Fraunhofer region of both beams, it is worth to note that the tails of the axial force curves in Fig. 4(a) exhibit non-monotonic behavior, which is seen linked to the spatial variations of cross-sectional intensity profiles of each beam in its Fresnel region (see Fig. 3).

 

Fig. 4. Typical dependence of the (a) axial and (b) transverse force components on the respective axial and transverse distances of the microsphere from the midpoint between the two beam waists for different waist separation distances.

Download Full Size | PDF

An alternative scheme may be found in plotting the dependence of the force constants f _x and f _z on s as shown in Fig. 5. Here we assume that, for radial distances from the midpoint (0,0, s 2) comparable with the particle radius, Q _xt ≈-f _x x _o and Q _zt ≈-f _z (z _o -s/2). For the pertinent range of s considered, our findings indicate that f _x decreases with s but remains positive-valued. This decay of transverse trap stiffness with s can be easily understood by recognizing the intensity profiles illuminating the microsphere near the trapping point. One finds that larger s would translate to smaller intensity gradient along the transverse direction as the two counterpropagating beams illuminate the particle with broader central lobe in their respective Fraunhofer regions (see Fig. 3). Note also that the same reason explains the slight increase in the transverse attractive range with s (see Fig. 4(b)). On the contrary, f _z can either be positive giving a stable trap when s>s _c or negative forming an unstable trap, when s<s _c (i.e. f _z (s _c )=0). Clearly, s _c may be accurately obtained as the zero-crossing of a f _z (s)-curve as exemplified by Fig. 5. Furthermore, the use of f _z (s)-curve to identify a maximum axial force constant f _z,max also suggests a procedure for choosing an optimal beam waist separation s _op , for example, by defining f _z (s _op )=f _z,max. In Fig. 5, we find s _op =30 µm which matches the estimated value applied in our previous experiments with the same set of parameters [6, 7]. Note however that the choice for s _op may also take into account the decrease of f _x with s and other parameters like size of the largest particle (for size-varied colloids) and desired dynamic axial range of manipulation as discussed below.

 

Fig. 5. Power-independent stiffness of a symmetric GPC-based CB trap as a function of beam waist separation.

Download Full Size | PDF

 

Fig. 6. Typical dependence of the (a) axial and (b) transverse force components on the respective axial and transverse distances of the microsphere from the midpoint between the two beam waists (s=40 µm) for different refractive indices n ₂. Spheres with radius a=1.5 µm and equal counterpropagating tophat beams of size R=1.5 µm are considered.

Download Full Size | PDF

At sufficient s>s _c , Figs. 6(a) and 6(b) illustrate the axial and transverse force curves for different particle indices n ₂. The result simply confirms that both force components increase with refractive index contrast between particle and suspending medium due to the expected increase in momentum change that light undergoes at the n ₁–n ₂ interface. Meanwhile, Figs. 7(a) and 7(b) provide the force curves for different sphere radii a. For sphere radii a≤R, a smaller sphere captures lesser beam power and thus translate to weaker forces particularly near the equilibrium point and thus giving weaker trap stiffness. In addition, the location of the leading axial force maximum is closer to the point of stable equilibrium for larger a as shown in Fig. 7(a). This is expected to diminish the axial range within which the particle can be stably positioned.

 

Fig. 7. Typical dependence of the (a) axial and (b) transverse force components on the respective axial and transverse distances of the microsphere from the midpoint between the two beam waists for different sphere radii a. Polystyrene spheres with refractive index n ₂=1.59 and equal counterpropagating tophat beams of size R=1.5 µm and separation s=40 µm are considered.

Download Full Size | PDF

Figures 8(a) and 8(b) show the axial and transverse force curves for varied differential powers s P ⁺- P ^-. In our simulation, we pattern the way P ⁺- P ^- is changed on how it is being efficiently done experimentally with a polarization scheme [6–9] that maintains constant total power P _t =P ⁺+P ^-. From zero differential power, we tune P ⁺- P ^- to a value that shifts the stable point of equilibrium to a cutoff axial offset Δz where the axial force constant decreases by 50%. The axial offset from -Δz to +Δz defines a competitively large dynamic range for axial position control of the particle. We also plot the axial force curve that gives another asymmetric potential with stable equilibrium point found at an axial offset of Δz/2. At xy-planes containing the axially displaced equilibrium points, the transverse force curves are determined. The manner in which the associated trap stiffness vary with the equilibrium position are shown in the respective insets.

The theoretical modeling presented above is expected to play a vital role for optical micromanipulation that utilizes GPC-technology. Envisioned colloidal and microbiological experiments that will benefit from the versatility of the GPC-based multiple-beam trapping system may be expected to handle polydisperse samples. Having a comprehensive theoretical understanding of how optical trap performance is affected by trapping parameters will aid in the actual implementation, optimization, calibration or troubleshooting of the system. In Fig. 9, we plot the f _z (s)-curves for microsphere radii a=1.0 µm, 1.5 µm, and 2.5 µm, say, to pre-analyze a trapping experiment with a size-polydisperse sample. The plot is indicative of the role that particle radius a plays in identifying an operational beam waist separation s _op >s _c . From the figure, one may anticipate the relative trap stiffness for the two smaller microspheres and the inability to (axially) trap the largest if the chosen s _op =30 µm (where f _z,max for microsphere with a=1.5 µm occurs). To simultaneously trap the three particles using three identical, symmetric CB traps having sufficient power, Fig. 9 suggest that s _op >35 µm must be implemented.

 

Fig. 8. Plot of the (a) axial and (b) transverse force curves for varied differential powers P ⁺-P ^-. Insets show the change in (a) axial and (b) transverse trap stiffness at three different points of stable equilibrium. The obtained dynamic range for axial position control is ±Δz corresponding to a depth of ~17.0 µm, in good comparison with experimentally obtained ~20.0 µm [6]. Polystyrene spheres with radius a=1.5 µm and index n ₂=1.59, and equal counterpropagating tophat beams of size R=1.5 µm and separation s=40 µm are considered.

Download Full Size | PDF

This is confirmed by the particular situation where s _op = 40 µm, which results only in a stable trap potential as shown by the force curves in Fig. 7 for the three different particle radii.

 

Fig. 9. Dependence of GPC-based CB trap’s axial stiffness on beam waist separation for different microsphere radii. The plot shows that a larger critical separation s _c (i.e. zero-crossing) is obtained for a larger particle.

Download Full Size | PDF

5. Conclusion

A rigorous analysis of the axial and transverse optical forces that refractive microspheres experience in GPC-based counterpropagating-beam traps was theoretically performed for the first time. One solid finding was that, by underpinning proper beam waist separation, GPC-based CB traps can offer stable optical potentials, which may be accurately considered harmonic over a significant range of the trapping volume. Analogous to those constructed with a pair of optical fibers, the GPC-CB traps were also found to be stiffer transversely than axially. Among the parameters affecting the trap forces and stiffness are beam waist separation, particle dimension and particle-host index contrast. Additionally, along with the inherently large transverse range for manipulation, the already large axial dynamic range can be further extended by selecting larger beam waist separation. The presented theoretical model plays a critical role in the accurate understanding of the 3D forces produced by GPCbased CB traps and is believed to advance the applications of the GPC trapping system (e.g. GPC-based traps with well-defined stiffness can simultaneously function as multiple force transducers/probes). The set of formulations outlined here will also serve as a framework for future modeling of CB traps with complex intensity and/or polarization spatial profiles.

It is well known that the model proposed here is valid for particles with radii much larger than the wavelength of the incident beams [4, 5]. Thus in the Mie regime, a≫λ, the geometrical optics approximation may offer accurate results. However, good agreement of results from geometrical optics approach with those obtained experimentally has also been achieved for particle sizes comparable to the wavelength [3]. For Rayleigh particles, a≪λ, calculation of the optical forces entails electromagnetic wave theory.