An analytical approach to sorting in periodic potentials

There has been a recent revolution in the ability to manipulate micrometer-sized objects on surfaces patterned by traps or obstacles of controllable configurations and shapes. One application of this technology is to separate particles driven across such a surface by an external force according to some particle characteristic such as size or index of refraction. The surface features cause the trajectories of particles driven across the surface to deviate from the direction of the force by an amount that depends on the particular characteristic, thus leading to sorting. While models of this behavior have provided a good understanding of these observations, the solutions have so far been primarily numerical. In this paper we provide analytic predictions for the dependence of the angle between the direction of motion and the external force on a number of model parameters for periodic as well as random surfaces. We test these predictions against exact numerical simulations.


I. INTRODUCTION
Transport of particles driven by external forces across modulated potential surfaces, and the associated ability to sort mixtures of particles according to some feature sensitive to this modulation, has attracted considerable recent interest [1,2,3,4,5,6,7,8,9,10,11,12]. The surface modulation may be periodic or random, and it may consist of traps, obstacles, or a combination of both. Typical sorting parameters include particle size, particle index of refraction, and particle mass. The underlying mechanism in all cases relies on the fact that while the particles attempt to follow the external force, the traps or obstacles induce systematic deviations in the trajectories so that there is a non-zero average angle between the trajectory and the external force. If this non-zero angle depends on some characteristic of the particle such as size, then particles that differ in this characteristic emerge at different angles and can therefore be sorted in a very effective way.
There are a number of experimental demonstrations of these capabilities [1,3,4,6,7,8], and there are also fairly extensive numerical simulations of various model formulations [2,5,10,11,12] that are quite successful in explaining the experimental phenomena. These simulations to some extent clarify the roles of a number of physical variables such as, for example, the direction of the external force, surface geometry, and ambient temperature, that affect the experimental outcomes. On the other hand, analytic results are relatively rare [9], and yet it is undoubtedly useful to have predictive analytic formulas, especially with a view to optimizing the sorting mechanism. In this paper we take a step in this direction, deriving approximate results that are shown to reproduce some of the important earlier numerical results. While our results are limited to certain parameter regimes, we believe that they may prove useful in designing experimental sorting potentials and in identifying likely regions of parameter space for more detailed numerical study.
In Sec. II we describe the model and define the quantities to be calculated, specifically, the angle between the average velocity of the particles and the external force. In our discussions we include periodic potentials as well as random potentials. In Sec. III (accompanied by an Appendix) we present approximate equations for the average velocity valid for sufficiently strong forces and/or high temperatures, and discuss how these approximations might be continued to further orders than those retained here. In Sec. IV we compare our theoretical results with those of numerical simulations for periodic potentials, and in the course of this comparison we refine our approximations so as to eliminate an unphysical outcome. With this adjustment we improve the theoretical predictions. In Sec. V we carry out a similar comparison for random potentials, and again find very good agreement between theory and numerical simulations. We summarize our results and indicate some future directions in Sec. VI.

II. THE MODEL
We consider identical non-interacting particles moving across a surface described by a two-dimensional potential V (x, y) of unit height or depth and unit period in both directions. The equation of motion for the particles is given by a Langevin equation for each component of the particle displacement [10], The dots denote time derivatives, T is the dimensionless temperature, and the thermal fluctuation terms ξ i (t) are Gaussian and δ-correlated, Note that the equations are in dimensionless form, with the scaled temperature T and force magnitude F given in terms of physical parameters as in Eq. (4) of [10].
The system is assumed to be overdamped so that inertial termsẍ andÿ have been dropped from the equations of motion. The constant external force vector is and all particles are assumed to start at the origin. Our goal is to calculate the direction and magnitude of the velocity of the particles averaged over the thermal fluctuations. The long-time limit of this velocity is denoted by v , and is decomposed into Cartesian components as The angle brackets denote averaging over the thermal noise. If the potential V (x, y) is random, then we also perform an average over realizations of the potential, denoting such an ensemble average by an overbar: v . To calculate the average velocity for a periodic potential, we write the Fokker-Planck equation for the concentration c(x, t) of the particles obeying the equations of motion (1), with initial condition c(x, 0) = δ(x). If this equation could be solved for the concentration c(x, t) at time t, then the desired average velocity vector follows from v = lim While an exact solution of Eq. (5) is not in general possible, an approximate solution for the concentration may be derived using a number of methods. We implement such an approximate solution in the next section, and also generalize the results to the case of a random potential by averaging over the disorder to obtain v . Of particular interest in sorting applications is the deflection angle α of the velocity from the external force direction, given by the relation [10] tan (If the potential is random, each velocity component v ⊥ , v x etc. in this formula should be replaced by its average over the disorder: v ⊥ , v x etc.). We will explore the behavior of this angle through that of each component v x and v y of the average velocity, and examine the accuracy of our approximations in capturing these behaviors.
To avoid overly complicated formulas, we adopt the convention throughout the remainder of this paper (except in the Appendix) that the notation v (or its components v x , v y , v ⊥ , v ) stands for the averaged quantity v (or, respectively, v x , v y , v ⊥ , v ). If the potential is random, v (or its components) stands for v (or its respective averaged components). Averages of all other quantities will be denoted explicitly using angle brackets and overbars as appropriate.

