Optical Sorting of Dielectric Spheres in a Standing Wave Field Calculated with Effective Velocities and Diffusion Constants

Zachary H. Levine
J. J. Curry

Recently, as part of a study of the ability of a Gaussian standing-wave interference field to sort small dielectric particles entrained in a moving fluid, we implemented a two-length-scale method to determine the trajectories of the particles.  In the first step, effective velocities and diffusion constants are determined as a function of field intensity by solving the Fokker-Planck equation in a force field determined by the Maxwell stress tensor.  In the second step, these effective constants are used to predict positions of particles having been swept across the interference field.  Here, we present the code which was used ito implement that method.

0. Overview

Quick start

The full program, prepared in Version 10.0, can be launched here.  The reader may prefer to execute the code in stages.  Files are written to the working directory so that the program will automatically use previously computed results.  The total run time for the default parameters is about 1 hour, of which about 98% is in the first section, solveFokkerPlanck.  Most of the remaining time is used by makeParticleTrace.  If the reader needs to adjust the plots, the feature of reading from saved data is particularly valuable.  The data from each section are written to files which are read by subsequent sections.

doJob = {True, True, True, True, True} ; (* Which of the parts below should be run *)
quickStart[] = (              (* The reader can evaluate the whole notebook here. *)
  If [ doJob[[1]]
  , Print[ plotFuchs[] ]
  ]
; If [ doJob[[2]]                 (* described in Section 1 (subsection proper) *)
  , solveFokkerPlanck[]         (* writes velDiffuse000.dat or similar *)
  ]
; If [ doJob[[3]]                 (* reads velDiffuse000.dat or similar *)
  , Print[ plotVelDiffuse[] ]   (* described in Section 2 (subsection plot) *)
  ]
; If [ doJob[[4]]
  , makeParticleTrace[]         (* described in Section 3 (subsection proper) *)
  ]
; If [ doJob[[5]]
    , Print[ plotParticleTrace[] ]  (* described in Section 3 (subsection plot) *)
  ]
) ;

Introduction

Recently, we examined the problem of sorting small dielectric spheres by radius.[1,2] The spheres have radii of about 50 nm to 300 nm by radius.  The scheme consists of introducing the spheres into air which flows diagonally across a standing-wave optical laser field as shown in Figure 1.  Particles with different radii experience different optical forces that lead to physical separation on the far side of the field.  The consequences have been worked out theoretically for a range of radii and field strengths.   In Ref. [1], we use a two-length-scale analysis to find the effective velocities and diffusion constants for these particles.  On the microscopic scale, we can consider the position of a sphere relative to a one-dimensional well representing the standing-wave laser field.  One dimension is sufficient because the macroscopic forces orthogonal to the optic axis are orders of magnitude weaker than those parallel to the optic axis.  Our mathematical approach is to construct the probability distribution of positions of the particle as a function of time.  The time rate of change of this distribution settles  to a constant value which we call an effective velocity.  A similar calculation to find the proportionality constant for the linear growth of the variance of the distribution in time yields an effective diffusion constant.   After a particle goes from one well to the next, it finds itself in essentially the same physical situation, since the laser field is periodic.   This allows us to be ignorant of the position of the particle on the scale length of an interference fringe, yet integrate these effective constants across an optical field with a slowly varying envelope  to find the distribution of particles at the exit.

The purpose of the present paper is to provide the codes used in Ref. [1] to calculate the effective velocities and diffusion constants as well as the sorting ability of a weakly-focused Gaussian beam.  The reader is encouraged to look at Ref. [1] for additional context and results.  In each section, we present the formalism, provide code to implement the formalism, and show some results.  The three sections describe the following:  (1) the evolution of a probability distribution; (2) the extraction of the effective velocity and diffusion constants from the probability distribution; and (3) the use of these constants to predict the capability of the system described above to sort small dielectric spheres.  Some additional constants and functions are given in the appendix.

By default, the code reproduces some of the results from Ref. [1].  The reader may change parameters in Sections 1 and 3 to explore other cases.

1. The Fokker-Planck Equation in a Fourier Basis

The one-dimensional probability distribution describing the position of a particle is taken to evolve according to a Fokker-Planck equation

((∂f(z,t))/(∂t))=D ((∂f(z,t))/(∂(z^2)))-μ(∂/(∂z))[F(z)f(z,t)]

where t is time, z is the position in one dimension (also the position along the optic axis), D is the diffusion constant for the sphere and is a function of the radius of the sphere and the visocity of the medium, μ is the mobility of the particle, and F(z) is a static force on the particle.  If the force is smooth and periodic, expanding f in a Fourier basis is efficient.  Let

f(z,t) = (1/(mΛ))(∑_(n=-∞)^∞)(Subscript[c, n])(t) (e^(i n k z))

where  Λ=2λ is the period of the interference field for a laser of wavelength λ, m is the number of periods of the interference field retained in the solution, k=2π/(mΛ), and the (Subscript[c, n])(t) are complex expansion coefficients.  In practice, m needs to be set sufficiently large so that f(mΛ/2,t)  is negligible.   For fixed m, this constraint leads to an upper limit on the range of validity of t because f(mΛ/2,t) eventually becomes non-negligible as a result of drift or diffuse or both.  Because f(z,t) is real, (Subscript[c, -n])=(Subscript[c, n]^( * )).   The specification of one particle in the solution domain of size mΛ at all times implies (Subscript[c, 0])=1.  Using these considerations, f may be simplified to

f(z,t) = (1/(nΛ))(1 + 2 Re(∑_(n=1)^∞)(Subscript[c, n])(t) (e^(i n k z))).

In the case of constant flow across a standing-wave interference field,

F(z) = (Subscript[F, 0])+(Subscript[F, 1])sin(2 π z/Λ)

where (Subscript[F, 0]) is a position-independent force due to convection and (Subscript[F, 1]) is the optical force.  If the background fluid has velocity (Subscript[v, d]), then (Subscript[F, 0])=(Subscript[v, d])/μ.  The optical force constant is more complicated, and we refer the reader to companion articles [1,2] for its determination.  Substituting the expansion of Eq. (2) into the Fokker-Planck Equation (1) leads to the system of ordinary differential equations

(d/(d t))(Subscript[c, n])(t) = -D (n^2)(k^2)(Subscript[c, n])(t) + i (Subscript[v, d])n k (Subscript[c, n])(t) -(μ/2)(Subscript[F, 1])n k (Subscript[c, n-m])(t)+(μ/2)(Subscript[F, 1])n k (Subscript[c, n+m])(t).

The fact that(Subscript[c, -n])=(Subscript[c, n]^( * )) means that we can, after the obvious substitution, solve only for (Subscript[c, n]) with n>0.  The equations are truncated by approximating (Subscript[c, n])=0 for sufficiently large n, i.e., for n>nTerm.   Eq. (5) decouples into m independent sets of equations.  In solid-state physics, such decoupling is known as Bloch’s theorem.  The decoupled equations are tridiagonal.  Hence, we can find solutions rapidly.  In practice, we simply set up the equations and pass them to NDSolve without explicitly decoupling the equations in the code.

