1. Introduction

Computation of the force of radiation on a given object, through evaluation of the electromagnetic field distribution according to Maxwell’s equations, followed by a direct application of the Lorentz law of force, has been described in Ref. [1]-[4]. In particular, for an isotropic, piecewise-homogeneous dielectric medium, the total force was shown to result from the force of the magnetic field acting on the induced bound current density, J _b × B= [(1 - 1/ε)∇ × H] × B, and from the force exerted by the electric component of the light field on the induced bound charge density at the interfaces between media of differing relative permittivity ε. The contribution of the E-field component of the Lorentz force is thus specified (Ref. [1]) by the force density ρ_b E = (ε ₀∇ ∙ E)E. (Note: As far as the total force and torque of radiation on a solid object are concerned, Barnett and Loudon Ref. [5] have recently shown that an alternative formulation of the Lorentz law - one that is often used in the radiation pressure literature - leads to exactly the same results. We discuss the relation between these two formulations in the Appendix, and extend the proof of their equivalence to the case of solid objects immersed in a liquid, which is a prime concern of the present paper.)

In a previous publication Ref. [6] we described the numerical implementation of the aforementioned approach to the computation of the Lorentz force, based on the Finite-Difference Time-Domain (FDTD) solution of Maxwell’s equations. (An alternative application of the FDTD method to problems of radiation pressure may be found in Ref. [7].) Our application of the FDTD method to the computation of the force exerted by a focused laser beam on a spherical dielectric particle immersed in a liquid showed that the forces experienced by the liquid layer immediately at the particle’s surface could impact the overall force experienced by the particle. A detailed discussion of the (conceptual) separation of the bound charges on the surface of the particle from the bound charges induced in the surrounding liquid at the solid-liquid interface is given in Ref. [2]. The analysis indicated that the contribution to the ρ_b E part of the Lorentz force by the component of the E-field perpendicular to the interface can be computed in different ways, depending on the assumptions made concerning the nature of the electromagnetic and hydrodynamic interactions between the solid particle and its liquid environment.

In the present study we apply the FDTD method to analyze the polarization dependence of the interaction between a focused laser beam and a small spherical particle trapped either in the air or in a liquid host medium (water). In the latter case, results from different methods of computing the contributions to the net force by bound charges at the solid-liquid interface will be compared, with the goal of quantifying the effects that might be possible to differentiate in experiments Ref. [8]-[9]. The polarization dependence of optical tweezers has been studied in the past, and the dependence of trap stiffness on polarization direction is well documented Ref. [10]-[13]. The goal of the present paper is not a re-evaluation of the existing models, but rather a demonstration of the applicability of our own new model Ref. [1]-[4] to the problem of trap stiffness anisotropy. Also, upon comparing our numerical results with experimental data, we gain further insight into the nature of the Lorentz force acting on solid objects in liquid environments, and identify a preferred interpretation of the proposed physical model. In a nutshell, there exists an ambiguity as to the nature of the effective force at the surface of the sphere, with at least three different models in contention. Carefully examining the trapping force on a dielectric sphere immersed in water provides enough information to establish one of the three competing models and rule out the other two.

The paper is organized as follows. Section 2 discusses example cases used to validate our numerical computations. In sections 3 and 4 we compare the computed stiffness of the trap for two orthogonal linear polarizations of a beam illuminating a micro-sphere suspended in air and in water, respectively. Summary and conclusions are presented in section 5.

 

Fig. 1. Surface force integrated over the left-half of the cylinder as function of the cylinder’s refractive index n_cyl . The incident plane-wave has vacuum wavelength λ₀ = 0.65μm, and the incidence medium’s refractive index is n_inc = 1.0 (solid lines) or n_inc = 1.3 (dashed lines). The exact results of the Lorentz-Mie theory (circles) are compared with those of our FDTD-based method (triangles and crosses). In the case of cylinder immersed in water, the charges induced at the solid-liquid interface were lumped together, then subjected to the average E-field at the interface; in other words, no attempts were made in these calculations to distinguish the force of the light’s E-field on the solid surface from that on the adjacent liquid. For each FDTD simulation the grid resolution Δ is indicated.

Download Full Size | PDF

2. Comparison to exact solutions

To validate our numerical procedures, we first consider the problem of computing the surface forces exerted by a plane-wave incident on a dielectric cylinder, with propagation direction being perpendicular to the cylinder axis. We consider in the YZ-plane of incidence a p-polarized plane-wave (E_y ,E_z ,H_x ), for which the electric-field component normal to the cylinder surface, E _⊥, is discontinuous; see Fig. 1. The free-space wavelength of the incident light is λ ₀ = 0.65μm and the radius of the cylindrical rod is r_cyl = 1 .0μm. The surface forces computed in our FDTD simulations can be compared with those obtained from the Lorentz-Mie theory of light scattering, Ref. [14]. In computing the Lorentz force of the E-field on the interfacial (bound) charges induced on the cylinder surface we used the average (Ref. [2]) of the E-field normal to the surface (the tangential component of the E-field is continuous across the interface); the surface-charge density is assumed to be proportional to the discontinuity of the E _⊥ field at the interface. While in the Lorentz-Mie theory E _⊥ at the surface is readily available, on the FDTD grid the cylinder surface is approximated as a discrete staircase, and the surface force density components ($F_y^{\mathit{\text{surf}}}$ ,$F_z^{\mathit{\text{surf}}}$ ) are evaluated along the coordinate axes, Ref. [6]. Figure 1 shows, for two different cases, the surface force integrated over one-half of the cylinder surface (180° ≤ θ ≤ 360°), as function of the refractive index n_cyl of the cylinder. The agreement between exact theory and numerical simulation is remarkable. Our numerical discretization error is of the order of 1-2% for n_cyl ranging from 1.5 to 3.4, consistent with the 30-60 points per wavelength discretization and first-order convergence due to staircase approximation of the geometry. For the high index-contrast case of n_inc /n_cyl = 1/3.4 the error (Fig. 1, Δ) was found to be dominated by the inadequate convergence to a time-harmonic state; the error was reduced (Fig. 1, ∇) when we increased the integration time.