III. APPROXIMATE SOLUTIONS
The experimentally interesting regimes involve an external force of magnitude sufficiently larger than the well depth (F > 1) so that particles do not easily become trapped. If a particle did become trapped, it would have to be extracted by a sufficiently large thermal fluctuation to continue moving across the surface. The particles of interest in a sorting context are those that do not become trapped during the course of the experiment. On the other hand, if the force is too large, it simply drags the particles along, the features of the potential become essentially invisible, and particle separation, which relies on the effect of the potential on the particle trajectories, does not occur. At the same time, the temperature of the system should not be so high as to obliterate the features of the potential. A temperature that is too high would again lead to an uninteresting situation in that particle separation would again not be observed. The regime of interest is thus that of an external force that is large but not too large, and a temperature that is as low as possible in theory and in practice.

A. First order approximation for periodic potentials
In the Appendix we show how a systematic large force and/or high temperature approximation may be obtained. The first-order approximation leads to the average velocity vector where for periodic potentials the functionQ(k) is defined byQ Here k 2 = k.k andV is the Fourier transform of the potential V (x) as defined in Eq. (A1). For later comparisons, it is useful to exhibit the result of this approximation for the simplest one-dimensional version of our problem, one with a simple cosinusoidal potential. This case has been studied extensively (see, for example, Chapter 11 of [13]), and exact solutions may be found for the average velocity. The equation of motion is We call the forcing U instead of F because for the simplest separable potential in two dimensions we will be able to identify the average speed v(U ) = ẋ with the components of the two-dimensional velocity via the re- A number of points are noteworthy. Firstly, a simple sinusoidal potential has minima (which we think of as traps) as well as maxima (which we think of as barriers) relative to a flat landscape [11]. Secondly, we observe that the net effect of the potential in one dimension is to slow down the particles. Thirdly, when the associations v x = v(F cos θ) and v y = v(F sin θ) are appropriate, we can immediately see that the deflection α is associated with the fact that the "slowing down" is different for each component for most angles θ.
The technique yielding the approximate formula (8) in the case of a periodic potential may readily be extended to random potentials, as simulated numerically in [14], for instance. Indeed, the method of deriving Eq. (8) was originally developed for the calculation of particle concentrations when advected by potentials which are random in space [15] and also in time [16,17]. The possibility of sorting on random potentials can thus be examined using analytic approximations.

B. First order approximation for random potentials
The appropriate kernel in formula (8) for a random potential isQ whereÊ(k) is the "energy spectrum" of the potential. The energy spectrum is defined as the Fourier transform of the (disorder-averaged) correlation function of the po- The disorder is assumed to be homogeneous, so that the correlation function depends only on the difference vector x; recall that the overbar denotes averaging over the ensemble of random potentials. It is again useful to exhibit the first order approximation for the simplest one-dimensional version of the problem. The equation of motion is where V (x) is now a zero-mean random modulation with spectrumÊ(k). The one-dimensional velocity averaged over realizations of the noise ξ(t) and of the random potential is then approximately given by A more explicit result requires the specification of the correlation function.
It is an interesting general result that regardless of the specific form of the correlation function, the average deflection angle is zero when the random potential in our two-dimensional scenario is isotropic, that is, when the correlation function V (x ′ ) V (x ′ + x) depends only on the magnitude |x| of the difference vector and not on its direction (see, for example, Eq. (14) of [14]). In this case the energy spectrumÊ(k) depends only on the magnitude k = |k| of its argument. Using (12), we consider the wave-vector integral in (8) by writing k in terms of a component k = k cos φ parallel to the force F, and a component k ⊥ = k sin φ perpendicular to F, where φ is the angle between k and F. With the two-dimensional integral written in polar coordinates, the average perpendicular velocity is found to be (15) The isotropy of the potential means that the energy spectrum is independent of the angle φ, and therefore the integration over φ may be performed in (15), giving the result v ⊥ = 0. We conclude that for isotropic random potentials the deflection angle is zero. This means that sorting by deflection angle is not possible in isotropic random potentials.

