Focusing and sorting of ellipsoidal magnetic particles in microchannels
Abstract
We present a simple method to control the position of ellipsoidal magnetic particles in microchannel Poiseuille flow at low Reynolds number using a static uniform magnetic field. The magnetic field is utilized to pin the particle orientation, and the hydrodynamic interactions between ellipsoids and channel walls allow control of the transverse position of the particles. We employ a farfield hydrodynamic theory and simulations using the boundary element method and Brownian dynamics to show how magnetic particles can be focussed and segregated by size and shape. This is of importance for particle manipulation in labonachip devices.
pacs:
Nowadays microscopic labonachip devices have become powerful tools to analyse, manipulate and control droplets Dreyfus et al. (2003); Stone et al. (2004); Link et al. (2004); Beatus et al. (2006, 2012), biological particles Pamme (2007); Bhagat et al. (2010); Sajeesh and Sen (2014) and active colloids Lindner (2014); Das et al. (2015); Bechinger et al. (2016). Different particle types can be separated in microfluidic channels where a steady Poiseuille flow is imposed Nott and Brady (1994); Squires and Quake (2005). In particular, positional control along the transverse direction of the channel is desirable in order to transport particles to outlets at different target positions. Under high Reynolds number flow, inertial forces lead to a migration of particles towards stable positions Segre and Silberberg (1961); Di Carlo (2009) which can be manipulated by feedbackcontrol Prohm and Stark (2014). In contrast, at low Reynolds number, which is the usual regime at the micron scale, spherical and elongated particles cannot achieve net transverse motion in the absence of external forces Bretherton (1962); Goldman et al. (1967); Pozrikidis (2005a); Guazzelli and Morris (2012). In this work, we focus on this low Reynolds number regime.
In labonachip devices, magnetic forces are commonly used to manipulate the position of microscopic particles Yan et al. (2012); Hejazian et al. (2015); Huang et al. (2017) or artificial microswimmers Peyer et al. (2013); Hamilton et al. (2017). For example, the segregation of different particle types can be realized by applying an external magnetic field gradient, which essentially acts as a body force Hejazian et al. (2015). Although this method can be used to segregate different types of particles, it can not be used to focus particles to a specific transverse target position. When a uniform magnetic field is applied instead of a gradient field, the particle will only experience a torque but no force, i.e. a uniform field is useful to change the orientation of the particle Erb et al. (2016), but it is not an intuitive way to achieve translation. Interestingly, Zhou et al. Zhou et al. (2017) recently showed that paramagnetic ellipsoidal particles can be focussed to the channel center by applying a static uniform magnetic field perpendicular to the flow. They managed to achieve net motion away from the wall by breaking the symmetry of cyclic updown motion Pozrikidis (2005a) of the ellipsoid.
Here we show that the particle position can be controlled not only to the channel center, but to arbitrary target channel positions by using a static uniform magnetic field to pin the orientation of the magnetic particles. Firstly, we show that the particle will continuously move either towards or away from the wall, purely by hydrodynamic particlewall interactions. Secondly, we demonstrate that the ellipsoidal particle can be focused to arbitrary transverse target positions just by a simple manipulation of the magnetic field.
We consider a permanent magnetic particle with prolate shape of volume , suspended in a Newtonian fluid of viscosity and density . The particle has magnetization , and it is assumed to be neutrally buoyant for simplicity. It has one semiaxis of length and two of length , where is the particle’s aspect ratio . The particle has a magnetic moment where is the magnetic moment parallel to the particle’s major axis and is the particle orientation angle [Fig. 1(a)(b)]. The particle is initially placed a distance away from a infinite plane wall located at , and it experiences a magnetic torque due to a uniform external field applied to the whole domain. We assume that is oriented in ()plane, , where is the strength and the orientation of the field which are both kept constant [Fig. 1(a)]. Note that we only consider inplane motion of the particle in this paper, because a strong magnetic field will orient the major axis of the particle inplane Almog and Frankel (1995). We introduce a nondimensional parameter that describes the strength of the magnetic torque compared to the hydrodynamic torque as
(1) 
where is the local shear rate of the flow around the particle. For example, when we assume that the particle magnetization where is the permeability of free space, particle size , water viscosity , water density , shear rate and magnetic field , the particle Reynolds number is and .
We use the boundary element method Pozrikidis (1992); Ishikawa et al. (2006); Mitchell and Spagnolie (2015) to solve for particle trajectory. When inertial effects are negligible, the flow field of a given point under Stokes flow can be described using a boundary integral formulation Pozrikidis (1992): where is the Green’s function, is the background flow, and is the viscous traction acting at a point on the particle surface. The Blake tensor Blake (1971) is used for the Green’s function to account for the walls. Integrating the traction force on the surface of ellipsoid gives the hydrodynamic force and torque acting on the particle. As the system is force and torquefree, these satisfy , , where is the hydrodynamic centre of the particle Kim and Karrila (1991). A given surface material point on the ellipsoid moves with a velocity , where , are the translational and rotational velocity of the particle, respectively. The surface of the ellipsoid is divided into triangular elements and nodes. The velocities are obtained by solving the dense matrix with a known vector and an unknown vector , and is a matrix with size based on equations above Ishikawa et al. (2006). The particle position is updated using the firstorder Euler method with a time step . The software is written in CUDA and all processes are parallelized Matsunaga et al. (2014).
First, we show that transverse motion can be manipulated by pinning the rotational motion of the particle in shear flow, (). The rotational motion of an ellipsoidal particle subjected to shear and a magnetic field was discussed in Ref. Almog and Frankel (1995), in the absence of a wall. The authors showed that the particle moves in the shear plane for sufficiently large and reaches a stable angle where it is pinned by the magnetic field. The general expression for the inplane rotational velocity is
(2) 
and is obtained by solving . Note that the first term of Eq. (2) is due to the magnetic torque aligning the particle towards the field orientation , and the second term is simply Jeffrey’s rotation of an ellipsoid in flow Jeffery (1922) with and Koenig (1975).
Figure 1(c) is a schematic of the motion of an ellipsoid in shear flow at different , now in the presence of a surface. In the absence of a magnetic field () the particle rotates and oscillates along the direction, but has no net displacement along Pozrikidis (2005a). However, when the magnetic field is strong enough to pin the orientation, the ellipsoid either continuously travels upwards or downwards. The transverse motion can be explained by hydrodynamic interactions between the pinned ellipsoid and the wall. The wall can be considered to act as an image stresslet Smart and Leighton (1991); Zhao et al. (2011); Nix et al. (2014), and the leading order contribution to the lift velocity arises from the stresslet component Blake (1971) evaluated for the stable angle :
(3) 
where , and to leading order it is sufficient to approximate by , which is its value in free space Kim and Karrila (1991), given by
(4) 
where , and are shape functions Kim and Karrila (1991); Sup () that are only a function of the eccentricity . Since [see inset of Fig. 1(d)], the stresslet changes its sign only for as shown in Fig. 1(d). Therefore, the particle moves away from the wall for , while it moves towards the wall for . Figure 1(e) shows simulation and theoretical results for the lift velocity under strong orientational pinning for different distances of the particle from the surface. Very good agreement is obtained for . Deviations occur close to the wall where higher order terms in Eq. (3) play a role. We also ignored the fact that the stresslet itself is modified due to the presence of the wall.
Next we show that the magnetic particle can be focused to an arbitrary transverse position under Poiseuille flow between two walls. This geometry is an approximation for high aspect ratio rectangular channels away from side walls. The background velocity profile is , with the shear rate at the wall and the distance between the two walls [Fig. 1(b)]. [Eq. (1)] can be locally defined as with describing the value at the wall. Since the rotational motion is much faster than the translational motion, the particle angle can be assumed to be quasistatic for a given position . Hence the farfield approximation of the positiondependent stable angle follows by simply solving [Eq. (2)] at each . As shown in Fig. 2(a)(c), the angle increases with because the local vorticity of the Poiseuille flow monotonically increases with while the magnetic contribution is constant throughout the space. Note that at the channel center since the local shear rate is zero, and hence the angle is only determined by the magnetic torque. Again, the farfield approximation to the transverse velocity is obtained by stresslet images (3), but with revised position factor to take into account the effect of two walls. The velocity is shown in Fig. 2(a)(c). Since is always positive, the sign of the stresslet alone determines the directional movement of the particle: equivalently, the stable angle determines the direction [see Fig. 1(d)].
A stable fixed point , which is determined by and , is required to focus the particle to a specific position. When the magnetic field is applied perpendicular to the flow direction () as shown in Fig. 2(a), the particles are focused to the channel center because the change in sign of to gives the conditions for a stable fixed point at [inset of Fig. 2(a)]. This relation is reversed when the magnetic field points in the flow direction () [Fig. 2(b)], and all particles move towards the walls. In general, a particle position satisfying is a stable fixed point, while a position with or is an unstable fixed point.
Hence by changing the direction of the magnetic field , it is possible to focus particles to arbitrary target positions. The stable fixed point can be predicted solving [Eq. (2)] for as
(5) 
where and are defined after Eq.(2). For example, when a field is applied in direction with [Fig. 2(c)], the stable fixed point of an ellipsoid () is and the particles are focused to this position. This we also confirm by performing boundary element simulations with different initial positions, shown in Fig. 2(d), and we observe that and obtained from simulations qualitatively agree with the farfield results. Note that the hydrodynamic contributions from the two walls are calculated considering the images of both walls Blake (1971) in the simulation, which is enough for relatively large channels , while higher order reflections are required for a much narrower channel Sup (); Mathijssen et al. (2016). Inset of Fig. 2(d) describes the stable fixed point [Eq. (5)] for ellipsoids () with . The figure shows shifts toward the bottom wall with increasing or , and the particles can be focused to arbitrary positions in the lower half of the channel (). By symmetry, the particles would be focused to the upper half of the channel for . Although the farfield theory predicts that the particles cannot cross the center line because , in reality they can do so because of their finite size, as confirmed by boundary element simulations [Fig. 2(d)].
To see the effect of side walls, at and , on the motion of the particles, we extended our simulation scheme to a rectangular channel geometry by using a triangular mesh both for the particles and the walls Sup (); Pozrikidis (2005b); Hu et al. (2012); Mortensen et al. (2005). When the particles are not too close to the walls we observe very similar trajectories, focusing points and focusing times as without side walls Sup () [Fig. 2(d)]. Moreover migration in direction is negligible.
Finally we show how a static magnetic field can be used to separate particles of different aspect ratio even in the presence of thermal fluctuations. The particles are initially uniformly distributed in the lower half of the channel Squires and Quake (2005), and we consider the same magnetic field and channel height as discussed above (, , ), wallshear rate , and viscosity of water ( ). We use Brownian Dynamics simulations at room temperature, solving the equations
(6)  
(7) 
for different particle size and aspect ratio . Here , and is the full 3D particle reorientation rate for the particle orientation with , Almog and Frankel (1995). is calculated from the translational diffusion tensor where is a symmetric 3x3 matrix Sup () and , where and are the respective longitudinal and transverse diffusion coefficients of an ellipsoid of aspect ratio with shape functions Kim and Karrila (1991); Han et al. (2006); Sup (); Cobb and Butler (2005); Corato et al. (2015). The rotational diffusion constant with the shape function Koenig (1975); Kim and Karrila (1991); Sup (). The random numbers and model Gaussian white noise with zero mean and ().
Distributions for 1000 particles of size for after s are shown in Fig. 2(e). Our results clearly show that particles of different shape can be separated to different target positions , given by Eq. (5), by applying a static magnetic field. 50% of the particles reach the target region in experimentally feasible times [7s () to 20s ()] and traveling distances ( mm). Note, the focusing times are even smaller for higher confinement Sup (). Efficient separation is only possible for particles of size , where the width of the steady state distribution Honerkamp (1994) is smaller than the distance between two peaks [see inset of Fig. 2(e)]. We find an approximate analytic expression for by linearizing the drift velocity around the fixed point , where only depends on the system parameters Sup (). We solve for the steady state distribution where we introduced a potential with which keeps the particle near its target position . Since and we obtain .
We have shown that the transverse position of magnetic ellipsoidal particles in microchannel Poiseuille flow can be controlled by a static magnetic field. This is due to the hydrodynamic interactions of the ellipsoids with the channel walls. Our method can be used to focus and segregate magnetic particles which is of importance for particle manipulation in labonachip devices.
Acknowledgements.
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 665440 and under the Marie SklodowskaCurie grant agreement No 653284.References
 Dreyfus et al. (2003) R. Dreyfus, P. Tabeling, and H. Willaime, Phys. Rev. Lett. 90, 144505 (2003).
 Stone et al. (2004) H. A. Stone, A. D. Stroock, and A. Ajdari, Annu. Rev. Fluid Mech. 36, 381 (2004).
 Link et al. (2004) D. R. Link, S. L. Anna, D. A. Weitz, and H. A. Stone, Phys. Rev. Lett. 92, 054503 (2004).
 Beatus et al. (2006) T. Beatus, T. Tlusty, and R. BarZiv, Nat. Phys. 2, 743 (2006).
 Beatus et al. (2012) T. Beatus, R. H. BarZiv, and T. Tlusty, Phys. Rep. 516, 103 (2012).
 Pamme (2007) N. Pamme, Lab Chip 7, 1644 (2007).
 Bhagat et al. (2010) A. A. S. Bhagat, H. Bow, H. W. Hou, S. J. Tan, J. Han, and C. T. Lim, Med. Biol. Eng. Comput. 48, 999 (2010).
 Sajeesh and Sen (2014) P. Sajeesh and A. K. Sen, Microfluid. Nanofluid. 17, 1 (2014).
 Lindner (2014) A. Lindner, Phys. Fluids 26, 101307 (2014).
 Das et al. (2015) S. Das, A. Garg, A. I. Campbell, J. Howse, A. Sen, D. Velegol, R. Golestanian, and S. J. Ebbens, Nat. Commun. 6, 8999 EP (2015).
 Bechinger et al. (2016) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Rev. Mod. Phys. 88, 045006 (2016).
 Nott and Brady (1994) P. R. Nott and J. F. Brady, J. Fluid Mech. 275, 157 (1994).
 Squires and Quake (2005) T. M. Squires and S. R. Quake, Rev. Mod. Phys. 77, 977 (2005).
 Segre and Silberberg (1961) G. Segre and A. Silberberg, Nature 189, 209 (1961).
 Di Carlo (2009) D. Di Carlo, Lab Chip 9, 3038 (2009).
 Prohm and Stark (2014) C. Prohm and H. Stark, Lab Chip 14, 2115 (2014).
 Bretherton (1962) F. P. Bretherton, J. Fluid Mech. 14, 284 (1962).
 Goldman et al. (1967) A. Goldman, R. Cox, and H. Brenner, Chem. Eng. Sci. 22, 653 (1967).
 Pozrikidis (2005a) C. Pozrikidis, J. Fluid Mech. 541, 105 (2005a).
 Guazzelli and Morris (2012) E. Guazzelli and J. F. Morris, A Physical Introduction to Suspension Dynamics (Cambridge University Press, 2012).
 Yan et al. (2012) J. Yan, M. Bloom, S. C. Bae, E. Luijten, and S. Granick, Nature 491, 578 (2012).
 Hejazian et al. (2015) M. Hejazian, W. Li, and N.T. Nguyen, Lab Chip 15, 959 (2015).
 Huang et al. (2017) W. Huang, F. Yang, L. Zhu, R. Qiao, and Y. Zhao, Soft Matter , (2017).
 Peyer et al. (2013) K. E. Peyer, L. Zhang, and B. J. Nelson, Nanoscale 5, 1259 (2013).
 Hamilton et al. (2017) J. K. Hamilton, P. G. Petrov, C. P. Winlove, A. D. Gilbert, M. T. Bryan, and F. Y. Ogrin, Sci. RepUK 7, 44142 EP (2017).
 Erb et al. (2016) R. M. Erb, J. J. Martin, R. Soheilian, C. Pan, and J. R. Barber, Adv. Funct. Mater. 26, 3859 (2016).
 Zhou et al. (2017) R. Zhou, F. Bai, and C. Wang, Lab Chip 17, 401 (2017).
 Almog and Frankel (1995) Y. Almog and I. Frankel, J. Fluid Mech. 289, 243 (1995).
 Pozrikidis (1992) C. Pozrikidis, Boundary integral and singularity methods for linearized viscous flow (Cambridge University Press, 1992).
 Ishikawa et al. (2006) T. Ishikawa, M. P. Simmonds, and T. J. Pedley, J. Fluid Mech. 568, 119 (2006).
 Mitchell and Spagnolie (2015) W. H. Mitchell and S. E. Spagnolie, J. Fluid Mech. 772, 600 (2015).
 Blake (1971) J. R. Blake, Math. Proc. Cambridge 70, 303 (1971).
 Kim and Karrila (1991) S. Kim and J. S. Karrila, Microhydrodynamics  Principles and Selected Applications (Dover Publications, Inc., 1991).
 Matsunaga et al. (2014) D. Matsunaga, Y. Imai, T. Omori, T. Ishikawa, and T. Yamaguchi, J. Biomech. Sci. Eng. 14 (2014).
 Jeffery (1922) G. B. Jeffery, P. Roy. Soc. AMath Phy. 102, 161 (1922).
 Koenig (1975) S. H. Koenig, Biopolymers 14, 2421 (1975).
 Smart and Leighton (1991) J. R. Smart and J. D. T. Leighton, Phys. Fluids 3, 21 (1991).
 Zhao et al. (2011) H. Zhao, A. P. Spann, and E. S. G. Shaqfeh, Phys. Fluids 23, 121901 (2011).
 Nix et al. (2014) S. Nix, Y. Imai, D. Matsunaga, T. Yamaguchi, and T. Ishikawa, Phys. Rev. E 90, 043009 (2014).
 (40) See Supplemental Material.
 Mathijssen et al. (2016) A. J. T. M. Mathijssen, A. Doostmohammadi, J. M. Yeomans, and T. N. Shendruk, J. Fluid Mech. 806, 35 (2016).
 Pozrikidis (2005b) C. Pozrikidis, J. Eng. Math. 53, 1 (2005b).
 Hu et al. (2012) X.Q. Hu, A.V. Salsac, and D. BarthèsBiesel, J. Fluid Mech. 705, 176 (2012).
 Mortensen et al. (2005) N. A. Mortensen, F. Okkels, and H. Bruus, Phys. Rev. E 71, 057301 (2005).
 Han et al. (2006) Y. Han, A. M. Alsayed, M. Nobili, J. Zhang, T. C. Lubensky, and A. G. Yodh, Science 314, 626 (2006).
 Cobb and Butler (2005) P. D. Cobb and J. E. Butler, J. Chem. Phys. 123, 054908 (2005).
 Corato et al. (2015) M. D. Corato, F. Greco, G. D’Avino, and P. L. Maffettone, J. Chem. Phys. 142, 194901 (2015).
 Honerkamp (1994) J. Honerkamp, Stochastic Dynamical Systems (VCH, New York, 1994).