In the second test problem we used the exact solutions for radiation pressure distribution in a solid dielectric prism illuminated by a Gaussian beam of light at Brewster’s incidence, Ref. [2].

Table 1. Definitions of the surface charge and surface force densities for methods I,II and III.

View Table | View all tables in this article

When the prism is immersed in water, the interfacial charges induced on the solid and liquid sides of the interface must be distinguished from each other. The interaction of these charged layers with the local E-field determines their contribution to the net force acting on the prism. Following the notation introduced in Ref. [2] we enumerate in Table 1 the three approaches to the computation of the surface force density acting on the solid side of the solid-liquid interface. Here ($F_{\parallel }^{\mathit{\text{surf}}}$,$F_{\perp }^{\mathit{\text{surf}}}$) are the surface force density components parallel and perpendicular to the solid’s surface; σ ₂ is the surface charge density belonging to the solid side; E _∥,E _1⊥,E _2⊥ are, respectively, the E-field components parallel to the interface, normal to the interface on the liquid side, and normal to the interface on the solid side; and ε ₀ E _g⊥ = ε ₀ ε ₁ E _1⊥ = ε0ε ₂ E _2⊥ is the electric displacement field D _∥ normal to the interface. Method III isolates the force acting on the solid side of the interface by introducing a small (artificial) gap between the solid and the liquid, then using the gap field E _g⊥ (derived from the continuity of the normal component of the displacement field D) to evaluate both the charge density σ ₂ and the average E-field that acts on the solid; this corresponds to Eq.(21) in Ref. [2]. While the conceptual introduction of a small gap at the solid-liquid interface is essential for the calculation of σ ₂, it gives rise to a mutual attractive force between the two charged layers thus pushed apart. This force, which is included in the total force experienced by the solid object as calculated by Method III, is probably inactive in practice and should be excluded. In Method II the contribution of this attractive force between the two charged layers (induced on the solid and liquid sides of the interface) is removed by using (E _1⊥ + E _2⊥)/2 as the effective perpendicular field acting on the surface charge density σ ₂ on the solid side of the interface; see Eq.(24) in Ref. [2]. In Method I the normal components of the E-field on the liquid and solid sides of the interface are used to compute the normal component of the average E-field as well as the induced charge density at the interface. In other words, the charges on the solid and liquid sides of the interface are lumped together, then subjected to the forces of the tangential E-field, which is continuous at the interface, and the perpendicular E-field, which is discontinuous at the interface (and, therefore, in need of averaging). Method I, which corresponds to Eq.(27) in Ref. [2], is also the method used in the case of a cylindrical rod immersed in water discussed earlier in conjunction with Fig. 1.

 

Fig. 2. Dielectric prism of refractive index n immersed in a host medium of refractive index n_inc .

Download Full Size | PDF

Table 2. Radiation force on the prism of Fig. 2, computed with the exact method and with the FDTD simulations. A p-polarized (H_x ,E_y ,E_z ) Gaussian beam illuminates the prism at Brewster’s angle θ_B . The net radiation force exerted on the prism is denoted by (F_y ,F_z ). The value of the mesh parameter Δ used in each simulation is indicated.

View Table | View all tables in this article

We use a Gaussian beam (λ ₀ = 0.65μm) having a full-width of 7.8μm at the 1/e point of the field amplitude, incident on the prism at Brewster’s angle θ_B =arctan(n/n_inc ); see Fig. 2. In the case of incidence from the free-space, n_inc = 1.0, the prism has refractive index n = 1.5, while for n_inc = 1.33 (water) the prism index is n = 1.995. The incident field amplitude’s peak value is E ₀ = 10³/n_inc V/m. Table 2 compares the total radiation force (consisting of the surface force density integrated over both surfaces, and the volume force density integrated over the prism volume) exerted on the prism, obtained with the exact method and with the FDTD simulations. As shown in Fig. 2, the total force has a component F_y directed along the bisector of the prism’s apex, and a second component F_z perpendicular to this bisector in the plane of incidence. (The exact value of F_z is zero in all cases.) For a prism immersed in water (n_inc = 1.33) the results obtained with the aforementioned methods I, II, and III of force computation are tabulated. The numerical results are generally in good agreement with the exact solutions; the largest disagreements occur in the case of Method I used in conjunction with the immersed prism, where (due to the specific set of parameters chosen for the simulation) the exact value of F_y happens to be so small (less than 10^-2 pN/m) that it falls below the numerical accuracy of our simulations. In the computations with method III the separation of the charges on the liquid and solid sides of the interface was simulated using an actual air-gap (width = 2Δ). The numerical and exact solutions in general are in good agreement with less than 1% error in the net force magnitude when Δ ≤ 20nm.

3. Single-beam optical trap in the air