This section is organized with two subsections:  preliminaries, where user-adjustable constants are given, and proper, which presents code to solve the equation.  The code relies on the Appendix to define (Subscript[F, 1]), μ, and D.

Code: the Fokker-Planck equation in a Fourier basis (preliminaries)

First, we set some constants which the reader may wish to overwrite in the second box below.  (A similar set of constants is given in Section 3.) The values given were used in Ref. [1] and can be plotted with the specific plot functions provided.   The reader is cautioned that changing most parameters requires some care.  For example, the Fuchs mobility and diffusion constants are specialized to air at 25 °C.  The viscocity is also only given for that temperature.  The index of refraction and the wavelength may be changed freely, although numerical issues could arise.  

The parameters powAreaAll and rPartAll are the set of powArea (i.e., optical intensity) and rPart (particle radius a) to solve for.

The reader must judge whether the Stokes force or the Fuchs relation is better in a particular case.  The optical force expressions are good for a particle radius up to a few hundred nanometers.[1]  Changing the numerical parameters requires monitoring the solution to ensure that it does not have spurious oscillations (requiring a greater value for nTerm) or significant amounts of probability at the ends of the domain (requiring larger m and larger nTerm).  The time step dt may need to be decreased if smaller particle radii are of interest.  On the other hand, dt may be increased for more efficiency when the particle radii are larger.

(* thermal parameters *)
temperatureC = 25. ;               (* temperature in Celsius *)
nu = 18.6 micro kg / (meter sec) ; (* viscosity of air at 25 degrees Celsius *)

(* optical parameters *)
indexOfRefr = 1.57 ;               (* index of refraction of polystyrene *)
lam = 1.064 micro meter ;          (* wavelength of light:  Nd:YAG (yttrium aluminum garnet) laser *)

(* cases to solve for *)
velAll     = {707} micro meter / sec ;                    (* velocity of particle in z direction *)
powAreaAll = Range[0, 400, 20] mega watt / meter^2 ;        (* optical power levels to solve for *)
rPartAll = Range[210, 270, 10] nano meter ;                   (* radii of particles to solve for *)

(* program options *)
flagFuchs = 1 ;                  (* 0 use Stokes' Law, 1 use Fuchs' expression for mobility *)
flagMie   = 1 ;                  (* 0 use Rayleigh limit, 1 use Mie expression Refs. [1,2] *)

(* program parameters *)
m = 32 ;                           (* number of periodic replicas   N in the formalism section *)
mExc = 3 ;                         (* number of periods to exclude when making moments *)
nTerm = 512 ;                      (* number of terms in the expansion to find *)
dt = 4 micro sec ;                 (* time step for solution of Fokker-Planck equation *)
nStep = 2^10  ;                    (* number of time steps in Fokker-Planck equation *)
nStepShow = 4 ;                    (* number of time steps to display (not counting t=0) *)
nPart = 2^5 ;         (* number of time steps for averaging to get velObs, diffuseObs *)
                      (* it is best if  nStep / nPart is an integer *)
imageSize = 125 mmToPoint ;        (* size of certain plots in millimeters *)
nPlotFokkerPlanck = 10 ;           (* maximum number of plots of Fokker-Planck solution *)

It is recommended that the reader leave the box above alone and overwrite modified values below.

(* reader modifications for Sections 1 and 2 go here *)
(*  for example run a single power and radius with moderate accuracy  *)
(*
powAreaAll = {100} mega watt / meter^2 ;
rPartAll = {210} nano meter ;
m = 16 ;
nStep = 2^9 ;
*)

We do a little more processing on these constants.  The dielectric constant of the medium is implicitly taken to be unity.

(* thermal parameters *)
temperatureK = (temperatureC + 273.15) kelvin ; (* temperature in Kelvin *)
kT = kBoltz temperatureK ;    (* thermal parameter   kelvin not defined but cancels here *)

(* derived logical variables *)
flagStokes = 1-flagFuchs ;            (* logical not on 0,1 Boolean variables *)
flagRayleigh = 1-flagMie ;

(* optical constants *)
bigLam = lam / 2 ;          (* period of interference fringe from wavelength of light *)
kWave = 2 Pi / lam ;        (* wavevector of light *)

(* solution domain constants *)
mInner = m-mExc ;      (* number of periods to retain for moments *)
                       (* of period of external potential *)
k = 2 * kWave / m ;    (* wave vector for periodic replicas *)
domainSize = m bigLam ;               (* size of simulation domain *)
innerDomainSize = mInner bigLam ;     (* size of useable simulation domain *)
nPeriodPlot = mInner / m ;            (* number of periods to plot *)
zMax = domainSize * nPeriodPlot / 2 ; (* highest co-ordinate in plot *)
sigma = bigLam / 8 ;                  (* standard deviation for initial Gaussian *)
nForAvg = Ceiling[nPart/2] ;   (* number of points to average over for effective velocity and diffusion constant *)
kPlotFokkerPlanck = 0 ;      (* counter of number of plots of Fokker-Planck distribution made *)
pFokkerPlanck = ConstantArray[{},nPlotFokkerPlanck]  (* plots of the Fokker-Planck distribution to be made *)

Global arrays are defined next.  The real Fourier basis basisCos and basisSin  will help represent the solution of the Fokker-Planck equation in Section 2.  We also find coefficients used in Section 2 to find the zeroth, first, and second moments of the distribution (momCoeff) which are a key step in determining (Subscript[v, eff]) and (Subscript[D, eff]).  We also define the set of dependent variables depVar, which correspond to those (Subscript[c, n])(t) (c[n] in the code) that we will solve for explicitly.

basisCos = Table [ Cos[n k z], {n,nTerm}] ;
basisSin = Table [ Sin[n k z], {n,nTerm}] ;
momCoeff = setMomCoeff[ k, nTerm, domainSize, mInner*bigLam ] ;
depVar = Table[ c[n], {n,nTerm}] ;

Code: the Fokker-Plank equation in a Fourier basis (proper)

The Fokker-Planck equation is solved for all desired cases by invoking solveFokkerPlanck.  It takes no arguments because it relies on global variables, set above.  For a particular optical intensity and particle radius (and fixed, global parameters) the solution is given as solFokkerPlanck.  The effective velocity and diffusion constants are extracted from it, as described in the next section.  Ultimately, these results are saved in a file because they can be used for a variety of macroscopic standing-wave optical field configurations without recomputing.  The output is a simple ASCII file with each case on its own line with a serial number such as velDiffuse000.dat.  The output file name is written to the file makeParticleTraceInputFileName.txt, which is read by both plotVelDiffuse and makeParticleTrace.  By default, these programs read the last file written, but makeParticleTraceInputFileName.txt can be modified by hand if the codes are invoked manually, i.e., not via quickStart.

