A flowing pair of particles in inertial microfluidics

A flowing pair of particles in inertial microfluidics gives important insights into understanding and controlling the collective dynamics of particles like cells or droplets in microfluidic devices. They are applied in medical cell analysis and engineering. We study the dynamics of a pair of solid particles flowing through a rectangular microchannel using lattice Boltzmann simulations. We determine the inertial lift force profiles as a function of the two particle positions, their axial distance, and the Reynolds number. Generally, the profiles strongly differ between particles leading and lagging in flow and the lift forces are enhanced due to the presence of a second particle. At small axial distances, they are determined by viscous forces, while inertial forces dominate at large separations. Depending on the initial conditions, the two-particle lift forces in combination with the Poiseuille flow give rise to three types of unbound particle trajectories, called moving-apart, passing, and swapping, and one type of bound trajectories, where the particles perform damped oscillations. The damping rate scales with Reynolds number squared, since inertial forces are responsible for driving the particles to their steady-state positions.


I. INTRODUCTION
The control of hydrodynamic flow fields on a microscopic scale are required in a variety of different applications in medicine, chemistry, and engineering [1]. Microfluidic labon-a-chip devices allow to sample [2] and sort cells [3][4][5], engineer flow patterns [6], and they can be used for fabricating metamaterials [7,8]. While a lot of research has been and is still being done in the field of low-Reynolds-number flows [9][10][11], especially industrial applications need a high throughput in the microchannels [12]. The necessarily increased flow velocities initiated the field of inertial microfluidics. Here, fluid inertia is no longer negligible and new phenomena arise [13]. One prominent example is the Segré-Silberberg effect [14,15], where rigid particles assemble in an annulus, halfway between channel center and wall, when pumped through a cylindrical channel. Its first observation in 1961 inspired many experimental works [16,17] as well as analytical calculations [18,19] and numerical simulations [20][21][22][23][24]. It can be rationalized by a lift-force profile, which a single particle experiences in the channel cross-section [17,25,26] and which can be used to implement an optimal-control scheme [27].
When the density of particles in the channel flow increases, they start to form microfluidic crystals or particle trains [28,29]. Here, the particles assemble in a linear or zig-zag pattern with a fixed axial distance typically ranging from 2.2 to 5 particle diameters [30,31]. A deeper understanding of these particle trains is important for cell analysis [32] and for understanding phonon excitations in microfluidic crystals [33,34]. As particle densities are still small, pair interactions of the particles can provide a first understanding. At vanishing Reynolds numbers pair interactions were studied by Batchelor and Green [9] in an unbounded shear flow. They found open and closed trajectories for a pair of particles. Similar trajectories also occur in Poiseuille flow [11]. Now, including inertia has a profound influence. In particular, * christian.schaaf@tu-berlin.de the flow field around a single particle in a linear shear flow changes noticeably by losing the fore-aft symmetry compared to low Reynolds numbers [35,36]. Applied to the trajectories of a particle pair, Kulkarni and Morris [37] showed that closed trajectories in linear shear flow are replaced by reversing and spiraling trajectories. In microfluidic channels flowing particle pairs, when staying together, perform damped oscillations at finite Reynolds numbers [29,38]. This observation was one motivation for the work reported in this article.
In the following we study the dynamics of a pair of two solid particles driven by Poiseuille flow through a rectangular microchannel. We perform lattice Boltzmann simulations of the Newtonian fluid and couple the particles to the fluid by the immersed boundary method. The lift force profiles of both particles are crucially influenced by their neighbors and strongly depend on their distance along the channel axis. We find strong differences of the profiles for the particles leading and lagging in flow. Furthermore, lift forces in general are larger, which should enhance inertial focusing. Interestingly, how they scale with the Reynolds number depends on the axial particle distance. A linear scaling at close distances reveals interactions determined by viscous forces while the quadratic scaling for larger distances shows the dominating inertial forces. Finally, we categorize the different types of trajectories, on which a particle pair moves, in terms of their initial lateral positions. When the particles stay together, damped oscillatory trajectories occur, which can be explained using the two-particle lift force profiles.
The article is organized as follows. In Section 2 we explain the set-up of our system, describe our implementation of the lattice-Boltzmann method, and how we couple the particles to the fluid. In Section 3 we present the results for the lift force profiles and the particle trajectories. We summarize and close with final remarks in Section 4.

