## Abstract

Using numerical solutions of Maxwell’s equations in conjunction with the Lorentz law of force, we compute the electromagnetic force distribution in and around a dielectric micro-sphere trapped by a focused laser beam. Dependence of the optical trap’s stiffness on the polarization state of the incident beam is analyzed for particles suspended in air or immersed in water, under conditions similar to those realized in practical optical tweezers. A comparison of the simulation results with available experimental data reveals the merit of one physical model relative to two competing models; the three models arise from different interpretations of the same physical picture.

© 2006 Optical Society of America

## 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** = (*ε*
_{0}∇ ∙ **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.

## 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} = 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].

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; *σ*
_{2} 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 *ε*
_{0}
*E*
_{g⊥} = *ε*
_{0}
*ε*
_{1}
*E*
_{1⊥} = *ε*0*ε*
_{2}
*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 *σ*
_{2} 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 *σ*
_{2}, 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 *σ*
_{2} 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.

We use a Gaussian beam (*λ*
_{0} = 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*
_{0} = 10^{3}/*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 Δ ≤ 20

*nm*.

## 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 *λ*
_{0} = 532*nm* (or 1064*nm*) plane-wave through a 0.97*NA*, 5*mm* 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*_{z}*dxdy* = 10*W*. 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* = 460*nm*, 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.0*fN*. 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 50

*nm*, 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 λ

_{0}= 1064

*nm*focused beam through a 0.9

*NA*objective.

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 = 50*nm* (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*.

## 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 λ_{0} = 532*nm* (or 1064*nm*) 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 λ_{0} = 532*nm*, the spherical particle has *d* = 460*nm*, *n* = 1.5, and the total power of the incident beam is *P* = 1.0*W*. 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.]

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.

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 λ_{0} = 1064*nm* laser beam focused through an oil-immersion ≈ 1.4*NA* 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* < 1400*nm*) 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* > 1400*nm*.

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 λ_{0} = 1064*nm*. 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* = 1260*nm* particle, but not so good for *d* = 1510*nm* and *d* = 1900*nm* 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.2*NA* 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.)

The agreement between theory and experiment is remarkable in this case, except for the *d* = 1660*nm* particle. (We repeated the simulation for the *d* = 1660*nm* 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* ≫ λ_{0}/*n*_{inc}
), the Lorentz force density (*ε*
_{0} ∇ ∙ **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*) = *α*
_{0}
**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 λ_{0} = 532*nm* and λ_{0} = 1064*nm* 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* = 460*nm*, *n* = 1.5, λ_{0} = 532*nm*) 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* = 460*nm*, *n* = 1.5, λ_{0} = 532*nm*) 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 λ_{0} = 1064*nm* 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 λ_{0}), 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:

In an isotropic, homogeneous medium of dielectric constant *ε* where the electric displacement field **D**(**r**) = *ε*
_{0}
**E**(**r**)+**P**(**r**) = *ε*0*ε*
**E**(**r**), one may write **P**(**r**) = *ε*
_{0}(ε - 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**
_{1}(**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:

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**
_{2}(**r**) = 1/2*P*
_{⊥} (**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}
= *ε*
_{0}(*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**
_{1}(**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*
_{0}
**x̂** = [*ε*
_{0}
*E*_{x}
(*x*) + *P*_{x}
(*x*)]**x̂** = *ε*
_{0}
*ε*(*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:

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

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*
_{0}
*r*
_{0}/*r*)*$\widehat{\theta}$*. In the first formulation, based on Eq.(A1), the bulk force is zero, and the surface force density is

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

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

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**
_{1}(**r**)and **F**
_{2}(**r**), written as

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**
_{1}(**r**)*dxdydz* = ∫*F*
_{2}(**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.]

Finally, we show that the total torque **T**
^{total} exerted on a given object is independent of whether **F**
_{1}(**r**)or **F**
_{2}(**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.

$$\phantom{\rule{2.5em}{0ex}}=\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}\text{}$$

$$\phantom{\rule{2.5em}{0ex}}=\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}$$

$$\phantom{\rule{2.5em}{0ex}}=-\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}$$

$$\phantom{\rule{2.5em}{0ex}}=\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).