The derivation of the scattering force and the gradient force on a spherical particle due to an electromagnetic wave often invokes the Clausius–Mossotti factor, based on an ad hoc physical model. In this article, we derive the expressions including the Clausius–Mossotti factor directly from the fundamental equations of classical electromagnetism. Starting from an analytic expression for the force on a spherical particle in a vacuum using the Maxwell stress tensor, as well as the Mie solution for the response of dielectric particles to an electromagnetic plane wave, we derive the scattering and gradient forces. In both cases, the Clausius–Mossotti factor arises rigorously from the derivation without any physical argumentation. The limits agree with expressions in the literature.

1. Overview

Recently, we made a theoretical study of a system to sort submicrometer dielectric spheres in the interference field of a laser in slowly flowing air [1, 2]. In the course of that project, we derived the scattering and gradient forces rigorously from Maxwell’s equations. The derivation was too detailed for that article, so we are presenting it here. The results agree with expressions from Harada and Asakura [3], as we show in the Appendix. In addition, we present the code we used to derive the force on a spherical particle used in [1, 2].

In this article, we develop code to generate the Mie scattering coefficients and the stress tensor formulas, and we combine them to form first the scattering force and then the gradient force. The scattering force comes first because it requires an incident plane wave, whereas the gradient force requires an incident standing wave that is a little more difficult to set up.

The solution for the response of a spherical dielectric particle in a vacuum was given more than 100 years ago [4]. The problem has been studied extensively by Bohren and Huffman [5]. Our formulation follows the textbook of Zangwill [6]. Since the problem is so well studied, we go directly to the solution. We omit an implicit time dependence in the plane wave traveling in the positive direction with the electric field linearly polarized in the direction. The electric field and magnetic field are thereby given by

where and are scalar functions related to the transverse electric and transverse magnetic parts of the solution, respectively, is the speed of light, and is a point in space. The International System of Units (SI) is used; in these units has the same dimensions as . Compared to the equations in Zangwill, the overall sign differs here, and we have inserted a factor of to simplify the implementation. We can obtain from by the substitution ; we will use this relationship later to simplify the calculation. The scalar function is given by

Here, and are the usual spherical coordinates, and the are associated Legendre polynomials. If , where is a spherical Bessel function, the expressions describe the incident fields for a plane wave in the case of an incident plane wave with spatial dependence . The sum of the incident and scattered fields is given by

for , where is a spherical Hankel function. For , we substitute and in the previous two expressions.

Explicit forms for the Mie coefficients and are given here for a particle with index of refraction and radius . We restrict attention to the case of a nonmagnetic dielectric sphere:

where and are Riccati–Bessel functions; a prime denotes differentiation with respect to the argument. The total scattering cross section is given by

written here in dimensionless form by incorporating the geometric cross section of the spherical particle; this form is commonly denoted by .

1.1 Mie Coefficients and Cross Sections

The function calculates the Mie coefficients and using equations (5) and (6); the table of those coefficients is . The functions and are and , respectively. The variable (with maximum value ) is an index used for Legendre coefficients in physics. The variable is , the product of the wavevector and the particle radius ; ranges from to with step size . The function returns a list pairing with the cross section normalized to the geometry. A simpler three-argument function calls a five-argument function. The parameter helps to tell how many terms to calculate to achieve convergence; 1.5 seems to work well, but the reader may wish to test this.

Next, we plot the Mie cross section, similar to the one found in Zangwill [6] and Bohren and Huffman [5]. This is slow: it took 442 seconds on a 3.7 GHz personal computer. The code could be written to run significantly faster, but it would become more cryptic. The main point here is to verify the correctness of the code and to clarify the exposition. The red line is the asymptotic value for large , a dimensionless parameter comparing the particle size to the wavelength. The fact that this value is exactly twice the geometric cross section is discussed in [3, 4].

Figure 1. Mie cross section for a sphere of radius and wavevector compared to the geometric cross section and its large limit, in red. For the physical significance of the sharp peaks, see [7].

1.2 Expansion of Mie Coefficients for Small

Here we develop the Taylor expansion of the Mie coefficients for small . This subsection confirms equations (8–11). The results are used in Section 5.

We start with a simplification.

The following definitions are motivated by equations (5) and (6).

Later, we will find we need these for . The index of refraction is .

