Continuous-feed optical sorting of aerosol particles

We consider the problem of sorting, by size, spherical particles of order 100 nm radius. The scheme we analyze consists of a heterogeneous stream of spherical particles flowing at an oblique angle across an optical Gaussian mode standing wave. Sorting is achieved by the combined spatial and size dependencies of the optical force. Particles of all sizes enter the flow at a point, but exit at different locations depending on size. Exiting particles may be detected optically or separated for further processing. The scheme has the advantages of accommodating a high throughput, producing a continuous stream of continuously dispersed particles, and exhibiting excellent size resolution. We performed detailed Monte Carlo simulations of particle trajectories through the optical field under the influence of convective air flow. We also developed a method for deriving effective velocities and diffusion constants from the Fokker-Planck equation that can generate equivalent results much more quickly. With an optical wavelength of 1064 nm, polystyrene particles with radii in the neighborhood of 275 nm, for which the optical force vanishes, may be sorted with a resolution below 1 nm. OCIS codes: (350.4855) Optical tweezers or optical manipulation; (350.4990) Particles. References and links 1. A. Ashkin, “Acceleration and trapping of particles by radiation pressure,” Phys. Rev. Lett. 24, 156 (1970). 2. A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, “Observation of a single-beam gradient force optical trap for dielectric particles,” Opt. Lett. 11, 288 (1986). 3. A. Jonás̆ and P. Zemánek, “Light at work: The use of optical forces for particle manipulation, sorting, and analysis,” Electrophoresis 29, 4813 (2008). 4. K. C. Neuman and S. M. Block, “Optical trapping,” Rev. Sci. Instrum. 75, 2787 (2004). 5. A. Ashkin, “History of optical trapping and manipulation of small-neutral particle, atoms, and molecules,” IEEE J. Sel. Top. Quantum Electron. 6, 841 (2000). 6. T. T. Perkins, “Optical traps for single molecule biophysics: a primer,” Laser Photon. Rev. 3, 203 (2009). 7. N. J. Alvarez, C. Jeppesen, K. Yvind, N. A. Mortensen, and O. Hassager, “The chromatic separation of particles using optical electric fields,” Lab Chip 13, 928 (2013). 8. M. Ploschner, T. C̆iz̆már, M. Mazilu, A. Di Falco, and K. Dholakia, “Bidirectional optical sorting of gold nanoparticles,” Nano Lett. 12, 1923 (2012). 9. Y. Pang and R. Gordon, “Optical trapping of 12 nm dielectric spheres using double nanoholes in a gold film,” Nano Lett. 11, 3763 (2011). 10. P. Jákl, T. C̆iz̆már, M. S̆erý, and P. Zemánek, “Static optical sorting in a laser interference field,” Appl. Phys. Lett. 92, 161110 (2008). 11. K. Grujic, O. G. Hellesø, J. P. Hole, and J. S. Wilkinson, “Sorting of polystyrene microspheres using a Y-branched optical waveguide,” Opt. Express 13, 1 (2005). 12. O. Brzobohatý, V. Karásek, M. S̆iler, L. Chvátal, T. C̆iz̆már, and P. Zemánek, “Experimental demonstration of optical transport, sorting, and self-arrangement using a ‘tractor beam’,” Nat. Phot. 7, 123 (2013). 13. A. Neild, T. W. Ng, and W. M. S. Yii, “Optical sorting of dielectric Rayleigh spherical particles with scattering and standing waves,” Opt. Express 17, 5321 (2009). #263305 Received 14 Apr 2016; revised 7 Jun 2016; accepted 9 Jun 2016; published 15 Jun 2016 (C) 2016 OSA 27 Jun 2016 | Vol. 24, No. 13 | DOI:10.1364/OE.24.014100 | OPTICS EXPRESS 14100 14. M. Guillon, O. Moine, and B. Stout, “Longitudinal optical binding of high optical contrast microdroplets in air,” Phys. Rev. Lett. 96, 143902 (2006). 15. T. C̆iz̆már, V. Garcés-chávez, K. Dholakia, and P. Zemánek, “Optical conveyor belt for delivery of submicron objects,” Appl. Phys. Lett. 86, 174101 (2005). 16. K. Ladavac, K. Kasza, and D. G. Grier, “Sorting mesoscopic objects with periodic potential landscapes: optical fractionation,” Phys. Rev. E 70, 010901(R) (2004). 17. S. Tatarkova, W. Sibbett, and K. Dholakia, “Brownian particle in an optical potential of the washboard type,” Phys. Rev. Lett. 91, 038101 (2003). 18. M. Horstmann, K. Probst, and C. Fallnich, “An integrated fiber-based optical trap for single airborne particles,” Appl. Phys. B 103, 35 (2011). 19. D. Rudd, C. López-Mariscal, M. Summers, A. Shahvisi, J. C. Gutiérrez-Vega, and D. McGloin, “Fiber based optical trapping of aerosols,” Opt. Express 16, 14550 (2008). 20. K. Xiao and D. G. Grier, “Sorting colloidal particles into multiple channels with optical forces: prismatic optical fractionation,” Phys. Rev. E 82, 051407 (2010). 21. D. R. Burnham and D. McGloin, “Holographic optical trapping of aerosol droplets,” Opt. Express 14, 4176 (2006). 22. O. Brzobohatý, M. S̆iler, J. Jez̆ek, P. Jákl, and P. Zemánek, “Optical manipulation of aerosol droplets using a holographic dual and single beam trap,” Optics Lett. 38, 4601 (2013). 23. A. E. Carruthers, J. P. Reid, and A. J. Orr-Ewing, “Longitudinal optical trapping and sizing of aerosol droplets,” Optics Express 18, 14238 (2010). 24. K. J. Knox, D. R. Burnham, L. I. McCann, S. L. Murphy, D. McGloin, and J. P. Reid, “Observation of bistability of trapping position in aerosol optical tweezers,” J. Opt. Soc. Am. B 27, 582 (2010). 25. D. McGloin, D. R. Burnham, M. D. Summers, D. Rudd, N. Dewar, and S. Anand, “Optical manipulation of airborne particles: techniques and applications,” Faraday Discuss. 137, 335 (2008). 26. K. Taji, M. Tachikawa, and K. Nagashima, “Laser trapping of ice crystals,” Appl. Phys. Lett. 88, 141111 (2006). 27. P. Jákl, A. V. Arzola, M. S̆iler, L. Chvátal, K. Volke-Sepúlveda, and P. Zemánek, “Optical sorting of nonspherical and living microobjects in moving interference structures,” Opt. Express 22, 29746 (2014). 28. M. Ploschner, M. Mazilu, T. C̆iz̆már, and K. Dholakia, “Numerical investigation of passive optical sorting of plasmon nanoparticles,” Opt. Express 19, 13922 (2011). 29. K. Xiao and D. G. Grier, “Multidimensional Optical Fractionation of Colloidal Particles with Holographic Verification,” Phys. Rev. Lett. 104, 028302 (2010). 30. R. F. Marchington, M. Mazilu, S. Kuriakose, V. Garcés-Chávez, P. J. Reece, T. F. Krause, M. Gu, and K. Dholakia, “Optical deflection and sorting of microparticles in a near-field optical geometry,” Opt. Express 16, 3712 (2008). 31. R. L. Smith, G. C. Spalding, K. Dholakia, and M. P. MacDonald, “Colloidal sorting in dynamic optical lattices,” J. Opt. A: Pure and Appl. Opt. 9, S134 (2007). 32. Z. H. Levine and J. J. Curry, “Scattering and Gradient Forces from the Electromagnetic Stress Tensor Acting on a Dielectric Sphere,” Mathematica J. (submitted) (2016). 33. P. Debye, “Der Lichtdruck auf Kugeln von belibigem Material, Ann. der Phyzik IV 335, 57 (1909). 34. J. P. Barton, D. R. Alexander, and S. A Schaub, “Theoretical determination of net radiation force and torque for a spherical particle illuminated by a focused laser beam,” J. Appl. Phys. 66, 4594 (1989). 35. A. Zangwill, Modern Electrodynamics, (Cambridge University, 2013) pp. 787–788. 36. R. N. C. Pfeifer, T. A. Nieminen, N. R. Heckenberg, and H. Rubinsztein-Dunlop, “Colloquium: Momentum of an electromagnetic wave in dielectric media,” Rev. Mod. Phys. 79, 1197 (2007). 37. Y. Harada and T. Asakura, “Radiation forces on a dielectric sphere in the Rayleigh scattering regime,” Opt. Commun. 124, 529 (1996). 38. Nicholas A. Krall and Alvin W. Trivelpiece, Principles of Plasma Physics, (San Francisco Press, Inc., 1986) p. 295. 39. S. Chandrasekhar, “Stochastic problems in physics and astronomy,” Rev. Mod. Phys. 15, 1 (1943). 40. N. A. Fuchs, The Mechanics of Aerosols, (Dover, 1989) p.27. 41. G. J. Borse, Numerical Methods with MATLAB, (PWS Publishing, 1997). 42. Z. H. Levine and J. J. Curry, (to be submitted for publication). 43. P. Morse, Thermal Physics 2nd ed. (Benjamin, 1969) pp. 228–236. 44. D. S. Dean and G. Oshanin, “Approach to asymptotically diffusive behavior for Brownian particles in periodic potentials: Extracting information from transients,” Phys. Rev. E 90, 022112 (2014). 45. P. P. Mathai, J. A. Liddle, and S. M. Stavis, “Optical tracking of nanoscale particles in microscale environments,” Appl. Phys. Rev. 3, 011105 (2016). 46. I. D. Nikolov and C. D. Ivanov, “Optical plastic refractive measurements in the visible and the near-infrared regions,” Appl. Opt.39, 2067 (2000). 47. P. Zemánek, A. Jonás̆, P. Jákl, J. Jez̆ek, M. S̆erý, and M. Lis̆ka, “Theoretical comparison of optical traps created by standing wave and single beam,” Opt. Commun. 220, 401 (2003). 48. P. T. Korda, M. B. Taylor, and D. G. Grier, “Kinetically locked-in colloidal transport in an array of optical tweezers,” Phys. Rev. Lett. 89, 128301 (2002). #263305 Received 14 Apr 2016; revised 7 Jun 2016; accepted 9 Jun 2016; published 15 Jun 2016 (C) 2016 OSA 27 Jun 2016 | Vol. 24, No. 13 | DOI:10.1364/OE.24.014100 | OPTICS EXPRESS 14101 49. S. Ahlawat, R. Dasgupta, R. S. Verma, V. N. Kumar, and P. K. Gupta, “Optical sorting in holographic trap arrays by tuning the inter-trap separation,” J. Opt. 14, 125501 (2012). 50. K. H. Lee, S. B. Kim, K. S. Lee, and H. J. Sung, “Enhancement by optical force of separation in pinched flow fractionation,” Lab Chip 11 354 (2011). 51. M. P. MacDonald, G. C. Spalding, and K. Dholakia, “Microfluidic sorting in an optical lattice,” Nature 426, 421 (2003). 52. H. A. Kramers, “Brownian motion in a field of force and the diffusion model of chemical reactions,” Physica 7, 284 (1940). 53. B. Derrida, “Velocity and Diffusion Constant of a Periodic One-Dimensional Hoping Model,” J. Stat. Phys. 31, 433 (1983).


