{'fraction_non_alphanumeric': 0.029947219459955633, 'fraction_numerical': 0.01028838063183661, 'mean_word_length': 4.623010752688172, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Separation of suspended particles in microfluidic systems by directional-locking in periodic fields

We investigate the transport and separation of overdamped particles under the action of a uniform external force in a two-dimensional periodic energy landscape. Exact results are obtained for the deterministic transport in a square lattice of parabolic, repulsive centers that correspond to a piecewise-continuous linear-force model. The trajectories are periodic and commensurate with the obstacle lattice and exhibit phase-locking behavior in that the particle moves at the same average migration angle for a range of orientation of the external force. The migration angle as a function of the orientation of the external force has a Devil's staircase structure. The first transition in the migration angle was analyzed in terms of a Poincare map, showing that it corresponds to a tangent bifurcation. Numerical results show that the limiting behavior for impenetrable obstacles is equivalent to the high Peclet number limit in the case of transport of particles in a periodic pattern of solid obstacles. Finally, we show how separation occurs in these systems depending on the properties of the particles.


I. INTRODUCTION
The transport of suspended particles in a two-dimensional (2D) periodic periodic lattice of obstacles has recently attracted significant attention as a means of separating particle species in microfluidic devices. A driving force induces particles to move through the system. Depending on the properties of the particles and the direction of the driving force θ with respect to the lattice orientation, the particles move through the lattice in commensurate phase-locked trajectories. That is, for a range of driving-force angles single species of particles will have an average trajectory in one of the lattice directions [p, q], whose angle is given by α = arctan(q/p). Further, that average trajectory angle will remain fixed over a range of driving-force angles. When the driving-force angle is outside that range, the particle trajectories will jump to a new fixed angle. The relation between the direction of motion α and the angle of the driving force θ depends on the properties of the particles. Consequently, the system has the potential to separate different species of particles, with the advantage that different particles migrate at different angles (vector chromatography), allowing for continuous fractionation [1,2].
These microfluidic systems can be broadly classified into two groups, depending on the nature of the energy landscape that is experienced by the suspended particles. In one group, the particles move through a periodic array of solid obstacles, such as in the separation devices based on deterministic hydrodynamics pioneered by Austin's group [3,4,5,6]. In terms of the corresponding energy landscape, the array of obstacles can be considered as a periodic array of hard-core repulsive potentials. In the second class of systems, particles move through smooth potential landscapes, with soft interactions between the particles and an external field. Examples include the optical fractionation methods pioneered by Grier's group, where colloidal particles are transported through an array of holographic optical tweezers [7,8,9,10,11].
Phase-locking behavior is common to transport through periodic structures in many systems [12,13,14,15,16] as well as to non-linear dynamical systems in general [17,18].
In the transport of particles in periodic systems, Lacasta et al. [19] studied a 2D periodic arrangement of wells (or traps) in a square lattice by means of numerical simulations of the corresponding Langevin equation. Numerical results showed the presence of periodic trajectories and lateral migration, in that the particles moved, on average, at an angle α different from the orientation angle of the driving force θ. They also showed the presence of phase locking in the α vs. θ curve, with clear plateaus at large Peclet numbers. In fact, the authors identify the observed migration as a deterministic phenomenon [20].
In a separate study, Lacasta et al. [21] investigated periodic landscapes that present either repulsive (obstacles) or attractive (traps) centers located on a square lattice. In both repulsive and attractive cases, they observe similar phase locking behavior, with the corresponding plateaus in the α vs. θ curve becoming evident at large Peclet numbers (Large Peclet numbers corresponds to low temperatures in Refs. [19] and [21]).
On the other hand, separable potentials, which can be written as a sum of periodic fields in each of the two principal directions of a square lattice [20], do not exhibit the same complex behavior, as clearly shown by Pelton et al. [11]. Experimental work using optical tweezers has also demonstrated the presence of periodic trajectories and locked-in states [7,11]. However, these experiments displayed a different type of locking in the presence of Brownian motion, in which the trajectories become commensurate only in a statistical sense and only some of the locked-in states are centered on commensurate directions (A hopping model for the observed dynamics is proposed by Gopinathan and Grier [9]).
Conversely, and in agreement with the discussed numerical results at high Peclet numbers, recent simulations and experiments on the motion of non-Brownian spherical particles through periodic arrays of obstacles clearly show that deterministic systems exhibit analogous phase-locking dynamics and periodic trajectories into commensurate lattice directions Phase-locking dynamics has also been observed in more complex transport systems, including the motion of an overdamped particle through a periodic potential under timeperiodic driving forces [12,15,22] and the dynamics of a particle driven through a deformable colloidal lattice [16].
In contrast to a significant number of experiments and simulations, there are few studies that present analytical results for the transport of particles in irreducible two-dimensional systems. Pelton et al. [11] studied the overdamped motion of a single particle in a periodic landscape. In the 1D case (linear fringes) they showed the existence of locked trajectories for driving forces that are not strong enough to drive the particles over the energy barrier.
They also showed that the behavior in separable (reducible) 2D potentials is analogous to the locking observed in 1D systems. Gleeson et al. [23] studied irreducible 2D periodic potentials and derived an iterative method for calculating the average velocity of the particles in inverse powers of the external force and particle's diffusivity.
In this work, we are interested in the deterministic limit (small diffusivity) for relatively small forces, i. e., driving forces for which the effect of the potential landscape is not negligible. Therefore, we consider both the deterministic evolution as well as the high Peclet number limit for finite forces. In particular, we consider a quadratic, continuous model for the 2D periodic potential that corresponds to a piecewise-continuous, linear-force model. The equation of motion for a Brownian particle traversing a periodic force field is the Langevin equation [24], which in the high-friction limit takes the form [25]: where F(x) is the periodic force field, F 0 is an external driving force and ξ(t) is the Langevin force describing the fluctuating force exerted by the fluid on the colloidal particle. The friction constant γ is given by 6πµa, where a is the radius of the colloid and µ is the viscosity of the fluid. The Langevin force is represented by a Gaussian distribution with zero mean, ξ(t) = 0, and δ correlation, ξ i (t)ξ j (s) = 2γkT δ(t − s)δ ij . Here, we shall focus on a spatially uniform external force, F 0 (x) ≡ F 0 and, for convenience, we choose a coordinate system with the x axis oriented along the direction of the force, F 0 (x) = F 0 e x (see figure 1).

B. Periodic Potential: Piecewise-linear model
We consider the case in which the periodic force can be derived from a potential field, Here, we shall model the two-dimensional, periodic landscape as a piecewise-smooth potential that is composed of repulsive centers of size R (obstacles) located in a square lattice with lattice spacing L > R. Specifically, we consider the periodic landscape shown in figure 1 which, in the unit cell, is given by, where the center of the coordinate system coincides with the center of the cell, r is the radial position, r 2 = x 2 + y 2 , and F max gives the magnitude of repulsive centers. Since we have chosen to align the x axis with the external force the obstacle lattice will be, in general, rotated with respect to the coordinate system. We shall refer to the rotation angle between the x axis of the coordinate system and the principal axis of the square lattice e 1 as the forcing angle θ (see figure 1). We shall also non-dimensionalize our variables using with f = F 0 /F max the ratio of the driving to the repulsive force.
We can include the external force in a modified potential field as V m (x, y) = V (x, y)+x f . This modified potential is also an inverted parabola, but with its center shifted to x 0 = (−f, 0). It is clear then that the family of curves perpendicular to the equipotential lines, i. e., the particle trajectories inside the circle, are also straight lines, with center at x 0 . The same result can be obtained from direct integration of Eqs. 3.
When f > 1 the x component of the velocity is always positive and all the trajectories are unbounded. Then, we can define the migration angle α of the particles as the asymptotic angle that the trajectory of the particles forms with the principal axis of the lattice e 1 , i. e., α = θ means that the particles move colinearly with the external force; α = 0 means that the particle moves parallel to the principal direction of the square lattice (see figs. 1 and 2).
Let us note that we are in the deterministic and high-friction limits, and therefore two independent trajectories cannot cross each other in real space. In addition, we consider two-dimensional trajectories in the plane, and thus for a given value of the forcing angle θ all the trajectories must have the same migration angle α. Therefore, we can obtain a single valued α vs. θ curve by determining the angle of a single trajectory as a function of the forcing angle.

A. Analytical solutions in the unit cell
In general, a global trajectory can be segmented into a number of successive collisions between the colloidal particle and the repulsive centers. In addition, we can classify each trajectory passing through a unit cell in terms of its incoming impact parameter b 0 which we define in the present coordinate system as the y coordinate of the particle when it enters the unit cell (see fig. 2a). For b 0 > 1 the trajectories are straight lines parallel to the x axis and do not interact with the obstacle at the center of the cell. For 0 < b 0 < 1, on the other hand, we have three different cases. For impact parameters larger than a maximum value, it can be shown from Eqs. 3 that the radial component of the velocity is negative outside the circle and positive inside it: In our piecewise approximation, this means that the particle will move around the circle with r = 1, until it separates at the top of the circle and then follows a straight line parallel to the x axis. A trajectory of this type is shown in fig. 2b: the incoming particle enters the circle at the point C 1 with impact parameter b 0 < −b m , and therefore leaves the circle at such that, for b c < b 0 < b m the particle enters the repulsive centers, but the separation from the circle still occurs at the top (or bottom for negative values of the impact parameter, see below). Such a characteristic trajectory is shown in fig. 2b: the particle enters the circle at the point A 1 and follows a straight line with center at (−f, 0). Then, the particle reaches the point A 2 on the circumference of the circle, which belongs to the region above b m . Thereafter, the trajectory becomes identical to those described for b 0 > b m , with the particle separating from the circle at the top (point A 3 ). We will refer to collisions with incoming parameter larger than b c as irreversible in that, independent of the exact value of the impact parameter, all the incoming trajectories collapse into a single outgoing trajectory, with outgoing impact parameter b f = 1.
Finally, for impact parameters smaller than the critical value, 0 < b 0 < b c , the particle enters the circle (e. g. point B 1 in figure 2b), moves in a straight line, and leaves the circle on the positive side of the x axis (point B 2 ). The outgoing impact parameter (or y coordinate) in this case is given by,

B. Poincare Map and Saddle-Point Bifurcation
In figure 2b we showed the different types of trajectories that a particle follows depending on the impact parameter. For each one of these trajectories, we also show the impact parameter that the particle will have in its next collision in a neighboring unit cell (see points A 1 , B 1 , and C 1 ). If the outgoing impact parameter b f ≡ b f (b 0 ) is defined as the y coordinate of the point where the particle leaves the cell, then the next incoming parameter will be b 0 = b f − β, with β = sin(θ), as shown in figure 2a. Therefore, for positive values of the incoming impact parameter, the impact parameter for the next collision is given by: where + = ( /2)(cos(θ) + sin(θ)) is the maximum possible value of the incoming impact parameter (see fig. 2b). The symmetric conditions apply to negative values of the impact parameter, with the only difference being that those trajectories leaving at the bottom of the cell come into the next cell from the top, that is: We can then investigate the global trajectories by studying the above transformation of the impact parameter, which is in fact a Poincare map of the impact parameter into itself b 0 = M (b 0 ) [26]. Our objective is to compute the migration angle α as a function of the forcing angle θ. We showed before that for a given θ all the trajectories have the same migration angle. Therefore we can investigate single trajectories using the Poincare map to determine the asymptotic direction of motion. This is particularly simple in the case of periodic trajectories.
Let us define a periodic trajectory to have period q if the particle goes through q different collisions before repeating its motion (Note that all periodic trajectories are commensurate with the obstacle lattice) [17]. This corresponds to q successive collisions with different incoming impact parameters before the sequence repeats itself. Periodic trajectories with period q correspond to fixed points b 0 of the q-times composed map of the impact parameter, The corresponding collision not only penetrates the parabolic regions but is also no longer irreversible. However, these trajectories are still locked into the α = 0 overall motion.
Finally, as the forcing angle increases, the map goes through a tangent bifurcation for the figure), when the stable and unstable fixed points meet at where the diagonal is tangent to the map (this bifurcation is sometimes also referred to as a saddle node bifurcation [27]).
Although in the vicinity of the bifurcation point the behavior is universal, the existence of irreversible collisions changes the global dynamics of the map, which in general is not chaotic.
For example, for angles slightly larger than θ b , the map exhibits the universal intermittent behavior associated with tangent bifurcations, with long-lived intervals of quasi-periodic motion of period q = 1, as the impact parameter goes across the b 0 ∼ b b region [27] (see ref. [28] for a discussion of the corresponding renormalization-group approach in a generic tangent bifurcation, and ref. [29] for a more recent perspective on tangent bifurcations in the context of Tsallis statistics).
However, the behavior outside this quasi-periodic region is also periodic due to repeated irreversible collisions. A typical trajectory with such intermittent behavior is shown in fig.   3b. The trajectory is quasi-periodic with q = 1 for b 0 ∼ b b . On the other hand, outside the near critical region, the trajectory is also periodic after two irreversible collisions, with period q = 9 as shown in figures 3b and 4.
In general, there are two types of irreversible collisions, which we may call positive and negative ones, with +1 and −1 as the outgoing impact parameter, respectively, independent of the value of the incoming parameter. Therefore, a trajectory becomes periodic whenever a second positive or negative collision occurs, which again highlights the commensurate nature of the periodic orbits.
We can still use the Poincare map to determine the phase locking behavior for all forcing angles. In the q = 9 periodic trajectory shown in fig. 3b, for example, the particle goes through the bottom-top periodic boundary condition only once, which means that its asymptotic migration angle will be arctan(1/8), as shown in figure 5. Figure 5 also shows the behavior near the critical angle θ c and the entire θ vs. α curve, which exhibits the typical Devil's staircase structure [30]. We also show some of the observed locking angles and their ordered structure, in which they form a Farey tree [31] .

