2326 J. Opt. Soc. Am. A/Vol. 23, No 9/September 2006 theorem of the Hankel function to each individual par The validation of the previous equations has been p ticle formed by comparing the predicted total field scattered by multiple particles with that obtained by the commercial Hh (elpe-P l)e,'p=> Jn(klPe-paDeinetg package CST MICROWAVE STUDIO. The comparison has been done with a number of particles, ranging from one to five submitted to one, two, and three incident plane waves. In XHn-n, (Alpp-Pale-i(n-n)pg.(9) all cases, the fields(E2, Hr, Hy) have been monitored and verified to agree well, both inside and outside the dielec In practice, the infinite sums are truncated at n tric cylinders. In addition, the force has been compared ax(n)(2+x-2+N), where N has to be properly chosen with the one obtained by the method presented in Ref 12 Ls a function of the number of particles in the system, in the way shown in Refs. 28 and 29, where bound cur- their size, permittivity, and separation, in order to achieve good convergence(values of N are specified in ticle. The agreement between the two calculation meth each case hereafter). The solution of the system of Eq (8) ods has been found to be excellent in all the situations requires the inversion of a square matrix of dimensions studied hereafter. a noticeable advantage of the Maxwell X(2N+1)=21L for N=10. For our purpose stress tensor approach, however, is that it requires fewer Gauss inversion is used, but more efficient algorithms can integration points, since the force is expressed as a line be implemented to improve the computation efficiency of integral instead of a surface integral. On the other hand the calculation. It should be noted that such matrix size is the method of bound currents and charges has the advan significantly smaller than the one obtained from a DDa tage of computing the local force density within the par method applied to particles of similar properties. Hence ticle, which is not accessible from the Maxwell stress ap- losing the ability of modeling arbitrary shapes and per- mittivities, we gain a considerable saving in computer memory and computational time. Once the w a) coefficients are solved for asia) are com- puted from Eq(6), and Eq) (p)is obtained from Eq (5). 3. EVIDENCE OF THE TRAPPING FORCE Finally, the total scattered field is expressed as Trapping of particles has been experimentally verified in Refs. 6 and 7 with a set of three incident gaussian beams E∞a(p)=∑Eq(p (10) The corresponding conclusions are here confirmed theo- retically on a simplified system of cylindrical particles submitted to three plane waves of identical amplitude The interior fields are obtained in a similar manner but Eo=1(V/m)propagating in the (xy) plane with angles are not required for the computation of the force on the f 7/2, 7/6, 11/6)(rad). The three incident plane waves cylinders in have a similar wavelength in free space of \=546 nm, and Knowing the scattered field, and therefore the total the particles are polystyrene cylinders of permittivity e fields (E, H)in outside all cylinders, we c 2.56E0 embedded in water(e=1.69Eo). The radius of the pute the time-averaged Lorentz force from the Maxwell's particles used in Ref 7 is 1 um, which is not a constraint stress tensor as in our case since size is not a limitation. However. the traps created by the three in-plane incident waves are F= del-Re((En)E*)+Re((Hn)H, maller than 1 um, as can be seen from the background pattern in Fig. 1, and therefore 1 um particles would feel the effect of multiple traps. To keep the particle size of the -E·E*n--(H·H*n same order as the traps, we resort to particles of radii =0.15A=819 nm and a=03A=1638mm. In Eq(11), Rel) denotes the real-part operator, n=i Y The trapping force on a single particle is shown in Fig where the background field pattern shows the interfer- he normal vector to the surface of the cylinder, and the ence of the three incident plane waves(only a relevant closed contour C surrounds the particle from the outside portion of the (xy) plane is shown, since the pattern is pe- thus requiring the knowledge of the fields outside the par- riodic in both directions) and the arrows show the force ticle only. In this case, the permittivity e and permeability acting on the particle as it sweeps the(xy) plane. The ar- u are those of the background medium. Equation(11)has rows in Fig. 1 show both the relative magnitude and di lready been extensively used in the literature but often rection of the force, where the starting point of each arrow applied to small particles. In such cases, gradient and corresponds to the location of the center of the cylinder at scattering forces can be separated, and a series of ap- which the computation is performed. It can immediately proximations can be used to simplify Eq (11)into expres- be seen that the trapping ability of the incident pattern sions with direct physical interpretations. However, in the significantly different for the two different sizes of pa case of large particles, or when many closely packed pal ticles: The high field intensity at(, y)=(0, 0)is a point of ticles are considered, no systematic approximation of Eq. stable equilibrium for a=0.15,, while it is a point of un 11)can be used. As we shall see hereafter, the results of stable equilibrium for a=0.3. It should be noted that a the force on large particles or on many particles do not al- similar behavior was found in Ref. 30 for a particle in the ways follow the intuition built from the study of small Rayleigh regime(a<X/20)as a function of its refractivetheorem of the Hankel function to each individual par￾ticle: Hn

1
k

−
pein
lp = n=− + Jn
k

−
qeinlq Hn−n