Introduction
Following demonstrations by Ashkin and co-workers [1,2], optical forces have been widely used to manipulate small particles with dimensions in the range of nanometers to micrometers [3,4,5].One of the more popular areas of application is biology where optical tweezers have been used to study a variety of particles in liquid, including cells, bacteria, viruses, lipid vesicles, and even DNA molecules [6].
Basic optical tweezers use a high numerical aperture (NA) microscope objective to focus a laser beam to a diffraction-limited spot.The strong spatial gradients in optical intensity produce trapping in three dimensions.Strong intensity gradients can also be produced by near-field structures [7,8,9,10,11], interference fringes [12,13,14,15,16,17], and optical fibers [18,19].Many of these techniques also utilize high NA optics.The trapping volume can be exceedingly small (< 1 μm 3 ) and located very close to a microscope objective.Micro-fluidics technology is used to deliver liquid-borne particles to the trapping area.Holographic arrays of tweezers have been used to obtain larger working volumes and to operate on more particles [20,21].
Optical forces have great potential for solving important measurement challenges with aerosol particles.It is possible to manipulate and measure aerosol particles without altering or perturbing the particles as is required in diffusion mobility analysis or in collection by mechanical filters.Measurement of aerosol particles arises in such disparate fields as climate science, nano-particle manufacturing, homeland security, healthcare, and forensics.Optical manipulation of aerosol particles has received some attention [13,14,21,22,23,24,25,26], but much less than manipulation of liquid-borne particles.
A desirable functional capability for any particle measurement is sorting, with size being a natural parameter to sort.Sorting makes it possible to measure particle size distributions, to manufacture mono-disperse particles, and to detect the presence of a specific type of particle.Some work has been done on sorting of particles in liquids [8,12,20,27,28,29,30,31], but almost none on sorting of aerosol particles [13].
High throughput is desirable for any sorting function.High throughput can be facilitated by parallel processing, which in turns requires a considerably larger working volume than the trapping volume of optical tweezers.Arrays of optical tweezers have been shown to provide parallelization [20,29].An alternative is to utilize low-NA optics to achieve a large sorting volume, through which a large number of particles can flow.
Here we propose and analyze a scheme for optically sorting airborne nanoparticles.It has the advantages of (1) utilizing a simple, low NA optical system with high particle throughput, (2) producing a continuous stream of continuously dispersed particles, and (3) exhibiting a size resolution better than 1 nm.
We have explored the performance of the sorting scheme using numerical simulations of optical and collisional forces.The reliability of the numerical results is enhanced by our use of two separate algorithms, which have been extensively compared against each other and against analytic results.One of the algorithms is a Monte Carlo solution to the equations of motion for individual particle trajectories through 6-dimensional phase space.The other approach uses the Fokker-Planck equation to determine effective velocities and diffusion constants obtained by averaging over the interference fringes of the optical field; these functions of optical intensity are then used to integrate the equations of motion over the large-scale structure of the field.
A key input to the simulations is the force of a standing wave (wavelength λ 0 ) on a particle of arbitrary radius a.We derive this force from first principles [32].It reproduces existing expressions for the force in the Rayleigh regime a λ 0 .
In this article we restrict our attention to particles in air, but the approach is qualitatively applicable and the results obtained are qualitatively similar to particles in liquids.In water, for example, the drag force is nearly two orders of magnitude larger, while Brownian motion is nearly one order of magnitude smaller.The former can be accommodated by adjusting controllable parameters.The latter is generally beneficial.

