## Abstract

Optical forces on a micro-bubble were computed using the Finite Difference Time Domain method. Non-paraxial Gaussian beam equation was used to represent the incident laser with high numerical aperture, common in optical tweezers. The electromagnetic field distribution around a micro-bubble was computed using FDTD method and the electromagnetic stress tensor on the surface of a micro-bubble was used to compute the optical forces. By the analysis of the computational results, interesting relations between the radius of the circular trapping ring and the corresponding stability of the trap were found.

© 2008 Optical Society of America

## 1. Introduction

Optical Tweezers (OT) use a tightly focused laser to trap and manipulate objects of size as small as several tens of nanometers to a size as large as several tens of micrometers. OT have been actively used in physics and biological experiments due to their unique non-invasive capability. OT are easily instrumented in optical microscopes that use a high Numerical Aperture (NA) objective lens. When a laser is focused by this high NA objective lens, the beam intensity gradient at the focus becomes large such that objects that have large refractive index are pushed in the direction of this gradient.

Since the first discovery of OT by Ashkin et al. [1] in 1986, OT have been widely used in various fields. OT have been used to trap and control DNA [2] and living cells [3]. Furthermore, OT have been used to construct micro-structures from micron-sized dielectric particles [4]. To understand the underlying physics of these applications and experiments, a theoretical computation is essential. The strength of an optical trap varies by the physical property of the target object and laser properties such as the mode, power, minimum waist size and NA. If we can find the most appropriate conditions of an OT for a given experiment, we can greatly enhance the performance or the effectiveness of the experiment. This can be achieved by simulations or physically based calculations. Simulations can not only be used to verify the experimental results but can be used to understand new physical behaviors and devise new extensions. Recently, there has been an increase in attempts to trap not only simple spherical particles but objects of various forms and compositions [5]. Simulations can be used to assess if a given object can be trapped or more difficultly, “how” we can trap objects that are not easily trapped using simple trapping methods.

There are basically three approaches used to compute the trapping forces of an OT based on its size compared to the wavelength of the laser. In the ray optics regime, the object is considered to be much smaller than the wavelength (D≪λ). D is the diameter of the object and λ is the laser wavelength. Conversely in the Rayleigh regime, the object is considered to be much larger than the wavelength (D≫λ). There are numerous relevant articles but we only note the work of Ashkin in the ray optics regime [6] and Harada et al. in the Rayleigh regime [7]. In practice however, we mostly deal with objects that are about the wavelength. And for these objects, one of the widely used computational methods is the FDTD (Finite-Difference Time- Domain) method [8]. The advantages of FDTD are that it is not constrained by object size, shape and compositions. In this study, we used FDTD to compute the electromagnetic field on a micro-bubble. A tightly focused laser is modeled by the non-paraxial Gaussian beam representation and the light-object momentum transfer is modeled by the Maxwell’s stress tensor.

## 2. Forces acting on objects

The force that acts on a target object immersed in a dielectric medium as illustrated in Fig. 1 due to an electromagnetic field can be represented using the Maxwell’s stress tensor [9] as

In Eq. (1), *ε* and *µ* are permittivity and permeability, **n** is the outward normal unit vector from the interior of surface **S** in Fig. 1, *c* is the speed of light in vacuum. **E** and **H** will be defined below. The surface and volume integration domain **S** and **V** represent the object surface and interior in Fig. 1, *da* and *dv* are infinitesimal area and volume. When the force in Eq. (1) is averaged for a monochromatic light in the time duration of a period *T*=*2π*/* λ* we get,

$$=\frac{1}{T}\sum _{T}{\int}_{s}\left[\epsilon \left(\mathbf{E}\xb7\mathbf{n}\right)\mathbf{E}+\mu \left(\mathbf{H}\xb7\mathbf{n}\right)\mathbf{H}-\frac{1}{2}\left(\epsilon {E}^{2}+\mu {H}^{2}\right)\mathbf{n}\right]\mathrm{da}.$$

The operator <> represents the time average. The last term vanished under steady-state assumption.

We modeled the tightly focused laser by the non-paraxial Gaussian beam [10]. When *ω* is the circular frequency and *λ* is the optical wavelength in medium, the Electric fields (**E**) and the Magnetic fields (**H**) can be shown in the Cartesian coordinate as

$$\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\times \mathrm{exp}\left[ik\left(\mathrm{xr}\mathrm{cos}\theta +\mathrm{yr}\mathrm{sin}\theta +z\sqrt{1-{r}^{2}}\right)\right]\mathrm{rdrd}\theta ,$$

$$\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\times \mathrm{exp}\left[ik\left(\mathrm{xr}\mathrm{cos}\theta +\mathrm{yr}\mathrm{sin}\theta +z\sqrt{1-{r}^{2}}\right)\right]\mathrm{rdrd}\theta .$$

In the above two equations, *z* is the beam propagation axis and the integrations were done in polar coordinates where *r is*
$\sqrt{{x}^{2}+{y}^{2}}$ and *θ* is the angle between *r* and the x axis. By defining *ω _{0x}* and

*ω*as the initial transversal Gaussian half widths in the

_{0y}*x*and

*y*axes, a as the ratio of initial amplitudes in x and

*y*directions,

*k*as the wave number (=

*2π*/