solveFokkerPlanck[] :=
  Module[{ outputFileName, output, vel, powArea, rPart, solFokkerPlanck
         , diffuse, velEff, diffuseEff, powAreaMega, rPartNano, velMicro }
  , outputFileName = getNextFileName["velDiffuse"]
  ; Export[outputFileName,{}]                                    (* empty file created to avoid name conflicts *)
  ; output = ConstantArray[{}, Length[rPartAll] Length[powAreaAll] Length[velAll] ]
  ; iCase = 0
  ; Do [
      vel = velAll[[iVel]]
    ; Do [
        rPart = rPartAll[[iRPart]]
      ; Do [
          powArea = powAreaAll[[iPowArea]]
        ; {solFokkerPlanck, diffuse} = solveFokkerPlanckA[ vel, powArea, rPart ]
        ; {powAreaMega, rPartNano, velMicro} = Round[ {powArea/mega, rPart/nano, vel/micro} ]
        ; {velEff,diffuseEff} = getVDeff[ solFokkerPlanck, powAreaMega, rPartNano, velMicro ]
        ; output[[++iCase]] = { flagFuchs, flagMie     (* flags *)
                            , powAreaMega            (* MW/m^2 *)
                            , rPartNano              (* nm *)
                            , velMicro               (* μm/s *)
                            , velEff/micro           (* μm/s *)
                            , diffuse / micro^2      (* μm^2/s *)
                            , diffuseEff / micro^2   (* μm^2/s *)
                            }
        , {iPowArea, Length[powAreaAll]}]        (* end loop over optical intensities *)
      , {iRPart, Length[rPartAll]}]            (* end loop over particle radii *)
    , {iVel, Length[velAll]}]        (* end loop over optical intensities *)
  ; export[ outputFileName, output ]
  ; Export["makeParticleTraceInputFileName.txt", outputFileName <> "\n"]
  ]

In a particular case, the Fokker-Planck equation is solved by solveFokkerPlanckA.

We begin somewhat abstractly by noting that a system of differential equations consists of the differential equations diffEq and their initial conditions initCond.  

The initial condition, taken to be a normal distribution with mean z=0 and standard deviation σ=Λ/8, is given in the Fourier basis for the inital values initVal.   Physically, this means the initial solution is almost entirely contained in the well centered at z=0, but extends over a substantial portion of that well.

The formal association of the initial values initVal with the parameters (Subscript[c, n])(0) or c[n][0]  is made with the variable initCond.   Only the values n=1, ..., nTerm are given by initVal.

Next, we set the optical force on the particle force1, which is (Subscript[F, 1]) in Eq. (5).  Optical forces have been discussed in by the authors elsewhere;[1,2] see also the Appendix.

We define a0 which is |(Subscript[μF, 1])/2|.  If (Subscript[F, 1])<0, physically, this means that the particle interchanges peaks and troughs.  By taking the absolute value, we ensure that the particle starts in a trough.

The differential equation [FokkerPlanckPlaneWave] is given in the first statement below.  The second statement encodes (Subscript[c, 0])(t)=1.  The third statement encodes (Subscript[c, -n])(t)=(Subscript[c, n]^( * ))(t).  The final statement truncates the series (Subscript[c, n])(t)  for n>nTerm.  The effect of the final three statements is to define the set of (Subscript[c, n])(t) which are not defined in initCond as 0, 1, or in terms of coefficients defined in initCond.  This closes the set of differential equations.

The equations are solved with NDSolve.

solveFokkerPlanckA[ vel_, powArea_, rPart_ ] :=
  Module[{system, initVal, initCond, force1, mobility, diffuse, a0, diffEq}
  , system = Flatten [ {diffEq, initCond} ]
;Print["solveFokkerPlanckA:  powArea ", powArea, " W (m^(-2))  rPart ", rPart, " m" ]
  ; initVal = Table [ Exp[ -k^2 n^2 sigma^2 / 2 ], {n,nTerm} ]
  ; initCond = Table [ c[n][0] == initVal[[n]], {n,nTerm} ]
  ; force1 =
      Which[
        flagRayleigh == 1
        , powArea getForceScaledRayleigh[rPart, kWave, indexOfRefr]
      , flagMie == 1
        , powArea getForceScaledMie[rPart, kWave, indexOfRefr]
      , True
        , Print["error:  force0 cannot be set"] ; Abort[]
      ] ;
  ; {mobility, diffuse} = 
      Which[
        flagFuchs == 1
        , {mobilityFuchs[rPart], diffuseFuchs[rPart]}
      , flagStokes == 1
        , {mobilityStokes[rPart], diffuseStokes[rPart]}   
      , True
        , Print["error:  mobility and diffuse cannot be set"] ; Abort[]
      ] ;
  ; a0 = Abs[ mobility force1 / 2 ]
  ; diffEq = Table [ c[n]'[t] == - diffuse  k^2 n^2 c[n][t] -  I vel k n c[n][t]
                                  + a0 k n c[n-m][t] - a0 k n c[n+m][t]
             , {n,nTerm} ]
  ; diffEq = diffEq /. {c[0][t]->1}                                      (* particle conservation *)
  ; diffEq = diffEq /. Table [ c[-n][t]->Conjugate[ c[n][t] ]           (* symmetry for real functions *)
                       , {n,m-1}]
  ; diffEq = diffEq /. Table [ c[n][t]->0, {n,nTerm+1,nTerm+m}] // Chop  (* truncation for large n *)
  ; { NDSolve [ system, depVar, {t, 0 sec, dt nStep} ], diffuse }
  ]

2. Effective velocities and diffusion constants from Fokker-Planck solutions

Our stategy here is to find the mean position ⟨z(t)⟩ and variance (σ^2)(t)  from the moments of the solution f(z,t) found in the previous section.  After an initial equilibration period, the first two moments grow linearly in time.  We will find (Subscript[v, eff])= (d/(dt))⟨z(t)⟩ and (Subscript[D, eff])=(1/2)(d/(dt))(σ^2)(t) at sufficiently large times such that these values may be approximated as constants.  By default, the first half of the period for which solutions are found is devoted to equilibration and the second half to obtaining the averages for the effective velocity and diffusion constants.

Code:  effective velocity and diffusion constants (preliminaries)

In order to find the zeroth, first, and second moments of the solution, we find corresponding values for each basis function in setMomCoeff.   This can be done globally.  The variables defined by setMomCoeff are used by setMeanAndVar which contains a contribution from the (Subscript[c, 0])=1 term fullRange^3/12.  The latter is related to the second moment of the uniform distribution.  Only the cosine basis functions contribute to the zeroth and second moments, and only the sine basis functions contribute to the first moment.