Sorting scheme
We consider a system designed to sort spherical particles by size.The basic configuration, given in Fig. 1, consists of a heterogeneous stream of aerosol particles, entrained in a global laminar air flow, moving across a Gaussian-mode standing wave at an oblique angle.A Cartesian coordinate system x, y, z (See Fig. 1) is chosen in which the z-axis coincides with the optical axis and the air flows parallel to the x-z plane.Particles of all sizes enter the optical field in a narrow input stream, but are physically separated by their interaction with the field.They exit the field at different locations depending on their size.
In the absence of an optical field, the particle stream flows in a straight line across the optical axis.A standing wave optical field presents a sinusoidal force directed along the optical axis.The large gradients in a standing wave make this force comparable in magnitude to that seen in optical tweezers.It is sufficient to impede long-range motion of particles in z, but does not effect motion in x and y.Sorting of particles arises from the size-dependence of the optical force combined with the dependence on distance from the optical axis.In the Rayleigh regime, for example, some particles are small enough that they are practically unaffected by the optical field and flow nearly undiverted in a straight line across the optical axis.Larger particles will be diverted into a purely x-directed motion when they get close to the optical axis where the optical field is largest.Particles of even larger size are diverted at lower field strengths that occur farther from the optical axis.Thus particles of different sizes are physically separated from each other.
The separation increases, by the same amount, as the particles exit the optical field and begin to flow, again, along streamlines.The net result is a continuous dispersion of particle sizes along z.Particle size distributions may then be measured in situ, or the different sizes may be separated permanently for retention or subsequent analysis, depending on the application.
This problem involves an optical force, a drag force from the air, and diffusion of particles due to random collisions with air molecules.The latter makes the problem three-dimensional, but for the conditions we explore, diffusion in the y-dimension is not so important.