*λ*), and lastly

*b*and

*c*as

**A ^{E}** and

**A**can be represented as

^{H}$${\mathbf{A}}^{H}(r,\theta )=\sqrt{\frac{\epsilon}{\mu}\pi}\{\left[\begin{array}{c}-{\omega}_{0y}^{2}\mathrm{exp}\left(-\frac{{\mathrm{cr}}^{2}}{2}\right)\\ -{\mathrm{ar}}^{2}\mathrm{sin}\theta \mathrm{cos}\theta {\omega}_{0x}^{2}\mathrm{exp}\left(-\frac{{\mathrm{br}}^{2}}{2}\right)\\ +{r}^{2}{\mathrm{cos}}^{2}\theta {\omega}_{0y}^{2}\mathrm{exp}\left(-\frac{{\mathrm{cr}}^{2}}{2}\right)\end{array}\right]\frac{\mathbf{i}}{\sqrt{1-{r}^{2}}}$$

$$+\left[\begin{array}{c}r\mathrm{cos}\theta {\omega}_{0y}^{2}\mathrm{exp}\left(-\frac{{\mathrm{cr}}^{2}}{2}\right)\\ -\mathrm{ar}\mathrm{sin}\theta {\omega}_{0x}^{2}\mathrm{exp}\left(-\frac{{\mathrm{br}}^{2}}{2}\right)\end{array}\right]\mathbf{k}\}.$$

The laser power can be obtained using a double integration. That is

$$\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}=\sqrt{\frac{\epsilon}{\mu}}\frac{{\pi}^{3}}{{\lambda}^{2}}{\int}_{0}^{1}\left[{a}^{2}{\omega}_{0x}^{4}\mathrm{exp}\left(-{\mathrm{br}}^{2}\right)+{\omega}_{0y}^{4}\mathrm{exp}\left(-{\mathrm{cr}}^{2}\right)\right]\frac{2r-{r}^{3}}{\sqrt{1-{r}^{2}}}dr,$$

where **S** is the Poynting vector that is equivalent to **E×H**.

## 3. Computation of the trapping force

There are two steps involved in the computation of a trapping force. In the first step, the electromagnetic field of a dielectric object immersed in a medium with an incident field is computed. To increase the accuracy, non-paraxial Gaussian beam was used to model the incident field. With these setups, the electromagnetic field was computed using the FDTD method. FDTD, first introduced by Yee in 1966 [11], is a numerical method for computing the Maxwell’s equations that are partial differential equations in both time and space.

FDTD requires that the space and time be discretized into units called grid distance (*dx*) and time step (*dt*). The time step is automatically determined by the stability condition when grid distance is determined [12]. Thus we only need to concentrate on the grid distance. It is generally difficult to determine the grid distance. The smaller the grid distance the better will be the accuracy. However, the grid distance is limited by the amount of available computer memory. As a rule of thumb, we have used one twentieth of the laser wavelength. After the grid distance is set, each grid is assigned with a physical value depending on whether the object or the medium is at that particular grid location. Subsequently, with the incident field properly set, numerical computation was performed until the field reached a steady-state. By then, Eq. (1) was used to compute the trapping force.