C. Limiting Behavior
In figure 6 we present the impact parameter map for different magnitudes of the external force. In the limit of large external forces the repulsive centers have negligible influence on On the other hand, for small forces, that is for forces f ∼ 1, all collisions become irreversible, in that b c → 0 for f → 1. Figure 6 shows that the map of the impact parameter tends to a piecewise-continuous map with only two regions: the region of unperturbed trajectories for |b 0 | > 1, and the region of irreversible collisions for |b 0 | < 1. This is a particularly interesting limit, in that it corresponds to impenetrable repulsive cores which could model, for example, the transport of a tracer particle through an array of impermeable obstacles in the limit of high Peclet numbers, as we investigate in more detail in the next section.

IV. STOCHASTIC TRANSPORT: HIGH PECLET NUMBER BEHAVIOR
In this section we investigate the stochastic transport of colloidal particles in the periodic landscape discussed in section II B, i. e., we consider diffusive transport in addition to the purely convective motion considered in the deterministic case. We will show that, at relatively high Peclet numbers, the average motion of the particles exhibits directional locking equivalent to that observed in the deterministic case.
In the presence of diffusive transport, the effective migration angle is given by the angle between the average velocity of the particles and the periodic lattice [32]. We first solve the Fokker-Planck equation for the probability density associated with the stochastic motion of the colloidal particles given by Eq. 1. The Fokker-Planck equation in non-dimensional variables reduces to: where the Peclet number is given by Pe = F max R/Dγ, with D the diffusivity of the colloidal particle. The asymptotic distribution of colloidal particles in the unit cell, P ∞ (x), corresponds to the steady state solution of equation 11 above inside the unit cell and using periodic boundary conditions [33]. We can then obtain the components of the average migration velocity, U i , as well as the migration angle α from their ratio [33], We obtained the stationary solution for the probability distribution using standard numerical methods. In figure 7 we present the results for the average migration angle as a migration and the forcing angle tends to a structure similar to those observed in the deterministic case, consisting of plateaus and steps. In fact, we show that the average effective angle exhibits, in the limit of high Peclet numbers, the same directional locking as that predicted in the deterministic case (see figure 7).
In figure 8 we present the asymptotic probability distribution in the unit cell for different Peclet numbers. We can see that, as the Peclet number increases and convective transport becomes dominant, the probability flux is dominated by a stream that leaves the obstacles from the bottom of the circle and is parallel to the force (x axis). This is consistent with our interpretation in the deterministic case. Note that the solutions presented here correspond to the limiting case f = 1.0 and, therefore, there is no penetration into the obstacles in the analogous deterministic case, as discussed in section III C.
The plot in figure 8 shows In the plots corresponding to Pe=100 and 500 we can see the re-entrance of the probability stream on the left boundary of the unit cell, as well as its collision with the obstacle. We can estimate the width of the probability peak to validate that it corresponds to diffusive spreading in the direction perpendicular to the convective motion. For this we follow the approximations done in the case of geometric ratchets [34], in which diffusion along the direction of the external force is neglected. In this case, where /2 is the distance between the streaming point and the cross-section measurement.
This estimate agrees well with the numerical results. For example, for Pe=500 the previous Eq. predicts ∆y ∼ 0.09 and the width of the peak is σ y ∼ 0.06.
Finally, in figure 9 we compare the results of our model to the case of impermeable solid obstacles. We discussed in the previous section that, in the limit of f ∼ 1 all collisions are irreversible (b c = 0 for f = 1) and do not penetrate the obstacles. In this case we expect the dynamics to be similar to that in the case of solid obstacles. This is only approximate, given that in the presence of diffusion, there is a non-zero probability of finding the particles inside the repulsive regions. Figure 9 shows that, in fact, the migration angles obtained in our model are very similar to the migration angles obtained in the case of solid obstacles for the same, relatively large, Peclet number.