We present computed results of the electromagnetic force exerted on a dielectric micro-sphere at and near the focus of a laser beam in the free space (n_inc = 1.0). The incident beam, obtained by focusing a λ ₀ = 532nm (or 1064nm) plane-wave through a 0.97NA, 5mm focal-length objective lens, propagates in the negative z-direction, as shown in Fig. 3. The plane-wave entering the objective’s pupil is linearly polarized along either the x-axis or the y-axis. The total power of the incident beam is P = ∫ S_zdxdy = 10W. Figure 4 shows the Poynting vector distribution (S-field) in the XZ-plane for a micro-sphere of refractive index n = 1.5 and diameter d = 460nm, offset by (x,y,z)-offset=(250,0,50)nm from the focal point of the y-polarized incident beam. Positive (negative) values of the z-offset represent the particle being displaced into the converging (diverging) half-space of the focused beam. Upon scattering from the micro-sphere, the large positive momentum acquired by the incident light along the x-axis, seen in Fig. 4, results in a net F_x force directed towards the beam axis. Figure 5 shows the computed (F_x ,F_z ) force components experienced by the micro-sphere as functions of the sphere center’s offset from the focal point of the lens (y-offset = 0). The top (bottom) row shows the case of an x-polarized (y-polarized) incident beam. Due to symmetry, y-offset is set to zero, and offsets along the positive half of the x-axis are the only ones considered. The y component of the force was found to be zero (within the numerical accuracy of the simulations) for these linearly-polarized beams. Figure 5 indicates that the particle is trapped laterally, but not vertically along the z-axis, since F_z < 0 for the entire range of offsets shown. (Note: A 460nm-diameter glass bead weighs approximately 1.0fN. If the focused laser beam shines on the bead from below, and if the laser power level is adjusted to a few microwatts, then the radiation pressure will work against the force of gravity to hold the bead in a stable trap along the z-axis.) The computed anisotropy of this trap in the lateral direction is s_l = 1 - (κ_x /κ_y ) = -0.15, where κ_x and κ_y are the trap stiffness coefficients along the x-axis given by ∂F_x /∂x for x- and y-polarized beams, respectively, Ref. [12]. The aforementioned s_l was computed at the x-offset value of 50nm, where z-offset ≈ 0μm is chosen to yield the maximum of F_x in the vicinity of the center of the small rectangles depicted in Fig. 5, where F_x is fairly insensitive to small variations of z. The computed stiffness anisotropy is plotted in Fig. 6 versus the particle diameter d for spherical beads having n = 1.5, trapped under a λ₀ = 1064nm focused beam through a 0.9NA objective.

 

Fig. 3. Computed distributions of the electric field intensity |E|² (units: [V ²/m ²]) in the xz- and yz-planes near the focus of a 0.9NA, f = 5.0mm, diffraction-limited objective. The λ₀ = 532nm plane-wave illuminating the entrance pupil of the lens is linearly polarized along the x-axis.

Download Full Size | PDF

 

Fig. 4. Distribution of the Poynting vector S (units: [W /m ²]) in and around a glass micro-sphere (n = 1.5, d = 460nm). The focused beam, obtained by sending a linearly-polarized plane-wave (polarization along y) through a 0.9NA objective, propagates along the negative z-axis. Sphere center offset from the focal point: (250,0,50)nm. The (S_x ,S_z ) vector-field is superimposed on the color-coded S_z plot on the right-hand side.

Download Full Size | PDF

 

Fig. 5. Plots of the net force components (F_x ,F_z ) experienced by a glass micro-sphere (d = 460nm, n = 1.5) versus the offset from the focal point in the xz-plane. The incidence medium is air, λ₀ = 532nm, the objective lens NA is 0.9, and the incident beam’s power is P = 1.0W. Top row: x-polarization, bottom row: y-polarization. The stiffness coefficients κ_x , κ_y are computed at the center of the small rectangles shown on the left-hand-side of the F_x plots.

Download Full Size | PDF

 

Fig. 6. Computed trap stiffness anisotropy s_l = 1 - (κ_x /κ_y ) versus particle diameter d, for micro-spheres of refractive index n = 1.5 trapped in the air with a λ₀ = 1064nm laser beam focused through a 0.9NA objective lens. The stiffness is computed at x-offset = 50nm, z-offset ≈ 0μm, where, for the chosen value of x-offset, the lateral trapping force F_x is at a maximum.

Download Full Size | PDF

For the offset ranges and particle diameters considered, the radiation force along the z-axis, F_z , was found to be negative (i.e., inverted traps are necessary to achieve stable trapping). The lateral trap anisotropy s_l , evaluated at x-offset = 50nm (with z-offset adjusted to yield maximum lateral trapping force F_x ), is seen in Fig. 6 to be positive for small particle diameters, negative when the particle size is comparable to the wavelength of the trap beam, and weakly positive for d > 1.4μm.

 

Fig. 7. (a) A linearly-polarized Gaussian beam having wavelength λ₀=532nm and 1/e (amplitude) radius r ₀= 4.0mm is focused through a 1.4NA oil immersion objective lens. The lens has focal length f=3.0mm and aperture radius (at the entrance pupil) R_a =2.85mm; the refractive index of the immersion oil is 1.47. A glass plate of the same index as the oil separates the oil from the water (n_water =1.33), where the focused spot is used to trap various dielectric beads. The marginal rays are lost by total internal reflection at the glass-water interface; there is also some degree of apodization due to Fresnel reflection at this interface, but the phase aberrations induced by the transition from oil/glass to water are ignored in our calculations. The neglect of the spherical aberrations thus induced is justifiable, so long as the trap is not too far from the oil/water interface, Ref. [12]. (b)-(d) Logarithmic plots of the intensity distribution for x-, y-, and z-components of the E-field at the focal plane. The relative peak intensities are |E_x |²:|E_y |²: |E_z |² ≈ 1000 : 9 : 200. The total intensity distribution at the focal plane, I(x,y) = |E_x |² +|E_y |² + |E_z |², (not shown) is elongated in the x-direction; the full-width at half-maximum intensity (FWHM) of the focused spot is 300nm along x and 196nm along y. The transmission efficiency of the oil/glass-to-water interface is 92.6%; that is, the overall loss of optical power to total-internal and Fresnel reflections at this interface is 7.4%.

Download Full Size | PDF

4. Single-beam optical trap in water