The Taylor series is taken next. A high order is necessary even for the small limit. In many cases, the spherical Neumann function enters the calculation, which leads to diverergence in this limit. The real and imaginary parts are represented by the suffixes 1 or 2 in the variable names. The series expansions are chosen so as to include terms to the lowest nonvanishing order.

The results are given here as mathematical formulas, which should agree with results in the cell above. Superscripts 1 and 2 refer to real and imaginary parts.

2. Maxwell Stress Tensor

For a harmonic electromagnetic field, the time-averaged electromagnetic force on a dielectric particle in a vacuum is given by [1, 2, 4]

where the angle brackets mean the time average, is the differential of the surface normal, and the integral is taken over any surface enclosing the particle. Thus the Maxwell stress tensor is given by

where

and where is the electric constant, is a dyadic, and is the identity matrix. The companion matrix is given by with . (The magnetic constant used in the references is related by .)

Since the particle is spherical, it is natural to do the integral on the surface of a sphere at radius . Since we are in a vacuum, the electromagnetic field exhibits no force there. Looking ahead, we will check to see that the result is independent of . The direction of is , so it is sufficient to calculate a mixed-basis tensor component , where refers to the component along and is a Cartesian direction. Rewriting equation (12) in spherical coordinates leads to

We accept the complexity of the mixed-bases tensor because the dot product leads to a single component in spherical coordinates, but the integral requires coordinates that are independent of the point of integration. In practice, these are Cartesian. Looking ahead, we will see that only the component will be nonzero, which further motivates the choice. The calculation will proceed by forming the electric field in spherical coordinates and making a row of , where the index implies spherical coordinates. The vector will be transformed to Cartesian coordinates by right-multiplying by a rotation matrix. (The same process is used for the magnetic field .)

The effect of the time average is the following: given , . The stress tensor will be computed from the second term. The first term is included later by adding the complex conjugate.

3. Scattering Force

For the scattering force, we take the incident electric field to be a plane wave traveling in the positive direction and linearly polarized in the direction. Explicitly,

Code: Symbolic Algebra

As discussed, we need to transform a vector from spherical to Cartesian coordinates. We give that function first.

The electric or magnetic stress tensor is computed from the following code. is a factor required for the time average.

gets the component, since the field components have the usual order (, , ); is the complex conjugate, but we need to help the kernel find simplifications. (The magnetic part will be found by symbolic substitution.) The conjugation is performed as follows.

is used for real expansion.

Following are some additional simplifications that we will use. We operate in the dimensionless radial variable called , which is (i.e. times ), where and do not enter the calculation independently. The assumptions are made because of the range of integration.

We introduced extra factors of into equations (1) and (2) because we implement rather than to enable the use of . The operation is named , the curl for a dimensionless radial variable.

Next we define a function to create the electric and magnetic fields (“em”), omitting the prefactor . Before the final answer for the force is obtained, the factor will be reintroduced as the variable .

This gives a simplified interface for a positive plane wave.

It is sufficient to set for both the scattering and gradient forces, since both occur in the limit of small particles. Setting does not capture all contributing terms to lowest order, and setting or higher does not cause the limit to change. For the scattering, we will have , which is for a plane wave going in the positive direction. For the gradient force where we have a standing wave, we decompose that into the sum of a plane wave going toward and one going toward . The backward wave is still polarized in the direction, so we produce it by mapping . Flipping is necessary because the direction of travel is ; if we wish to preserve , it is necessary to change the sign of . We will always set and ; these may be set to zero to suppress the external field or the induced field, respectively.

We next determine the electric and magnetic fields. Although we could increase , the run time increases dramatically. For example, the case of from [1, 2] was run overnight.

For , use 1 for a quick test, 2 to get the main results, and 3 or more to confirm that a higher-order expansion does not change the limit. The run time increases rapidly as is increased. Run times can vary, but typically it is best to start the code and come back to it after a few minutes to a couple of hours. The parameter is used both here and for the gradient force in Section 4.

Next, we make rows from the electric part and magnetic part (multiplied by ) of the stress tensor in Cartesian coordinates.

We do the integrals over next, yielding the three Cartesian components of the electric part of the stress tensor, after azimuthal integration.

The azimuthal integrals of the and components of the electric part of the stress tensor are zero, as they should be. The component is nonzero in general and will be used later.