II. METHODS
We study a pair of solid spherical particles moving in a microfluidic channel flow at moderate Reynolds numbers. In the L 2w 2h 0 x z ∆z x lead x lag a Figure 1. Left: A schematic of the microfluidic channel. We use a rectangular channel with width 2w, height 2h, and length L. Only particle motion in the x, z-plane (gray color) is considered. Right: Detailed view in the x, z-plane with Poiseuille flow along the z-axis. Two particles start with axial distance ∆z and lateral coordinates x lead , x lag (measured against the center line) for particles leading and lagging in flow, respectively.
following we shortly explain the microfluidic setup and the lattice-Boltzmann method to simulate the hydrodynamic flow.

A. Microfluidic setup in the simulations
The channel of length L has a rectangular cross section (width 2w and height 2h) with an aspect ratio w/h = 0.5 (see Fig. 1). The channel is filled with a Newtonian fluid with density ρ and kinematic viscosity ν. To drive a Poiseuille flow in the rectangular channel, we apply a constant pressure force [39]. The flow is characterized by the channel Reynolds number Re = 2wu max /ν with the maximal flow velocity u max . It is directed along the z-direction, which we call axial direction, while movements perpendicular to the flow direction are referred to as lateral movements.
In this Poiseuille flow we place two spherical particles with the same radius a and neutral buoyancy. At moderate Reynolds numbers fluid inertia becomes relevant and both particles experience lateral lift forces f lift . They push the particles into the x, z plane containing the shorter cross-sectional axis [23] and ultimately cause inertial focusing onto a specific position. To determine these lift forces in our simulations, we fix the particles' positions on the cross-sectional x axis and measure the forces, which the fluid exerts on them. The lift forces are crucially influenced by the presence of the second particle and we will illustrate how they depend on the axial particle distance ∆z in Sec. III A. The particles flow with different velocity along the channel axis depending on their positions in the channel cross section. So, when we measure the lift-force profiles, we let the particles move with their center-of-mass velocity and keep ∆z constant. This means that we effectively act with an axial force along the flow direction on each particle, resulting in small changes of the lift forces according to the Saffman effect [40], as we demonstrate below. In Sec. III B we analyze trajectories of the particle pair. Here, they can evolve freely without constraints.
Finally, along the flow direction we use periodic boundary conditions. To ensure that the particles do not interact with their mirror images, we use a channel length of L = 30a + ∆z.

B. Lattice-Boltzmann method
To solve the Navier-Stokes equations, we use the lattice Boltzmann method (LBM) in 3D based on 19 different velocities vectors (D3Q19) [41] and rely on the Bhatnagar-Gross-Krook (BGK) collision operator [42]. In the LBM the fluid is modeled by a one-particle probability distribution f i ( x,t), which is determined on a cubic lattice with lattice spacing ∆x. The distribution function depends on the lattice vectors x and the index i stands for the 19 discretized velocity vectors c i pointing to the edges and the faces of a cube and the zero velocity. Now, f i ( x,t) evolves during time ∆t according to two alternating steps: streaming: where f eq i is a second-order expansion of the Maxwell-Boltzmann distribution in the mean velocity and τ is the relaxation time of the BGK model.
Macroscopic quantities like the density ρ and the momentum density ρ u are defined via the zeroth and first moments of the distribution function: Typically, in the LBM the density of the fluid is set to 1. The viscosity is related to the relaxation time τ [43], where c 2 s = 1/3 is the speed of sound measured in LBM units. To ensure incompressibility of the fluid, simulations have to be performed at small Mach numbers Ma = u max /c s . One additional constraint arises from the immersed-boundary method, which we use to implement the two particles (see Sec. II C). It gives best accuracy for relaxation times τ ≤ 1 or ν ≤ 1/6 [44]. In order to vary the Reynolds number, we fixed the viscosity to ν = 1/6 and modified the maximum flow velocity ensuring that Ma < 0.1, which corresponds to density variations of less than 1%.
The channel flow was driven by a constant body force according to the Guo-force scheme [45] and we used regularized boundary conditions at the walls [46]. The lattice-Boltzmann simulations were performed with the code provided by the Palabos project [47], which we supplemented by the implementation of particles using an immersed boundary method. Finally, we discretized the width of the channel along the x axis by 75 lattice cells.
We use the same simulation code as in our previous publication [23]. We only added the event-based Euler method to prevent overlap between particles, which we describe at the end of the following section.