In computing the optical trap properties in a liquid host medium (refractive index n_inc = 1.33), a collimated laser beam having λ₀ = 532nm (or 1064nm) was focused through an NA ≈ 1.4 immersion objective designed for diffraction-limited focusing within an immersion liquid of refractive index n_oil = 1.47; see Figs. 7, 8. In our first set of simulations corresponding to λ₀ = 532nm, the spherical particle has d = 460nm, n = 1.5, and the total power of the incident beam is P = 1.0W. For Methods I, II, and III of computing the surface-force discussed in section 2, Figs. 9–11 show the net force components (F_x ,F_z ) exerted on the micro-sphere illuminated by x- or y-polarized light. For most of the offset range considered in Fig. 9, the F_x force component computed with method I for x-polarized light is opposite in direction to the F_x computed for y-polarized light, indicating the impossibility of lateral trapping with x-polarization. Such a marked difference in the behavior of F_x for different polarization states exhibited by method I, contradicts the experimental observations which show trapping (with similar strengths) for both polarization states. We conclude, therefore, that Method I is unphysical and must be abandoned. [The root of the problem with Method I can be traced to the lumping together of the solid and liquid charges at the interface, which weakens the negative contribution of $F_x^{\mathit{\text{surf}}}$ , thus allowing the positive contribution by the magnetic part of the Lorentz force (acting on the micro-sphere volume) to push the particle away from the focal point. In contrast, for y-polarized light, the contribution of the magnetic Lorentz force to F_x is negative, thus enabling lateral trapping.]

 

Fig. 8. Computed distributions of the electric field intensity |E|² (units: [V ²/m ²]) in the xz-and xz-planes near the focus of a λ₀ = 532nm beam in water. The focused spot is obtained by sending an x-polarized plane-wave through an oil-immersion ≈ 1.4NA objective. The oil and water are separated by a thin glass slide, index-matched to the immersion oil, Fig. 7.

Download Full Size | PDF

The force distributions computed with Method II (Fig. 10) and Method III (Fig. 11) show qualitatively similar behavior for the two polarization states, and indicate trapping along both x-and z-directions. F_x is strongest at lateral offset values on the order of one micro-sphere radius. While both methods II and III result in trapping of the micro-bead, the ratio of the maximum lateral restoring forces F_x (sampled inside the dashed rectangles in Figs. 10 and 11) for x- and y-polarized beams is found to be 0.92 for method II and 1.2 for method III. Thus, although Figs. 10 and 11 exhibit qualitatively similar behavior, their quantitative estimates of the restoring forces are sufficiently different to enable one to distinguish Method II from Method III based on experimentally determined values of stiffness anisotropy.

 

Fig. 9. Plots of the net force components (F_x ,F_z ), computed with Method I, for a glass micro-bead (d = 460nm, n = 1.5) versus the offset from the focal point. The host medium is water (n_inc = 1.33), λ₀ = 532nm, the objective lens NA is ≈ 1.4, and the incident beam’s power is P = 1.0W. Top row: x-polarization, bottom row: y-polarization. The non-trapping behavior of x-polarized beam, which is contrary to experimental observations, indicates the invalidity of Method I used in these calculations.

Download Full Size | PDF

Figure 12 shows the dependence of the trap stiffness anisotropy s_l = 1 -(κ_x /κ_y ) on bead diameter d, computed with Method II for polystyrene micro-beads (n = 1.57) trapped in water (n_inc = 1.33) under a λ₀ = 1064nm laser beam focused through an oil-immersion ≈ 1.4NA objective lens (with a glass slide used to separate oil from water, as shown in Fig. 7). For d < 850nm, the lateral trap stiffness is found to be smaller for the particle offsets along the polarization direction, that is, s_l > 0. For larger particles (850 < d < 1400nm) the trap stiffness in the x-direction becomes greater for x-polarization than that for y-polarization, i.e., s_l < 0. The trap stiffness anisotropy reaches a minimum before returning to positive values for d > 1400nm.

Superimposed on Fig. 12 are two sets of experimental data. The green triangles correspond to measurements carried out with a system similar to that depicted in Fig. 7, operating at λ₀ = 1064nm. In practice, this system suffered from chromatic (and possibly spherical) aberrations; defects that are not taken into account in our computer simulations. The agreement with theoretical calculations is good for the d = 1260nm particle, but not so good for d = 1510nm and d = 1900nm particles. The second set of experimental data (solid blue circles in Fig. 12) is from Table II of Rohrbach Ref. [9]. These were obtained with a 1.2NA water immersion objective, corrected for all aberrations and, therefore, presumably operating in the diffraction limit. All other experimental parameters were the same as those used in our simulations. (Considering the loss of marginal rays at the oil/water interface, apodization due to Fresnel reflection losses, and the relatively small beam diameter at the entrance pupil of the simulated lens, our focused spot should be fairly close to that used in Rohrbach’s experiments.)

 

Fig. 10. Plots of the net force components (F_x ,F_z ), computed with Method II, for a glass micro-bead (d = 460nm, n = 1.5) versus the offset from the focal point. The host medium is water (n_inc = 1.33), λ₀ = 532nm, the objective lens NA is ≈ 1.4, and the incident beam’s power is P = 1.0W. Top row: x-polarization, bottom row: y-polarization. The lateral trapping force F_x at the center of the small rectangle in the case of x-polarization is weaker than the corresponding F_x in the case of y-polarization (ratio = 0.92).

Download Full Size | PDF

 

Fig. 11. Plots of the net force components (F_x ,F_z ), computed with Method III, for a glass micro-bead (d = 460nm, n = 1.5) versus the offset from the focal point. The host medium is water (n_inc = 1.33), λ₀ = 532nm, the objective lens NA is ≈ 1.4, and the incident beam’s power is P = 1.0W. Top row: x-polarization, bottom row: y-polarization. The lateral trapping force F_x at the center of the small rectangle in the case of x-polarization is stronger than the corresponding F_x in the case of y-polarization (ratio = 1.2).