V. RELEVANCE TO MICROFLUIDIC DEVICES
In experimental work, it has been shown the possibility to separate particles using differences in the value of the bifurcation angle. In dimensional variables, the bifurcation angle is given by L sin(θ b ) = R/f , which shows its dependence on the properties of the particles.
Specifically, the properties of the particles could come into play through the force ratio f or through the effective size of the repulsive centers R.
In the case of optical lattices, both the magnitude of the repulsive force as well as the characteristic size of the repulsive centers depend on the particle size [11]. In fact, our results for the bifurcation angle are analogous to the calculation by Pelton et al. [11] for the critical angle to escape a single barrier. Our results can then be used to investigate the behavior at angles above the bifurcation.
Finally, we can investigate the transport of suspended particles in a pattern of solid obstacles using the limiting case f ∼ 1 discussed in III C.
In this case, sin(θ b ) = R/L, which depends on the size of the particles through the effective size of the obstacles, R = R 0 + a, where R 0 is the size of the obstacles and a is the radius of the particles. In this trivial approximation the bifurcation angle can also be obtained by straightforward geometrical considerations. The size dependency implies that larger particles will become unlocked from the α = 0 direction at larger angles of the driving force, as observed in experiments.
Consider, for example, the separation of particles of two different sizes, a 1 = 4µm and a 2 = 6µm, in a lattice of (solid) cylindrical obstacles of radius R 0 = 5µm and with a lattice constant L = 25µm. Let us assume that the particles are driven at an average velocity of U ≈ 10µm/s, which corresponds to a large value of the particles Peclet number, Pe∼ 10 3 . Then, considering only hard-sphere interactions (excluded volume effects) between the spheres and the obstacles the corresponding effective size of the obstacles is R 1 = 9µm and R 2 = 11µm for particles of radius a 1 and a 2 , respectively. Similarly, the dimensionless length depends on the size of the particles. Specifically, 1 = 25/9 and 2 = 25/11, which results in different bifurcation angles for the two type of particles, θ b1 = 21.1 • and θ b2 = 26.1 • for sizes a 1 and a 2 , respectively. Therefore, for a driving force oriented at any angle θ relative to the obstacle lattice such that θ b1 < θ < θ b2 , the particles will separate. In particular, the small particles will move in the lattice direction [1,2] with α 1 = 26.56 • and the large particles will be locked at α 2 = 0 • .

VI. SUMMARY
We have shown that the transport of particles in a periodic lattice of repulsive centers exhibits analogous behavior to that observed in microfluidic systems. The simplicity of the parabolic repulsive potentials allowed us to integrate the trajectories explicitly and showed the existence of periodic trajectories that are commensurate with the obstacle lattice. We also showed that the motion can be determined by means of a Poincare map of the incoming impact parameter into itself, which shows that there is a tangent bifurcation at the critical forcing angle for which the locking becomes different from α = 0. The entire migration-angle vs. forcing-angle curve shows the characteristic Devil's staircase type of structure common to phase-locking systems. Finally, we showed that the limiting behavior for impenetrable obstacles is equivalent to the high Peclet number limit in the case of transport of particles in a periodic pattern of solid obstacles. Therefore, our previous results provide insight into the separation problem in the case of periodic potential landscapes as well as in the case of periodic patterns of solid obstacles. In fact, we discuss a straightforward application of our results to calculating the bifurcation angle in both solid and optical lattices.