The magnetic terms are similar.

The integral over is done next. Although these variables contain the word , they lack some constant factors to be forces, hence the prefix . These factors will be included after a few manipulations. The notation ending with is for the component. The successive variables lack the constant, but the final result is the force.

Combining the electric and magnetic terms into a total force at this point leads to some cancellations, so it is convenient to do it next.

Up to this point, we have done the manipulations without telling the kernel that and are, in fact, spherical Bessel functions and . Moreover, we use (a complex function of real arguments) because we will be obtaining real and imaginary parts shortly.

As discussed previously, we have calculated only one of two terms in the force. Next, we add the complex conjugate.

We introduce real and imaginary parts for the Mie coefficients.

The expression is still quite complicated. However, knows (for ) that the Wronskian of the spherical Bessel functions [8]. We will see that the force does not depend on the radius of the integration sphere, as long as . We expect this on physical grounds: there should be no force on the vacuum region surrounding the physical sphere. Mathematically, the Wronskian plays a key role in achieving this condition.

We next include some constant prefactors, namely , and . The factors arise from the definition of the stress tensor. The factor comes from the fact that when we performed the two angular integrals over the sphere, we did not yet include the dimensional constant associated with the area element. This extra factor of is the only exception to the radius appearing in the dimensioness variable . The substitution allows powers of to be canceled.

The real and imaginary parts of the Mie coefficient are and , respectively, and similarly, for they are and .

4. Gradient Force

In the previous section, we used the general form for the electric and magnetic fields for the Mie expansion as input to the Maxwell stress tensor to derive the force on a particle. We do the same in this section, borrowing from the previous section to the extent possible. However, the spherical particle is in a standing wave field instead of a plane wave traveling in the direction. The standing wave is a linear superposition of two plane waves going toward and toward , so we need a solution from the latter field. We obtain this by symmetry from the existing solution. Sending does what is needed: the polarization of is unchanged, but the direction of propagation changes sign. This is implemented with the wrapper to the function described in Section 3.

Although the standing wave is the sum of two plane waves, more precisely, a phased sum is required. Although we do not show it here, the answer is proportional to . This is because the dielectric sphere is located at the coordinate origin, and we need to have a maximum of the electric field there. We do this as follows. (In the variable names, ▽ stands for gradient.)

From here, the manipulation is the same as for the scattering force, so the code is written without further explanation. Certain terms integrate to zero.

First, the electric terms in the stress tensor are calculated. As previously, the and components
are 0.

Second, the magnetic terms in the stress tensor are calculated. Again, the and components are 0 and the component is nonzero in general.

Finally, the electric and magnetic terms are combined and the force is found.

5. Forces in the Limit of Small

Our next task is to determine the scattering and gradient forces on spheres up to the leading order in . Physically, these are spheres that are small compared to a wavelength. We combine the results of Section 3 for the scattering force and Section 4 for the gradient force with those of Section 2 for the Mie coefficients. In Section 1.2, we showed that , , , , , , and using superscripts and for the real and imaginary parts, respectively. Therefore, the scattering force, to lowest order in , retains only the terms in , and for the gradient force, the lowest-order coefficient is . The terms were given in equation (9). The key point is that the Clausius–Mossotti factor appears with the second and first powers in the two terms. The result falls out of a Taylor expansion without any appeal to physical arguments about the response of dipoles. All of the terms are proportional to , leading to the physically required result that if the sphere in fact contains a vacuum , there is no electromagnetic response and thus no force. However, the denominator is characteristic only of the terms giving the lowest-order response. Higher terms have different dependencies on the index of refraction , such as , as seen previously. The gradient and scattering forces properly exist only in the limit of small .

Having selected the lowest-order terms for the two force expansions, these are the forces to lowest order in .

These formulas agree with those given by Harada and Asakura [3], as shown in the Appendix.

Summary

Our goal was to derive the scattering force and the gradient force rigorously from the Mie solution and the Maxwell stress tensor. We began by presenting the Mie solution for a plane wave incident on a dielectric sphere in a vacuum and showing that our implementation matches a figure from a textbook. We then presented the formula for the force in terms of a surface integral of the Maxwell stress tensor, which we take on a sphere of arbitrary radius centered on and including the whole dielectric sphere. We analyzed the scattering force first, giving a formula for the force in terms of the Mie coefficients and then taking the limit as the radius of the sphere tends to zero. This yields agreement with the usual formula for the scattering force in a vacuum. The forces were reformulated for a standing wave, and a similar program was carried out, leading to agreement for the gradient force.