There are several notes worth mentioning. Firstly, we used the UPML (Uni-axial perfectly matched payer) to represent the absorbing boundary layer [12]. Secondly, trapping forces were computed at the dashed circle in Fig. 2, located slightly offset from the object boundary. Thirdly, for computing the forces at different locations, we translated the incident beam instead of translating the object. For example in Fig. 2, the object is centered at the computational domain and the focused Gaussian beam is moved relative to the object. Fourthly, accuracy of the incident field was achieved by using scattered-field formulation at each grid point. Scattered-field formulation requires substantial amount of computational time. To reduce this, we used a look-up table to cache the incident field values in both time and space. Scattered-field formulation can be formulated as

$$\{\begin{array}{c}{\mathbf{E}}_{\mathrm{inc}}=\mathrm{electric}\phantom{\rule{.2em}{0ex}}\mathrm{incident}\phantom{\rule{.2em}{0ex}}\mathrm{field}\\ {\mathbf{E}}_{\mathrm{scat}}=\mathrm{electric}\mathrm{scattered}\phantom{\rule{.2em}{0ex}}\mathrm{field}\\ {\mathbf{H}}_{\mathrm{scat}}=\mathrm{magnetic}\phantom{\rule{.2em}{0ex}}\mathrm{scattered}\phantom{\rule{.2em}{0ex}}\mathrm{field}\end{array}.$$

For our case the conductivity (*σ*) was ‘0’, thus we only needed to consider region where (*ε*-*ε _{0}*) was not ‘0’. This fact was leveraged and the cached incident fields were saved only at the grid points that were near the object.

Fifthly, when assigning the physical property at grid locations that were immediately adjacent to the object surface boundary, we used a value that was the weighted sum of the physical properties of the object and the medium. The weight was determined based on the interior (object) / exterior (medium) volume ratio of the cell whose center is the grid point. Lastly when the forces were computed using Eq. (2), the integration was performed on polar coordinates (*0*≤*θ*≤*2π*, *0*≤*φ*≤*π*), with the radius of the object. The field values were obtained by linearly interpolating field values at adjacent grid points.

## 4. Computation of the trapping force on a micro-bubble

When an object with volume (*V*), specific gravity (*γ*), density (*ρ*) is placed inside a solution, the buoyant force can be calculated as F_{b}=* _{γ}V*=

*, where*

_{ρg}V*g*is the gravitational acceleration constant. For a micro-bubble with a diameter (d

_{b}) of 2um placed in water of specific gravity

*γ*of 9211.9 [kg/m

^{2}s

^{2}], the buoyant force is 0.0407 pN. This is why when micro-bubbles are placed in a water chamber, they float onto the top surface. When they are put under a focused laser, they are pushed away from the focal point due to their low refractive index compared to water. Because of this strange behavior, trapping of micro-bubbles require special techniques. One technique scans a focused Gaussian laser beam around a micro-bubble forming a Circular Trapping Ring (CTR). The result is that the micro-bubble becomes caged inside this ring and by moving this ring, the micro-bubble can be manipulated [13]. An interesting thing is that in this case, the object is trapped through the balance of the trapping force and buoyant force. Other forces such as the gravitational force and the forces due to the thermal motion of water molecules can be neglected [13,14]. For example, the buoyant force of a 2µm micro-bubble is 100 times greater than the weight. There are also other ways for trapping micro-bubbles. For example, Mohanty et al. [15] have used dual line optical tweezers. Other methods were also surveyed in this literature.

The trapping forces of a CTR on a micro-bubble can be decomposed into axial and transversal components. Axial force equilibrium of a micro-bubble is achieved when the axial trapping force is balanced by the buoyant force. However, the transverse force equilibrium is achieved by the circular symmetry of the optical forces on a micro-bubble. In this paper, we will calculate the trapping force imposed by a laser at one location on the CTR. This configuration is illustrated in Fig. 3. And this knowledge will be used to find the optimal way of stably trapping a micro-bubble under a CTR.

Figure 4 illustrates the radial and axial forces on a micro-bubble whose diameter d_{b} is 2.0µm and the refractive index is 1.0001. The simulated environment consists of an objective lens placed in an inverted geometry with the laser propagating in +z direction and the surrounding medium is water whose refractive index is 1.33. The focused Gaussian laser beam was modeled with non-paraxial Gaussian beam equation. Following numbers were used. Wavelength (*λ*)=1064nm, NA=1.2155, Power=10mW. The grid spacing in FDTD calculations were *dx* of 53.2nm and *dt* of 8.872×10^{-17}sec.

#### 4.1 Axial force

To understand the combined effect of both optical force and buoyant force, we plotted the summation of axial force and buoyant force for inverted (Fig. 5) and upright (Fig. 6) geometry. It is where the force is ‘0’ where the trap is in equilibrium. When a micro-bubble is disturbed from an equilibrium shown at inset **B** in Fig. 5, it will move away from the equilibrium. This is an unstable equilibrium. However when disturbed at inset **A** in Fig. 5, it will experience a force that opposes its direction of motion. This is a stable equilibrium.

From Fig. 5, we can deduce an interesting relation between the trap diameter, d_{tr} and the axial trapping stiffness. From the plots of various axial forces as a function of a trap diameter, we can see that the trapping stiffness becomes larger for smaller trap diameter. The slope at the inset **A** denotes the trapping stiffness. However, for a trap diameter below certain value, a critical diameter, trapping zone was non-existence. For example in Fig. 5, trapping is possible when d_{tr} is 2.5µm but not possible for 2µm. Thus we can find when micro-bubbles are trapped by a CTR, there is a limiting trap diameter for a successful trap. By the calculations done for a micro-bubble of diameter d_{b} of 2µm, the trap diameter (d_{tr}) must be greater than 2.5µm.

Figure 6 shows similar results for different objective lens geometry, an upright case. We can see similar results for the relation of a trap diameter d_{tr} to the axial force. This result also verifies the result shown in the literature where Jones et al. [13] have experimentally shown that strongest trap was achieved when d_{tr} is about 1.7 times the d_{b}, the diameter of the micro-bubble. In this case, **B** is stable and **A** is unstable. Other than this, the trapping condition was similar to the inverted case. Smaller trap diameter resulted stiffer trap. For both Inverted and upright geometry, the trapping zone is near the focus, within the depth of focus. Thus micro-bubbles can be easily monitored during trapping.

#### 4.2 Radial force

Figure 7 plots the radial forces at certain z locations near the focus for a stationary trap (not a CTR). From this plot, the effect of a CTR can be deduced by placing the radial force along the circumference of a CTR. For example, a CTR of d_{tr}=2 µm can be thought of as the circumferential distribution of a radial force on this plot when x is half of d_{tr}. That is 1 µm.

We can see that the radial force suddenly drops near x=1.5µm. It was experimentally shown [13,14] that when trap diameter d_{tr} is large, the radial force becomes weak and micro-bubbles could not be moved quickly along the transversal direction. Figure 7 supports these experiments. Now let’s try to analyze the behavior of a micro-bubble under a CTR. When the micro-bubble is located at the center of a CTR, it will not experience any transversal force. This is due to the symmetry of radial forces that originates along the circumference of the CTR. Now assume a micro-bubble is disturbed and it moves slightly along a direction we call **d⃗**. The force symmetry breaks down and the most prominent force change would be along **d⃗**. Now let’s think about only the radial force along **d⃗**. Before being disturbed, both ends along **d⃗** would have been pushing with a same force magnitude. However when disturbed, one end would be closer to the center of the micro-bubble and the other end would be further away. This result in one end pushing stronger and the other and pushing weaker. For example in Fig. 7, when the CTR trap diameter is 1µm and z location is at the origin, the radial force density acting on the micro-bubble would be about 10/*π*d_{tr} [pN/m, inset A]. However when disturbed by 0.2µm from the center, the end that becomes closer would push at about 5.5/*π*d_{tr} [pN/m, inset B] and the opposite end will push at about 11.5/*π*d_{tr} [pN/m, inset C]. Under this circumstance, the external forces act in favor of the disturbance and trapping is not possible. Generally, in Fig. 7, the region that shows positive force derivative is unstable trapping zone and negative force derivative region are stable trapping zone.

Figure 8 plots the sum of optical force and buoyant force using a scale of (F/|F|)**×**|F|^{0.3}. This scale was selected to amplify the change in directions. Figure 8 (a) and (b) are for the inverted and upright, respectively. From this, we can see that for trapping micro-bubbles, we should first place the laser focus at the same depth as the micro-bubble and gradually decrease the trap diameter just above the critical diameter. When the trap diameter becomes smaller than the critical diameter, the micro-bubble will be pushed into +z (inverted) or -z (upright) direction. The region towards the right of the red border drawn in Fig. 8 shows the trapping zone. When micro-bubbles are placed inside this region they will eventually reach axial force equilibrium. However, when they are outside this region, they will not be trapped and escape from the trap.

During the micro-bubble trapping, the required amount of laser power is also a concern. Figures 5 and 6 were calculated with a laser power of 10mW. Optical forces were significantly larger than the buoyant force of 0.0407 pN. Since optical forces are directly proportional to the laser power, for smaller laser power the trapping forces and trappable region becomes smaller. At certain limit, when the optical force is weaker than the buoyant force, the trapping will not be feasible. Figure 9 illustrates the axial force that is the sum of the trapping forces and buoyant force under different laser power. We can see that the axial force is proportional to the laser power. When the laser power is small such that the optical axial force is smaller than the buoyant force, axial trapping will not possible.

## 5. Conclusions

In this study, we have sown a method for the computation of optical forces on a micro-bubble under a tightly focused Gaussian beam. FDTD method was used with non-paraxial Gaussian beam representation. We have shown the relations between optical forces on a micro-bubble and the trap diameter of a CTR. We have concluded that there exists an optimal trap diameter that can result in the strongest trapping forces. These calculations can be used to better understand the physics behind the trapping and manipulation of micro-bubbles using optical tweezers.

## Acknowledgement

This work was supported by the Korea Science and Engineering Foundation (KOSEF) grant (No. R01-2007-000-11650-0) funded by the Korea government (MOST).

## References and links

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

**2. **Steven Chu, “Laser manipulation and atoms and particles,” Science **253**, 861–866 (1991). [CrossRef] [PubMed]

**3. **A. Ashkin, J. M. Dziedzic, and T. Yamane, “Optical trapping and manipulation of single cells using infrared laser beams,” Nature **330**, 769–771 (1987). [CrossRef] [PubMed]

**4. **In-Yong Park, Seung-Yong Sung, Jong-Hyun Lee, and Yong-Gu Lee, “Manufacturing micro-scale structures by an optical tweezers system controlled by five finger tips,” J. Micromech. Microeng. **17**, 82–89 (2007). [CrossRef]

**5. **Peter John Rodrigo, Lóránd Kelemen, Carlo Amadeo Alonzo, Ivan R. Perch-Nielsen, Jeppe Seidelin Dam, Pál Ormos, and Jesper Glückstad, “2D optical manipulation and assembly of shape complementary planar microstructures,” Opt. Express **15**, 9009–9014 (2007). [CrossRef] [PubMed]

**6. **A. Ashkin, “Forces of a single-beam gradient laser trap on a dielectric sphere in the ray optics regime,” Biophys. J. **61**, 569–582 (1992). [CrossRef] [PubMed]

**7. **Yasuhiro Harada and Toshimitsu Asakura, “Radiation forces on a dielectric sphere in the Rayleigh scattering regime,” Opt. Commun. **124**, 529–541 (1996). [CrossRef]

**8. **Wei Sun, Shi Pan, and Yuchi Jiang, “Computation of the optical trapping force on small particles illuminated with a focused light beam using a FDTD method,” J. Mod. Opt. **53**, 2691–2700, (2006). [CrossRef]

**9. **J. A. Stratton, *Electromagnetic Theory* (McGraw-Hill Inc., New York and London, 1941), Chap. 2.

**10. **Guoquan Zhou, Ruipin Chen, and Junlang Chen, “Propagation of non-paraxial nonsymmetrical vector Gaussian beam,” Opt. Commun. **259**, 32–39 (2006). [CrossRef]

**11. **Kane S. Yee, “Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media,” IEEE Trans. Antennas Propag. **AP-14**, 302–307 (1966).

**12. **Allen Taflove and Suans C. Hangess, *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, Third Edition (Artech House Inc., Norwood, MA, 2005).

**13. **P. H. Jones, E. Stride, and N. Saffari, “Trapping and manipulation of microscopic bubbles with a scanning optical tweezer,” Appl. Phys. Lett. **89**, 081113 (2006). [CrossRef]

**14. **N. B. Simpson, D. McGloin, K. Dholakia, L. Allen, and M. J. Padgett, “Optical tweezers with increased axial trapping efficiency,” J. Mod. Opt. **45**, 1943–1949 (1998). [CrossRef]

**15. **S. K. Mohanty, R. S. Verma, and P. K. Gupta, “Trapping and controlled rotation of low-refractive-index particles using dual line optical tweezers,” Appl. Phys. B87, 211–216 (2007). [CrossRef]