C. Higher order approximations
It is useful to explore the behavior and magnitude of the second and higher order approximations for the velocity. We do not do this in full generality, but only for the simplest one-dimensional problems (10) and (13).
Following the procedure outlined in the Appendix, a second iteration is somewhat cumbersome but straightforward, and leads in the periodic case to the approximate velocity . (16) The last term is the second order correction. In a subsequent section, where we make comparisons with exact results, we comment on some higher order terms in this series. Here we simply point out that for T = 0 the full infinite series sums to the known exact solution [13] In the case of a one-dimensional random potential, continuing the iteration to the next order extends Eq. (14) to Finally, while we recognize that the approximation methods based on the truncation of a series precludes a straightforward determination of the regimes of validity, one might expect that high temperature and/or strong external forcing (T ≫ 1 and/or F ≫ 1) would be sufficient to make the approximations presented in this section reasonably accurate. As we shall see later, such an assertion requires some caveats, but numerical simulation results help resolve the issues which arise.

IV. RESULTS FOR PERIODIC POTENTIALS
In this section we test the approximate results against exact simulation results for two-dimensional periodic potentials, explain why large F is not necessarily sufficient for agreement at the orders developed in the last section, and modify our approximations nonperturbatively so as to correct for the source of disagreement.
If the potential is 1-periodic and even in both x and y, such as the one used in [10], it may be expanded in a Fourier series as where M may be infinite, although in our calculations we retain only a finite number of modes in each direction.
After some algebra, Eq. (8) yields explicit expressions for the average velocity in the x and y directions as sums over Fourier modes, where the coefficients a nm depend on the specific potential, and Note that to this order of approximation the Fourier coefficients a nm appear only as a 2 , so that the resulting velocities are the same regardless of the overall sign of the potential, e.g., whether the potential consists of traps or of obstacles.
To examine the usefulness of the approximate formulas (19), we choose the periodic potential used in [10], with g(x, y) = 5 [cos(2πx) + cos(2πy) − 2B], and with various values of the parameter B (associated with different particle sizes [10]). The negative sign in the numerator indicates that we are considering traps (but, as noted above, this sign does not matter for the assessment of the validity of the approximation). For each value of B, the potential is first expressed in terms of its Fourier series (18) to determine the Fourier coefficients a nm . Figure 1 shows the angle tan α defined in Eq. (7) plotted against tan θ for F = 8 and T = 0.1, and with B = 0.9. Note the degree of agreement between the exact and approximate results, especially the positive and negative regions of deflection angle which reflect the "terrace phenomenon" or "locked-in" states observed in experiments [1,3,4,6,7,8]. Note also the accord in the order of magnitude of the deflection. These indications point to the qualitatively satisfactory performance of the approximate formula (8)  (21) is well approximated with M = 6; however our results with M as low as 1 show that the correct order of magnitude of the deflection can be predicted using only the lowest Fourier harmonics of the potential.
The Fourier decomposition provides an opportunity to understand the role of various symmetry contributions to the sorting capability of the surface. For this purpose, we examine two special cases with M = 1. The potential V a (x, y) = cos(2πx) + cos(2πy) corresponds to having a 01 = a 10 = 1, with all other a nm being zero. This separable potential was used in [11] and recently in [18], and we discuss it further subsequently. For the parameters used in Fig. 1, the potential V a yields deflection angles which are negative for all forcing directions. In contrast, the potential corresponding to a 11 = 1, with a nm = 0 otherwise, gives positive deflection angles for all forcing directions. It is the (appropriately weighted) combination of the potentials V a and V b in the Fourier series of (21) which generates the crossover from negative to positive deflection angle in the M = 1 curve of Fig. 1.
In Fig. 2 we show results of numerical simulations of sorting in the separable potential V a at temperature T = 0.1 and for various magnitudes of the forcing F . As noted following Eq. (10), separable potentials allow the application of one-dimensional results to the two-dimensional case by setting v x = v(F cos θ), v y = v(F sin θ). The dashed curves in Fig. 2 show the deflection angle predicted using the first order approximation (11). While the agreement is quite good at the higher angles θ, it is clearly not quantitatively satisfactory at low angles. The second order approximation from Eq. (16) yields the solid curves, and further improves the agreement at the higher values of θ and pushes this agreement toward lower θ, but does not greatly improve the low-θ situation. The difficulty can easily be traced in the case of the separable potential and, by inference, for more complex potentials, by recalling that v y = v(F sin θ). Clearly, at low angles the argument of v y is necessarily small no matter the magnitude of F . Equations (11) and (17) make it clear that at sufficiently small argument the one-dimensional approximate velocity becomes negative. This unphysical result is reflected in the misbehavior of v y at small θ and leads to the disagreements observed at small angles.
To test this hypothesis, we introduce an enhanced approximation, which consists of calculating v x and v y from the first-order or second-order iteration approximation, as we have done so far, but then replacing v y by max[v y , 0] and v x by max[v x , 0]. This eliminates the unphysical negative values of the components, replacing them by zero, without affecting positive outcomes. We shall call this the adjusted truncation, and its effect on the predictions is to cut off the curves to the left of the dotted line tan α = − tan θ in Fig. 2. The agreement between the adjusted truncation and the numerical results is clearly very good , with the remaining differences arising from the abrupt replacement of negative velocities by zero velocities in the adjusted truncation (and, for the second order case, a very small region near θ = 0, where a positive deflection angle is erroneously predicted).
We end this section by repeating our earlier assertion that the perturbation theory becomes more accurate with increasing temperature. In fact, the series expansion developed in the Appendix is a perturbation series about the limit of infinite forcing, and non-zero temperature acts to regularize this series at finite forcing. The accuracy of truncations of the series thus improves with increasing force and with increasing temperature. Although the sorting capability of the system decreases with increasing temperature (and with increasing force), it is useful to display some higher temperature results explicitly. We do so in the one-dimensional case. By continuing the iteration method introduced in the Appendix, we write successive approximations to the average velocity as with .
Here we writeŨ = U/(2π) to keep the formulas simple; note that v 1 and v 2 have already appeared in Eq. (16). Figure 3 show v/U vs T at forcing value U = 3 as the number of terms retained in the series is increased. Note that at high temperatures, even the simple first order (N = 1) truncation already gives accurate results. At lower temperatures, retaining higher order terms in the approximations gives results which are closer to the exact values. As the temperature approaches zero for this U value, the approximations all fail; this effect depends on the forcing U , and at higher values of U all the approximations remain accurate even at T = 0, as noted following Eq. (16).