C. Modeling of solid particles
We modeled the colloids by an immersed boundary and couple them to the fluid by the method proposed by Inamuro [48]. This immersed boundary method ensures the no-slip boundary condition at the particle surface by iteratively refining the body force g i acting on the surrounding fluid nodes. To interpolate the fluid velocity at a position between the nodes, we follow Peskin [49] and use an interpolation kernel, which considers all neighboring nodes in a sphere with radius 2∆x. For further details on these methods we refer the reader to the original publications and our previous work [23].
Furthermore, the approach assumes that the particles are filled with a Newtonian fluid which is unphysical. In order to compensate for this, we also apply Feng's rigid body approximation [50] and add an additional force acting on the particle so that it moves like a solid particle. With all contributions the equations of motion for the colloids are given by where i is the particle index, r, v, and ω are, respectively, the position, velocity, and angular velocity. Finally, M and I stand for the mass and moment of inertia. For some trajectories, which we show in Sec. III B, the two particles touch each other. To avoid overlap, we use an eventbased Euler step for the colloids. When the particles are so close that they would overlap in the next time step, we reduce the time step to ∆t so that the particles just touch, perform the collision between the particles, and finish the remaining time step with length ∆t − ∆t with the new values for the particle velocities and angular velocities. To realize the collision, we follow Ref. [51] and consider two rough hard spheres [51] so that during collision also angular momentum is exchanged.

A. Lift force profiles
The dynamics of a particle pair in inertial microfluidics is best captured by the particle-particle lift force profiles. They quantify the lift forces, the particles experience in the presence of the other particle either leading or lagging in flow [see Fig.  1 for the geometry]. Zero forces correspond to fixed points or equilibrium positions in the channel cross section and the magnitude of the force indicates how fast the particles become focused on their equilibrium positions. When both particles are mirrored at the channel axis the lift force reverses sign. More importantly, when the flow direction is reversed such that the leading particle becomes lagging and vice versa, the lift force profiles change since due to secondary flow in the inertial regime the leading and lagging particles experience different flow fields.
In previous work we already analyzed the lift-force profile for a single rigid particle [20,23]. This force profile changes sign when mirrored at the channel center and scales with the Reynolds number ∝ Re 2 (Fig. 2). Typically, one finds an unstable fixed point in the channel center and stable off-centered fixed points or equilibrium positions (indicated by a negative slope) along symmetry axes in the channel cross section. In this work we will focus on the two-particle lift-force profiles as a function of the axial particle distance ∆z (measured along the flow direction), the lateral coordinates x lead , x lag and the channel Reynolds number Re as indicated in Fig. 1. In what follows we concentrate on particles with radius a/w = 0.4.