Download Full Size | PDF

 

Fig. 12. Trap stiffness anisotropy s_l = 1 - (κ_x /κ_y ) versus particle diameter d, computed with method II for polystyrene micro-beads (n = 1.57) trapped in water (n_inc = 1.33) under a λ₀ = 1064nm laser beam focused through an oil-immersion ≈ 14NA objective lens. Numerically, κ_x and κ_y were evaluated with a forward finite-difference approximation at x-offset = 50 - 100nm, offset discretization Δx = 50nm, and z-offset = z ₀ where F_z (x = 0,z ₀) = 0. The computed value of z ₀ for each particle is listed on the right-hand side. The solid triangles (green) represent experimental data obtained with a system similar to that depicted in Fig. 7, operating at λ₀ = 1064nm, and suffering from chromatic (and possibly spherical) aberrations. The solid circles (blue) represent experimental data from Table II of Rohrbach Ref. [9]. Although Rohrbach uses a 1.2NA water immersion objective in his experiments, the comparison is warranted here because our simulated 1.4NA focused beam, upon transmission from oil to water, loses its marginal rays and acquires characteristics that are not too far from those of a 1.2NA diffraction-limited focused spot. (Our experimental methodology is similar to that of Rohrbach, the only difference being that we use the oil-immersion objective depicted in Fig. 7(a), whereas Rohrbach uses a water-immersion objective.)

Download Full Size | PDF

The agreement between theory and experiment is remarkable in this case, except for the d = 1660nm particle. (We repeated the simulation for the d = 1660nm particle using the exact NA used in Rohrbach’s experiments; the resulting value of s_l , however, did not change very much.) While it is possible that the larger beads used in the aforementioned experiments have been described inaccurately (i.e., n and d deviating from the specifications), and while it is true that the aberrations of our own measurement system need to be understood and properly modeled, we cannot rule out the possibility that the radiation forces acting on the surrounding liquid (and ignored in the simulations) could have impacted the results of these measurements.

We mention in passing that, in the ”Rayleigh particle” limit (d ≫ λ₀/n_inc ), the Lorentz force density (ε ₀ ∇ ∙ E)E + J _b × B reduces, in the dipole approximation, to (p ∙ ∇)E + ∂ p/∂t × B, assuming linear dependence of the polarization p on the local E-field, namely, p(r,t) = α ₀ E(r,t), and neglecting the small terms due to the inhomogeneity of the magnetic field. For the smaller particles, the experimental data was found in Ref. [9] to be in good agreement with the computed results based on the Rayleigh-Gans approximation.

5. Concluding remarks

We have computed the radiation force of λ₀ = 532nm and λ₀ = 1064nm linearly polarized laser beams focused through air (NA = 0.9) and through water (NA ≈ 1.4) on various spherical micro-particles. For example, the restoring force acting on a small glass bead (d = 460nm, n = 1.5, λ₀ = 532nm) trapped in the air was found to be roughly 10-20% stronger when the polarization was aligned with, rather than perpendicular to, the particle offset direction. This and related predictions (described in Section 3) now await experimental verification to confirm the validity of the assumptions underlying our theoretical model.

The radiation force on a small glass bead (d = 460nm, n = 1.5, λ₀ = 532nm) immersed in water was evaluated using three different models for the partitioning of the radiation force between the solid particle and its liquid environment. Method I, which lumps together the induced bound charges on the solid and liquid sides of the interface, was ruled out on account of its unsubstantiated prediction of a strong polarization-dependence of the trap behavior. While methods II and III both predict the trapping of the bead under x- and y-polarized beams, the ratio of the maximum lateral restoring force for polarization directions parallel and perpendicular to the particle offset direction was found to be 0.92 for method II and 1.2 for method III. The dependence of the trap stiffness anisotropy on the physical model underlying the computation is thus seen to be large enough to enable one to accept or reject specific models. In the case of trapping polystyrene beads (n = 1.57) in water under a λ₀ = 1064nm focused beam, comparison with available experimental data strongly suggests that Method II is the correct method of computing the force of radiation on objects immersed in a liquid environment. (Method III yielded stiffness anisotropies of 0.07, -0.004, -0.1, -0.25, -0.34, -0.06, -0.007 for particle diameters d = 532nm, 690nm, 850nm, 1030nm, 1280nm, 1500nm, 1660nm, respectively.) For small particles (d less than or equal to λ₀), we found excellent agreement between Method II and experiment; see Fig. 12. The disagreement between theory and experiment for larger particles may indicate the need for more accurate measurements in this regime. Alternatively, the deviation of the observed stiffness anisotropy of larger particles from model calculations may be a hint that the radiation forces acting on the surrounding liquid and/or the attendant hydrodynamic effects should not be ignored in such calculations.

Appendix. Equivalence of total force (and torque) for two formulations of the Lorentz law

The Lorentz law of force, F = q(E+V×B), may be written in two different ways for a medium in which the macroscopic version of Maxwell’s equations are satisfied. The two formulations are:

(A1)$${\mathbf{F}}_1\left(r\right)=-(\nabla ·\mathbf{P})\mathbf{E}+\left(\frac{\partial \mathbf{P}}{\partial t}\right)\times \mathbf{B}$$

(A2)$${\mathbf{F}}_2\left(r\right)=(\mathbf{P}·\nabla )\mathbf{E}+\left(\frac{\partial \mathbf{P}}{\partial t}\right)\times \mathbf{B}$$