V. RESULTS FOR RANDOM POTENTIALS
We next test our approach for random potentials, for which the first order approximation to the velocity is given in general in Eq. (8) and Eq. (12), and in one dimension in Eq. (14). The second order iteration for a one-dimensional random potential is shown in Eq. (17). Recall that sorting is not possible in an isotropic random potential. We must therefore choose an anisotropic potential. In particular, we consider a two-dimensional (separable) potential generated by adding together independent, zero-mean modulations in the x and y directions, with correlation functions given by The energy spectrum of this potential has the form whereÊ 1 andÊ 2 are the one-dimensional spectra given by the Fourier transforms of V 1 and V 2 . Using spectrum (28) in the formula (8) yields the average velocity components We note that these formulas imply that the average velocity component perpendicular to F is non-zero in general (even ifÊ 1 =Ê 2 ): The separable potential (26) is special because it reduces the two-dimensional sorting problem to motion in uncoupled one-dimensional potentials in the x and y directions. The results (29) for the two-dimensional case follow from Eq. (14) by noting that v x = v(F cos θ) and v y = v(F sin θ).
A comparison between the first-order adjusted truncation results and numerical simulations for a twodimensional separable potential with correlation functions is presented in Fig. 4. While the agreement is not as good as in the periodic case, the theory clearly captures the numerical behavior rather well, and the effect of the order of the approximation is again manifest.