Parameter study
In Fig. 3 we demonstrate how the presence of another particle influences the lift force profiles and the equilibrium positions. We keep the axial distance of the particle pair fixed at ∆z = 3a and plot force profiles of the leading particle for different lateral positions x lag of the lagging particle [ Fig. 3a)] and vice versa [ Fig. 3b)]. Overall one recognizes that the profile is drastically influenced by an adjacent particle and lift forces generally are larger compared to the single-particle case. Thus, inertial focusing is enhanced.
For the leading particle [see Fig. 3(a)] we only find stable fixed points in the channel side opposite to the location of the lagging particle, the other fixed points have disappeared. However, the new equilibrium positions are closer to the channel center compared to the single-particle case (black line) and, ultimately, for x lag = 0 the stable fixed point is in the channel center. In contrast, when the leading particle resides in the upper half of the channel, the fixed point of the lagging particle in the other channel side [ Fig. 3b)] becomes unstable and stable equilibrium points only exist very close to the upper channel wall for sufficiently large x lead . Interestingly, the configuration with x lead = x lag = 0 is not stable. Finally, note that the lift-force profiles of the leading and lagging particles differ from each other due to secondary flow as stated in the beginning. Figure 3. Lateral lift force profiles along the short axis for the leading (a) and lagging (b) particle for Re = 5, axial distance ∆z/a = 3, and particle radius a/w = 0.4. The curve parameters are the positions of the lagging (a) or leading (b) particle, respectively. The black dotted line corresponds to the single-particle force profile.
When we increase the distance ∆z of the two particles along the flow direction, the lift-force profiles are more similar in shape to the profile of a single particle, however shifted downwards (leading particle) or upwards (lagging particle) (Fig. 4). In particular, for the cases where the other particle is relatively close the the channel center (x/w < 0.3), two stable equilibrium positions in the two respective sides of the channel are still present. In addition, when the neighboring particle is in the channel center (x/w = 0), the force profile agrees with the single-particle case close to the channel center (blue lines in Fig. 4) but the stable equilibrium positions are located closer to the channel walls. Finally, by increasing the axial distance between the particle pair, the strength of the lift forces decreases compared to Fig. 3.
The existence of stable fixed points in the lateral force profiles of both particles does not necessarily define a stable particle configuration, since particles closer to the channel center move faster than particles near the channel walls. For a stable pair configuration the fixed points of both particles have to be at the same distance from the channel center. From Fig. 4 we observe that this might be possible for x lag /w = −x lead /w Figure 4. Lateral lift force profiles along the short axis for the leading (a) and lagging (b) particle for Re = 5, particle radius a/w = 0.4, and at the larger axial distance ∆z/a = 5 compared to Fig. 3. The curve parameters are the positions of the lagging (a) or leading (b) particle, respectively. The black dotted lines correspond to the single-particle force profile. Figure 5. Lateral forces for a particle pair as a function of the axial distance ∆z with the leading particle at x lead /w = −0.4 and the lagging particle at x lag /w = 0.2. The dotted lines show the single particle lift forces at x lead and x lag , respectively. around 0.4, which we will indeed confirm further below in Sec. III B. In contrast, when the particles are close together, for example at ∆z/a = 3 as in Fig. 3, such a stable pair configuration is not possible.
In Fig. 5 we fix the lateral positions of the leading and lagging particles and show how the lift forces change with the axial distance ∆z. Below ∆z/a ≈ 3 the lift forces for both particles are positive so that they are pushed in positive x direction towards the upper channel wall. Thus the leading particle at x lead /w = −0.4 moves closer to the channel center while the lagging particle at x lag /w = 0.2 approaches the upper wall. So, the distance of the particles grows and effectively they are repelled from each other. For an axial distance ∆z/a > 3 the lift force on the lagging particle at x lag /w = 0.2 becomes negative. Thus, both particles are pushed together, which can be described as an effective attraction. Increasing ∆z further initiates further sign changes. Finally, for long distances both particles approach the lift forces of the respective single particle (dotted line). The small offset between the pair force and the single particle force can be explained by the Saffman effect. To measure the lift force, we have to fix the axial distance of these particles. In doing so we effectively accelerate the lagging particle and decelerate the leading particle. This additional axial force along the flow direction leads to the small contribution in the lateral direction and explains the offset.
Further, we study how a variation of the Reynolds number influences the two-particle force profiles. Again, we plot them at an axial distance ∆z/a = 3 as in Fig. 3 but now for Re = 20 instead of Re = 5 (see Fig. 6). We immediately recognize that in contrast to Fig. 3 the force profiles are similar in shape to the one-particle profile but shifted upwards (lagging particle) or downwards (leading particle) with increasing lateral distance. We saw a similar behavior already in Fig. 4 for Re = 5 at larger axial distance ∆z/a = 5. In both cases the strength of the lift forces are similar to the single-particle forces, while in Fig. 3 the two-particle induced forces are considerably larger than the inertial forces on a single particle. In addition, the lift forces in Fig. 6 rescaled by ρν 2 Re 2 are smaller than in Fig. 3, which suggest that the usual scaling with Re 2 does not apply. We study this in more detail in the next paragraph.