Optical force on a dielectric sphere in a standing wave
Expressions for the force on a spherical particle in terms of the Mie coefficients have been known since the classic work of Debye [33]; a general formula has been given by Barton et al. [34].We consider the force exerted on a spherical dielectric particle of arbitrary radius by a standing wave.
The force exerted by an arbitrary electromagnetic field on an arbitrary object can be obtained by generalizing the point charge in the Lorentz force equation to a charge density and a current density and then integrating over the object volume.Gauss's law and Ampere's law make it possible to describe the charge and current densities in terms of material properties, i.e., dielectric permittivity ε and magnetic permeability μ.The Maxwell-Faraday equation and Gauss's law for magnetism can then be used to obtain the electromagnetic force on an object [35] in terms of the divergence of the Maxwell stress tensor and the time rate-of-change of the Poynting vector S. The integral is over the volume of the object.The brackets indicate a time-average over optical frequencies, Ī is the identity matrix and E and B are the electric and magnetic field vectors, respectively.These fields must be obtained from the incident fields by applying the appropriate boundary conditions at the object surface.EE and BB are dyadics.We will only consider situations in which the electromagnetic field has a single optical frequency time-dependence, so the last term in Eq. ( 1) does not contribute to the time-averaged force.In a material medium, the electromagnetic stress tensor must be accompanied by corresponding terms describing the momentum carried by the medium [36].In the case of a particle in vacuum, such terms vanish.We will proceed from this point as if the electromagnetic properties of air are identical to those in vacuum.
Considering a spherical particle, it is easier to evaluate the force when the volume integral is converted to a surface integral using the Divergence Theorem where the surface of integration must be outside the particle, i.e., the fields must be those external to the particle [32].Both incident and scattered fields are included in E and B.
We are concerned with the force generated by a low numerical aperture standing wave, so we consider an infinite plane standing wave with electric field where ω/(2π) is the optical frequency, k 0 = 2π/λ 0 is the wave number of the optical field, λ 0 is the free-space wavelength, and z is the coordinate along which the waves propagate.This field creates a time-averaged intensity where is the power per unit area in each plane wave.
We have obtained the force on a spherical particle of arbitrary radius in a plane standing wave by numerically evaluating Eq. ( 3).The validity of the results is ensured by our use of two different approaches to the numerical evaluation, which were implemented in isolation.In one approach, the incident and scattered fields were derived analytically (in terms of Bessel and Hankel functions) and then evaluated numerically at the surface of the particle.In the second approach, the fields inside the particle were derived analytically (in terms of Bessel functions only).The fields on the external surface of the particle were then obtained from the continuity of the tangential components of the electric and magnetic fields and the continuity of the normal components of the electric displacement and magnetic induction.The force obtained from numerical evaluation of these fields is in agreement everywhere with the force obtained by the first approach.
The five-term expansion of Eq. ( 3) about k 0 a = 0 is an excellent analytic approximation to the time-averaged force of a standing wave on a spherical particle in vacuum for the range of parameters we consider.The Mie coefficients a and b are given by Zangwill [35].Values for the coefficients c i are given in Table 1.The first term in Eq. ( 7) represents a product of the incident electric field and the scattered field, whereas the following two terms represent the product of the scattered field with itself.For k 0 a → 0, Eq. ( 7) reduces to  2. Comparison of different expressions for the optical force versus a/λ 0 for λ 0 = 1064 nm, n p = 1.57, and k 0 z = π/4.The force is normalized by the factor I 0 πa 2 , so that the plotted functions give the force per unit of optical power incident on the sphere.The dashed black line is the Rayleigh expression Eq. ( 8).The blue curve is the 5-term expression of Eq. (7).The red curve is the actual force obtained from high-order numerical evaluation of Eq. (3).which agrees with the expression for the gradient force in the Rayleigh limit [37] when the optical field is a plane standing wave.The Clausius-Mossotti term (n 2 p − 1)/(n 2 p + 2), frequently introduced in the literature by an ad hoc argument involving scooped spheres and dipoles, is seen here to originate naturally as the lowest order consequence of the electromagnetic boundary conditions at the surface of a dielectric sphere, in particular, through the Mie coefficient b 1 [32].
Figure 2 compares three different expressions for the force on a spherical dielectric particle in a standing wave: the Rayleigh expression of Eq. ( 8), our analytic approximation Eq. ( 7), and the essentially exact force obtained from evaluating Eq. (3) to order = 30.The force is normalized by the factor I 0 πa 2 , so that the plotted functions give the force per unit of optical power incident on the sphere.The force is everywhere proportional to sin 2k 0 z, with z = λ 0 /8 in the figure .The Rayleigh expression diverges qualitatively from the exact solution for a/λ 0 > 0.1, while our approximate expression is in excellent agreement for the range of a/λ 0 shown in Fig. 2, including the irregular behavior.The force vanishes for particles of certain sizes [47].This is to be expected when the particle size has a certain relationship to an integral number of fringes of the optical field.It is perhaps surprising that the force changes sign even though the particle center remains fixed relative to the standing wave.When the force is negative, a particle with index of refraction greater than that of the medium prefers to sit in the intensity maximum.When the force is positive, the particle prefers to have its center in the intensity minimum.For our scheme, this small shift of position has no practical significance.
The vanishing of the optical force for certain discrete radii is potentially useful for tying particle sizes to an SI length standard because the radii at which the zero-crossings occur are proportional to the optical wavelength.The latter can be measured accurately to high precision.
There is also a dependence on the particle index of refraction, but this parameter can also be measured independently.
The preceding results have been derived for a plane standing wave, but we can use them directly in simulating the force from a Gaussian TEM 00 mode standing wave because the forces arising from the finite radial extent of the beam are small enough to be ignored relative to the axial force from the interference fringes.Considering a Rayleigh length much greater than the spatial extent of the sorter, the TEM 00 mode looks like a collimated beam and the force can be approximated by where is the force per unit irradiance exerted by a plane standing wave when the particle's position is such that 2k 0 z = π/2, I 0 is the intensity of each plane wave at the waist, and w 0 is the e −2 intensity radius.The Rayleigh length corresponding to a given beam waist radius is

Collisional forces
In this section, we describe the general theory of the transport of spherical particles in a fluid, then discuss two models for calculating particle trajectories under the combined influences of optical forces and collisional forces.One model solves the three-dimensional equations of motion numerically, using a Monte Carlo algorithm to account for the random nature of collisional forces.Here, collisions are incorporated as random changes in particle velocity, so the six-dimensional phase-space probability distribution is utilized.A second model greatly simplifies the picture of collisions, by assuming that the velocity-space distribution equilibrates with the background gas in a time much shorter than the time required for a particle to move a significant distance.This and the quasi-1D symmetry of the optical field allows the particle dynamics along the optical axis to be described with a one-dimensional spatial probability distribution function.Effective velocities and diffusion coefficients describe the modification of trajectories produced by the one-dimensional fringes of a given magnitude.The particle trajectories perpendicular to the optical axis are described by integrating the equations of motion across the radially-inhomogeneous optical field using the effective fluid parameters appropriate to each location.The use of two different numerical approaches adds confidence to the correctness of results.The second method is quite a bit faster computationally.We begin this section with background common to both models, Section 5 discusses the Monte Carlo specific issues, and Section 6 discusses the use of effective velocities and diffusion constants.
The trajectory of a particle under the influence of collisions with air molecules is stochastic and can be described by the time evolution of the single-particle phase-space probability distribution function f such that gives the probability that a particle will, at time t, be found within the volume d 3 r about r and within the velocity range d 3 v about v after having an initial position r i and initial velocity v i at time t i , hereafter taken to be 0. The probability distribution function is normalized such that integration over all final positions and velocities yields a probability of one.The time evolution of f is described by the Boltzmann equation in which ∇ v is the gradient in velocity-space, F represents the optical forces on the particle, m p is the particle mass, and d f dt c represents the stochastic forces of collisions between the particle and surrounding molecules.We will not consider situations in which the particle density is large enough that interactions between particles are significant.
The momentum exchange between particle and air molecule in a given collision is such a small fraction of the particle momentum that the particle trajectory through phase space can be described as a continuous process using the Fokker-Planck formalism [38].In the present case, the collision cross section is independent of velocity and both the collisions and medium are isotropic.This allows us to write the collision operator as where the vector C is the coefficient of convection and the scalar D v is the coefficient of velocity-space diffusion.Realizing that the net flux in velocity-space must vanish when the particle probability distribution f has equilibrated with the air molecules, and equating the convection term with the drag force, it is possible to express the Fokker-Planck coefficients in terms of the particle mobility μ as where v 0 is the velocity of the medium, k B is the Boltzmann constant, and T is the absolute temperature of the medium.Chandrasekhar [39] has given the fundamental solution to Eq. ( 13), which generalized to a moving medium is where is the time constant for equilibration of the particle kinetic temperature with the air temperature.
For a particle with a = 100 nm in air at standard temperature and pressure, τ ≈ 50 ns.The Einstein-Smoluchowski equation relates the diffusion constant D to mobility by and thus to