In an isotropic, homogeneous medium of dielectric constant ε where the electric displacement field D(r) = ε ₀ E(r)+P(r) = ε0ε E(r), one may write P(r) = ε ₀(ε - 1)E(r). In the absence of free charges and free currents in such a medium, Maxwell’s first equation, ∇ ∙ D(r) = 0, implies that the volume density of bound charges within the medium is zero, that is, ρ_b = - ∇ ∙ P(r) = 0. This leaves for the Lorentz force density in the bulk of the medium, according to Eq.(A1), only the magnetic contribution, namely, F ₁(r) = (∂ P/∂t) × B. Equation (A2), however, leads to an entirely different result. Using Maxwell’s equations in conjunction with the Lorentz law as expressed in Eq.(A2) yields:

(A3)$${\mathbf{F}}_2\left(r\right)=\frac{1}{4}{\epsilon }_0\left(\epsilon -1\right)\nabla \left({\mid E_x\mid }^2+{\mid E_y\mid }^2+{\mid E_z\mid }^2\right).$$

In other words, the Lorentz force density in the bulk of a homogeneous medium according to the second formulation is proportional to the gradient of the E-field intensity, irrespective of the state of polarization of the optical field.

At the boundaries of isotropic media, the two formulations again lead to substantially different force densities. The contribution of the B-field to the surface force is negligible as, for non-magnetic media (μ = 1), the B-field is continuous at the boundary; also, the discontinuity of ∂ P/∂t, if any, is finite. In contrast, the perpendicular component of the E-field at the boundary, E _⊥, has a sharp discontinuity, which results in a Dirac δ-function behavior for the E-field gradient ∇E _⊥(r). When evaluating the integral of (P ∙ ∇)E across the boundary between the free space and a medium of dielectric constant ε, one finds a force density (per unit interfacial area) F ₂(r) = 1/2P _⊥ (E _a⊥ - E _b⊥); here P _⊥ is the normal component of the medium’s polarization density at the surface, E _a⊥ is the perpendicular E-field in the free space region just outside the medium, and E _b⊥ is the perpendicular E-field within the medium just beneath the surface, Ref. [5]. Note that P _⊥, being identical to the bound surface charge density σ_b , may be derived from the discontinuity in the magnitude of E _⊥ at the surface, namely, P _⊥ = σ_b = ε ₀(E _a⊥ - E _b⊥).

In the first formulation based on Eq.(A1), the surface charge density σ_b derived from – (∇ ∙ P) is equal to P _⊥, as mentioned before. This σ_b , however, must be multiplied by the effective E-field at the boundary to yield the surface force density (per unit area) F ₁(r) = σ_b E. The tangential component E _∥ of the E-field is continuous at the boundary and, therefore, defined unambiguously. The perpendicular component, however, can be shown to be exactly equal to the average E _⊥ at the boundary, namely, E _⊥ = 1/2(E _a⊥ + E _b⊥).

We give two examples from electrostatics to demonstrate the difference between the two formulations of the Lorentz law as concerns the E-field contribution to the force within isotropic media. In the first example, shown on the left-hand side of Fig. A1, the medium is heterogeneous with a dielectric constant ε (x). Let the D-field be constant and oriented along the x-axis, that is, D(x) = D ₀ x̂ = [ε ₀ E_x (x) + P_x (x)]x̂ = ε ₀ ε(x)E_x (x)x̂. (Clearly ∇ ∙ D = 0 and ∇ × E = 0, as required by Maxwell’s electrostatic equations.) From the first formulation, Eq.(A1), we find:

(A4)$$F_{1x}^{\mathit{surface}}(x=0)=-\frac{1}{2}{P}_x(0)\left[E_0+E_x\left(0\right)\right]=-\left(\frac{D_0}{2{\epsilon }_0}\right)\left[1+\frac{1}{\epsilon \left(0\right)}\right]P_x\left(0\right);x=0$$

(A5)$$F_{1x}^{\mathit{bulk}}\left(x\right)=-\frac{E_x\left(x\right)d{P}_x\left(x\right)}{\mathit{dx}}=\left(\frac{1}{2{\epsilon }_0}\right)d\frac{{\left[D_0-P_x\left(x\right)\right]}^2}{\mathit{dx}};0<x<L$$

(A6)$$F_{1x}^{\mathit{surface}}\left(x=L\right)=\left(\frac{D_0}{2{\epsilon }_0}\right)\left[1+\frac{1}{\epsilon \left(L\right)}\right]P_x\left(L\right);x=L$$

The second formulation, Eq.(A2), yields:

(A7)$$F_{2x}^{\mathit{surface}}\left(x=0\right)=\frac{1}{2}{P}_x(0)\left[E_x\left(0\right)-E_0\right]=-\left(\frac{D_0}{2{\epsilon }_0}\right)\left[1-\frac{1}{\epsilon \left(0\right)}\right]P_x\left(0\right);x=0$$

(A8)$$F_{2x}^{surface}(x)=P_x(x)d{E}_x(x)/dx=-(1/2{\epsilon }_0)d{P}_x^2(x)/dx;0<x<L$$

(A9)$$F_{2x}^{\mathit{surface}}\left(x=L\right)=\left(\frac{D_0}{2{\epsilon }_0}\right)\left[1-\frac{1}{\epsilon \left(L\right)}\right]P_x\left(L\right);x=L$$

Clearly the two formulations predict different distributions for the force density, both at the surface and in the bulk. However, the total force - obtained by adding the surface contributions to the integral of the force density within the bulk- turns out to be exactly the same in the two formulations (F ^total = 0).

In the second example, shown on the right-hand-side of Fig. A1, the medium is homogeneous (dielectric constant=ε), and the D-field profile is D(r) = (D ₀ r ₀/r)$\widehat{\theta }$. In the first formulation, based on Eq.(A1), the bulk force is zero, and the surface force density is

(A10)$$F_2^{\mathit{surface}}\left(r,\theta ={\theta }_0,1\right)=\pm \left(\frac{1}{2{\epsilon }_0}\right)\left(1-\frac{1}{{\epsilon }^2}\right){\left(\frac{D_0{r}_0}{r}\right)}^2\hat{\theta };\theta ={\theta }_0,{\theta }_1$$