Scaling of the lift force with Re
We now take a closer look on how the lift force scales with the Reynolds number Re. We already realized that for small particle distances the two-particle lift forces no longer scale with Re 2 as in the single-particle case. However, it is also clear that for large distances this scaling has to be recovered since the influence of the two particles on each other strongly decreases. To analyze this aspect in more detail, we fix the leading particle at x lead /w = 0.3 and vary x lag . We determine the maximum value of the magnitude of the lift force profile for the lagging particle and plot it versus Reynolds number for several particle distances [52]. Figure 7 shows the results in double-logarithmic scale. One clearly recognizes a power-law scaling with exponent α: f max lift ∝ Re α . In the inset, we plot α versus ∆z for both the leading and lagging particle. Indeed, Figure 6. Lateral lift force profiles along the short axis for the leading (a) and lagging (b) particle for Re = 20, ∆z/a = 3, and a/w = 0.4. The curve parameters are the positions of the lagging (a) or leading (b) particle, respectively. The black dotted lines correspond to the single-particle force profile. we find α = 2 for ∆z/a > 7. When the particles approach each other, the scaling exponent decreases to almost α = 1 for ∆z = 3a. This scaling helps to further understand the character of the lift force, in particular, when two particles interact. A particle disturbs the fluid flow, which then influences the motion of nearby particles through a viscous coupling. This is the dominant contribution to the lift force at small distances as indicated by the linear scaling of the lift force with Re. The inertial contribution takes over at large distances, where the disturbance flow from the neighboring particle is weak, and one recovers the typical scaling for the inertial lift force, f lift ∝ Re 2 . So, our analysis confirms the picture of Ref. [29], which explicitly speaks about a viscous disturbance flow.
So, both Figs. 5 and 7 indicate that beyond the distance ∆z/a ≈ 7 the particles essentially do not interact. In Ref. [53] it is argued that hydrodynamic interactions in a microchannel are screened on distances larger than the width of the channel cross section. In our case taking a particle radius of a = 0.4w, a distance of 7a corresponds to 2.8w, which is close to the channel width of 2w. This explains our observation. Figure 7. Maximum value of the lift force of the lagging particle plotted versus Re for different axial distances ∆z. The leading particle is fixed at x lead /w = 0.3 while x lag is varied. The particle radius is a/w = 0.4. The dashed lines indicate scaling laws ∝ Re and ∝ Re 2 . Inset: Scaling exponent α from f max lift ∝ Re α plotted versus ∆z for both particles. Always the leading particle is fixed at x lead /w = 0.3, while the position of the lagging particle is varied. Figure 8. Color-coded lift force profiles in a two-dimensional representation plotted versus x lead and x lag for the leading (left) and lagging (right) particle for Re = 5 and ∆z/a = 4.

Contour plots
In Sec. III B we analyze possible trajectories for a pair of solid particles moving under the influence of the lateral lift forces. To rationalize these trajectories, it is instructive to use a two-dimensional representation of the respective liftforce profiles of the leading and lagging particles (see Fig. 8). Again, we clearly recognize the asymmetry of the profiles between the leading and the lagging particle, while each profile is symmetric under reflection at the channel center. The white lines indicate zero crossings of the lift force, so stable and unstable equilibrium points.