Stokes Fuchs
Fig. 3.The mobility for a spherical particle in air at standard temperature and pressure according to Fuchs [40] compared to the g → 0 limit.
we treat collisions of the particle with air molecules as velocity-space diffusion.A particle's random velocity walk then leads deterministically to a diffusion in real space.In our model based on effective velocities and diffusion constants, collisions are treated directly as a diffusion in real space.Both approaches lead to the same result.Fuchs [40] gives for the mobility of a particle in a gas in terms of the gas viscosity η and molecular mean free path g .The constants A = 1.25, Q = 0.42, and b = 0.87 were obtained empirically and g = 68 nm at standard temperature and pressure.In the large particle limit (a/ g 1), the Fuchs expression results in Stokes' Law for the drag force on a particle For small particles (a/ g < 10), the numerator in the Fuchs expression accounts for the transition from fluid to molecular dynamics in which the mobility μ varies as a −2 .The Fuchs mobility is compared, in Fig. 3, to mobility in the limit g → 0.

Monte Carlo simulation of particle trajectories
Multiplying the Boltzmann equation by v and integrating over all phase-space and multiplying the Boltzmann equation by r and integrating over all phase-space yields the equations of motion for a particle where and the force in the term F(r,t) m p is due entirely to optical forces.The trajectory of a particle is determined numerically from this system of equations using a fourth-order Runge-Kutta discretization method [41], given an initial velocity and position.The collision term is dealt with by randomly selecting a change in velocity Δv ≡ v final − v initial from the time-dependent probability distribution function Eq. ( 15) given a discrete time interval Δt.Then The interval Δt = 10 ns was found to achieve satisfactory results.It is much longer than the mean collision time of order 10 fs, yet small enough that a representative particle does not have time to change its position significantly in the external (optical) force field.A particle with radius a = 40 nm and moving at thermal velocity moves only about 1 nm in 10 ns.Larger particles move even shorter distances.
Although numerically intensive, the Monte Carlo approach is fully three-dimensional and can easily accommodate an arbitrary optical field and other forces.

Particle trajectories using effective velocities and diffusion constants
Here, we describe an alternate approach to calculating particle trajectories.We work on two length scales: (a) the submicrometer length scale of an individual optical fringe, and (b) the millimeter length scale of the system as a whole.We solve the Fokker-Planck equation for the spatial probability distribution (in contrast to the velocity-space probability distribution in the previous section) of a particle as a function of time in a 1D optical interference pattern.By considering the time-evolution of the probability distribution of a particle, it is possible to find an effective velocity and effective spatial diffusion constant in a regime in which the position and second moment both change linearly in time.After crossing a fringe, the particle finds itself facing a symmetry-equivalent optical landscape, so the macroscopic motion of the particle can be defined by the microscopic constants.Integration of these constants through the experimental setup leads to predictions which are comparable to the Monte Carlo results.Software for this section is available [42].

Fokker-Planck equation in a plane wave basis
The mathematical content of the "slowly varying" assumption is that we model the mean motion and its variance in one dimension and that we neglect the small forces in the x direction.Because the particle is large compared to the air molecules, it undergoes overdamped Brownian motion and satisfies the Fokker-Planck equation [43] ∂ The diffusion constant D of the particle is tied to the mobility μ by Eq. ( 17), F(z) is the static force on the particle, and f (z,t) is the probability distribution for the location of a particle at (It is the 1D analogue of the distribution used in the previous section, integrated over all velocity space.)Unlike the previous section, here we assume that the particle achieves thermal equilibrium before it moves a significant distance.Physically, f represents a localized function.However, we describe it as an artificially periodic function for mathematical convenience, with the artificial period chosen to be larger than the extent of f .
Since the external force is smooth and periodic, we anticipate a rapidly converging expansion for the probability distribution f (z,t) in a Fourier basis: where Λ = λ 0 /2, the c n (t) are complex expansion coefficients, and k = 2π/(NΛ), where N is the number of interference fringes between the artificial periodic replicas of the solution. Because which implies c 0 = 1.Putting these facts together, We consider the case of a standing wave optical field and a convective background with the form where F 1 is the optical force from Eq. ( 9).The influence of the background fluid with velocity v 0 is implemented by setting the spatially invariant force term to F 0 = v 0 /μ using Eq. ( 18).Substituting Eq. ( 25) and Eq. ( 28) into the partial differential Eq. ( 24) leads to the system of ordinary differential equations The relation c −n = c * n means that it is sufficient to solve only for n > 0. Terms with n < 0 may be eliminated from the equation in favor of terms with positive n.In the special case v 0 = 0, the solution is symmetric and the c n are real.Moreover, the sum is truncated after a finite number of terms N term , i.e., we assume c n = 0 for n > N term .These relations suffice to limit the number of c n which need to be found to N term .Because F(z) has period Λ, by Bloch's theorem the system Eq.( 29) splits into N independent sets of coupled equations.Furthermore, each independent system of equations is tri-diagonal.The net result is that the computational burden is very low, with a few seconds required to solve the equations in a particular case.A typical solution to the Fokker-Planck equation is shown in Fig. 4.An initial Gaussian quickly relaxes within its well, shown by the shift between the blue and green curves, and then more slowly populates the first well downwind as the central well is depopulated.A very small population of the first well upwind is also visible.

Analytic limits and program tests
The Fokker-Planck equation has certain recognizable limiting cases.If D = 0 and F 1 = 0, Eq. ( 24) reduces to the advection equation (a one-way wave equation).If F(z) = 0, we obtain the diffusion equation at zero velocity and the drift-diffusion equation at finite velocity.The case of neglecting diffusion is given in Appendix A and that of neglecting drift is given in Appendix B. Setting any one parameter among D, v 0 , and F 1 to zero leads to an analytic result for Eq. ( 29).We verified all three cases numerically for Gaussian starting distributions.