Supplemental Material for “Focusing and sorting of ellipsoidal magnetic particles in microchannels”
Appendix A Table of shape functions
a.1 Stresslet
Shape functions Kim and Karrila (1991) appear in Eq. (4) are functions of eccentricity as
(8)  
(9)  
(10)  
(11)  
(12) 
a.2 Brownian dynamics

matrix M in translational diffusion tensor
(13) 
shape coefficients Kim and Karrila (1991) for rotational diffusion
(14) (15) (16)
Appendix B Validation of boundary element simulation
In this section, we validate the accuracy of our software by calculating the translational velocity and rotational velocity of a sphere close to a wall where analytic expressions exist Goldman et al. (1967). The problem setup is the same as the schematic Fig.1(a) in the main text, but a sphere (radius ) with no magnetic moment is used for the validation.
Appendix C Theory for steady state distribution
In the following we want to calculate how strongly a particle is trapped near a stable point, and compare it with thermal fluctuations. Therefore we linearize Eq. (2) of the main text around and by taking and and find the linearized stable solutions . Retaining only the linear terms we obtain the relation
(17) 
where we have substituted by (note that drops out here). It can easily be checked that . Now we can calculate the velocity near the fixed point as , and linearize it which results in with
(18) 
where we skipped the dependence on of the parameters ,, and . Note that and has dimension . We then write the effective Langevin equation for the motion of the ellipsoidal particle in the direction near the fixed point as
(19) 
Since the stable orientation is in the direction the particle fluctuations in this direction are around its longitudinal axis, hence we use here its longitudinal friction coefficient . We also introduced a quadratic potential
(20) 
defined via . The stationary solution of Eq. (19) is a Gaussian distribution Honerkamp (1994)
(21) 
with standard deviation . Since , we get the scaling , see inset of Fig. 2(f) of the main text.
Appendix D Effect of confinement
Figure 4 shows the focusing velocity between two meshed walls, which take into account the higher order reflections. If the distance between the two walls takes the value used in the main text, the farfield theory matches the simulation well. Hence this also validates the method of using two meshed walls. When the distance is decreased, the velocity starts to deviate from the farfield theory, and the true velocity is faster than the prediction. Also note that the particle can easily cross the centerline for higher confinement. The focusing point shifts slightly toward the channel center , but this effect is not large, in particular for .
Method
The numerical method used here is based on previous studies Hu et al. (2012); Pozrikidis2005, with the confinement effect taken into account by using two meshed walls. Contrary to the mesh on the particle, noslip boundary conditions are applied for the walls. The dimensions of the two walls (  ) and the mesh size ( ) were changed depending on the confinement distance .
Appendix E Effect of side walls: focusing in a rectangular channel
Finally, we show particle focusing in a rectangular channel (width , height ). Figure 5 shows crosssectional stream lines for two different channel aspect ratios ( and ). In the case of two infinite planar walls, the particle is focused to as discussed in the main text. If the channel is rectangular, most particles still focus to . However particles that are located distances less than away from a side wall reach a curved focusing region at or the top () or bottom () walls. This effect can be controlled by increasing the channel aspect ratio .
Method
We conducted full 3D simulations based on a method used in previous studies Hu et al. (2012); Pozrikidis2005. The boundary integral formulation in this system is
(22) 
where the notations describe integrals over the ellipsoidal particle (P), side walls (W) and outlet cap (O) respectively. is the Poiseuille velocity profile inside a rectangular channel Mortensen et al. (2005), is the normal vector pointing into the channel, and is a additional pressure drop due to the presence of the particle:
(23) 
where is the flow rate. The characteristic shear rate is defined at the wall in the twowall system, while the characteristic shear rate at is used in the rectangular channel. The channel length and mesh size are used in the simulation. In order to avoid error arising from the coarse mesh, we kept the distance between particle and walls larger than . The accuracy of the simulation was validated in a cylindrical channel by comparing with the analytical solution (data not shown).