setMomCoeff[ k_, nTerm_, domainSize_, fullRange_ ] :=                   (* momCoeff defined *)
  Module[{a, i0, i1, i2, mom0Coeff, mom1Coeff, mom2Coeff}
  , i0 =  2 Integrate [     Cos[ a z ], {z,-fullRange/2,fullRange/2}]
  ; i1 = -2 Integrate [ z   Sin[ a z ], {z,-fullRange/2,fullRange/2}]
  ; i2 =  2 Integrate [ z^2 Cos[ a z ], {z,-fullRange/2,fullRange/2}]
  ; mom0Coeff = i0 /. {a->k Range[nTerm]}
  ; mom1Coeff = i1 /. {a->k Range[nTerm]}
  ; mom2Coeff = i2 /. {a->k Range[nTerm]}
  ; {domainSize,fullRange,mom0Coeff,mom1Coeff,mom2Coeff}
  ]

setMeanAndVar[ coeff_, momCoeff_] :=                                     (* momCoeff used *)
  Module[{domainSize,fullRange,mom0Coeff,mom1Coeff,mom2Coeff,norm
         ,zSqMean,zMean,zVar,probInFullRange}
  , {domainSize,fullRange,mom0Coeff,mom1Coeff,mom2Coeff} = momCoeff   (* unpack *)
  ; norm = fullRange + Re[coeff] . mom0Coeff                         (* divide by norm to get probability norm *)
  ; zSqMean = (( (fullRange^3)/12 + Re[coeff] . mom2Coeff )) / norm   (* <z^2> *)
  ; zMean   = Im[coeff] . mom1Coeff / norm                            (* <z> *)
  ; zVar    = zSqMean - zMean^2                                       (* <z^2>-z^2 *)
  ; probInFullRange = norm / domainSize                               (* should be 1 or very close to it *)
  ; {probInFullRange,zMean,zVar}
  ]

Code: effective velocity and diffusion constants (proper)

We are expecting a single solution to the Fokker-Planck equation.  Typical results are shown in Figure 2.  It is unpacked from solFokkerPlanck to solRule.  The dependent variables depVar are described in the section above.  The numbers in the solution are extracted from their form as rules and placed into solVec.  There must be nTerm of the (Subscript[c, n]) defined.  The solution is repackaged into a function of time by solOfT, which in turn is used to define ftnOfT which is f(z,t) from Eq. (3), omitting the normalization constant 1/mΛ.  The array tStat is a set of times for which the moments probInFullRange, zMean, and zVar are found after the numerical values of (Subscript[c, n])(t) are placed into coeff.  The mean positions are grouped into time intervals in zMeanA.  With each time interval, the position as a function of time is fit with a line.  Finally, the velocity observed velObs in each time interval is found by taking the time derivative.  The variance is processed in an analogous manner yielding an observed diffusion constant diffuseObs in each time interval.

The observed velocities and diffusion constants are not quite the effective constants.  If we have set the solution parameters properly, the solution relaxes to a situation in which the mean position and the variance both grow linearly in time.   (If the solution does not run for enough time, intrawell relaxation modifies the behavior.  However, if the solution runs too long, then the probability unphysically wraps from the right end of the solution domain to the left end.)  The values of the observed velocity and diffusion constant for the times in which they are nearly constant are averaged to make the effective velocity and diffusion constant.

getVDeff[ solFokkerPlanck_, powAreaMega_, rPartNano_, velMicro_ ] :=
  Module[{ solRule, solVec, solOfT, ftnOfT, tStat, zMeanA, coeff, probInFullRange, zMean, zVar
         , zMeanFitTab, velObs, velEff, iT, zVarA, zVarFitTab, diffuseObs, diffuseEff}
  , If [ Length[ solFokkerPlanck ] ≠ 1
    , Print["getVDeff:  more than one solution to the Fokker-Planck equation ... unexpected"]
    ; Abort[]
    ]
  ; solRule = First[ solFokkerPlanck ]
  ; solVec = depVar //. solRule
  ; If [ Length[solVec] != nTerm
    , Print["getVDeff:  bad length for solVec ", Length[solVec], " expected ", nTerm ]
    ; Abort[]
    ]
  ; solOfT[t_] := Table [ solVec[[n]][t], {n,nTerm}]
  ; ftnOfT[t_] := 1 + 2 ( Re [ solOfT[t] ] . basisCos
                        - Im [ solOfT[t] ] . basisSin )
(* -------------------- optional plot ftnOfT here is desired ---------------- *)
  ; If [ kPlotFokkerPlanck < nPlotFokkerPlanck
    , tShow = (nStep/nStepShow) Range[0,nStepShow] dt
    ; ftnOfZ = Table[ ftnOfT[ tShow[[i]] ], {i,Length[tShow]}]
    ; ftnOfZnm = ftnOfZ /. {z->zNM nano}
    ; pFokkerPlanck[[++kPlotFokkerPlanck]] = 
        Plot[ ftnOfZnm
        , {zNM,-bigLam /nano, 9/2 bigLam / nano}
        , FrameStyle->Thick
        , FrameLabel->{ { "Particles per domain size", ""}
                      , {  "z (nm)"
                        , Row[{ "(Subscript[I, 0])=" <> IntegerString[powAreaMega] <> " MW/"
                              , Superscript["m","2"]
                              , " a=" <> IntegerString[rPartNano] <> " nm"
                              , " (Subscript[v, d])=" <> IntegerString[velMicro] <> " μm/s"}]
                        }
                      }
        , ImagePadding->{{50,10},{100,100}}
        , LabelStyle->{Directive[Bold,Medium],FontFamily->"Helvetica"}
        , PlotStyle->{Green, Orange, Blue, Brown, Black, Purple}
        ]
 ; Print[ pFokkerPlanck[[kPlotFokkerPlanck]] ]
    ]       (* -----------------end optional plot ftnOfT here is desired ---------------- *)
  ; tStat = Range[nStep] dt                                               (* statistics collected at these times *)
  ; coeff = Table [ solOfT[ tStat[[i]] ], {i,Length[tStat]}]              (* coefficients c_n(t) *)
  ; {probInFullRange,zMean,zVar} = setMeanAndVar[ coeff, momCoeff ]
  ; zMeanA = Partition [ zMean, Floor[nStep/nPart] ]                      (* positions yield effective velocity *)
  ; zMeanFitTab = Table [ Normal[ LinearModelFit [ zMeanA[[i]], iT, iT ] ]
                 , {i,nPart} ]
  ; velObs = D [ zMeanFitTab, iT ] / dt
  ; velEff = Mean[ Take[ velObs, -nForAvg ] ]
  ; zVarA = Partition [ zVar, Floor[nStep/nPart] ]             (* variances yield effective diffusion constants  *)
  ; zVarFitTab = Table [ Normal[ LinearModelFit [ zVarA[[i]], iT, iT ] ]
                 , {i,nPart} ]
  ; diffuseObs = 0.5 D [ zVarFitTab, iT ] / dt
  ; diffuseEff = Mean[ Take[diffuseObs,-nForAvg] ]
  ; {velEff, diffuseEff}
  ]

Code: effective velocities and diffusion constants (plot)

Here, we provide code to plot the velocities and diffusion constants given default input.   The reader will almost certainly need to adapt these functions for other cases.