Extraction of effective velocities and diffusion constants
The goal of this section is to obtain two functions v eff and D eff characterizing the linear increase in the mean and variance rate of f (z,t) as a function of particle diameter and field intensity.The system of differential equations is started with a Gaussian distribution of width σ = Λ/8 centered in one of the wells.The times need to be large enough so that there is intra-well equilibration.Numerically converged results are achieved with the parameters in Table 2.
The mean and variance are calculated as a function of time and locally fit to a linear function.The parameters of Table 2 were chosen so that the results yield nearly constant values for v eff and D eff in the second half (in time) of the simulation after the intra-well equilibration.Typical results are shown in Fig. 5 for the solution corresponding to Fig. 4. Dean and Oshanin have recently observed similar behavior and analyzed the transient term [44].
Effective velocities for several particle radii from 40 nm to 140 nm are shown in Fig. 6.They have been normalized to the velocity at zero optical force.Small particles are relatively unaffected by the optical field, but the effect grows with particle radius.The differences in these Fig. 5. Time evolution for the velocity in z, v obs , which is the average over 1/16 of the total time range, of a spatial distribution of particles of radius a = 100 nm in a standing wave with I 0 = 0.4 GW/m 2 and a flow velocity v 0 = 0.707 mm/s.A second average is taken in the shaded interval to determine v eff , the large-time limit of v obs .
effective velocities permit optical sorting, as we see below.

Use of effective constants of motion
In principle, the effective velocities and effective diffusion constants could be used as input to Monte Carlo modeling where the steps need to be small compared to changes in the fringeaveraged optical intensity which is typically 1 mm.By comparison, a direct application of Monte Carlo to the drift-diffusion equation requires steps which are small compared to Λ, equal to 532 nm in our example.However, the simplicity of the geometry allows a more direct approach.We take the xcomponent of velocity to be a constant v x = v 0 sin θ with θ defined by Fig. 1 We will also use v z = v 0 cos θ .Because the system is translationally invariant in z, the z motion does not affect the accumulation of either displacement or variance in z position.In the y direction, we need to make an additional physical assumption that the total diffusion in y is small compared to the scale length of the optical field.We test this assumption at the end.
The mean displacement in z compared to displacement without an optical field is given by where a particle starts at position (x 0 , 0, z 0 ) and ends at (x 1 , 0, z 1 ), making use of the uniform motion in x.We use the subscript "eff" to denote an average over multiple fringes of the optical field.The variance in z 1 (σ 2 z ) has two contributions.First, we have the effect of 1 dimensional diffusion in z during the mean crossing time.This is given by σ The important contributions to Eqs. ( 30) and ( 31) come from different regions.For the velocity integration, the differences in the central region drive the whole phenomenon.For the diffusion integration, since most of the path is through a relatively less intense field, and since the results add in quadrature, the reduction (or even a small increase) in the effective diffusion constant at high intensity has only a few percent effect on the result compared to simply using the zero-field value for the duration of the transit.
A second contribution to σ 2 z arises from the variations in crossing times.After the time t 1 = (x 1 − x 0 )/v x , the distribution of particles acquires a variance in x given by σ 2 x = 2Dt 1 .Neglecting a term of order σ x /(x 1 − x 0 ) throughout this argument, the standard deviation in the crossing times is given by σ t = σ x /v x .The drift in z acquires a corresponding standard deviation where the average is the integral in Eq. ( 30).These two processes -1D diffusion and variable time t 0 drift -are independent, so the total variance in x is given by The variance of position in the y direction is given by the diffusion constant times the crossing time; symbolically, If we pick the highest value of D used in this study, namely 2.2 × 10 −10 m 2 /s, picking the parameters from Table 2 yields σ y = 43 μm which is small compared to the beam waist.This justifies the assumption that the diffusion in y may be neglected.Methods of tracking nanoparticles optically has been reviewed recently [45].It is possible in principle to observe v eff and D eff directly by optical means.