B. Two-particle trajectories
We now present the possible trajectories, which two particles traverse as a result of the coupled lift-force profiles presented in Sec. III A and advection in the Poiseuille flow. The different types occur depending on the starting lateral positions and the axial distance. Thus, in Fig. 9 we categorize them in a diagram for the starting lateral positions x lag and x lead , while keeping the starting axial distance fixed. We identified four different kinds of coupled particle movements, which we term moving apart, passing, swapping, and damped oscillations.
The first three types of trajectories we name unbound as their particles drift apart and reach their equilibrium lateral positions at large axial distances, where they do not influence each other anymore. We will analyze these trajectories in more detail further below. We also observe bound trajectories, where the two particles ultimately perform damped oscillations about their equilibrium lateral positions. They occur in the narrow red region in Fig. 9, where the particles occupy opposing channel sides with the lagging particle only little faster than the leading such that they can stay together. We start with describing the bound trajectories.  ber [11], here the amplitudes decrease in time indicating that damping is an inertial effect. Such a damped motion is not possible in Stokes flow as it would violate kinetic reversibility of the Stokes equations. Interestingly, the damped oscillatory two-particle trajectories were also observed in experiments by Lee et al. [29]. However, the authors report that the particles move apart symmetrically ("symmetric repulsive interactions"), while there is an asymmetry when approaching each other ("asymmetric attractive interactions"), which we do not observe [compare Fig. 10 b)]. Maximum and minimum displacements are always in phase. Ultimately, the particles reach their final lateral equilibrium positions, which agree with the positions of single particles.
In the following we analyze how oscillation frequency Ω and damping rate γ behave as a function of the Reynolds number (see Fig. 11). We determined Ω by measuring the time between maximal displacements and γ by an exponential fit for the amplitudes decaying in time. For the oscillation frequency Ω (inset of Fig. 11) we find a linear scaling with the Reynolds number, which indicates that the oscillations are due to the viscous coupling between the particles. In contrast, the damping rate scales quadratically with the Reynolds number since inertial lift forces drive them to their equilibrium positions. According to Fig. 7 these inertial forces act here as a pertubation. Note that our findings on the damped oscillations are in full agreement with Ref. [54]. Figure 11. Damping rate γ of the particle-particle distance and oscillation frequency Ω (inset) plotted versus Re. The initial conditions are x lag /w = −0.2, x lead /w = 0.24, and ∆z 0 /a = 5. Linear and quadratic fits in Re are indicated, respectively.
The dynamics of the oscillating particle pair, which we discussed in Fig. 10, can be nicely illustrated using lift-force contour plots similar to the one we determined in Sec. III A (see Fig. 8 and video in SI) but now for Re = 10. We start with the initial conditions x lag /w = −0.2, x lead /w = 0.24, and ∆z/a = 5 .
We can now follow the particles in the lift-force contour plots in video to understand their trajectories in the x, z plane. In the beginning the lagging particle is faster as it is closer to the center. The signs of both lift forces are such that they push the particles towards the walls. In this phase, since the lagging particle is faster the axial distance decreases. The leading particle turns around and moves away from the upper wall, while the lagging particle still moves towards the lower wall. Thus, both forces are negative. Ultimately, the lagging particle is closer to the wall than the leading particle. It moves slower and the axial distance increases. By following the trajectories further, the signs of the lift forces always indicate the lateral direction of the moving particles. In the end they show damped oscillations about the zero lines of the contour plots in agreement with the spiraling motion in the x, z plane. Finally, after a few oscillations the particles reach their stable equilibrium positions, where the lift forces are zero.
Interestingly, we find that all bound particle pairs performing damped oscillations assemble at an axial distance of ∆z/a ≈ 4.1 independent of their initial conditions or the Reynolds number, which we varied between 2 and 20. The value of this axial distance is in good agreement with experimental and theoretical results [28,30,31,54]. The scaling of the lift force with Re (cf. Fig. 7) indicates that at this equilibrium distance the particle interactions are dominated by a viscous disturbance flow as already mentioned above and in Refs. [29,54]. The shape of this flow does not depend on the Reynolds number, which explains why the equilibrium distance is independent of Re. In this article we explain the formation of stable particle pairs by the fact that the lift forces acting on them become zero and that both particles have the same distance from the channel center so that they drift with the same velocity. An alternative and intuitive explanation for the formation of crossstreamline pairs is given in Ref. [30]. A particle creates a viscous disturbance flow, which contains eddies or vortices on the opposite side of the channel as indicated by the streamlines in the co-moving frame in Fig. 12. The second particle then occupies the center of an eddy, where it does not move relative to the first particle. Since the viscous disturbance flow is independent of Re, the position of the eddy does not change with Re in agreement with our argument in the previous paragraph. The advantage of the lift force profiles introduced in this article is that they not only describe equilibrium positions but that they also determine the full dynamics of a particle pair. We demonstrated this before when describing the oscillatory motion of the bound particle pair.
In addition to the cross-streamline pairs Hood and Roper [54] also formulated a theory that predicts stable samestreamline pairs. However, in our analysis particles moving on the same side of the channel never form bound states. This result implies that two-particle interactions are not sufficient to explain particle lattices, where all particle assemble on one side of the channel.

