A pair of particles in inertial microfluidics: eﬀect of shape, softness, and position † deeper understanding of the inertial Soft Matter

Lab-on-a-chip devices based on inertial microfluidics have emerged as a promising technique to manipulate particles in a precise way. Inertial microfluidics exploits internal hydrodynamic forces and the mechanical structure of particles to achieve separation and focusing. The article focuses on the hydrodynamic interaction of two particles. This will help to develop an understanding of the dynamics of particle trains in inertial microfluidics, which are typical structures in multi-particle systems. We perform three-dimensional lattice Boltzmann simulations combined with the immersed boundary method to unravel the dynamics of various mono- and bi-dispersed pairs in inertial microfluidics. We study the influence of diﬀerent starting positions for mono- and bi-dispersed pairs. We also change their deformability from relatively soft to rigid and choose spherical and biconcave particle shapes. The observed two-particle motions in the present work can be categorized into four types: stable pair, stable pair with damped oscillations, stable pair with bounded oscillations, and unstable pair. We show that stable pairs become unstable when increasing the particle stiﬀness. Furthermore, a pair with both capsules in the same channel half is more prone to become unstable than a pair with capsules in opposite channel halves.


Introduction
Biomedical studies of biological cells play a crucial role in diagnosing several fatal diseases 1 such as malaria, 2 cancer, 3,4 and the human immunodeficiency virus (HIV) infection, 5 to name a few. In recent years, microfluidic 6 lab-on-a-chip devices 7 have emerged as a promising technique to precisely manipulate and control biological cells needed on a commercial level. 8 Lab-on-a-chip devices for separating particles and biological cells rely on external force fields or internal hydrodynamic forces. 9 Lab-on-a-chip techniques using external force fields include optical tweezers, 10 dielectrophoresis, 11 magnetophoresis, 12 and acoustophoresis. 13 On the other hand, deterministic lateral displacements at low Reynolds number using appropriately placed pillars 14 and inertial microfluidics [15][16][17][18][19][20] exploit internal hydrodynamic forces to achieve particle separation. In contrast to common microfluidic lab-on-a-chip devices in which fluid inertia is negligible, inertial microfluidics operates in an intermediate range between Stokes and turbulent regimes, where flow is still laminar. The possibility to attain high throughput makes inertial microfluidics an attractive option among other microfluidic particle-separation techniques. In this article, we investigate how different factors such as shape, softness, and position influence the motion of a pair of particles in a microfluidic channel.
Inertial focussing was first reported by Segré and Silberberg for solid particles with the well-known tubular pinch effect 21,22 and then spurred numerous theoretical, 23-26 computational, 19,27-30 and experimental [31][32][33][34] studies to gain a deeper understanding of this phenomenon. Inertial migration of a neutrally buoyant solid particle is essentially the result of a balance between shear-gradient lift force, directed towards the channel wall, 35 and a wall repulsion force, directed towards the channel center. 36 Applying an external force to the solid particle along the channel axis, the resulting Saffman force 37 also contributes significantly to the force balance in the cross-stream direction. For almost five decades from the observation of the Segré-Silberberg effect, inertial migration was mainly a blue-skiesresearch problem. In 2007 Di Carlo et al. 38 demonstrated ordering and separation of particles in a microchannel using inertial focussing, which gave rise to the field of inertial microfluidics. Subsequently, Hur and co-workers 39 experimentally showed inertia-driven separation and enrichment of deformable capsules. Soft capsules experience an additional lift force induced by the deformability of the cell, which drives them towards the channel center and which thereby competes with the inertial lift force. 40 For biomedical applications of inertial microfluidics it is crucial to strive for a deeper understanding of the inertial migration of deformable capsules. The dynamics of soft capsules is widely studied at low Reynolds numbers in canonical flows [41][42][43][44][45] and cellular blood flow simulations. [46][47][48][49] Only recently, inertial microfluidics of deformable cells has caught attention. 20,[50][51][52][53][54][55][56][57][58] Using two-dimensional numerical simulations, Shin and Sung 51 showed that the final equilibrium position for a given value of capsule elasticity varies with the Reynolds number. In their following work 54 they identified different types of capsule dynamics such as tank-treading, swinging, and tumbling. According to Kilimnik et al. 50 the steady-state position of a deformable capsule in the lateral direction is independent of the flow rate and it only depends on the softness. This was later verified by Schaaf and Stark 20 using 3D lattice Boltzmann simulations. They also proposed that the inertial lift force in the vicinity of the channel center scales as a cube of capsule radius. More recently, Alghalibi and colleagues 58 investigated cross-stream inertial migration of a soft particle in pipe flow. Interestingly, they found that irrespective of the particle's softness and the flow rate, the particle always ends up at the channel center for an initial particle diameter of 0.2 times the pipe diameter.
Practical applications require manipulation of several soft particles at a time and, hence, it is important to study how particles interact in flow. There are a few systematic studies in literature to understand the interaction of solid and soft particles in inertial microfluidics. Yan et al. 59 investigated the dynamics of two solid particles in linear shear flow generated by two plates moving against each other. They observed that particles either settle down on the centerline where flow velocity is zero or exhibit oscillatory motion when released from symmetric initial positions. Schaaf et al. 60 reported that for two closely placed solid particles in a channel flow, the inertial lift forces acting on the leading and lagging particles differ significantly. In their subsequent work 61 they extended the understanding of a particle pair to study the behavior of particle trains 62 in inertial microfluidics. At low Reynolds numbers, Lac et al. 63 and Aouane et al. 64 investigated the interaction of two soft capsules in shear and Poiseuille flows, respectively. The inertial counterparts of these two problems were also studied. Doddi and Bagchi 65 observed a transition from irregular self-diffusion 66 to different types of spiraling motions (outward, fixed orbit, and inward) with increasing Reynolds number. Lan and Khismatullin 67 reported that a pair of deformable cells in a rectangular microchannel undergoes either swapping or passing trajectories, while deformable cells in a series experience damped oscillations. Like Lan and Khismatullin, 67 Schaaf et al. 60 also observed these swapping and passing trajectories for a pair of solid particles. Krüger et al. 68 looked into the suspension of deformable capsules in Poiseuille flow at finite Reynolds numbers. They discovered that particleparticle interactions enhance the migration towards the channel center for sufficiently large Reynolds and capillary numbers.
In the present work we perform three-dimensional lattice-Boltzmann simulations combined with the immersed boundary method 20,60,69 to investigate the dynamics of mono-and bi-dispersed particle pairs. They flow through a square microchannel in the inertial regime. We systematically investigate the effect of different starting positions for mono-and bi-dispersed pairs. We also vary their deformability from relatively soft to rigid and choose spherical and biconcave particle shapes. Reynolds number, particle radius, and axial particle distance are kept the same in all the cases. The observed dynamics of the different mono-and bi-dispersed particle pairs can be categorized into four types: stable pair, stable pair with damped oscillations, stable pair with bounded oscillations, and unstable pair.
The outline of the article is as follows. In Section 2 we discuss the computational setup and simulation parameters of the problem and explain the lattice Boltzmann simulations. Results from the numerical investigations are presented in Section 3 and we close with concluding remarks in Section 4.   1 shows the schematic of our microfluidic setup. We study the hydrodynamic interaction and inertial migration of a flowing pair of soft particles in a quadratic microchannel. We start from a fully-developed Poiseuille flow and place two particles in the midplane at y = 0. We apply periodic boundary conditions in the z-direction along the channel axis and the noslip boundary condition at the channel walls in the remaining directions. The length of the microchannel is large enough so that particles do not interact with their periodic images unless the axial separation becomes as large as the channel length. The pressure gradient or force density a required to pump the fluid through the channel is calculated from the analytical In order to obtain the macroscopic picture of the flow field, fluid density r and momentum density ru at each lattice node can be calculated using the zeroth and first moments of the populations f i : We enforce the no-slip condition at channel walls using the regularized method of Latt and Chopard. 75,76 The numerical accuracy of our simulations is maintained by setting the relaxation time t r 1. Please refer to ref. 20 for further details. The numerical solution of the lattice Boltzmann equation is obtained using an open-source parallel lattice Boltzmann solver from the Palabos project, 77 while implementation of the immersed boundary method and particle dynamics uses in-house libraries.
2.2.2 Modeling of soft capsules. We employ the finite element method to discretize the surface of the particle with P 1 elements, 78 i.e., linear triangular elements. Tessellation of the particle surface is achieved by dividing each triangle of an icosahedron into four triangles (for more details refer to ref. 69) and the newly created nodes are then radially shifted to the circumsphere. The process is repeated until the distance between two neighboring nodes approaches the lattice spacing. This procedure of generating mesh from a highly symmetric solid (icosahedron in the present case) enables us to achieve superior mesh quality in terms of isotropy and homogeneity. 69 The same procedure is used for the solid particles as well.
The capsule dynamics is determined by a combination of strain, bending, and volume energies. For the strain energy we rely on the Skalak model 79 valid for an isotropic and homogeneous capsule. It is formulated as a function of the strain invariants of the Green strain tensor and given as where k s is the shear modulus and k a is the area dilation modulus. The strain invariants I 1 and I 2 can be expressed in terms of in-plane principal extension ratios l 1 and l 2 : The discretized version of the bending energy proposed by Helfrich 80,81 takes the form where k b is the bending modulus, y ij is the angle between the area normal vectors of two neighboring surface elements, and y 0 ij is the reference angle from the undeformed (initial) state. Implementation of the bending energy prevents cusp formation and buckling of soft capsules. In this work, we consider soft particles with an impermeable membrane. To implement the constraint of constant volume, at least approximately, we choose a standard method and incorporate the volume energy 82 where k v , V and V 0 are the volume modulus, actual capsule volume, and reference volume, respectively. To keep the volume changes below 1%, we choose k v = 75 000rn 2 /W 2 . Note, these changes do not affect the dynamics of the particle.
Following Krüger et al. 68 and Schaaf and Stark, 20 we set k a /k s = 2 and k b /(k s a 2 ) = 2.87 Â 10 À3 . We can compare our parameters with the ones for red blood cells. To mimic their the mechanical behavior in blood flow simulations, the ratio k b /(k s a 2 ) is set to 2.36 Â 10 À3 , 47 which is very close to the value in present simulations. However, the ratio k a /k s in current simulations is significantly smaller than that used for a healthy red blood cell, which allows our spherical capsules to deform significantly.
The total force acting on node i can be obtained from the principle of virtual work; 83 given as The derivatives of strain, bending, and volume energies were derived analytically and then inserted directly into the solver to save computational costs. The analytical form of the derivatives can be found in ref. 84.