1
k
p −
qe−i
n−n
pq.
9 In practice, the infinite sums are truncated at N =max
n
− +→−N +N, where N has to be properly chosen as a function of the number of particles in the system, their size, permittivity, and separation, in order to achieve good convergence (values of N are specified in each case hereafter). The solution of the system of Eq. (8) requires the inversion of a square matrix of dimensions L
2N+1=21L for N=10. For our purpose, a direct Gauss inversion is used, but more efficient algorithms can be implemented to improve the computation efficiency of the calculation. It should be noted that such matrix size is significantly smaller than the one obtained from a DDA method applied to particles of similar properties. Hence, losing the ability of modeling arbitrary shapes and per￾mittivities, we gain a considerable saving in computer memory and computational time. Once the wn
q coefficients are solved for, an s
q are com￾puted from Eq. (6), and Escat
q

is obtained from Eq. (5). Finally, the total scattered field is expressed as Escat

=q=1 L Escat
q

.
10 The interior fields are obtained in a similar manner but are not required for the computation of the force on the cylinders in our case. Knowing the scattered field, and therefore the total fields
E,H in the region outside all cylinders, we com￾pute the time-averaged Lorentz force from the Maxwell’s stress tensor as10,27 F =  C d

2 Re
E · nˆ E * +  2 Re
H · nˆ Hˆ * −
4
E · E * nˆ −  4
H · H * nˆ .
11 In Eq. (11), Re{ } denotes the real-part operator, nˆ =
ˆ is the normal vector to the surface of the cylinder, and the closed contour C surrounds the particle from the outside, thus requiring the knowledge of the fields outside the par￾ticle only. In this case, the permittivity
and permeability  are those of the background medium. Equation (11) has already been extensively used in the literature but often applied to small particles. In such cases, gradient and scattering forces can be separated, and a series of ap￾proximations can be used to simplify Eq. (11) into expres￾sions with direct physical interpretations. However, in the case of large particles, or when many closely packed par￾ticles are considered, no systematic approximation of Eq. (11) can be used. As we shall see hereafter, the results of the force on large particles or on many particles do not al￾ways follow the intuition built from the study of small particles. The validation of the previous equations has been per￾formed by comparing the predicted total field scattered by multiple particles with that obtained by the commercial package CST MICROWAVE STUDIO. The comparison has been done with a number of particles, ranging from one to five, submitted to one, two, and three incident plane waves. In all cases, the fields
Ez ,Hx ,Hy have been monitored and verified to agree well, both inside and outside the dielec￾tric cylinders. In addition, the force has been compared with the one obtained by the method presented in Ref. 12 in the way shown in Refs. 28 and 29, where bound cur￾rents and charges are used within the volume of the par￾ticle. The agreement between the two calculation meth￾ods has been found to be excellent in all the situations studied hereafter. A noticeable advantage of the Maxwell stress tensor approach, however, is that it requires fewer integration points, since the force is expressed as a line integral instead of a surface integral. On the other hand, the method of bound currents and charges has the advan￾tage of computing the local force density within the par￾ticle, which is not accessible from the Maxwell stress ap￾proach. 3. EVIDENCE OF THE TRAPPING FORCE Trapping of particles has been experimentally verified in Refs. 6 and 7 with a set of three incident Gaussian beams. The corresponding conclusions are here confirmed theo￾retically on a simplified system of cylindrical particles submitted to three plane waves of identical amplitude E0=1
V/m propagating in the
xy plane with angles /2,7/6,11/6 (rad). The three incident plane waves have a similar wavelength in free space of =546 nm, and the particles are polystyrene cylinders of permittivity
c =2.56
0 embedded in water

=1.69
0. The radius of the particles used in Ref. 7 is 1 m, which is not a constraint in our case since size is not a limitation. However, the traps created by the three in-plane incident waves are smaller than 1 m, as can be seen from the background pattern in Fig. 1, and therefore 1 m particles would feel the effect of multiple traps. To keep the particle size of the same order as the traps, we resort to particles of radii a =0.15 =81.9 nm and a=0.3 =163.8 nm. The trapping force on a single particle is shown in Fig. 1 where the background field pattern shows the interfer￾ence of the three incident plane waves (only a relevant portion of the
xy plane is shown, since the pattern is pe￾riodic in both directions) and the arrows show the force acting on the particle as it sweeps the
xy plane. The ar￾rows in Fig. 1 show both the relative magnitude and di￾rection of the force, where the starting point of each arrow corresponds to the location of the center of the cylinder at which the computation is performed. It can immediately be seen that the trapping ability of the incident pattern is significantly different for the two different sizes of par￾ticles: The high field intensity at
x,y=
0,0 is a point of stable equilibrium for a=0.15 , while it is a point of un￾stable equilibrium for a=0.3 . It should be noted that a similar behavior was found in Ref. 30 for a particle in the Rayleigh regime
a /20 as a function of its refractive index. 2326 J. Opt. Soc. Am. A/Vol. 23, No. 9/September 2006 Grzegorczyk et al.