In the second formulation, Eq.(A2) yields:

(A11)$$F_2^{\mathit{surface}}\left(r,\theta ={\theta }_0,1\right)=\pm \left(\frac{1}{2{\epsilon }_0}\right){\left(1-\frac{1}{\epsilon }\right)}^2{\left(\frac{D_0{r}_0}{r}\right)}^2\hat{\theta };\theta ={\theta }_0,{\theta }_1$$

(A12)$$F_{2x}^{\mathit{bulk}}(r,\theta )=-\left(\frac{1}{{\epsilon }_0\epsilon }\right)\left(1-\frac{1}{\epsilon }\right)\left(\frac{D_0^2{r}_0^2}{r^3}\right)\hat{\mathit{r}};{\theta }_0<\theta <{\theta }_1$$

 

Fig. A1. Examples of electrostatic field distributions in two dielectric media. On the left-hand side, the medium’s dielectric constant ε(x) varies along the x-axis while the electric displacement D = D ₀ x̂ remains constant. On the right-hand side,ε is constant within the crescent-shaped medium (r ₀ < r < r ₁, θ ₀ < θ < θ ₁), while the azimuthally oriented D(r) = (D ₀ r ₀/r)$\widehat{\theta }$ decreases with the inverse of the radius r.

Download Full Size | PDF

Once again, the force distributions in the two formulations are quite different, but the total force is exactly the same, that is,

(A13)$${\mathbf{F}}^{\mathit{total}}=-\left(\frac{1}{{\epsilon }_0}\right)\left(1-\frac{1}{{\epsilon }^2}\right){\left(\frac{D_0{r}_0}{r}\right)}^2\sin \left[\frac{\left({\theta }_1-{\theta }_0\right)}{2}\right]\hat{\mathbf{r}}.$$

Next we prove the generality of the above results by showing that the total force obtained from Eqs. (A1) and (A2) is the same under all circumstances. The proof of equivalence between the two formulations with regard to total force (and also total torque) is due to Barnett and Loudon Ref. [5], who also clarified the role of surface forces in the second formulation. In what follows we outline this proof of equivalence along the same lines as originally suggested by Barnett and Loudon Ref. [5], then extend their results to objects that are immersed in a liquid (or surrounded by an isotropic, homogeneous medium of a differing index of refraction).

Consider the difference between F ₁(r)and F ₂(r), written as

(A14)$$\mathbf{\Delta F}\left(\mathbf{r}\right)={\mathbf{F}}_2\left(\mathbf{r}\right)-{\mathbf{F}}_1\left(\mathbf{r}\right)=\left(\mathbf{P}\bullet \nabla \right)\mathbf{E}+\left(\nabla \bullet \mathbf{P}\right)\mathbf{E}=\partial \frac{\left(P_x\mathbf{E}\right)}{\partial x}+\partial \frac{\left(P_y\mathbf{E}\right)}{\partial y}+\partial \frac{\left(P_z\mathbf{E}\right)}{\partial z}$$

When ∫ΔF(r)dxdydz over the volume of an object is being evaluated, individual terms on the right-hand-side of Eq.(A14), being complete differentials, produce expressions such as (P_x E) |_{x_max} - (-P_x E) |_{x_min} after a single integration over the proper variable. These expressions then reduce to zero because x_min ,x_max are either in the free-space region outside the object, where P = 0, or they are outside the beam’s boundary, where E = P = 0. Either way, the volume integral of ΔF(r)evaluates to zero, yielding ∫F ₁(r)dxdydz = ∫F ₂(r)dxdydz.

If the object of interest happens to be immersed in (or surrounded by) a medium other than the free space, where, for example, the condition (P_x E) |_{x_min} = 0 cannot be ascertained, one must proceed to assume the existence of a narrow gap between the object and its surroundings, as shown in Fig. A2. Under such circumstances, the radiation force exerted on the boundary of the object should be computed using E _g(r), the gap field, rather than E _a(r), the field in the immersion medium just outside the object. However, the introduction of the gap creates oppositely charged layers (across the gap), and the force of attraction between these proximate charges accounts for the difference between F ^surface(r) computed using E _g(r) and E _b(r) on the one hand, and that computed using E _a(r) and E _b(r) on the other hand. Once the effect of this attractive force on the object is discounted, the remaining force may be computed by ignoring the gap field and assuming, for the first formulation, that the effective E-field is 1/2(E _a + E _b), while, for the second formulation, the E-field discontinuity is E _a - E _b. [Note: In the first formulation, using 1/2(E _a + E _b) for the effective field leads to Method II discussed in Section 2, whereas using 1/2(E _a + E _b) leads to Method III.]

 

Fig. A2. When the object is surrounded by a medium other than the free space, a narrow gap may be imagined to exist between the object and its surroundings, so that, within the gap, the polarization density P(r) may be set to zero. The integration boundary is then placed within the gap to ensure that the integrated force densities corresponding to Eqs. (A1) and (A2) lead to the same total force.

Download Full Size | PDF

Finally, we show that the total torque T ^total exerted on a given object is independent of whether F ₁(r)or F ₂(r) is used to compute the torque. In the following derivation, the fact that P(r) = 0 when r is outside the proper boundaries of the object is used in both steps. In the first step, terms containing the integrals of complete differentials are omitted. (As before, such integrals reduce to expressions that vanish outside the object’s boundaries.) In the second step, integration by parts is used to simplify the integrands.