2.2.3
Modeling of solid particles. The motion of a neutrally buoyant solid particle is governed by Newton's law and its equivalent for the angular momentum: where, S, M, I, U, and x are the particle position, mass, moment of inertia, linear velocity, and angular velocity, respectively. Force F fluid and torque T fluid exerted by the fluid on the particle are computed using the iterative procedure of Inamuro, 85 which ensures the no-slip condition at the boundary of the solid particle. Suzuki and Inamuro 86 demonstrated that for simulating the dynamics of solid particles at finite Reynolds numbers correctly, the effect of a fictitious fluid inside the solid particle must be taken into account when using the immersed boundary method. Therefore, one needs to introduce Feng's rigid body approximation 87 with the help of additional forces and torques, F Feng and T Feng , respectively. 2.2.4 Immersed boundary method. The hydrodynamics of soft and solid particles in inertial microfluidics is essentially a fluid-structure interaction (FSI) problem on the micron scale. In the present case, the motion of the bulk fluid in a microchannel is responsible for the movement of particles; in turn, particles exert forces on the fluid via their surfaces. We implement this interaction between fluid and particles using the immersed boundary method (IBM) of Peskin. 88,89 The immersed boundary method is one of the fluid-structure interaction methods that belong to the mixed Eulerian-Lagrangian framework. There also exist fully Eulerian and Lagrangian FSI methods (see Fig. 1 of ref. 90 for more details). The classical IBM treats vertices of the P 1 elements of the discretized particle surfaces as Lagrangian markers. The advection velocity : x i of Lagrangian markers is obtained using an interpolation of the velocity values u from neighboring Eulerian points, which in our case are the lattice nodes of the lattice Boltzmann method: In contrast, the body force density f exerted by Lagrangian markers on neighboring Eulerian points is computed using the spreading operation: Following ref. 88 the discretized delta function d(X À x i (t)) in eqn (15) and (16) are represented as d( The interaction between soft capsules and the bulk fluid is resolved using the classical IBM. However, for solid particles we employ interpolation and spreading operations in an iterative manner 85,86 to obtain the force F fluid and torque T fluid directly from the flow field. Then, advection velocity U is obtained using eqn (13). This particular approach is known as the multi-direct-forcing variation of the immersed boundary method. 91,92

Simulation parameters
In our microfluidic setup we simulate two types of particle pairs, mono-and bi-dispersed (see Fig. 2). Mono-dispersed pairs consist of spherical particles of the same stiffness. Bi-dispersed pairs are classified into two subtypes; first, a pair with particles of different shapes but the same stiffness, and second, a pair with particles of varying stiffness but the same shape. We consider the combination of a spherical (leading) and biconcave (lagging) particle for the first type, and for the second type, two spherical particles. The initial orientation of a biconcave particle in the flow field is shown in a magnified view in Fig. 2. We apply the transformation x bicon: ¼ 1:2x sphere ; y bicon: ¼ 1:2y sphere ; to the vertices of a spherical particle of radius a to model the biconcave shape. 93,94 In all the cases we set the non-dimensional particle radius a and the initial axial separation Dz equal to 0.4 and 1, respectively. For the biconcave particle the radius characterizes the spherical particle from which it is generated using eqn (18). The value of the Reynolds number Re = 2WU max /n is fixed to 10 in all the simulations. Here, U max is the maximum fluid velocity in the channel center and n is the kinematic viscosity. Note, in a typical inertial microfluidics setup, the Reynolds number varies from 1 to 100 17 . Following Schaaf and Stark, 20 we use the Laplace number to define the softness of the particle; it is given as In our simulations we consider cases with relatively soft (La = 10) and stiff (La = 100) particles (please refer to Fig. 2). We also include a mono-dispersed pair of solid particles (case with La = N in Fig. 2) to compare the dynamics of soft and solid particles of spherical shape. For comparison we note that for a red blood cell suspended in plasma, one has a = 4 Â 10 À6 m and k s = 5.3 Â 10 À6 N m À1 , 47 and the Laplace number turns out to be 0.0212. In practice, there might exist a viscosity contrast between the fluid within soft particles (cytoplasm in biological cells) and the bulk liquid (Z cyto /Z bulk 4 1). For example, the cytoplasmic viscosity of the red blood cell is five times the plasma viscosity Z plasma = 1 mPa s. 47 Moreover, tumour cells such as A549 can have cytoplasmic viscosity in the range of 144.8 to 1390.7 Pa s, 95 which is several orders of magnitude higher than the plasma viscosity. We shortly summarize what influence this will have. An increase in the viscosity of the interior fluid or cytoplasm attenuates the deformation of soft particles and biological cells. 50 Thus, while for a very soft particle a viscosity contrast close to one should not significantly alter the migration dynamics, a viscosity contrast significantly larger than one makes particles less deformable and modifies their equilibrium positions. In the case of relatively stiff particles, the viscosity contrast will have a negligible influence on the particle dynamics. Moreover, for biconcave soft particles in inertial flows with Re B O(1), there exists a critical value of viscosity contrast above which the transition from tank-treading to tumbling motion takes place. 52 Later, we will see that biconcave particles in the present work already exhibit tumbling motion. With this in mind, we set the density and viscosity of the soft particle identical to the values of the bulk fluid.
In each simulation we fix the leading particle position in the x-direction to 0.2 and vary the position of the lagging particle from À0.4 to 0.4 in steps of 0.2. The leading particle's starting position is kept sufficiently far from the stable equilibrium point of isolated particles (see Fig. 5 and 9) so that disturbances from the lagging particle can have a significant effect on the motion of the leading particle. The implemented LBFEIBM solver has been employed in a series of works by our group 19,20,60,61 (see Appendix A for benchmarking). According to our previous work we discretize the channel width with 90 lattice nodes for the case of soft particles and 75 lattice nodes for solid particles. In our simulations, the necessity of a high grid resolution arises to ensure the proper coupling between particles and the bulk fluid. Moreover, relatively more lattice nodes are required in the case of soft particles to capture the particle deformation and lift force accurately. Lastly, the surface of soft and solid particles is discretized using 20 480 P 1 elements. Fig. 2 Classification of different types of particle pairs with radius a. Each pair type is simulated for the same initial distance Dz, for different values of softness La lead , La lag , and five sets of different initial positions x lead , x lag . Alphanumeric labels M1-M3 and B1-B4 are the identification tags given to each case. Note that the mono-dispersed pair with La = N corresponds to the case of solid particles.

Results and discussion
In this section, we first explain how the overall dynamics of a particle pair is categorized in the present work for the different types of particle pairs as illustrated in Fig. 2. A detailed explanation of the hydrodynamic behavior is provided later in Sections 3.1 and 3.2 for the different mono-and bi-dispersed pairs, respectively.
In Poiseuille flow the inertial lift force migrates a particle towards a specific position between channel center and wall, where softer particles are located closer to the center. Depending on the direction of migration, the particle either accelerates or slows down. For two closely placed particles the lift force acting on each particle is modified compared to the single-particle case 60 and the modified cross-streamline motion causes particles to approach or move away from each other. Notably, the initial interaction and positions of particles determine their dynamics at a later stage.
In Fig. 3 we use the time evolution of the axial particle distance Dz to identify four different types in the dynamics of a flowing particle pair: stable pair, stable pair with damped oscillations, stable pair with bounded oscillations, and unstable pair. In our notation the stable particle pair either directly approaches a constant particle separation at long times (see video M1) or the particle distance oscillates. In the second case we find both, damped oscillations ( see video M2) and ongoing bounded oscillations (see video B1). The source of damped and bounded oscillations is different as we discuss below in Sections 3.1 and 3.2. Finally, particles of an unstable pair continuously drift apart until the particle distance approaches the channel length ( see video B3).
The diagram in Fig. 4 summarizes the dynamics of the particle pairs, which we observed in our simulations for different mono-and bi-dispersed pairs when the initial position of the lagging particle is varied. The mono-dispersed pair M1 with soft capsules is always stable irrespective of the starting position x lag . It is situated in the channel center as we will see below.
On increasing the particle stiffness from La = 10 to La = 100 (label M2), the pair is stable when the leading and lagging particles initially occupy different channel halves and becomes unstable when the particles start in the same channel half. The transition between both cases occurs through damped oscillations. For the mono-dispersed pair M3 composed of rigid particles (La = N), we observe that only the symmetric pair with (x lead , x lag ) = (0.2, À0.2) is stable, reminiscent of the observations in ref. 60.
The dynamics of the bi-dispersed pairs B1 and B2 is qualitatively the same as that of pairs M1 and M2, however, with the difference that now the B1 pair performs bounded oscillations. Lastly, the symmetric placement of leading and lagging particles in the B3 and B4 pairs gives a stable pair, which we also observe for the B3 pair with x lag = 0. In summary, irrespective of particle-pair type when the leading and lagging particles start from the symmetric positions (x lead , x lag ) = (0.2, À0.2) the pair is stable. Also, stable pairs become unstable when increasing the particle stiffness. Finally, a pair with both particles in the same channel half is more prone to become unstable than a pair with particles starting in the opposite channel halves.
Before proceeding, we add a comment. Prohm and Stark 19 demonstrated that a freely placed solid particle in Poiseuille flow might either migrate to the main axis or the diagonal depending on the particle size and cross section of a microchannel. Furthermore, these zero lift-force positions can either be stable, saddle points, or unstable. Particles located on unstable positions migrate to stable positions in response to a small perturbation. In the work of Schaaf and Stark, 20 it was observed that the majority of implemented deformable capsules migrated to the channel diagonal. In our case we initialize all particles in the midplane at y = 0. We observe that in all cases they essentially stay in the midplane during migration since for symmetry reasons the lift force component normal to the midplane is zero. Fig. 3 Characterization of the dynamics of a particle pair based on the temporal evolution of the particle distance Dz. For all stable pairs the axial separation eventually equilibrates to a constant value or undergoes bounded oscillations. In contrast, the leading and lagging particles of an unstable pair drift apart until Dz approaches the channel length.

Mono-dispersed pairs
3.1.1 Single particles. To understand similarities and differences in the dynamics of a mono-dispersed pair compared to a single particle, we briefly look at single particles with the same softness as in the pairs M1-M3. Fig. 5 shows the inertial migration of a spherical particle for different values of softness. Particle radius and Reynolds number are fixed to 0.4 and 10, respectively. To make position x and time t non-dimensional in Fig. 5, the channel half-width W and viscous time scale W 2 /n are used, respectively. While rigid particles show the typical inertial focusing, the particle size in the present case is large enough so that the softness corresponding to La = 10 is sufficient for the lift forces due to deformability to dominate over the inertial lift forces. Hence, the particle settles at the channel center independent of the starting position. Note that smaller particles need smaller La or larger deformability to occupy the channel center. 20 We also show in the inset of Fig. 5 how the particle with rigidity La = 10 transforms to a bullet-shaped particle in steadystate, while the particle with rigidity La = 100 stays nearly spherical.
3.1.2 Particle pairs. For the mono-dispersed pairs M1-M3 the lateral positions of the leading and lagging particles and their axial distance are plotted versus time in Fig. 6. The softest particles (pair M1) migrate to the channel center (see video M1), however, the time needed by the leading and lagging particles to reach the center varies with the initial position of the lagging particle. For x lag = 0 and À0.2 the leading and lagging particles in the M1 pair migrate with nearly constant axial distance in time. Furthermore, when the leading and lagging particles initially move on the same streamline, x lead = x lag = 0.2, the leading particle experiences a very strong push towards the channel center, which significantly enhances its axial migration velocity. However, once it has come close to the channel center in a very short time, it takes a longer time to fully reach the channel center. Finally, we note that, initially, the axial particle distance for the M1 pair with x lag = 0.4 and À0.4 grows at the same rate. Nevertheless, the steady-state value of Dz for x lag = 0.4 is twice as large as that for x lag = À0.4. To quantify the hydrodynamic interaction of soft particles in the M1 pair, we compare the migration time t m the leading and lagging particles need to reach the channel centerline with the corresponding time of a single particle. Such a comparison allows to demonstrate how strongly the total lift force is influenced by the neighboring particle in a pair. Thus, in Fig. 7, we plot the ratio t m pair /t m single for different values of the starting position x lag . We find that the presence of the lagging particle with x lag = AE0.4 slows down the migration of the leading particle to the channel center, while for the remaining three starting positions the migration of the leading particle is enhanced. Moreover, the lagging particle migrates faster when both particles are in different channel halves initially and vice versa.
We now focus our attention on the dynamics of particles in the M2 pair. For x lag = 0.2, similar to the M1 pair, the lagging particle initially pushes the leading particle towards the channel center. Thus, leading and lagging particles drift at a different axial velocity, which continuously increases the axial distance and makes the pair unstable. Note that we stop the simulations once the axial particle distance Dz is beyond 10 and therefore close to the channel length, where a particle interacts with the image of the other particle. Continuing the simulations of the leading particle from the momentary position at x = 0.055 without the other particle present shows that it migrates towards the center instead of drifting back to the known equilibrium position at x E 0.3. Thus, we identify a second type of stable position at x = 0, which we attribute to the strong wall repulsion for the large particles, we use in our simulations, and the fact that the inertial lift force is small close to the center. The M2 pair with x lag = 0.4 also turns out to be unstable since already initially, they drift with very different axial velocities. Interestingly, for x lag = À0.4 the particle pair is stable although with a larger axial distance compared to the corresponding M1 pair with La = 10, while for the symmetric initial positions with x lag = À0.2 the axial distance stays close to the starting value similar to the case La = 10.
A new feature is observed for the initial position x lag = 0. The M2 pair is stable but reaches a constant axial distance close to the initial value only after extended damped oscillations in both the lateral positions and Dz (see video M2), as the insets in Fig. 6 show. Similar oscillations in the axial distance were reported earlier in numerical simulations of a pair of solid particles 60 and a train of soft particles 67 in a rectangular microchannel. Experimental evidence on damped oscillations of solid particles in a microchannel can be found in Lee et al. 96 Recently, based on theoretical calculations, Hood and Roper 97 suggested that in such oscillations, inertia and viscosity switch their roles, i.e., viscous interactions of a particle with the channel wall and other particles initiate oscillations while inertial forces dampen them.
We further investigated the nature of these damped oscillations for x lag = À0.05, À0.1, and À0. 15. We find that as the initial position of the lagging particle moves away from the channel center, the maximum amplitude of oscillations decreases [see Fig. 8(a)]. We assign this to the reduced viscous coupling, when the particles are further away from each other and an increased damping due to the (inertial) lift forces, which force the particle to a specific lateral position. 60 In all cases, particles equilibrate with an axial separation very close to its initial value. In the inset of Fig. 8(b) we visualize the trajectories of the leading and lagging particles in the center-of-mass frame for x lag = 0. They reveal the inward spiraling motion of the leading and lagging particles. For decreasing x lag the extent of the spirals shrink because of the reduction in the value of max(Dz). As the main plot in Fig. 8(b) shows, the damping rate of the oscillations in Dz first does not vary significantly for decreasing x lag but then increases sharply.
Finally, we discuss the case of a pair of solid particles (M3) in Fig. 6. Only the symmetric placement with x lag = À0.2 and x lead = 0.2 provides a stable flowing pair, reminiscent of the findings in ref. 60. Unlike soft particles with x lag = À0.2, the initial distance between two solid particles no longer remains constant over time. It grows until the particles have reached their final lateral positions. Earlier, we saw that in monodispersed pairs M1 and M2 with x lag = AE0.4, the axial separation grows at the same rate only at early times. Now, the temporal increase of the axial distance Dz is the same for the whole time interval for these two cases since the lagging particles already start close to their equilibrium positions (see Fig. 5) and only the leading particle has to relax.

Bi-dispersed pairs
In this subsection we concentrate on the dynamics of different bi-dispersed pairs. We employ an approach similar to that in the previous discussion on mono-dispersed pairs. We first look at the migration of a biconcave particle so that later it can be compared with bi-dispersed pairs.
3.2.1 Single particles. Fig. 9 shows the temporal evolution of the lateral position of a biconcave particle with rigidities La = 10 and 100. It migrates to a stable equilibrium point very close to that shown in Fig. 5 for a spherical particle of the same rigidity. However, for La = 10 the inertial migration is significantly slower compared to the spherical particle of identical La value.
At low Reynolds numbers a biconcave particle in shear flow undergoes either tank-treading, tumbling, or swinging motion depending on the shear rate, softness, and viscosity ratio. 98 Here, for more rigid particles at La = 100 we observe a solidbody like tumbling motion in the y = 0 plane, while for the softer biconcave particle at La = 10 a deformation-driven tumbling motion is observed as illustrated in Fig. 10. Initially [ Fig. 10(b and c)], the particle rotates because the lower part being closer to the center moves faster. Later, when the particle aligns itself with the flow direction [ Fig. 10(d)], it transforms into a parachute-like shape [ Fig. 10(f)]. Eventually, the transition to a convex disc-like shape [ Fig. 10(h)] completes one tumbling cycle, during which the particle has moved closer to the channel center. This also increases the period of the tumbling motion. Note that the biconcave particle with La = 10 does not fully reach the channel center at x = 0 and also leaves the midplane at y = 0 [inset of Fig. 9]. Therefore, it does not have a symmetric shape and performs a tumbling motion. We call it wobbling since the tumble axis varies in time, as we explain further below in Fig. 12. This is qualitatively similar to the observation of Fedosov et al. 99 for the tumbling-to-slipper transition of a red blood cell in cylindrical microchannels. A slipperlike form of the biconcave particle occurs in Poiseuille flow at low Reynolds numbers, 100 which in steady state undergoes tank-treading motion. Thus, we observe a transition from deformationdriven tumbling to wobbling when the particle approaches the center. It also occurs in the bi-dispersed pair B1 as discussed below.
3.2.2 Particle pairs. In Fig. 11 for the bi-dispersed pairs B1 and B2 we plot the lateral positions of the leading and lagging particles and their axial distance versus time. Similar to the M1 pair with La = 10 also the B1 pair with the same softness moves to the channel center. The axial distance in steady state is larger compared to the M1 pair. This is mainly because the biconcave particle migrates slower than spherical particles (compare Fig. 5 and 9). Due to the tumbling of the lagging biconcave particle, bounded oscillations in the axial distance and x lead occur. Except for x lag = À 0.2, all the other biconcave particles in Fig. 7 Migration time t m of particles in the mono-dispersed pair M1 relative to a single particle plotted versus the initial position x lag for fixed x lead . t m is the time needed to reach the channel centerline with a distance smaller than 10 À3 . Note that it is not possible to include the point x lag = 0 for the lagging particle since it is already at the equilibrium position for La = 10. the B1 pairs exhibit the transition from deformation-driven tumbling to wobbling close to the channel center as described earlier. We illustrate this transition for x lag = 0.2 in Fig. 12 using the trajectory of a Lagrangian marker initialized on the periphery of the biconcave particle in the midplane at y = 0. Initially, the deformation-driven tumbling is indicated by the zig-zag trajectory, which then transitions to a helical trajectory due to wobbling (see also video B1). The different tumbling motions are highlighted using blue and red planes in Fig. 12. An exception is observed for the biconcave particle with initial position at x lag = À0.2. It migrates very quickly towards the channel center in contrast to the single particle in Fig. 9 and performs directly the wobbling motion.
At stiffness La = 100 stable particle pairs only occur for x lag = À0.2 and À0.4. The pronounced bounded oscillations have vanished although modulations with small amplitudes are visible. All other B2 pairs are unstable and axially move away from each other while reaching their equilibrium positions with finite distance from the center. An exception is the case with x lag = 0. The leading particle is strongly dragged to the center by the biconcave particle and ultimately reaches the center, which we already identified as an additional stable lateral position when discussing the dynamics of the M2 pairs.
We now discuss the bi-dispersed pairs B3 and B4 where leading and lagging particles are spherical and have different stiffnesses. In Fig. 11 we plot lateral positions and axial distances versus time. When the leading particle is softer (La = 10) than the lagging particle (La = 100), the B3 pairs with x lag = 0 and À0.2 form stable pairs. The leading soft particle succeeds to drag the stiffer particle to the channel center, where it occupies the additional stable position so that the pair can stay bounded. In the remaining cases the lagging particle migrates towards the off-center stable position, which ultimately makes the B3 pairs unstable.
In the B4 pairs the rigidity values are switched between the particles. As a result both leading and lagging particles move towards the channel center. A special case occurs when both particles start from the symmetric lateral positions x lead = Àx lag = 0.2. They obstruct each other from aligning at the channel center, which results in oscillations in the lateral positions. However, the axial distance between the particles only fluctuates little about its initial value and the pair is stable. In the rest of the cases the leading and lagging particles successfully align at the channel center. They drift apart at different axial velocities due to their different equilibrium shapes (see Fig. 5), which directly affect the drag acting on each particle in the axial direction. Hence, even if the leading and lagging particles have the same lateral position, their axial distance does not remain bounded. As shown in Fig. 11, for the unstable pairs the axial distance increases at the same rate once both particles are at the channel center. We did not simulate the unstable pairs for longer times because it is computationally expensive as the axial distance grows very slowly. On extrapolating the trend of the curve for x lag = 0.2, one would need E400 time units for Dz to reach Dz = 10. We also simulated the B4 pairs with the starting position of the leading particle x lead = 0. However, we did not find any stable configurations in these cases.

Concluding remarks
Inertia driven separation and focusing of soft particles in microfluidics has direct relevance for biomedical applications. Inertial microfluidics has recently emerged as a robust and precise technique for particle manipulation since it relies entirely on internal hydrodynamic forces. Moreover, inertial microfluidics functions at finite Reynolds numbers, which allows faster processing of samples. In the present article we investigated the hydrodynamics of a particle pair in inertial microfluidics using 3D numerical simulations. We combined the lattice Boltzmann method to compute the fluid flow with the finite element method to model particle dynamics. To couple fluid flow and particle motion, the immersed boundary method was used.  In this work we considered mono-(M1-M3) and bi-dispersed (B1-B4) particle pairs of different softness and particle shape (see Fig. 2). For each pair we investigated the influence of different starting positions of the lagging particle and identified four types: stable pair, stable pair with damped oscillations, stable pair with bounded oscillations, and unstable pair.
We found that the starting position of the lagging particle does not affect the stability of the mono-dispersed pairs with stiffness La = 10 (M1). On increasing the stiffness to La = 100 (M2), the mono-dispersed pairs become unstable when the leading and lagging particles start in the same channel half. Moreover, when the lagging particle starts at x lag = 0, the pair undergoes damped oscillations. The oscillations arise due to hydrodynamic interactions between particles, while inertial forces damp them. The maximum amplitude of the damped oscillations decreases as we shift the starting position of the lagging particle away from the leading particle. Finally, from mono-dispersed pairs of perfectly solid particles (M3) only the pair with initially symmetric lateral positions about the centerline retains its stability.
When we change the shape of the lagging particle from spherical to biconcave, bi-dispersed pairs with particle stiffness Fig. 11 Lateral positions of the leading and lagging particles and their axial distances Dz plotted versus time for different initial positions x lag and for fixed x lead = 0.2. The bi-dispersed pairs B1 and B2 with different rigidity as well as B3 and B4 consisting of leading and lagging particles with different rigidity La moving in the channel midplane at y = 0 are considered. Inset: The lift-off in the y-direction for the biconcave particle in the B1 pair is shown. Fig. 12 A transition from the deformation-driven tumbling (blue plane) to wobbling (red plane) motion occurs when a biconcave particle approaches the channel center. The particle has rigidity La = 10 and moves in the B1 pair with x lag = 0.2. The distance from the midplane is shown.

Soft Matter Paper
Open Access Article. Published on 19 April 2021. Downloaded on 11/22/2021 9:52:27 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

View Article Online
La = 10 are always stable. They show bounded oscillations in the axial distance due to the tumbling motion of the biconcave particle. For La = 100 (B2) the biconcave shape does not change the qualitative behavior when compared to the mono-dispersed pairs (M2) of the same particle stiffness. Finally, the soft leading (La = 10) and stiff lagging (La = 100) particles in the bi-dispersed pair B3 assemble at the channel center. The pair is stable when the lagging particle starts at x lag = 0 or À0.2, while the remaining starting positions result in unstable pairs. When the stiffness of the leading and lagging particles are interchanged (pair B4), both particles migrate towards the channel center irrespective of the starting position and they ultimately drift apart due to their different shapes. Only for symmetric starting positions about the channel center a stable pair with lateral oscillations occur.
Our work explores and classifies different two-particle motions, which should be observable in the flow of monoand bi-dispersed particle suspensions in the regime of inertial microfluidic. To explicitly test the different types of pair dynamics, one could use microchannels with inlets, which place both particles at specific lateral positions.
Finally, in inertial microfluidics the multi-particle states consist of linear trains either moving on the same streamline or in a staggered configuration as observed in experiments. 96,101 As we demonstrated in ref. 61 for rigid particles, knowing the bounded two-particle states one can directly infer the lattice constant of the staggered trains, while the lattice constant of linear trains on the same streamline typically have a value determined by the range of the two-particle interaction. In this context, the steady-state value of the axial distance between the leading and lagging particles is summarized in Fig. 13 for mono-and bi-dispersed pairs, when they are stable. Thus, the present study will help to understand the multi-particle train configurations occurring in inertial microfluidics for soft particles and their mixtures, which will be more involved compared to the case of ref. 61, where only rigid particles of the same size were studied.

Conflicts of interest
There are no conflicts of interest to declare. Fig. 13 The equilibrium axial distance Dz eq between the leading and lagging particle is plotted versus x lag for stable particle pairs. discretize the channel width. The surface of the solid particles is approximated using 20480 P 1 elements. The temporal evolution of the ratio Dz/a is plotted in Fig. 15, along with its steady-state value reported by Schaaf et al. 60