Unbound trajectories
We already introduced the unbound trajectories, where no stable pair configuration exists (green,blue, and organge patches in Fig. 9). In the moving-apart trajectories (green patch with |x lead | < |x lag |) the leading particle is faster than the lagging one. The distance between both simply increases and they independently migrate towards their equilibrium positions due to inertial focusing as given by the single-particle lift force profile. Interestingly, even when the particles start on their single-particle equilibrium positions in the same half of the channel, they do not keep their distance fixed but move on a moving-apart trajectory (green line in Fig. 13). The rea- son is that these single-particle fixed points are not stable at a small particle distance, which also shows again the asymmetry in the force profiles. The lift forces push the leading particle closer to the center and the lagging particle closer to the walls. This causes a non-zero relative velocity and the particles move apart ending at a larger axial distance. Finally, in the narrow green stripe with x lead > x lag (see Fig. 9) the particles initially approach each other. However, the lagging particle also drifts towards the wall so that it becomes slower than the leading particle and they both just move apart.
Passing trajectories occur since the lagging particle is closer to the channel center and therefore faster than the leading particle. So, they change the order in axial direction (blue line in Fig. 13). During this overtaking the displacement of the two particles is asymmetric. The particle closer to the channel center is displaced much stronger and the offset is clearly visible after the passing event. Then the axial distance increases and the particles assume their single-particle positions due to inertial focusing. We note that the particles do touch during the overtaking as we did not implement any lubrication approximation. However, due to our event-based Euler step they do not overlap.
When moving on swapping trajectories (orange line in Fig. 13), the faster lagging particle does no succeed to overtake the leading particle. Instead, the particles come close to each other and then swap the lateral position, which makes the leading particle the faster particle so that they keep their axial order. When the particles move apart, they have interchanged their distances to the channel center. For example, for the orange trajectories in Fig. 13 one finds x after lead ≈ −x before lag and vice versa. Note similar trajectories in linear shear flow were called reversing trajectory [37] In general we see very similar types of trajectories also at low Reynolds numbers [11] indicating that they are governed by the viscous particle coupling and the Poiseuille flow profile. Inertial forces are responsible for focusing the particles on positions determined either by the two-particle or single-particle lift force profiles. Although the two-particle trajectories studied in this article are often unbound, in preliminary results we find that they are the fundamental building blocks in the formation of multi-particle lattices. Our goal is to explain the formation of particle lattices using also these trajectories in a future work.

IV. CONCLUSION
Understanding pair interactions of two particles in inertial microfluidics is an important step for understanding collective dynamics such as the formation of particle trains.
In this work we studied the lift force profiles and the trajectories of a pair of two solid particles driven by Poiseuille flow through a rectangular microchannel. The lift force profiles of both particles are strongly influenced by their neighbors and depend on the particle distance along the channel axis. They clearly differ between the leading and lagging particles and the lift forces are stronger compared to a single particle. The increased lift force should enhance particle focusing by driving them faster towards their equilibrium positions. At close distance the lift force profiles differ strongly from the profile of a single particle and do not allow for stable pair configurations. However, when increasing the axial distance or the channel Reynolds number, the profiles appear similar in shape but are shifted by constant forces. Interestingly, at small axial distances below ∆z/a = 4 the strength of the lift forces scales with Re indicating that hydrodynamic interactions between the particles are dominated by viscous forces, while for distances ∆z/a = 10 and larger scaling is quadratic in Re showing the importance of inertial forces. In between, the scaling follows Re α with the exponent varying smoothly from 1 to 2 while increasing ∆z. Finally, we presented the lift force profiles of leading and lagging particles in a two-dimensional representation as a function of both lateral particle positions.
These two-dimensional plots determine the coupled dy-namics and the trajectories of two floating particles. We identified four types of particle trajectories depending on the initial lateral position of the leading and lagging particles. Three of them are unbound, where the particle distance ultimately increases until both particles reach their single-particle equilibrium positions. In the moving-apart trajectories the leading particle is mostly faster than the lagging particle and the pair drifts apart. If the lagging particle is much faster, it overtakes and thereby changes axial order with the leading particle in what we call passing trajectories. If the lagging particle is not much faster, it only approaches the leading particle but then they exchange their lateral positions and move apart again. Thus, they move on swapping trajectories. Finally, bound trajectories occur for x lag ≈ −x lead , where axial distance and lateral positions of the particles perform damped oscillations while reaching their equilibrium values. As such a damping does not occur in Stokes flow, it is clearly an inertial effect. Consequently, the damping rate scales with the Reynolds number squared, while the oscillation frequency increases linearly in Re. Interestingly, for the specific particle radius studied here all oscillating trajectories ultimately end at an axial distance of ∆z/a ≈ 4.1 independent of the initial conditions and the Reynolds number. With our investigations we hope to shed further light on the collective dynamics and ordering of particles flowing through microfluidic channels at moderate Reynolds numbers. It should be useful in designing microfluidic crystal structures as well as developing and improving particle separation and sorting techniques in inertial microfluidics.