Transport in randomly-fluctuating spatially-periodic potentials
The motion of overdamped particles in a one-dimensional spatially-periodic potential is considered. The potential is also randomly-fluctuating in time, due to multiplicative colored noise terms, and has a deterministic tilt. Numerical aimulations show two distinct parameter regimes, corresponding to free-running near-deterministic particles, and particles which are trapped in local minima of the potential with intermittent escape flights. Perturbation and asymptotic methods are developed to understand the drift velocity and diffusion coefficient in each parameter regime.
The subject of this paper is the statistical characterization of the motion of particles whose (one-dimensional) position is governed by the Langevin equation
where and are random functions of time. The positive constant is the deterministic particle velocity, i.e., the speed at which particles would travel in the absence of the stochastic terms, and is a spatial frequency which governs how the noise sources affects the particle at position . This equation may be viewed as describing overdamped one-dimensional motion in a spatially-periodic potential which fluctuates randomly in time, with an additional tilt. Denoting the potential as , the equation of motion is
with the potential given by
Here is interpreted as the constant tilt of the potential.
Equations of motion of this form arise in a number of diverse applications. In microelectronic circuit design applications, for example, it is important to accurately predict the effect of external noise sources upon the phase angle of nonlinear self-sustaining oscillators [1, 2]. Under some assumptions, the equation of motion for has been shown [3, 4] to be of the form (1), with the tilt representing the designed oscillator frequency in the absence of noise, and noise sources (voltage or current sources) being denoted by and . The noise sources cause a diffusive drift of the phase of the oscillator, known as phase noise. Noise of a finite bandwidth (colored noise) may also induce a shift in the average oscillator frequency . It has recently been demonstrated  that a stochastic perturbation method may be applied to accurately predict the effects of the noise on the oscillator dynamics in the limit where is much larger than the noise intensity, this being the limit of interest for engineering applications.
Further examples of models of the form (1) occur in the study of weakly coupled phase oscillators , where the noisy terms and describe fluctuations in the mean field of the oscillators , and in an agent-based model of stock market price movements . For all these examples, information on the qualitative behavior of the first and second moments of the process described by equation (1) is very desirable. In this paper we therefore examine the mean velocity and diffusion coefficient for an ensemble of trajectories in all regimes of parameter space, using extensive numerical simulation, perturbation approaches, and appropriate non-perturbative methods.
The random functions and appearing in equation (3) are independent, stationary, zero-mean Gaussian processes, which are characterized by their variance and correlation function :
(Angle brackets are used throughout the paper to denote ensemble averages). We briefly consider the case where the stochastic terms and are white noise sources (so that is a delta function), but most of our analysis is focused on coloured noise, where the presence of a non-zero correlation timescale for the noise has several nontrivial effects. Note that we do not add an independent white noise term to the right hand side of equation (2), as used to model thermal diffusion in  for example. Neglecting such diffusion effects at this stage allow us to clearly demonstrate the role of the fluctuations in the potential; moreover, it is expected that low levels of additive white noise will not qualitatively affect the results determined here.
Figure 1 provides an illustration of the effects of various parameters in the problem. These parameters are discussed in detail in section 3 below, but the basic features are clear from the three trajectories plotted in Figure 1. If the potential tilt is sufficiently large to dominate the fluctuating part of the potential, then the particle moves essentially according to , with small corrections due to the influence of the fluctuating potential. The solid red line illustrates a realization of such a case, which we class as being in the running regime. On the other hand, if the tilt is relatively weak compared to the fluctuating-in-time periodic-in-space potential, the motion of the particle is dominated by the latter. Indeed, for sufficiently slow fluctuations, we expect the particle to become trapped in a local minimum of the potential, with infrequent escapes to other nearby energy wells. The trajectory plotted as the black dashed line illustrates an example of this so-called trapping regime. Note that the contours and shading represent the gradient of the fluctuating part of the potential (shading and dotted contours denote negative values, while solid contours and white background indicate positive values), and the trajectory closely follows the position of a local minimum (zero gradient contour). We defer a full explanation until appropriate nondimensional parameters are introduced in section 3, but note that our main result is the theoretical analysis of motion in the two asymptotic parameter regimes—running and trapping—introduced here, and depicted schematically on the parameter plane of Figure 2.
The remainder of this paper is structured as follows. In section 2 we examine the white-noise limit for the random processes and in equation (3). More interesting results emerge when the noise processes have a non-zero correlation time ; the analysis of such colored noise cases begins with the introduction of nondimensional parameters in section 3. Methods of numerical simulation are detailed in section 4 and theoretical analysis of the two main parameter regimes is undertaken in sections 5 and 6. Section 7 concludes with a discussion of results.
2 White noise
If the noise sources and in equation (1) are white (i.e., in equation (4)), then the probability density function giving the particle concentration at position at time is the solution of a Fokker-Planck equation . Some care is needed when interpreting this white-noise limit, as the noise sources appear in equation (1) in a multiplicative fashion. It is typical in such problems that the Itô and Stratonovich interpretations of the white noise terms in equation (3) lead to different results. Rather remarkably however, in this particular example the same Fokker-Planck equation is found using both interpretations, as the symmetry of the two noise sources causes the so-called spurious drift term to vanish. Thus either the Itô or Stratonovich interpretation of the white noise limit of equation (1) leads to the following equation for :
It is straightforward to determine from this equation that the mean value of grows at a constant rate ; indeed, multiplying equation (5) by and integrating over space yields the exact relation
so that the mean velocity is just . Similarly, the diffusion coefficient of equals the coefficient of the final term of equation (5), i.e. . Note this coefficient shows no dependence on the value of the potential tilt parameter . As we show below, qualitatively very different results (including diffusion coefficients which depend on ) are found when the noise sources have a non-zero correlation timescale. We therefore proceed to consider the case of colored noise in the next section.
3 Colored noise and nondimensional parameters
A physically more realistic model of random fluctuations than the infinitely-fast variation of white noise is provided by colored noise, for which the correlation function defined in equation (4) decays from its peak value to zero over some characteristic timescale . One example of such a correlation function is
This correlation function is particularly convenient for numerical simulation, and also has the property which is important for some of our theoretical analysis below. This smoothness of at zero argument is related to the finiteness of the variance of the processes and , but it is not always guaranteed: if and are Ornstein-Uhlenbeck processes, for example, then this smoothness condition is violated. Nevertheless, for the purposes of this initial examination of the properties of the equation (2) we assume the smoothness property when needed, and use the correlation function (7) in all numerical simulations.
Having determined a characteristic timescale , it proves convenient to introduce dimensionless variables (temporarily denoted by tildes) by measuring time in units of and space in units of . The variance of the noise sources has dimensions of (length/time), while the tilt has dimensions length/time. The dimensionless equation of motion for may therefore be written as
where the dimensionless tilt is
and the Kubo number  is defined as
Note the dimensionless noise sources and in equation (8) have unit variance and unit (dimensionless) correlation time. The physical meaning of the Kubo number can be understood best in the case of zero tilt, i.e., with in equation (1). Then particles moving in the fluctuating potential will travel a (dimensional) distance of order during the correlation time of the noise processes; the Kubo number is the ratio of this lengthscale to the reciprocal spatial frequency of the potential. Thus the limit, for example, corresponds to a potential which fluctuates rapidly in time, so that a particle moving in it travels a very short distance (compared to the period of the potential) during the time in which values of and change significantly. On the other hand, taking corresponds to a frozen (quenched) limit, where particles feel the spatial variation of the potential much sooner than the temporal variation.
For the remainder of the paper we concentrate on the dimensionless form of the equation (8) and henceforth drop the tilde notation. We focus particularly on the range of behaviors given by varying the two positive parameters and . Results will be presented for the large-time mean velocity and the diffusion constant of an ensemble of particles. The mean velocity is defined as
or equivalently by the large-time asymptotic form of the mean position:
Assuming diffusive behavior over sufficiently long timescales, the diffusion constant is defined by the large-time asymptotic relation
where is the variance of the distribution of particle positions at time :
Alternatively, we may use the limit definition
4 Numerical simulations
Gaussian random functions may be constructed using a combination of a large number of Fourier modes, as follows:
The amplitudes and are random numbers from a Gaussian distribution of mean zero and unit variance. The random frequencies are from a distribution shaped as the Fourier transform of — for the correlation function (7) this means the are also Gaussian distributed, with mean zero and (nondimensional) variance one. The function constructed in this way is Gaussian in the limit [11, 12]; we use in the numerical simulations reported here.
In each realization, the independent noise functions and are generated as above. The ordinary differential equation (8) is then solved numerically, with initial condition , to determine . By averaging over a large ensemble (typically ) of realizations, we determine the time-dependent moments and , and the variance . A linear fit to the average position and variance at large (we use nondimensional end-times of to for this fitting) gives, respectively, the mean velocity and the diffusion constant from the slope of the fitted lines, in accord with the definitions (12) and (13).
Results from these numerical simulation are plotted as symbols in Figures 3 and 4. Note the same symbol types are used as in the schematic parameter plane of Figure 2, with the constant- results shown in Figure 3, and the constant- results in Figure 4. We choose to plot all results as functions of in order to highlight several important scalings. Panel (a) of each figure shows the normalized mean velocity , while the diffusion coefficient appears in the panel (b). Note all scales are logarithmic. It appears that two distinct regimes of simple scalings arise: broadly speaking these regimes appear to correspond to being respectively much smaller than, and much larger than, unity. However, as we show in sections 5 and 6 below the detailed picture is somewhat more complicated, with the border between the running regime (section 5) and the trapping regime (section 6) having the shape shown in the schematic Figure 2. Analytical results are derived below for limiting behaviour in each regime and these are shown as the various curves in Figures 3 and 4.
5 Running regime and perturbation method
When the mean velocity is sufficiently large, or the noise parameter is sufficiently small, particles moving in the fluctuating potential according to equation (8) may be treated as slightly perturbed from their deterministic paths. In this parameter regime (to be defined in detail below) particles typically move with a speed close to and move swiftly over the potential landscape (i.e., covering many spatial periods during the fluctuation timescale of the potential); we therefore call this the running regime of parameters. This regime is of particular interest to designers of oscillator circuits and some of the results of this section have been previously derived (using different methods) in .
To study this regime more carefully, it is helpful to move with a reference frame traveling at the deterministic mean velocity . By introducing the new variable defined as
equation (8) may be rewritten as
We combine the noise sources and time-oscillatory terms to create two new stationary Gaussian random processes and , defined by
From these definitions, it is straightforward to verify that and each have mean zero, and correlation function
However, unlike and , the new source terms and are not independent, since their cross-correlation is given by
In terms of the new noise processes, the equation of motion (18) is
Next, we assume a regular perturbation expansion for in powers of the parameter . The results of such an approach will be useful if the distance traveled by the particle in the characteristic time for the potential to decorrelate is significantly less than the spatial period of the potential. Provided this condition holds, the errors in approximating and by their series expansions about do not accumulate over time. In terms of the dimensionless parameters introduced in section 3, the spatial period of the potential is of order one. The decorrelation time of the potential felt by the particle in the moving reference frame depends on ; from (20) we can estimate this as . The distance traveled by the particle during this time is of order and the condition for the perturbation expansion to hold is therefore expressible as , or equivalently as
We adopt this condition as defining the running regime of parameter values.
We proceed to assume a regular perturbation expansion for for the limit of the form
where each is a stochastic process. Substituting this into equation (22), expanding the trigonometric terms, and gathering powers of yields the following system of equations for the unknown functions , , and :
5.1 Mean velocity in running regime
The solution of (25) with initial condition is
and substitution into (26) yields
where the substitution has been used. This result implies that in the long-time limit , we obtain a finite mean value for at leading order , given by
The mean velocity defined in equation (11) is therefore given to leading order (as ) by
5.1.1 Large- asymptotics of mean velocity
To investigate the asymptotic behavior of (32) at large we integrate by parts:
For sufficiently smooth correlation functions , this integration by parts may be repeated in an iterative fashion, to obtain an asymptotic series is inverse powers of . To leading order in , we obtain
and we conclude that deviations of the mean velocity from are of order in this regime of parameter values.
5.2 Diffusion coefficient in running regime
The variance of (which equals the variance of ) is given to leading order by
and from (28) we obtain
The slope of is given by its derivative with respect to :
and as this reaches the finite limit giving from equation (15):
This leading-order result for is plotted with solid lines in panel (b) of Figures 3 and 4. While it matches numerical results well for low and low cases, it is noticeable that even well within the running regime (with ) it is very inaccurate (not even appearing on the scale of the Figures) for (Fig. 3(b)) and for (Fig. 4(b)). The reason for this loss of accuracy is examined in detail in the following section.
5.2.1 Large- asymptotics of diffusion coefficient
We consider the asymptotic behavior of the diffusion constant at large values of . For correlation functions which are sufficiently smooth at , the order integral in equation (38) can decrease rapidly with increasing . For example, using the correlation function of equation (7) in equation (38) leads to the prediction that the diffusion constant is exponentially small at large values:
Because the first non-vanishing term in the perturbation expansion for may be exponentially small, we proceed to calculate to higher orders. The contribution to at order is given by
and after extensive algebraic manipulation this can be shown to yield a contribution to the diffusion constant of
The large- asymptotic behavior of this term may be calculated (in a similar fashion to the mean velocity case above) as
For the correlation function of equation (7) we therefore have the following small-, large- asymptotic behavior:
This corrected form for is plotted with dashed lines in Fig. 3(b) and Fig. 4(b) and shows excellent agreement with numerics within the running regime. We conclude that the deviation from the leading-order result (38) noted above is due to the order- effects dominating the order- term when is sufficiently large.
5.3 Summary of running regime
In brief, condition (23) defines the running regime of parameter values and within this regime the mean velocity and diffusion coefficient may be determined using a stochastic perturbation expansion. To order the mean velocity is given in equation (32) and the diffusion coefficient in equation (38). Asymptotic analysis for large reveals that the velocity scales as as , while can have an order term which is exponentially small, see equation (39). In this case the order- contribution of equation (42) will dominate.
6 Trapping regime
In section 5 we analyzed particle motion under the condition that , which corresponds to the running regime. At the other extreme, i.e. when , the dynamics of particles may be understood in terms of trapping at low-velocity positions, with intermittent releases and flights between trapping points. To analyze the dynamics in this trapping regime, it is useful to rewrite equation (8) as
The phase and the amplitude are each random functions of time, and are given in terms of the original noise processes and as
conversely, the noise sources may be expressed as and .
6.1 Diffusion coefficient in trapping regime
Taking the limits and in equation (45) allows us to derive some useful results. In this limit, the particle trajectories are locked to the phase variable: (up to the addition of multiples of ), and in particular the velocity of the particle is given explicitly as a function of time by . The phase velocity in systems of this type has been extensively studied [14, 6], and known results may be applied to find the diffusion coefficient in this limit . The mean value of is zero, and the phase velocity correlation function may be expressed in terms of the correlation function of the noise sources as
In this trapping limit the effective diffusion coefficient for the particles may therefore be approximated by setting and proceeding as follows:
Given the noise correlation function, this integral may be calculated numerically after applying equation (47); for the case this yields the value . This value is plotted with a dot-dashed line in Figures 3(b) and 4(b) and agrees with the numerical results for parameters well within the trapping regime .
6.2 Mean velocity in trapping regime
In the strong trapping limit of the mean velocity of the particles is necessarily zero. Numerical results (e.g. Fig. 4(a)) indicate that if is non-zero but small then the mean velocity scales linearly with , and we are motivated to study the behavior of equation (45) as a small- perturbation about the strong-trapping case considered above. As we shall see, it is a challenging problem to describe the dynamics in this case, because of short (in time) but large (in space) intermittent escapes of the particle from the strong-trapping solution. As a consequence, we will revert to a simpler toy problem in an attempt to understand the observed numerical scaling of the mean velocity in this regime. We begin by giving the exact equation for the deviation of the particle position from its strong-trapping limit :
Analysis of this equation is significantly complicated by the dynamics of the non-Gaussian phase velocity . For example, noting that the phase velocity may be expressed as
and that and are independent Gaussian processes, we see that the one-time distribution of is the same as that of , where is a zero-mean and unit variance Gaussian random function and is the amplitude given by equation (46) above. This fact may be used to find the distribution of in closed form , but here we use it to develop an analytically-tractable toy model which demonstrates similar behavior to the full dynamics.
We consider splitting time into intervals of length , and in each interval replacing the time-varying functions and with random constants chosen from the appropriate one-time distributions. This frozen-noise approximation thus neglects all variation of the phase and amplitude variables on timescales shorter than , but it allows the dynamics of the corresponding to be found be solving the constant-coefficient equation
within each -interval. Here is given by and is approximated by , with , , and being independent Gaussian random numbers with mean zero and unit variance, chosen anew for each new time interval.
We focus on the mean velocity , which is equal to since has zero mean. The ensemble-averaged -velocity within the toy model is calculated by solving (51) within each -interval, and then averaging over the possible values of and . There are essentially two cases to consider. If the random numbers and satisfy the inequality
then a steady state solution of (51) exists for , and so after a transient the particle reaches this steady state and the velocity decays to zero. For sufficiently long time intervals , we can neglect the transient and so consider the effective velocity to be zero in cases satisfying (52). If the inequality is not satisfied, however, then no steady state solution of equation (51) exists, and these running solutions have a time-averaged velocity calculated (see the appendix) to be
To determine the ensemble mean velocity it remains only to average these running velocities over the distribution of and values satisfying the inequality (52). After some manipulation (see appendix) it is possible to show that the resulting mean velocity is linear in , taking the form
Using this form for the mean velocity gives (after normalizing by ) the dot-dashed curves in Figures 3(a) and 4(a). It is apparent from Fig. 4(a) that the linear dependence on in (54) qualitatively matches the numerical results, but the prefactor is too small by a factor of approximately 2. Given the approximations made in the simplification to the toy problem used here, this level of agreement is considered rather satisfactory.
To summarize: we have examined the behavior of first and second moments of the solution of the equation (1) (nondimensionalized to equation (8)) with independent random forcing terms and . Numerical simulation results are supported by asymptotic analysis in the two main regimes of the – parameter plane, termed the running regime (section 5) and the trapping regime (section 6)—see also Figure 2 and definition (23).
The smoothness of the noise correlation function at was assumed at the beginning of section 3. If the smoothness property is violated by, for example, using Ornstein-Uhlenbeck processes for and , then the analysis of the trapping regime which relies on and being finite-variance processes is not valid. The perturbation analysis in the running regime is unchanged, but note the leading-order diffusion term (38) typically scales as for small and large if is non-smooth and so the order- term (43) will not dominate as it does when is smooth.
It is interesting to compare the results for our model with related work on noisy transport in periodic potentials. The non-monotone dependence of the diffusion coefficient on , seen in Fig. 3(b) for and in Fig. 4(b) for , is reminiscent of Figure 1 of . The authors of  consider overdamped motion in a tilted, spatially-periodic potential which (unlike our randomly-fluctuating case) is static in time, but to which white noise effects (not included by us) are added. They find the diffusion coefficient in their model is greatly increased at a critical value of the tilt parameter corresponding to the onset of deterministically running solutions. We find the peak value of in our model occurs for ratios of order unity (actually approximately 0.5 to 0.7), which corresponds to the border between the running and trapping regimes at high values of or . Because particles may be either running or trapped in any given realization, the overall variance of across the ensemble is particularly large for , leading to a diffusion coefficient which exceeds that found in cases where each single realization is more typical of the ensemble average.
Our perturbation expansion method for the running regime bears some resemblance to that used in  to examine Stokes drift in a fluid due to the presence of periodic waves. Substitution of the Fourier expansions (16) for and in equation (8) allows the equation of motion to be written as
with the sum admitting the interpretation as a superposition of waves of various amplitudes (and directions) and random frequencies . Our result (32) may therefore be compared (by taking ) to the drift velocity determined by a perturbation exapansion in . Unlike , however, the randomness of our waves means there is no preferred wave direction, and so we find a zero drift velocity (we also do not consider the additive white noise diffusion effects used in  and ). The directional symmetry of the random wave field is broken when is non-zero and in this case the deviation of from the deterministic velocity is of order (from equation (32)), which (like in ) is the square of the wave amplitude. It would be of considerable interest to extend our methods to include the effects on the particle trajectories of additive white noise terms as in , , , in addition to the multiplicative colored noise terms considered here.
Finally, we note the possibility of using spatially-periodic potentials for sorting of microparticles, which has attracted considerable recent attention , , , . If experimental periodic potentials can also be made to fluctuate randomly in time (perhaps by modulation of the laser intensities generating optical lattices as in , ) then the results of this paper would provide a useful basis for predicting mean velocities and spreads of particles being sorted within the device.
This work is funded by Science Foundation Ireland under Investigator Award 06/IN.1/I366 and Research Frontiers Programme 05/RFP/MAT0016. The author benefited from helpful discussions with Fergal O’Doherty and Grant Lythe.
Appendix: Mean velocity in toy model for trapping regime
In section 6.2 we introduced a simplified toy model for the trapping regime within which we can gain some analytical insight into the behavior of the mean velocity. Recall equation (51) generates a non-zero contribution to the mean velocity in a given -interval if the parameters and do not satisfy inequality (52)—here , and are unit variance, mean zero normal random numbers in each -interval. Indeed, in this case equation (51) immediately yields an integral solution for the time taken to travel the spatial period :
which gives the equation (53) in the main text for the running velocity in that particular -interval (assuming is sufficiently large). In all cases the sign of is the same as the sign of the quantity .
Next, we consider the averaging of (53) over the possible values of and . The values of are chosen from the range with density , while the values of are from the range with density . When averaging over all values of and satisfying the inequality
we must consider three regions of the – parameter plane, as follows:
Region 1: to , with to . In this region is positive.
Region 2: to , with to . In this region is positive.
Region 3: to , with to . In this region is negative.
Because we are interested only in the leading-order behavior of at small deterministic velocity , we focus on the linear term in an expansion in small :
It is easy to show that Region 2 does not contribute to the mean velocity at this order, while Regions 1 and 3 each contribute an amount equal to
After the change of variable we obtain
where the final expression has been obtained by taking in the integrand. Adding the equal contributions from Regions 1 and 3 then gives the asymptotic expression of equation (54) of the text for the averaged velocity.
-  A. Demir, A. Mehrotra, and J. Roychowdhury, “Phase noise in oscillators: a unifying theory and numerical methods for characterization,” IEEE Trans. Circuits and Systems I, 47, 655-674 (2000).
-  A. Demir, “Phase noise and timing jitter in oscillators with colored-noise sources,” IEEE Trans. Circuits and Systems I, 49, 1782-1791 (2002).
-  J.P. Gleeson,“Phase diffusion due to low-frequency colored noise,” IEEE Trans. Circuits and Systems II, 53, 183-186 (2006).
-  F. O’Doherty and J.P. Gleeson, “Phase diffusion coefficient for oscillators perturbed by colored noise,”IEEE Trans. Circuits and Systems II, 54, 435-439 (2007).
-  Y. Kuramoto, Chemical oscillations, waves, and turbulence, Springer, Berlin (1984).
-  J.P. Gleeson, “The mean field of weakly coupled oscillators exhibits non-smooth phase noise,” Europhys. Lett., 73, 328-334 (2006).
-  J.P. Gleeson, “Passive motion in dynamical disorder as a model for stock market prices,” Physica A, 351, 523-550 (2005).
-  P. Reimann, C. Van den Broeck, H. Linke, P. Hänggi, J. M. Rubi, and A. Pérez-Madrid, “Giant acceleration of free diffusion by use of tilted periodic potentials,” Phys. Rev. Lett., 87, 010602 (2001).
-  H. Risken, The Fokker-Planck equation: methods of solution and applications, 2nd ed., Springer-Verlag, Berlin (1996).
-  N.G. Van Kampen, Stochastic processes in physics and chemistry, 3rd ed., North-Holland (2007).
-  P.R. Kramer, O. Kurbanmuradov, and K. Sabelfeld, “Comparative analysis of multiscale Gaussian random field simulation algorithms,” J. Comp. Phys., 226, 897-924 (2007).
-  R.H. Kraichnan, “Diffusion by a random velocity field,” Phys. Fluids, 13, 22 (1970).
-  C.M. Bender and S.A. Orszag, Advanced mathematical methods for scientists and engineers, Springer, New York (1999).
-  J.P. Gleeson, “Exactly solvable model of continuous stationary 1/f noise,” Phys. Rev. E, 72, 011106 (2005).
-  K. M. Jansons and G. D. Lythe, “Stochastic Stokes drift,” Phys. Rev. Lett., 81, 3136-3139 (1998).
-  C. Van den Broeck, “Stokes’ drift: an exact result,” Europhys. Lett., 46, 1-5 (1999).
-  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).
-  M.P. MacDonald, G. C. Spalding and K. Dholakia, “Microfluidic sorting in an optical lattice,” Nature, 426, 421 (2003).
-  A. Gopinathan and D. G. Grier, “Statistically locked-in transport through periodic potential landscapes,” Phys. Rev. Lett., 92, 130602 (2004).
-  J. P. Gleeson, J. M. Sancho, A. M. Lacasta, and K. Lindenberg, “Analytical approach to sorting in periodic and random potentials,” Phys. Rev. E, 73, 041102 (2006).