plotVelDiffuse[] :=
  Module[{ inputFileName, datRaw, flagFuchs, flagMie, powAreaMega
         , rPartNano, velMicro, velEffMicro
         , diffuseMicroSq, diffuseEffMicroSq, velRel, rPartNanoUniq, nRPartNano
         , powAreaMegaUniq, nPowAreaMega, velMicroUniq, nVelMicro
         , velRelB, p1a, diffuseEffMicroSqB, p1b}
  , inputFileName = Import["makeParticleTraceInputFileName.txt"]
  ; datRaw = Import[inputFileName]
  ; { flagFuchs, flagMie, powAreaMega, rPartNano, velMicro, velEffMicro
    , diffuseMicroSq, diffuseEffMicroSq} = Transpose[datRaw]
  ; velRel = velEffMicro / velMicro
  ; rPartNanoUniq = Union[rPartNano]
  ; nRPartNano = Length[rPartNanoUniq]
  ; powAreaMegaUniq = Union[powAreaMega]
  ; nPowAreaMega = Length[powAreaMegaUniq]
  ; velMicroUniq = Union[velMicro]
  ; nVelMicro = Length[velMicroUniq]
  ; velRelB = Partition[ pair[powAreaMega,velRel], nPowAreaMega ]
  ; p1a = 
      ListPlot [ velRelB
      , FrameLabel->{ Row[{"Intensity (MW/", Superscript["m","2"], ")"}]
                    , Row[{Subscript["v","eff"], "/", Subscript["v","d"]}]
                    , Row[{Subscript["v","d"], " =", IntegerString[ First[velMicroUniq] ], " μm/s"}]
                    , ""
                    }
      , FrameStyle->Thick
      , ImageSize->imageSize
      , Joined->True
      , PlotStyle->{Black, Red, Green, Blue, Orange, Cyan, Purple}
      , Epilog->{
          inset["210",  20, 0.35]
        , inset["260", 200, 0.25]
        , inset["270", 300, 0.70]
        }
      ]
  ; diffuseEffMicroSqB = Partition[ pair[ powAreaMega, diffuseEffMicroSq ], nPowAreaMega ]
  ; p1b =
      ListPlot [ diffuseEffMicroSqB
      , FrameLabel->{ Row[{"Intensity (MW/", Superscript["m","2"], ")"}]
                    , Row[{ Subscript["D","eff"], " (",
                            Superscript["μm","2"], "/s)"}]
                    }
     , FrameStyle->Thick
     , ImageSize->imageSize
     , Joined->True
     , PlotStyle->{Black, Red, Green, Blue, Orange, Cyan, Purple}
     , Epilog->{
         inset["210",  37, 40]
       , inset["260", 245, 20]
       , inset["270", 300, 60]
       }
     ]
  ; GraphicsRow[{p1a,p1b}]
  ]

The default input parameters which include the radii 210 nm, 220 nm, ..., 270 nm lead to Figure 3.   If the reader changes the input parameters, at a minimum, it will be necessary to change the positions of the labels.  The labels represent the particle radius in nanometers.  For the default parameters, the force vanishes for a radius of 275 nm.  As this radius is approached, the particles are affected less by the standing wave optical field in both drift velocity and diffusion constant.

3. Application of effective velocities and diffusion constants to optical sorting

In principle, (Subscript[v, eff]) and (Subscript[D, eff])can be used in a variety of optical landscapes based on slowly varying interference fields.  Here, we restrict attention to the situation discussed in Ref. [1] where the optical intensity envelope varies transversely as a Gaussian but does not vary along the optic axis z.  For weakly focused beams, this can be a good approximation.  As in Ref. [1], spherical dielectric particles are entrained in a gentle wind and blown across a standing-wave optical interference field.  We are interested in finding the ability of this system to sort the particles by radius.  Thus, we are interested in the mean and standard deviation of the position at exit.  We assume the particles move uniformly transverse to the optic axis (in the x direction) and remain in the y=0 plane.  The mean position in z compared to where the particles would arrive without the optical field is given by

(Subscript[z, 1]) = (Subscript[v, x]^(-1))(∫_(Subscript[x, 0])^(Subscript[x, 1]))((Subscript[v, eff])(x) - (Subscript[v, d]))dx.

Here, the particle is released at the position ((Subscript[x, 0]),0,(Subscript[z, 0])) and travels to ((Subscript[x, 1]),0,(Subscript[z, 1])). Typically, (Subscript[z, 1])<0, (Subscript[x, 0])=-(Subscript[x, 1]), and the beam is focused at (0,0,0).

As discussed in Ref. [1], there are two contributions to the variance of the distribution describing the particle on exit.  Considering the particle as being in uniform motion, the effective diffusion constant leads to a contribution

((((Subscript[σ, z]^((D))) ))^2)= 2(Subscript[v, x]^(-1))(∫_(Subscript[x, 0])^(Subscript[x, 1]))(Subscript[D, eff])(x) dx.

A second contribution arises due to the diffusion in x, which leads to variations in the crossing times (Subscript[σ, t]).  When this variation is small compared to the crossing time itself, the drift in z acquires a contribution to its standard deviation given by

(Subscript[σ, z]^((v))) = ⟨(Subscript[v, eff])⟩(Subscript[σ, t])

where the effective velocity is averaged over the crossing.  The two processes are independent, so the total variance in the z position of the sphere on output is given by

(Subscript[σ, z]^2)= ((((Subscript[σ, z]^((D))) ))^2)+((((Subscript[σ, z]^((v))) ))^2).

Code: application of effective velocities and diffusion constants to optical sorting (preliminaries)

w0 = 0.5 milli meter ;      (* Gaussian beam width at 2σ of intensity ... used in Section 3. *)
dxdz = 1 ;                  (* slope of fluid velocity *)
x1 = 3 w0 ;                 (* stopping plane *)
x0 = -x1 ;                  (* starting plane *)
z1 = x1  / dxdz ;           (* stopping z without optical field *)
z0 = x0  / dxdz ;           (* starting z *)
powAreaPlotMin = -1. ;      (* the reader can set these constants *)
powAreaPlotMax = -1. ;      (* negative values lead to the default behavior *)
deltaPowAreaPlot = -1. ;

It is recommended that the reader leave the box above alone and overwrite modified values below.

(* reader modifications for Section 3 go here *)

The optical intensity profile has the form

I(x) = (Subscript[I, 0])exp( - ((x^2)/(2 ((((Subscript[w, 0])/2))^2))) ) .

The expression is intentionally not simplified to emphasize that (Subscript[w, 0])/2 is the standard deviation of the intensity.

pow[x_,i0_] := i0 Exp [-x^2/((2 (w0/2)^2))] ;

Code: application of effective velocities and diffusion constants to optical sorting (proper)

We read the effective constants which we calculated in Sections 1-2 into datRaw.  After reading, the variables are unpacked.   Variables which end in SI prefixes such as powAreaMega need to be multiplied by the appropriate number (e.g., mega=(10^6)) to be returned to the SI.   