(A15)$$\Delta {\mathbf{T}}^{\mathit{total}}=\iiint \mathbf{r}\times \mathbf{\Delta F}\left(r\right)\mathit{dxdydz}$$
$$=\iiint \mathbf{r}\times [\partial \left(P_x\mathbf{E}\right)+\partial x+\frac{\partial \left(P_y\mathbf{E}\right)}{\partial y}+\frac{\partial \left(P_z\mathbf{E}\right)}{\partial z}]\mathit{dxdydz}$$
$$=\iiint \{[y\partial \frac{\left(P_y{E}_z\right)}{\partial y}-z\partial \frac{\left(P_z{E}_y\right)}{\partial z}]\hat{\mathbf{x}}+[z\partial \frac{\left(P_z{E}_x\right)}{\partial z}-x\partial \frac{\left(P_x{E}_z\right)}{\partial x}]\hat{\mathbf{y}}+[x\partial \frac{\left(P_x{E}_y\right)}{\partial x}-y\partial \frac{\left(P_y{E}_y\right)}{\partial z}]\hat{\mathbf{z}}\mathit{dxdydz}$$
$$=-\iiint [\left(P_y{E}_z-P_z{E}_z\right)\hat{\mathbf{x}}+\left(P_z{E}_x-P_x{E}_z\right)\hat{\mathbf{y}}+\left(P_x{E}_y-P_y{E}_x\right)\hat{\mathbf{z}}]\mathit{dxdydz}$$
$$=\iiint (\mathbf{E}\times P)\mathit{dxdydz}$$

In an isotropic medium E and P will be parallel to each other and, therefore, their cross-product appearing on the right-hand-side of Eq.(A15) will vanish. Consequently ΔT ^total = 0, completing the proof that both formulations of the Lorentz law yield the same overall torque on the object.

It is remarkable that the two formulations in Eqs. (A1) and (A2) yield the same total force (and torque) exerted on a given object under quite general circumstances. In fact, until recently, the present authors had relied solely on the first formulation, assuming (falsely) that the second formulation is incapable of producing consistent answers for certain well-defined problems Ref. [1]-[4]. This situation has now changed with the equivalence proof provided by Barnett and Loudon Ref. [5]. It appears, therefore, that the classical electromagnetic theory (with its reliance on the macroscopic properties of matter, such as polarization density P, displacement field D, dielectric constant ε, etc.) cannot decide which formula, if any, provides the correct force distribution throughout the object (even as both formulas yield the correct values for total force and torque). A resolution of this interesting problem may have to await the verdict of future experiments.

Acknowledgments

We are grateful to Rodney Loudon and Stephen Barnett for illuminating discussions and, in particular, for sharing the preprint of their latest manuscript, Ref. [5], with us. Thanks are also due to Alexander Rohrbach for his extensive comments on an early draft of this paper, and to Koen Visscher for granting us access to his laboratory. This work has been supported by the Air Force Office of Scientific Research (AFOSR) contracts FA9550-04-1-0213, FA9550-04-1-0355, and F49620-03-1-0194.

References and links

1. M. Mansuripur, “Radiation Pressure and the linear momentum of the electromagnetic field,” Opt. Express 12, 5375–5401 (2004), http://www.opticsexpress.org/abstract.cfm?id=81636. [CrossRef]   [PubMed]  

2. M. Mansuripur, A.R. Zakharian, and J.V. Moloney, “Radiation Pressure on a dielectric wedge,” Opt. Express 13, 2064–2074 (2005), http://www.opticsexpress.org/abstract.cfm?id=83011. [CrossRef]   [PubMed]  

3. M. Mansuripur, “Radiation Pressure and the linear momentum of light in dispersive dielectric media,” Opt. Express 13, 2245–2250 (2005), http://www.opticsexpress.org/abstract.cfm?id=83032. [CrossRef]   [PubMed]  

4. M. Mansuripur, “Angular momentum of circularly polarized light in dielectric media,” Opt. Express 13, 5315–5324 (2005), http://www.opticsexpress.org/abstract.cfm?id=84895. [CrossRef]   [PubMed]  

5. S. M. Barnett and R. Loudon, “On the electromagnetic force on a dielectric medium,” submitted to J. Phys. B: At. Mol. Phys. (January 2006).

6. A.R. Zakharian, M. Mansuripur, and J.V. Moloney, “Radiation Pressure and the distribution of the electromagnetic force in dielectric media,” Opt. Express 13, 2321–2336 (2005), http://www.opticsexpress.org/abstract.cfm?id=83272. [CrossRef]   [PubMed]  

7. R. Gauthier, “Computation of the optical trapping force using an FDTD based technique,” Opt. Express 13, 3707–3718 (2005), http://www.opticsexpress.org/abstract.cfm?id=83817. [CrossRef]   [PubMed]  

8. A. Rohrbach and E.H.K. Stelzer, “Three-dimensional position detection of optically trapped dielectric particles,” J. Appl. Phys. 91, 5474–5488 (2002). [CrossRef]  

9. A. Rohrbach, “Stiffness of Optical Traps: Quantitative agreement between experiment and electromagnetic theory,” Phys. Rev. Lett. 95, 168102 (2005). [CrossRef]   [PubMed]  

10. W.H. Wright, G.J. Sonek, and M.W. Berns, “Radiation trapping forces on microspheres with optical tweezers,” Appl. Phys. Lett. 63, 715–717 (1993). [CrossRef]  

11. W.H. Wright, G.J. Sonek, and M.W. Berns, “Parametric study of the forces on microspheres held by optical tweezers,” Appl. Opt. 33, 1735–1748 (1994). [CrossRef]   [PubMed]  

12. A. Rohrbach and E. H. K. Stelzer, “Trapping forces, force constants, and potential depths for dielectric spheres in the presence of spherical aberrations,” Appl. Opt. 41, 2494–2507 (2002). [CrossRef]   [PubMed]  

13. 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?id=80240. [CrossRef]   [PubMed]  

14. P.W. Barber and S.C. Hill, Light Scattering by Particles: Computational Methods (World Scientific Publishing Co. 1990).