VI. CONCLUSION
In summary, we have derived approximate analytic or quadrature forms for the average velocity of particles moving in a periodic potential or random potential. The average velocity in general deviates from that of the applied constant external force, the deviation depending on some particle characteristic such as size. This can then be used to sort particles that differ in this characteristic. A systematic perturbation series valid for large external forces and/or high temperatures is shown to capture the behavior observed in experiments and numerical simulations in physically interesting parameter regimes when the angle between the external force and the crystallographic x axis (or y axis) is relatively large, e.g., when the force lies near the crystal diagonal, even when the perturbation series is truncated at low orders. The results are not nearly as good when the force lies near one of the crystal axes, because the truncated perturbation series can then lead to unphysical negative velocity components. We have proposed an adjustment to the simple perturbation expansion whereby negative velocity components are set to zero, and have shown that this adjusted truncation scheme leads to very good agreement with numerical simulation results even for low temperatures.
Further directions in this work are plentiful, and here we list just a few. A first direction would be a generalization to reduced-symmetry potentials. For example, in [12] we considered a potential of the form which has different length scales in the x and y directions. Another extension would be to the approximate calculation of the diffusion tensor (see, e.g., Eq. (4.24) of [15]), in particular to compare values for diffusion in the transverse and parallel directions to the direction of transport [12]. Reimann et al. [19] supply an exact expression for the effective diffusion in onedimensional problems, which should also be applicable to the two-dimensional case with a separable potential. These and other related investigations are in progress.

APPENDIX A: ITERATION
In this Appendix we outline a systematic approximation scheme that leads to our main results in Sec. III, starting from the Fokker-Planck equation (5). We call u(x) = −∇V (x) for convenience. The case where u is a Gaussian random field perturbing the strong external bias F has been examined in, for example, Sec. 4.2.2 of [15]. In our case the field u may be random or periodic, and we mimic the derivation of the systematic expansion of [15] to determine our results for the long-time average velocity.
Defining the (spatial) Fourier transformĉ of the concentration field c bŷ the Fokker-Planck equation is transformed to the equation ∂ĉ ∂t +ik·Fĉ+ i (2π) 2 dp k·û(p)ĉ(k−p, t)+k 2 Tĉ(k, t) = 0, (A2) with initial conditionĉ(k, 0) = 1. Note that the Fourier transform of u isû (p) = −ipV (p). (A3) As a consequence we havê which we will use below. We also perform a Laplace transform in time, defininḡ After Laplace transforming Eq. (A2) and using the initial condition, our goal becomes the solution of the integral equation where the propagator (or "free Green function" [15]) is defined as The following limiting behavior of the propagator will be important later: If the exact solutionc(k, s) of Eq. (A6) could be found, then the thermal-averaged velocity defined in Eq. (6) can be determined using standard limiting theorems for Laplace transforms as v = lim s→0 is 2 ∂c ∂k k=0 . (A9) However, as Eq. (A6) is not in general exactly solvable, we seek instead an approximate solution forc(k, s).
A standard approach to approximating the solution of an integral equation is by iterating; i.e., first settinḡ c(k, s) = P s (k) (neglecting the second term on the right hand of Eq. (A6)), then substituting this approximate solution into (A6) to obtain an updated solution, etc. After two iterations this procedure yields: c(k, s) = P s (k) − i (2π) 2 P s (k) dp k ·û(p)P s (k − p) − 1 (2π) 4 P s (k) dp k ·û(p)P s (k − p) with the dots signifying the further terms in the formal iteration series. Successive powers of the propagator appear at each iteration, and therefore the approximation is assumed to be better for large F and large T .
Consider now the velocity v in Eq. (A9) which results from using the approximation that terminates the series (A10) by neglecting the dots. First, we differentiatec with respect to the component k j and evaluate at k = 0: ∂c ∂k j k=0 = −i s 2 F j − i (2π) 2 P s (0) dpû j (p)P s (−p) − 1 (2π) 4 P s (0) dpû j (p)P s (−p) × dq (−p) ·û(q)P s (−p − q). (A11) Following (A9), we now multiply this equation by is 2 and examine the limit as s → 0, and find that the first term of (A11) yields F j , i.e. the unmodified influence of the external force on the velocity. The second term of (A11) reduces in the limit to 1 (2π) 2 dpû j (p) lim = lim s→0 s s + p 2 T − ip · F s + (p + q) 2 T − i(p + q) · F = 0, unless p = 0 or p + q = 0.
Together with the contribution from the first term of (A11), we thus have the approximation v ≈ F + i (2π) 4 dpû(p) 1 T p 2 − ip · F p ·û(−p).

(A15)
Equation (8) follows on replacingû from equation (A3); the random potential case then requires a further ensemble average over the disorder to yield equation (12).