The file created by solveFokkerPlanck allows for several velocities.  Each velocity is treated independently.  For now, only one velocity is allowed in makeParticleTrace.

makeParticleTrace[] :=
  Module[{ inputFileName, datRaw, flagFuchs, flagMie
         , powAreaMega, rPartNano, velMicro, velEffMicro
         , diffuseMicroSq, diffuseEffMicroSq, powAreaMegaUniq, nPowArea, rPartNanoUniq
         , nRPart, velMicroUniq, nVel, partition, powArea, rPart, vel, velEff
         , diffuse, diffuseEff, nCase, output, i0, velX, velInterp, diffuseInterp
         , tCross, diffuseX, sigmaX, sigmaT, i, serialOutput, outputFileName}
  , inputFileName = Import["makeParticleTraceInputFileName.txt"]
  ; datRaw = Import[inputFileName]
  ; { flagFuchs, flagMie, powAreaMega, rPartNano, velMicro, velEffMicro
    , diffuseMicroSq, diffuseEffMicroSq } = Transpose[datRaw]
  ; velEffMicro       = Map[ Max[#,0] &, velEffMicro,    {-1} ]     (* negative values are replaced by zero *)
  ; diffuseEffMicroSq = Map[ Max[#,0] &, diffuseEffMicroSq, {-1} ]  (* usually have a very small magnitude *)
  ; powAreaMegaUniq = Union[powAreaMega]
  ; nPowArea = Length[powAreaMegaUniq]
  ; rPartNanoUniq = Union[rPartNano]
  ; nRPart = Length[rPartNanoUniq]
  ; velMicroUniq = Union[velMicro]
  ; nVel = Length[velMicroUniq]
  ; If [ nVel != 1
    , Print["makeParticleTrace:  only set up for one velocity at a time"]
    ; Abort[]
    ]
  ; partition[expr_] := N[ Partition[expr,nPowArea] ]
  ; powArea    = powAreaMega mega // partition
  ; rPart      = rPartNano nano // partition
  ; vel        = velMicro micro // partition
  ; velEff     = velEffMicro micro // partition
  ; diffuse    = diffuseMicroSq micro^2 // partition
  ; diffuseEff = diffuseEffMicroSq micro^2 // partition
  ;  If [ powAreaPlotMin < 0
     , powAreaPlotMin = Min [ powAreaMegaUniq ] ( mega watt / meter^2 )
     ]
  ;  If [ powAreaPlotMax < 0
     , powAreaPlotMax = Max [ powAreaMegaUniq ] ( mega watt / meter^2 )
     ]
  ;  If [ deltaPowAreaPlot < 0
     , deltaPowAreaPlot = (powAreaPlotMax-powAreaPlotMin)/(nPowArea-1)
     ]
  ; nCase = Length[ powArea ]
  ; output = ConstantArray[{}, nCase]
  ; Do [
      velX = vel[[i,1]] dxdz        (* velocity in z determines velocity in x *)
    ; velInterp     = Interpolation [ pair[ powArea[[i]], velEff[[i]]     ] ]
    ; diffuseInterp = Interpolation [ pair[ powArea[[i]], diffuseEff[[i]] ] ]
    ; tCross        = (x1-x0) / vel[[i,1]]
    ; diffuseX      = diffuse[[i,1]]
    ; sigmaX        = Sqrt[ 2 diffuseX tCross ]
    ; sigmaT        = sigmaX / velX
    ; output[[i]] =                                      (* columns in particleTrace000.dat *)
        Table [ 
          { Round[ i0 / ( mega watt / meter^2 ) ]                           (* intensity (MW/(m^2)) *)
          , Round[ rPart[[i,1]] / ( nano meter ) ]                         (* particle radius (nm) *)
          , Round[ vel[[i,1]] / ( micro meter / sec ) ]                (* drift velocity (Subscript[v, d]) (μm/s) *)
          , (( (z1-z0) -
               (1/velX) NIntegrate[ velInterp[ pow[x, i0] ], {x,x0,x1} ]        (* (Subscript[δ, z]) (μm) *)
            )) / (micro meter)
          , Sqrt[ (2/velX)                                                      (* (Subscript[σ, z]) (μm) *)
               NIntegrate[ Max[ diffuseInterp[ pow[x, i0] ], 0], {x,x0,x1} ]
            + ( NIntegrate[ velInterp[ pow[x, i0] ], {x,x0,x1} ] / (x1-x0)
                * sigmaT )^2
          ] / (micro meter)
          }
       , {i0, powAreaPlotMin, powAreaPlotMax, deltaPowAreaPlot}
       ]
    , {i,nCase}]
  ; serialOutput = StringReplace[inputFileName, Except[DigitCharacter]->""]  (* extract serial number *)
  ; outputFileName = "particleTrace" <> serialOutput <> ".dat"
  ; export[outputFileName, Flatten[output,1] ]
  ; Export["plotParticleTraceInputFileName.txt", outputFileName <> "\n"]
  ; Print["makeParticleTrace:  wrote ", outputFileName]
  ]

Code: application of effective velocities and diffusion constants to optical sorting (plot)

Needs["ErrorBarPlots`"]
plotParticleTrace[] :=
  Module[{ inputFileName, datRaw, powAreaMega, rPartNano, velMicro
         , deltaZMicro, sigmaZMicro
         , powAreaMegaUniq, nPowArea, rPartNanoUniq, nRPart, velMicroUniq, nVel
         , dZPlotData
         }
  , inputFileName = Import["plotParticleTraceInputFileName.txt"]
  ; datRaw = Import[inputFileName]
  ; { powAreaMega, rPartNano, velMicro, deltaZMicro, sigmaZMicro } = Transpose[datRaw]
  ; powAreaMegaUniq = Union[powAreaMega]
  ; nPowArea = Length[powAreaMegaUniq]
  ; rPartNanoUniq = Union[rPartNano]
  ; nRPart = Length[rPartNanoUniq]
  ; velMicroUniq = Union[velMicro]
  ; nVel = Length[velMicroUniq]
  ; If [ nVel != 1
    , Print["plotParticleTrace:  only set up for one velocity at a time"]
    ; Abort[]
    ]
  ; dZPlotData =
      Partition [ pair [ pair [ powAreaMega, deltaZMicro ], Map [ ErrorBar, sigmaZMicro ] ]
      , nPowArea ]
  ; ErrorListPlot[ dZPlotData
     , FrameLabel->{ Row[{"Central Light Intensity (MW/",Superscript["m","2"], ")"} ]
                   , Row[{"(Subscript[δ, z]) (μm)"} ]
                   }
     , FrameStyle->Thick
     , Joined->True
     , LabelStyle->{Directive[Bold,Medium],FontFamily->"Helvetica"}
     , PlotStyle->{Black, Red, Green, Blue, Orange, Cyan, Purple}
     , Epilog -> {
         inset[ "standard error", 10, 1000]
       , inset[ "210", 110, 800]
       , inset[ "260", 220, 420]
       , inset[ "270", 320, 100]
       , inset[ Row[{ Subscript["v","x"], "="
                    , Subscript["v","z"],  "=1.414 mm/s"}]
         , 1150, -0.32]
       }
     ]
  ]

The default conditions lead to the plot in Figure 4.  The labels refer to the particle radius.  Notice how the ability to sort particles improves dramatically as the radius approaches 275 nm for which the particle does not feel the optical force and passes through the system undeflected.  The standard deviation of position is nearly constant because most of the value is accumulated in the low field regions on entry to or exit from the center of the beam which are the same for all conditions.

Summary

We present code to calculate effective velocities and diffusion constants for a spherical dielectric particle in a standing-wave laser field as reported in a conventional journal.[1]  A second, smaller code (included in the notebook) uses these effective quantities to predict the motion of dielectric spheres across an optical field which varies on a length scale large compared to the interference fringes.

Conclusion

The code presented here was used to make predictions about the ability to sort dielectric spheres in a system consisting of a standing wave laser interference pattern.[1]  By presenting this code, anyone can verify our results.  Readers can vary the parameters, and, with somewhat more work, change the system geometry to explore related optical landscapes.  While the notebook has been tested for the results given in Ref. [1] readers who wish to calculate new situations should be alert to numerical issues.

Appendix: Additional constants and functions

Constants

These constants should not be changed by the reader.

mega  = 1*^+6 ;         (* some International System of Units (SI) prefixes *)
centi = 1*^-2 ;
milli = 1*^-3 ;
micro = 1*^-6 ;
nano  = 1*^-9 ;
               (* SI units are like comments since they are all unity.
                  However, they force conversion to floating point. *)
sec   = 1. ;              (* time unit is 1 second *)
meter = 1. ;              (* length unit is 1 meter *)
kg    = 1. ;              (* mass unit is 1 kg *)

newton = kg meter / sec^2 ;  (* SI derived units *)
joule = newton meter ;
watt = joule / sec ;
                             (* fundamental constants *)
kBoltz = 1.38064852*^-23 joule / kelvin ;    (* 2014 CODATA value for Boltzmann constant NIST SP 959 *)
cLight = 299792458 meter / sec ;             (* exact value for speed of light *)

Mobility and diffusion formulas of Stokes and Fuchs

The diffusion constant and the mobility are related by the Einstein-Smoluchowski equation

D = μ (Subscript[k, B])T

where (Subscript[k, B]) is the Boltzmann constant and T is the temperature.   In general, the mobility is given by μ=(Subscript[v, d])/F.   Stokes’ formula F=6π ν a (Subscript[v, d]), relating the force and velocity of a sphere of radius a moving in a fluid with viscocity ν,  in effect defines the mobility.  The Stokes mobility  is given next, along with the corresponding diffusion constant from the Einstein-Smoluchowski relation, Eq. (11).

mobilityStokes[rPart_] :=  (6 Pi nu rPart)^-1             (* general for "large" spheres *)
diffuseStokes[rPart_] := kT mobilityStokes[rPart]         (* from Einstein-Smoluchowski relation *)

For small particles in air this relationship breaks down.[3]  We use a different relation, shown below.  These are plotted in Ref. [1].  The diffusion constants are from Ref. [3].  The interpolation is made in the logarithm of both the independent variable, the particle radius, and the dependent variable, the diffusion constant.

diffuseFuchs[rPart_] := 
  Module[{radius, dCM2secInv, diffuse, logDiffuse, logRadius, diffuseFuchsInterp}
  , radius =
      { 1., 2., 5., 10., 20., 50., 100.,
        200., 500., 1000., 2000., 5000., 10000.}  nano meter
  ; dCM2secInv =   (* constants from Fuchs[3]*)
      { 1.28*^-2, 3.23*^-3, 5.24*^-4, 1.35*^-4, 3.59*^-5, 6.82*^-6, 2.21*^-6
      , 8.32*^-7, 2.74*^-7, 1.27*^-7, 6.10*^-8, 2.38*^-8, 1.38*^-8
      }
  ; diffuse = dCM2secInv (centi meter)^2 / sec   (* diffusion constant in SI units *)

  ; logDiffuse = Log[diffuse]
  ; logRadius  = Log[radius]
  ; diffuseFuchsInterp = Interpolation[ pair[logRadius,logDiffuse] ]
  ; Exp[ diffuseFuchsInterp[ Log[rPart] ] ]
  ]

mobilityFuchs[rPart_] := diffuseFuchs[rPart] / kT        (* from Einstein-Smoluchowski relation *)

The reader can choose to make a plot next by setting doPlot[[iFuchsPlot]] to 1 in the code subsection of Section 1.  It is also given in Figure 5. The plots show that the Fuchs diffusion constant approaches the Stokes value for large radius.  Also, the Einstein-Smoluchowski relation of Eq. (11) makes these constants proportional to each other.

plotFuchs[] :=
  Module[{pFuchsPlotA, pFuchsPlotB}
  , pFuchsPlotA = 
      LogLogPlot [ {diffuseFuchs[a], diffuseStokes[a]}
      , {a,1.*^-9,1.*^-5}
      , FrameLabel->{ "sphere radius a (m)"
                    , Row[{ "diffusion constant D ("
                           , Superscript["m","2"]
                           , Superscript["s","-1"]
                           , ")"
                      }]
                    }
      , ImageSize->125 mmToPoint
      , PlotRange->{{1*^-9,1*^-5},{1*^-12,1*^-6}}
      , PlotStyle->{Green,Blue}
      , Epilog->{
          insetLogLog["Fuchs",6*^-9,1*^-7]
        , insetLogLog["Stokes",1*^-8,2*^-10]
        }
      ]
  ; pFuchsPlotB =
      LogLogPlot [ {mobilityFuchs[a], mobilityStokes[a]}
      , {a,1.*^-9,1.*^-5}
      , FrameLabel->{ "sphere radius a (m)"
                    , Row[{ "mobility μ (m "
                          , Superscript["s","-1"]
                          , Superscript["N","-1"]
                          , ")"
                      }]
                    }
      , ImageSize->125 mmToPoint
      , PlotRange->{{1*^-9,1*^-5},{1*^8,1*^14}}
      , PlotStyle->{Green,Blue}
      , Epilog->{
          insetLogLog["Fuchs",8*^-9,1*^13]
        , insetLogLog["Stokes",1*^-8,5*^10]
        }
      ]
  ; GraphicsRow[{pFuchsPlotA,pFuchsPlotB}]
]

Optical forces

We just give code for the optical forces which are explained in Ref. [2]; see also Ref. [4].

jBox[el_,x_] := Module[{y}, D [ y SphericalBesselJ[ el, y ],  y ] /. {y->x} ]
hBox[el_,x_] := Module[{y}, D [ y SphericalHankelH1[ el, y ], y ] /. {y->x} ]
jj[el_,x_] := SphericalBesselJ[ el, x ]
hh[el_,x_] := SphericalHankelH1[ el, x ]
       (* Mie a and b coefficients from Zangwill [4], Eq. 21.43 *)
   (* el is angular momentum index *)
   (* indexOfRefr is index or refraction *)
   (* ka is wave vector of light times radius of particle *)
abCoeff[el_, indexOfRefr_, ka_] :=
   Module[{nka}
   , nka = indexOfRefr ka ;
   ; aCoeff = - ( jj[el,nka] jBox[el,ka]  -               jj[el,ka]  jBox[el,nka] ) /
                ( jj[el,nka] hBox[el,ka]  -               hh[el,ka]  jBox[el,nka] )
   ; bCoeff = - ( jj[el,ka]  jBox[el,nka] - indexOfRefr^2 jj[el,nka] jBox[el,ka] ) /
                ( hh[el,ka]  jBox[el,nka] - indexOfRefr^2 jj[el,nka] hBox[el,ka] )
   ; {aCoeff,bCoeff}
   ]

Next, we provide the function that actually calculates the force on the particle of radius rPart and index of refraction indexOfRefr in an optical field with wavevector kWave.  The validity limits are discussed in Ref. [1].  The long, final expression was computed using the code in Ref. [2].  The force is scaled to an intensity of 1 W/(m^2).

getForceScaledMie[rPart_, kWave_, indexOfRefr_] :=
  Module[{elMax, aTab, bTab, ka}
  , elMax = 5 ;
  ; {aTab,bTab} = abCoeff[ Range[elMax], indexOfRefr, kWave*rPart ]     (* calculate Mie coefficients *)
  ; a1[el_] := Re[ aTab[[el]] ]                                (* functions in last line are table lookup *)
  ; a2[el_] := Im[ aTab[[el]] ]
  ; b1[el_] := Re[ bTab[[el]] ]
  ; b2[el_] := Im[ bTab[[el]] ]
  ; (4*Pi*(1470*a2[3] + 1120*a1[2]*a2[3] + 1575*a1[4]*a2[3] -           (* use Mie coefficients to find the *)
     1890*a2[4] - 1575*a1[3]*a2[4] - 2016*a1[5]*a2[4] + 2310*a2[5] +    (* optical force on a dielectric sphere *)
     2016*a1[4]*a2[5] +                                                 (* Ref. [1], Table 1 *)
     630*a2[1]*(1 + a1[2] + b1[1]) - 70*a2[2]*(15 + 9*a1[1] + 16*a1[3] +
     5*b1[2]) + 245*a2[3]*b1[3] - 189*a2[4]*b1[4] + 154*a2[5]*b1[5] -
     630*b2[1] - 630*a1[1]*b2[1] - 630*b1[2]*b2[1] +
     1050*b2[2] + 350*a1[2]*b2[2] + 630*b1[1]*b2[2] + 1120*b1[3]*b2[2] -
     1470*b2[3] - 245*a1[3]*b2[3] - 1120*b1[2]*b2[3] - 1575*b1[4]*b2[3] +
     1890*b2[4] + 189*a1[4]*b2[4] + 1575*b1[3]*b2[4] + 2016*b1[5]*b2[4] -
     2310*b2[5] - 154*a1[5]*b2[5] - 2016*b1[4]*b2[5]))/(105*cLight*kWave^2)
  ]

For particles of small radius, the Rayleigh limit may be used instead.

getForceScaledRayleigh[rPart_, kWave_, indexOfRefr_] :=                 (* Ref. [1], Eq. 8 *)
  (8 Pi kWave rPart^3 / cLight) (indexOfRefr^2-1)/(indexOfRefr^2+2)

Files

The three functions getNextFileName create a file name that has the base name fileBase and the extension fileExt (default “dat”).  There is also a serial number with a default of 3 digits.  The program passes “velDiffuse”, so it returns the first file name of “velDiffuse000.dat”, “velDiffuse001.dat”, ... which does not already exist.   The purpose is to avoid overwriting potentially valuable files if the code is rerun.  However, if all file names are taken, the last one, e.g., “velDiffuse999.dat” is returned.   The code below is also used to make a second output file “particleTrace000.dat”, or “particleTrace001.dat”, etc.

getNextFileName[ fileBase_ ] := getNextFileName[ fileBase, "dat" ]
getNextFileName[ fileBase_, fileExt_] := getNextFileName[ fileBase, fileExt, 3 ]
getNextFileName[ fileBase_, fileExt_, nDigit_] :=
  Module[{nTot}
  , nTot = 10^nDigit
  ; Do[
      nextFileName = fileBase <> IntegerString[ iSerial, 10, nDigit ] <> "." <> fileExt
    ; If [ Not[ TrueQ[ FileExistsQ[ nextFileName ] ] ]
      , Break[]                          (* The loop stops if a file name is not found on the disk. *)
      ]
    , {iSerial,0,nTot-1}]
  ; nextFileName
  ]

The program directs communication between sections by writing “makeParticleTraceInputFileName.txt” which tells makeParticleTrace which file to read, and “plotParticleTraceInputFileName.txt” which tells plotParticleTrace which file to read.  These files are written by default, but the reader can change the names to work with a data set written on an earlier run of part of the notebook.  Setting the appropriate elements in doJob at the top of the notebook to False is necessary for the reader-given names to take effect.

We wrap Export because when it writes data files, it does not produce a terminating newline character; export does.

export[ fileName_, expr_ ] := Export[ fileName, Append[ expr, {} ], "Table" ]

Plots

mmToPoint = 72/25.4 ;                  (* allows plot sizes in millimeters rather than points *)
inset[string_, posX_, posY_] :=                                          (* wrapper for Inset *)
  Inset[ Style [ string, Directive[Bold,Medium],FontFamily->"Helvetica"]
  , {posX, posY}, Left
  ]
insetLogLog[string_, posX_, posY_] := inset[ string, Log[posX], Log[posY] ]  (* inset for log-log plots *)

The function pair is best illustrated by example: if x={x1,x2,x3}, y={y1,y2,y3}, then pair[x,y] is {{x1,y1},{x2,y2},{x3,y3}}.   In addition to being useful for preparing input to functions in the  Plot family, it is also used for Interpolate.

pair[x_,y_] := Transpose[{x,y}]

Tips

To restart cleanly, select Evaluation>”Quit Kernel”, then Cell>”Delete All Output”.  In addition to Control/Enter, the numeric Enter key also launches an evaluation.

Acknowledgments

The mention of commercial products does not imply endorsement by NIST.

References

[1] J. J. Curry and Z. H. Levine, “Continuous feed sorting of optical particles,” Optics Express 24, 14100-14123 (2016).
DOI: 10.1364/OE.24.014100

[2] Z. H. Levine and J. J. Curry, “Scattering and Gradient Forces from the Electromagnetic Stress Tensor Acting on a Dielectric Sphere,” The Mathematica Journal, in press.

[3] N. A. Fuchs, “The Mechanics of Aerosols”, Dover, 1989, p. xiv.

[4] A. Zangwill, Modern Electrodynamics, Cambridge:Cambridge, 2013, p. 788.