Conclusion

The Claussius–Mossotti term that appears in expressions for the scattering force and the gradient force is seen to be implicit in the rigorous Mie solutions to Maxwell’s equations. By finding the Maxwell stress tensor for a plane wave or a standing wave acting on a dielectric sphere, we are able to show the response in lowest order is in agreement with a widely used formula for the scattering force and the gradient force, respectively. Since part of the derivation includes a tenth-order Taylor expansion of special functions, it is difficult to see how the result could be obtained without computer-assisted algebra.

Appendix: Reconciling the Harada–Asakura Formulas with the Present Results

First, we wish to match equation (17) to the scattering force as given in equation (12) of [3], hereafter called equation (HA12). Variables with superscript (HA) are from the reference. We consider only the case of the external medium being a vacuum, so . This implies , the index in the sphere, so that the Clausius–Mossotti factor is present in both our equation (17) and equation (HA12). The same holds for equation (18) and equation (HA16). Next, [3] considers a Gaussian beam profile, so we simply pick the point in the middle, setting . The factor is the intensity at the center of the beam . The beam intensity is related to the electrical field by . Given these expressions, all the factors can be matched by inspection.

Next, we wish to match equation (18) to the gradient force given in equation (HA16). The one additional equation to note is that . The factor of 2 occurs in the numerator because, by our definition, represents the field of one of the two interfering beams. The 2 in the denominator appears in the intensity-field conversion equation. The twos in are due to the physical fact that a standing plane wave interference pattern is periodic with half the wavelength of the electric fields of the plane waves of which the interference pattern is composed. Again, a match can be made by inspection.

Acknowledgments

Eric L. Shirley provided a key step in the derivation.

References

	J. J. Curry and Z. H. Levine, “Continuous-Feed Optical Sorting of Aerosol Particles,” Optics Express, 24(13), 2016 pp. 14100–14123. doi:10.1364/OE.24.014100.
	Z. H. Levine and J. J. Curry, “OptSortSph: Optical Sorting in a Standing Wave Field Calculated with Effective Velocities and Diffusion Constants,” Journal of Research of the National Institute of Standards and Technology, 121, 2016 pp. 420–421. doi:10.6028/jres.121.020.
	Y. Harada and T. Asakura, “Radiation Forces on a Dielectric Sphere in the Rayleigh Scattering Regime,” Optics Communications, 124(5-6), 1996 pp. 529–541. doi:10.1016/0030-4018(95)00753-9.
	G. Mie, “Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen,” Annalen der physik, 330(3), 1908 pp. 377–445. onlinelibrary.wiley.com/doi/10.1002/andp.19083300302/pdf.
	C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles, New York: Wiley, 1998.
	A. Zangwill, Modern Electrodynamics, Cambridge: Cambridge University Press, 2013.
	P. Chylek, J. T. Kiehl and M. K. W. Ko, “Optical Levitation and Partial-Wave Resonances,” Physical Review A, 18(5), 1978 pp. 2229–2233. doi:10.1103/PhysRevA.18.2229.
	M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, 10th printing, New York: Dover, 1972, Equation (10.1.6).
Z. H. Levine and J. J. Curry, “Scattering and Gradient Forces from the Electromagnetic Stress Tensor Acting on a Dielectric Sphere,” The Mathematica Journal, 2017. dx.doi.org/doi:10.3888/tmj.19-1.

About the Authors

Zachary H. Levine is a physicist at NIST and a Fellow of the American Physical Society. He has concentrated on various aspects of the interaction of light and matter, including the calculation of dielectric constants, second harmonic coefficients and higher optical tensors.

J. J. Curry is an electrical engineer at the National Institute of Standards and Technology. His current interests include manipulation and measurement of particles using optical forces. In recent years, he has worked on problems of interest to the lighting industry.

Zachary H. Levine
The National Institute of Standards and Technology
Gaithersburg, MD 20899-8410
zlevine@nist.gov

J. J. Curry
The National Institute of Standards and Technology
Gaithersburg, MD 20899-8441
john.curry@nist.gov