Results
We show optical sorting of aerosol particles by simulating conditions readily implemented in an actual device.Table 3 gives numerical values for simulated parameters.High-power frequency-stabilized lasers with wavelength λ 0 = 1064 nm are commonly available.A Fabry-Perot cavity of modest finesse (≈ 1000) and small numerical aperture (≈ 10 −3 ) can be used to establish a standing wave with intensities great enough to achieve the same magnitude of forces commonly used in optical tweezers.The larger beam waist of the present approach requires more power than optical tweezers, but it renders diffusive transport in y irrelevant and makes a high throughput possible.We assume spherical particles of polystyrene, with a real index of refraction n p = 1.57[46].There is no optical absorption by the particles and we do not consider interactions between the particles.The small spatial scale and low flow rate considered ensure laminar flow with Reynolds number Re ≈ 0.2.
The particles are entrained in the 1 mm/s laminar flow of air at standard temperature and pressure from a point source.Although an experimental implementation would utilize a finite diameter particle source instead of a point source, predictions for the former can be estimated by convolving our results over any size input aperture.In doing so, it should be remembered that the flux density from an aperture will be concentrated in the center of the aperture with a parabolic profile because the flow velocity vanishes at the edge of the aperture.The release point is x 1 = −3w 0 .The air velocity vector is in the x-z plane and oriented at 45 • to the optical axis.
Figure 7 shows mean trajectories, from the Monte Carlo model, for several populations of particles characterized by radii from 40 nm to 110 nm travelling across an optical field with I 0 = 2 GW/m 2 .These simulated trajectories reflect what is shown schematically in Fig. 1, with the dispersion of different particle sizes along the z coordinate in the output plane.The smallest particles travel in nearly a straight line across the optical field, intersecting the output plane near z 1 = 0. Larger particles are deflected by the optical field and intersect the output plane at negative values of z 1 .All particle trajectories have been projected onto the x-z plane in this figure, but the simulation was fully three-dimensional.
Figure 8 shows, for the same data, spatial distributions of the particle populations in the output plane.The solid lines are Gaussian fits and provide a reasonably accurate, though not detailed, description of the distributions.This simulation included 1600 individual particles of each size.It follows from Fig. 8 that a continuous uniform distribution of particle sizes going into the sorter will produce a continuous distribution in the output plane.Any given point, or interval, in the output plane corresponds to a distribution of particle sizes whose mean value is a function of z.The width of each size distribution is also a function of z and is about 10 nm (FWHM≈ 2σ a ) for radii near a = 70 nm.
Table 4 shows a comparison of our two methods for calculating particle trajectories.Here, we compare the positions and one standard deviation widths of the distributions shown in Fig. 8 (Monte Carlo) with the corresponding values obtained from the effective velocity method.
The widths of the spatial distributions are determined by diffusion during the time the particles transit the optical field from the point source to the output plane.Although the optical field suppresses diffusion for some part of that transit time, depending on particle size, this  suppression is not very important in determining distribution widths.Most of the diffusive expansion of the particle stream occurs immediately after release from the point source because diffusive expansion is proportional to t 1/2 .Figure 9 compares the distribution widths in the output plane for I 0 = 0.4 GW/m 2 with those obtained when no optical field is present (I 0 = 0).The difference in distribution width at the output is due almost entirely to the difference in particle size.
We use the effective constants of motion approach to explore particle sorting over a larger range of particle sizes and a range of optical intensities.First, we look at the mean deflection z 1 of particles in the output plane as the result of interaction with the optical field.Figure 10 gives |z 1 | as a function of particle radius for a number of values of I 0 .(All values of z 1 are negative or zero.)In general, the deflection increases like the absolute magnitude of the force (Fig. 2), though the increase is not linear with either radius or intensity.The radius a = 275 nm corresponds to the first zero-crossing of the force (Fig. 2).The deflection also vanishes here, but unlike the force, it does not reverse sign for larger a.The function of the standing wave in this sorting configuration is simply to provide a barrier to particle movement along z and a reversal in the sign of the force simply corresponds to a change of π in the phase of the fringes, but the nature of the barrier remains essentially the same.Thus, as a increases toward the zerocrossing, deflection decreases.It vanishes at the zero-crossing, but then increases again as the magnitude of the force increases with further increases in a.
The vanishing of the optical force imposed by a standing wave on particles of certain sizes is known [47], but its potential exploitation for determination of absolute particle sizes has not been realized.Particles with radii corresponding to the zero-crossings exit the optical field at the edge of the sorted particle stream with both larger and smaller particles deflected to the left in Fig. 7.This enables an easy identification of the desired particles because one need only find the edge of the sorted stream, instead of an absolute position.The radius corresponding to the zero-crossing is a function only of the wavelength of the optical field and the index of refraction of the particle (and the medium when it differs from unity).The wavelength can be adjusted to nearly any desired value and measured absolutely to any practically useful precision by common methods.The index of refraction can also be determined by other well-known means, although not with the same level of precision.Nevertheless, the precision is likely good enough to obtain particle size values to better than 0.1 % absolute accuracy.
The dispersion of particle sizes along the output plane is the derivative dz 1 /da of the deflection.In Fig. 10, one can see ranges of particle radius that are potentially most attractive for sorting particles by size.The Rayleigh regime is one such region, using relatively large intensities.Even better are the regions on either side of the zero-crossing (a = 275 nm) where the deflection is the most sensitive to particle radius.The neighborhood of maximum deflection (a = 205 nm), on the other hand, is less desirable for sorting.Although the deflection is large, a wide range of particle sizes end up at the same position in the output plane.
More generally, values of a for which the force has an extremum or a zero-crossing also correspond to zero dispersion.There are two such values for the range covered in Fig. 10, but there are many more such values as a increases because the force continues to alternate between positive and negative values.Although zero-dispersion is generally to be avoided, these points, particularly the zero-crossing points, offer the possibility of creating bimodal size distributions.For example, an input distribution covering the range a = 250 nm to a = 300 nm, will result in bimodal distributions across the output.At z 1 = 29 μm, the distribution will have modes with mean values a = 272 nm and a = 278 nm.Unimodal size distributions, including monodisperse distributions, can be created by utilizing regions of the deflection curve farther from the extrema and zero-crossing points.
The sorting resolution depends on, in addition to dispersion dz 1 /da, the width of the spatial distributions created by each particle size, the distributions illustrated in Fig. 8.More specifically, the resolution is the standard deviation σ a of the size distribution as a function of particle radius (as opposed to position) at a given location in the output plane.Resolution has the units of length.In addition, we define resolving power R as the ratio where σ a and a are the standard deviation and mean of the relevant size distribution.The resolving power as a function of particle radius for I 0 = 0.4 GW/m 2 is shown in Fig. 11.
Resolving powers larger than 100 are indicated in regions near the zero-crossing point.There are much wider ranges of particle size for which the resolving power exceeds 10.In these regions, resolving powers in excess of 100 can be obtained by combining two sorting stages in series.A resolving power of 100 is sufficient for many applications.Similar resolving power has been reported for an array of optical tweezers [16].
We have discussed the deflection and resolving power as though the features in those curves are fixed with respect to particle size.In fact, high resolving power regions can be shifted relative to particle size by changing the optical wavelength.This allows for considerable flexibility in creating and measuring size distributions.
Finally, we consider particle throughput.There are two factors affecting throughput that are relevant to this specific geometry for sorting particles.One is perturbation of the optical field by the particles being sorted and another is the maximum flux of particles that can be put through a point-like source.
We imagine obtaining a standing wave of sufficient intensity inside of a resonant optical cavity.The power loss from the cavity by scattering from particles must be small compared to mirror loss so that the power in the cavity is not dependent on the number of particles.Power loss in an empty cavity is P c (1 − R) where P c is the circulating power and R is the mirror reflectivity.Power loss by optical scattering from N particles with scattering cross section σ is Fig. 11.Resolving power versus particle radius for I 0 = 0.4 GW/m 2 .In general, the size distribution is bimodal because two particle sizes, one smaller than the zero-force size (a ≈ 274 nm) and one larger, will be deflected to the same position.The blue curve is obtained using only one mode of the distribution.The green curve includes both.The inset shows the probability that a particle observed at z 1 = −29 μm has radius a, assuming a uniform distribution of input radii.
Assuming a Gaussian mode with beam waist w 0 and particles distributed uniformly across the beam between x = − √ 2w 0 and x = + √ 2w 0 , where we have used erf(2) ≈ 1. Making scattering losses negligible in comparison to mirror losses, gives the constraint Particle throughput Γ is the product of N and the velocity V c / √ 2 across the optical field divided by the distance, which we have taken as 2 √ 2w 0 .Thus, Using the values R = 0.999, a = 50 nm giving σ = 1.7 × 10 −17 m 2 , V c = 1 × 10 −3 m/s, and w 0 = 0.5 mm gives the constraint Γ 3 × 10 7 s −1 .Commercially-available aerosol generators are generally limited to particle densities of 10 14 m −3 .This limits the flux of particles that can be generated by a source of finite size.Considering a tube source with cross-sectional radius r = 18 μm and a cross-sectionally averaged flow velocity of 1 mm/s, the corresponding flux is Γ = 10 2 s −1 .This is not a fundamental limit as there are a number of changes that can be made to increase this limit.Particle sources, as many as 100, can be placed in parallel along the length of a cavity, as long as sufficient separation is maintained.The finesse of the cavity can be increased by an order of magnitude, allowing the flow velocity and flux to be increased proportionally.It may also be possible to increase the particle density, by as much as an order of magnitude, for the short times during which the particles are being sorted.
A throughput of 10 4 s −1 was referenced by Ladavac, Kasza, and Grier [16] for sorting of colloids by an array of optical tweezers, but the work cited [48] did not discuss a throughput higher than a few colloids per second.Recent efforts on sorting a continuous stream of colloids have achieved throughputs of tens of particles per second [7,12,29,49,50,51].Schemes involving arrays of optical tweezers are inherently limited in throughput.Particles being sorted in a tweezer array are more likely to be found in a trap than between traps [48] and trap volumes are not large enough to interact with more than 1 or 2 particles simultaneously.

