## Abstract

We theoretically investigate the three-dimensional (3D) trapping force acting on a microsphere held in a pair of counterpropagating beams produced by the generalized phase contrast (GPC) method. In the case of opposing beams of equal power, we identify the range of beam waist separation *s* that results in a stable 3D optical potential-well by assessing the dependence of the axial and transverse force curves on *s*. We also examine how the force curves are influenced by other parameters such as size and refractive index of the microsphere. Aside from force curves of beam tandems with equal powers, we also numerically calculate force curves for cases of beam pairs having disparate relative strengths. These calculations enable us to elucidate the large dynamic range for axial position control of microparticles in GPC-based counterpropagting-beam traps.

©2006 Optical Society of America

## 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=2*a* 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 d**F**^{+}
_{‖} and the perpendicular force d**F**^{+}
_{⊥}, may be resolved when a ray is incident to an infinitesimal region d*S* on the sphere’s surface:

where *n*
_{1} is the refractive index of the host medium, d*P* 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 d*P* through d*S* 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 d**F**^{+}
_{‖} d**F**^{+}
_{⊥} with unit vectors *x̂*, *ŷ*, and *ẑ* over the illuminated sphere surface.

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}
+*a*sin*θ*cos*φ*, *y*
_{p}
=*y*
_{o}
+*a*sin*θ*sin*φ*, and *z*
_{p}
=*z*
_{o}
-*a*cos*θ* 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^{2}φ or cos^{2}
*φ* 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*
_{1}.

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*
^{2} at *z*=0) occurs at *z*~3.8 µm which is consistent with the predicted position *z*
_{I,max}=*n*
_{1}
*R*
^{2}/*λ*, when considering that the beam (*λ*=830 nm) propagates in water - commonly used particle host medium with index *n*
_{1}=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.

## 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*
_{2}=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).

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.

At sufficient *s*>*s*
_{c}
, Figs. 6(a) and 6(b) illustrate the axial and transverse force curves for different particle indices *n*
_{2}. 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*
_{1}–*n*
_{2} 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.

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.

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.

## 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.

## Acknowledgments

This work has been funded by the European Science Foundation through the Eurocores-SONS program and through the European-FP6-NEST project, ATOM3D, the Danish Technical Research Council (FTP) and partially by an internal grant awarded by Risø National Laboratory.

## References and links

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

**2. **A. Constable, J. Kim, J. Mervis, F. Zarinetchi, and M. Prentiss, “Demonstration of a fiber-optical light-force trap,” Opt. Lett. **18**, 1867–1869 (1993). [CrossRef] [PubMed]

**3. **E. R. Lyons and G. J. Sonek, “Confinement and bistability in a tapered hemispherically lensed optical fiber trap,” Appl. Phys. Lett. **66**, 1584–1586 (1995). [CrossRef]

**4. **G. Roosen and C. Imbert, “Optical levitation by means of two horizontal laser beams: a theoretical and experimental study,” Phys. Lett. A **59**, 6–8 (1976). [CrossRef]

**5. **E. Sidick, S. D. Collins, and A. Knoesen, “Trapping forces in a multiple-beam fiber-optic trap,” Appl. Opt. **36**, 6423–6433 (1997). [CrossRef]

**6. **P. J. Rodrigo, V. R. Daria, and J. Glückstad, “Real-time three-dimensional optical micromanipulation of multiple particles and living cells,” Opt. Lett. **29**, 2270–2272 (2004). [CrossRef] [PubMed]

**7. **P. J. Rodrigo, V. R. Daria, and J. Glückstad, “Four-dimensional optical manipulation of colloidal particles,” Appl. Phys. Lett. **86**, 074103 (2005). [CrossRef]

**8. **I. R. Perch-Nielsen, P. J. Rodrigo, and J. Glückstad, “Real-time interactive 3D manipulation of particles viewed in two orthogonal observation planes,” Opt. Express **18**, 2852–2857 (2005). [CrossRef]

**9. **P. J. Rodrigo, L. Gammelgaard, P. Bøggild, I. R. Perch-Nielsen, and J. Glückstad, “Actuation of microfabricated tools using multiple GPC-based counterpropagating-beam traps,” Opt. Express **13**, 6899–6904 (2005). [CrossRef] [PubMed]

**10. **J. Glückstad, “Phase contrast image synthesis,” Opt. Commun. **130**, 225–230 (1996). [CrossRef]

**11. **J. Glückstad and P. C. Mogensen, “Optimal phase contrast in common-path interferometry,” Appl. Opt. **40**, 268–282 (2001). [CrossRef]

**12. **J. W. Goodman, *Introduction to Fourier Optics*, Second Edition (McGraw-Hill, New York, 1996).