Conclusions
The trajectories of aerosol particles in air flowing through an optical field are determined by the interplay of optical forces, convection, and diffusion.The latter two are produced by random collisions of the particles with air molecules.
We have derived, from first principles, the optical force on spherical dielectric particle of arbitrary radius in a standing wave.For a/λ 0 → 0, our evaluation of the Maxwell stress tensor is equivalent to the expression denoted as the gradient force by Harada and Asakura [32,37].For larger radii, the force alternates between positive and negative values with decreasing amplitude, a result that has been observed previously [47].
We have evaluated collisional forces using two distinct approaches, implemented in isolation from each other, in order to validate our results.The first method integrates the microscopic equations of motion to determine individual particle trajectories for large populations of particles.Brownian motion is incorporated with Monte Carlo methods.The second method solves the macroscopic fluid equations of motion.Computational efficiency is achieved by averaging velocity and the diffusion coefficient over the fine-scale features of the optical field before integrating the macroscopic equations of motion.These methods have enabled us to explore particle trajectories in air flowing through a standing wave Gaussian mode optical field.
We have presented and demonstrated a relatively simple optical scheme capable of sorting particles with high throughput and high resolution resolving power.This scheme diverges from the paradigm of high numerical aperture optics based on optical tweezers.

A. Fokker-Planck equation without diffusion
If we neglect diffusion, or equivalently, consider T = 0, we may obtain an analytic result for the effective velocity of a point particle traversing the interference fringe.If we fix the drift velocity v 0 > 0, there is some sufficiently small value of the optical potential I 0 for which the optical well no longer traps the particle.In this case, the effect of the optical potential is to modify the average velocity of travel through the potential, but not to bring it to zero.The average is that of a point particle moving through a complete spatial period Λ.We introduce v(z) to be the instantaneous speed at position z.The time t 1 to cross an interference fringe is given by We use Eq. ( 40) to obtain the effective velocity If there is no optical force, i.e., F 1 = 0, then v 0 = μF 0 is recovered.As the optical force F 1 increases, the effect is to reduce the average speed.Ultimately, the average speed falls to zero with a square root singularity as F 1 → F 0 from below.In this limit, the particle is marginally trapped.The greatest sensitivity of the effective velocity v eff = Λ/t 1 on small changes in the optical force F 1 occurs near marginal trapping.Our program reproduces the analytic result of Eq. (41).

B. Fokker-Planck equation without drift
If we have no drift but have diffusion, and consider the limit of large well depth, Kramers [52] has given an analytic solution.The transition rate is given by where the potential energy surface is given by near the quadratic minimum at z = z A and near the quadratic maximum at z = z C .The exact form of the potential surface does not enter the quoted result.For this numerical test, we use Stokes' Law, Eq. ( 19), and the Rayleigh force law, Eq. ( 8).The energy of a particle in a light field is given by the integral of the 1D optical force From Eq. ( 45), we may give the energy difference between the top and the bottom of the well as Similarly, by a second order Taylor expansion of the cosine  Fig. 12.The escape rate from the well of a standing wave versus intensity I 0 .Kramers' analytic formula valid for small escape rates is shown in blue.The escape rate found numerically from the effective fluid constant model is in red.The ratio of the two is in the inset for the same range of intensities.
As a validity check, we compare our solutions to those obtained from Kramers' formula for transition rates of a particle over a barrier due to Brownian motion [52].To make this comparison, we need to connect the transition rates to the diffusion constant.Derrida [53] gives the diffusion constant and the effective velocity as and where the dimensional constant Λ is added here.Equation (47) may also be obtained by comparing at large times the Gauss-Weierstrass solution to the heat equation for the continuous case to the discrete Poisson distribution arising from hopping.The k (hop) ± are the hopping rates with (+) and against (−) the convective current.For simplicity, we confine the check to the case without a convective fluid, i.e., v eff = 0 so k . The transition rates found from the Kramers' rate equation discussed above and the Fokker-Planck equation are shown in Fig. 12.The numerical solution is seen to approach Kramers' asymptotic formula for large values of the trapping, which is a key test of our program.

# 263305 Received 14 Fig. 1 .
Fig. 1.The scheme for sorting particles by size consists of a narrow stream of aerosol particles flowing across a Gaussian mode standing wave at an oblique angle.The figure is not to scale; in particular, the optical fringes and particles are much finer than shown.

# 263305 Received 14 Fig. 4 .
Fig. 4. The solution of the Fokker-Planck equation for a = 100 nm, Λ = 532 nm, I 0 = 0.4 GW/m 2 , and v d = 707 μm/s.The solutions occur starting from the initial condition in blue, then in four steps of 64 μs represented sequentially by green, red, gold, and black curves.

Fig. 6 .
Fig. 6.Normalized effective velocity across interference fringes versus optical intensity I 0 for different particle radii in nm (shown next to each curve). z

Fig. 8 .
Fig. 8. Spatial distributions in the output plane (x = 1.5 mm) for each particle size.

# 263305 Received 14 2 Fig. 9 .Fig. 10 .
Fig. 9. Distribution width in the output plane versus particle size with and without an optical field.

Table 1 .
(7)ct values of the c i coefficients in Eq.(7) 2 p ), demonstrating that diffusion in velocity space and in real space are manifestations of the same underlying phenomenon.In our Monte Carlo model

Table 2 .
Parameters used for converged solution of Eq. (29).Small a means 40 nm ≤ a ≤ 110 nm; medium a means 120 nm ≤ a ≤ 170 nm; large a means 180 nm ≤ a ≤ 300 nm

Table 3 .
Parameter values used in simulations