3.1. Lateral forces at low and high Reynolds numbers

As representative examples, the lateral forces exerted on a particle are plotted in figure 3(a,b) when the particle is placed on the $y$ -axis or on the diagonal, respectively, for $Re=144$ and $514$ . In figure 3(b), $r\;(=\sqrt{y^{2}+z^{2}})$ represents the distance from the channel centre and $F_{r}$ denotes the lateral force along the diagonal. In both figures, the outward forces toward the channel walls act on the particle when it is located close to the channel centre, whereas significant wall repulsion forces act on the particle when it is placed near the channel walls. In all cases computed, one equilibrium position exists at a location between the channel centre and the channel wall, which is determined as the intersection of the force curve with the horizontal axis. The intersection on the $y$ -axis, shown in figure 3(a), represents the channel face equilibrium position, and the intersection on the diagonal, shown in figure 3(b), represents the channel corner equilibrium position.

Miura et al. (Reference Miura, Itano and Sugihara-Seki2014, Reference Miura, Harada, Itano and Sugihara-Seki2015) observed the distributions of particle centres at downstream channel cross-sections, and determined the locations of equilibrium positions as local maximum positions of the probability density function of the particles. At $Re=144$ , the equilibrium positions thus determined are present only at the centre of the channel faces, while at $Re=514$ , both the channel face and channel corner equilibrium positions are present. In figure 3(a,b), experimentally obtained equilibrium positions are indicated by an asterisk for $Re=144$ and a cross for $Re=514$ , respectively. The locations of these equilibrium positions showed excellent agreement with the present numerical results, although there are some deviations in the experimental measurements (see figure 8 a,b). At $Re=144$ , the present numerical simulation predicts the presence of channel corner equilibrium positions, whereas they are not observed experimentally. We think that the channel corner equilibrium positions are unstable for $Re$ smaller than a certain critical value, $Re_{c}$ , and so they cannot be observed experimentally for $Re<Re_{c}$ . This is discussed below.

Figure 4.

(a) Map of the lateral forces $(F_{y},F_{z})$ and trajectories of the particle centre projected over the channel cross-section at $Re=144$ . Each arrow at the lattice represents the direction of the lateral force, and its magnitude is indicated by the colour. The triangle represents the channel face equilibrium position, which is stable, and the square shows the channel corner equilibrium position, which is unstable. The thick red line represents the heteroclinic orbit that joins the channel face equilibrium position and the corner equilibrium position. (b) Particle distribution obtained experimentally at $Re=144$ . (Miura et al. Reference Miura, Itano and Sugihara-Seki2014) and the heteroclinic orbit obtained numerically. Each dot represents the location of a particle centre observed in the cross-section.

Figure 4(a,b) show a map of the lateral forces $(F_{y},F_{z})$ in the first quadrant of the channel cross-section at $Re=144\;({<}Re_{c})$ , and the corresponding experimental particle distribution results (Miura et al. Reference Miura, Itano and Sugihara-Seki2014). In figure 4(a), each arrow represents the direction of the lateral force acting on a particle located at the origin of the arrow, with the magnitude expressed by the colour. Several representative trajectories of the particle centre projected over the cross-section are also plotted with black lines. In figure 4(b), each dot represents the location of a particle centre observed in the cross-section; the data are summed over the region of the azimuthal angle between $0$ and ${\rm\pi}/4$ based on symmetry considerations.

In figure 4(a), the red square represents the channel corner equilibrium position on the diagonal, the location of which is determined from the force curve shown in figure 3(b). Trajectories starting from the neighbourhood of the corner equilibrium position lead away from this location, indicating that the corner equilibrium position is unstable. In addition, trajectories starting from various initial positions approach the channel face equilibrium position (red triangle in figure 4 a), indicating that the channel face equilibrium position is stable. This is the reason why no particles are observed near the channel corners and most particles are placed near the centre of the channel face in figure 4(b).

From the trajectories illustrated in figure 4(a), we notice that the particles migrate in two stages. In the first stage, they move in nearly radial directions until they reach a single line that represents a heteroclinic orbit joining the channel face equilibrium position and the corner equilibrium position (red line in figure 4 a). In the second stage, the particles travel along the heteroclinic orbit toward the channel face equilibrium position. In accordance with this prediction, figure 4(b) shows that the particles observed experimentally are almost aligned along the heteroclinic orbit, with most particles focused near the channel face equilibrium position. The first stage is a relatively fast process, while the second stage is slower. We will show in § 3.3 that the time scale, or the travelling distance in the flow direction during each stage, depends on not only the particle position but also the Reynolds number.

The present results shown in figure 4(a) at $Re=144$ are consistent with previous studies for microchannel flows for $1\leqslant Re\leqslant 120$ . For the same size ratio with the present study ( $d/D=0.11$ ), Abbas et al. (Reference Abbas, Magaud, Gao and Geoffroy2014) made microscope measurements of particles suspended in square microchannel flows and numerical analyses of the trajectories of a single particle based on the force coupling method. The trajectories at $Re=120$ shown in figure 8 of their paper are comparable to those shown in figure 4(a), although the locations of the channel face equilibrium positions obtained in the present study are closer to their experimental results rather than their numerical results (see our figure 8 a).

At lower $Re\;(=40)$ and for larger particles $(d/D=0.22)$ , the distribution of the lateral forces in the cross-section was calculated by a finite element computation (Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009) and by a perturbation analysis (Hood, Lee & Roper Reference Hood, Lee and Roper2015). While not shown here, we have also plotted the map of the lateral forces and computed the trajectories for $Re=50$ and $d/D=0.11$ (Nakagawa et al. Reference Nakagawa, Miura, Otamo, Makino and Sugihara-Seki2014), which were found to have similar features to figure 4(a); outward forces are exerted on the particle in the inner region of the heteroclinic orbit, whereas in the outer region significant repulsion forces act, with the magnitude being larger nearer to the centre of the channel face and weaker nearer to the corner. There is also a stable equilibrium position at the centre of each channel face, and an unstable equilibrium position at each corner. These features are also observed in figure 1 of Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009), despite the difference in size ratio $d/D$ . Hood et al. (Reference Hood, Lee and Roper2015) estimated the lateral forces using a regular perturbation method from the Stokes flow regime. Although the trajectories and the lateral force distribution shown in figure 8 of their paper are qualitatively similar to the present results, there seem some discrepancies in the region near the channel wall and around the channel corner equilibrium positions; the lateral forces in that region are rather weak compared to Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009) and the present study. Large repulsion forces on the particles near the channel wall were reported experimentally by Choi et al. (Reference Choi, Seo and Lee2011) in the range of $2.4\leqslant Re\leqslant 60$ . Since the effect of the channel wall was calculated by the method of reflections in Hood et al. (Reference Hood, Lee and Roper2015), a finite number of reflections may lead to the rather weak repulsion force when a large particle is placed close to the channel wall.

Figure 5.

(a) Map of the lateral forces $(F_{y},F_{z})$ and trajectories of the particle centre projected over the channel cross-section at $Re=514$ . At this $Re$ , both the channel face equilibrium position (triangle) and the channel corner equilibrium position (square) are stable. The thick red solid line represents the heteroclinic orbit and the thick red dashed line represents the separatrix. (b) Particle distribution obtained experimentally at $Re=514$ (Miura et al. Reference Miura, Harada, Itano and Sugihara-Seki2015) and the heteroclinic orbit obtained numerically. For additional explanation, see figure 4.

The map of the lateral forces computed at $Re=514\;({>}Re_{c})$ and the corresponding experimental results of the particle distribution (Miura et al. Reference Miura, Harada, Itano and Sugihara-Seki2015) are shown in figure 5(a,b), respectively. Several trajectories of the particle centre are plotted with black lines in figure 5(a). The trajectories indicate that at this Reynolds number, both the channel face and channel corner equilibrium positions are stable, and separatrix (red dashed line) as well as a heteroclinic orbit (red solid line) appear in figure 5(a). The intersection of these two lines represents a saddle point. The trajectories in figure 5(a) also suggest a two stage migration of the particles: migration up to the heteroclinic orbit in the first stage and migration along the heteroclinic orbit toward either the channel face or corner equilibrium position in the second stage. This property is exemplified in figure 5(b), which shows the alignments of particles along the heteroclinic orbit near the channel face and channel corner equilibrium positions.

To the authors’ knowledge, there are no comparable studies investigating the lateral forces on a particle suspended in square channel flows at this $Re$ , except Chun & Ladd (Reference Chun and Ladd2006). They reported from their numerical simulation, using the lattice Boltzmann method, that only the channel corner equilibrium positions are stable and the channel face equilibrium positions are unstable above $Re=500$ . This result does not match our experiments (Miura et al. Reference Miura, Itano and Sugihara-Seki2014, Reference Miura, Harada, Itano and Sugihara-Seki2015) and the present numerical simulation shown in figure 5(a). For low $Re\;({\approx}100)$ , Chun & Ladd (Reference Chun and Ladd2006) found stable equilibrium positions in corners and faces while in the following numerical studies (Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Abbas et al. Reference Abbas, Magaud, Gao and Geoffroy2014), including the present study, and in experimental studies (Kim & Yoo Reference Kim and Yoo2008; Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Choi et al. Reference Choi, Seo and Lee2011; Abbas et al. Reference Abbas, Magaud, Gao and Geoffroy2014; Miura et al. Reference Miura, Itano and Sugihara-Seki2014), the corner corner equilibrium positions were indicated to be unstable. With regard to the discrepancy at low $Re\;({\approx}100)$ , Abbas et al. (Reference Abbas, Magaud, Gao and Geoffroy2014) reported from their numerical simulations, adopting a very short computational domain, that the particle can be pushed towards the channel corner instead of the centre of the channel face at $Re=120$ , suggesting that the inter-particle interaction in the flow direction may result in the presence of the channel corner equilibrium position at low $Re$ . Chun & Ladd (Reference Chun and Ladd2006) employed a periodic channel length of $l/d=5$ for their simulations at $Re=100{-}1000$ . Thus, the discrepancy observed at high $Re\;({\approx}500)$ may also originate from the hydrodynamic interaction between particles, since the wake of a particle could be larger and longer as $Re$ increases. We will examine the effect of the channel length on the particle behaviour at high $Re$ in a future study.

Figure 6.

(a) Map of the lateral forces $(F_{y},F_{z})$ and trajectories of the particle centre projected over the channel cross-section at $Re=260$ . At this $Re$ , an additional stable equilibrium position (red circle) appears on the heteroclinic orbit (thick red solid line). (b) Particle distributions obtained experimentally at $Re=260$ (Miura et al. Reference Miura, Itano and Sugihara-Seki2014) and the heteroclinic orbit obtained numerically. For additional explanation, see figure 4.

3.2. Lateral forces near the critical Reynolds number

The present numerical study showed that both the channel face and channel corner equilibrium positions always exist in the range of $Re$ computed, but in fact the channel corner equilibrium positions are unstable until the Reynolds number exceeds a certain critical value, $Re_{c}\;({\approx}260)$ . Figure 6(a,b) show a map of the lateral forces computed at $Re=260$ and the corresponding experimental particle distribution (Miura et al. Reference Miura, Itano and Sugihara-Seki2014), respectively. Near this $Re$ , it was shown experimentally that additional equilibrium positions appear at azimuthal angles of ${\sim}0.5{-}0.7$  rad, as well as the channel face equilibrium positions (Miura et al. Reference Miura, Itano and Sugihara-Seki2014, Reference Miura, Harada, Itano and Sugihara-Seki2015). Corresponding to this observation, the present computation showed the presence of another stable equilibrium position on the heteroclinic orbit at $Re=260$ , which is indicated by a red circle in figure 6(a). This location, of azimuthal angle $0.59$  rad, is close to the equilibrium position observed experimentally. One of the characteristics at this $Re$ is that the magnitudes of the lateral forces are very small along the heteroclinic orbit. This property may be demonstrated in figure 6(a), in which the trajectories (black lines) enter the heteroclinic orbit (red line) at almost right angles. The particles in the second stage move very slowly along the heteroclinic orbit, taking a long time to reach the stable equilibrium position (see figure 7). The magnitudes of the lateral forces are so small that other equilibrium positions may be present on the orbit. The slow movement in the second stage may be reflected in the experimental observation that the particles shown in figure 6(b) are aligned along almost entire circumference of the heteroclinic orbit and the focusing near the equilibrium positions is weak.

In order to examine the stability of the channel corner equilibrium position from the lateral force profile, we computed the azimuthal component of the lateral force, $F_{{\it\theta}}$ , along a short circular arc passing the corner equilibrium position at $Re=144$ , 260 and 514. The computations in the range of ${\it\theta}$ from ${\rm\pi}/4$ to ${\rm\pi}/4-0.025$  rad indicated that, as ${\it\theta}$ decreases from ${\rm\pi}/4$ , $F_{{\it\theta}}$ decreases monotonically from zero at $Re=144$ , whereas $F_{{\it\theta}}$ increases monotonically from zero at $Re=514$ . Since the radial component of the lateral force, $F_{r}$ , vanishes at the corner equilibrium position, varying from positive to negative values with $r$ along the diagonal, this result implies that the corner equilibrium position is a saddle point at $Re=144$ and it is a stable equilibrium point at $Re=514$ . The magnitudes of $F_{{\it\theta}}/{\it\rho}U^{2}d^{2}$ at ${\it\theta}={\rm\pi}/4-0.025$  rad are about $2\times 10^{-4}$ at $Re=144$ and 514. On the other hand, at $Re=260$ the maximum value of $F_{{\it\theta}}/{\it\rho}U^{2}d^{2}$ in this range of ${\it\theta}$ is about one-tenth of this value, i.e.  $F_{{\it\theta}}$ is extremely small, indicating that this Reynolds number is close to the critical value, $Re_{c}$ . The similar computations at $Re=300$ also yielded small values of $F_{{\it\theta}}$ in the same range of ${\it\theta}$ . As such subtle variations in the lateral forces may be beyond the accuracy of the current numerical method, further analyses of higher precision are necessary to determine the critical Reynolds number.

Figure 7.

Particle trajectories in (a) the $(x,y)$ plane and (b) the $(x,z)$ plane, starting from $(y/(D/2),z/(D/2))=(0.5,0.375)$ , at $Re=144$ (solid line), 260 (dashed line) and 514 (dash-dotted line). The horizontal dotted lines denote the asymptotic coordinates the trajectories approach. The dashed thin lines represent the trajectories, starting from $(y/(D/2),z/(D/2))=(0.2,0.1)$ , at $Re=260$ . Note that the $x$ -scale for this case is shown on the upper horizontal axis.

3.3. Entry length for focusing

The trajectories shown in figures 4(a)–6(a) clearly demonstrate the two stage properties of the particle migration. As a representative example of the particle behaviour at each $Re$ , the projections of the trajectories onto the $(x,y)$ and $(x,z)$ planes, starting from $(y/(D/2),z/(D/2))=(0.5,0.375)$ in the cross-section, are plotted at $Re=144,260$ and 514 in figure 7. The horizontal dotted lines represent the coordinates of the asymptotic position where the trajectory approaches: the channel face equilibrium position at $Re=144$ (triangle in figure 4 a), the stable equilibrium position on the heteroclinic orbit at $Re=260$ (circle in figure 6 a) and the channel corner equilibrium position at $Re=514$ (square in figure 5 a). The first stage continues until the particle reaches the heteroclinic orbit. The travelling distance in the flow direction during the first stage can be estimated from figure 7 as $L_{1}/D\approx 50,35$ and 30 at $Re=144,260$ and 514, respectively. The value of $L_{1}/D\approx 50$ at $Re=144$ agrees with the result for $Re=120$ reported by a numerical simulation of Abbas et al. (Reference Abbas, Magaud, Gao and Geoffroy2014), although the initial position of the trajectory is different.

After the first stage, the particle translates much slower along the heteroclinic orbit towards the stable equilibrium position in the second stage. Figure 7 shows that the second stage continues until the travelling distance $L_{2}/D\approx 600$ and 100 from the channel inlet, at $Re=144$ and 514, respectively, indicating that $L_{2}$ is about ten times and several times larger than $L_{1}$ at $Re=144$ and 514, respectively. These values of $L_{2}$ were estimated from the axial distance where the centre of particle approaches within $0.005\times (D/2)$ to the asymptotic position. The value of $L_{2}/D\approx 600$ at $Re=144$ is about a half of the result of Abbas et al. (Reference Abbas, Magaud, Gao and Geoffroy2014) for $Re=120$ , which is presumably due to the difference in the initial position as well as the Reynolds number, or the method of estimate for $L_{2}$ . They reported that $L_{2}/D\approx 1000$ at both $Re=60$ and $120$ , and figure 9 of their paper showed that, at $Re=12$ , the particle did not approach within a close distance to the asymptotic position even at $x/D=3000$ , indicating that $L_{2}/D>3000$ at $Re=12$ . The present study indicated that the axial distance of the second stage is, in general, longer for smaller $Re$ , which seems consistent with the results of Abbas et al. (Reference Abbas, Magaud, Gao and Geoffroy2014).

Near the critical Reynolds number, however, the lateral velocity of the particle moving along the heteroclinic orbit is so small that the particle keeps migrating in much longer axial distance in the second stage, as seen from the curve for $Re=260$ in figure 7. As another example for this Reynolds number, a trajectory starting from $(y/(D/2),z/(D/2))=(0.2,0.1)$ is also plotted by thin dashed lines in figure 7(a,b). Note that the horizontal axis for these curves is shrunk by one-third (the scale is shown on the upper axis). In this case also, the particle approaches asymptotically the stable equilibrium position on the heteroclinic orbit (circle in figure 6 a). The travelling distance before reaching the stable equilibrium position is estimated to be $L_{2}/D\approx 1100$ from figure 7. Thus, it is demonstrated that, at $Re\approx Re_{c}$ , $L_{2}$ is much larger than $L_{1}$ , possibly by two orders of magnitude. Therefore, we conclude from the present results that an increase in $Re$ decreases $L_{1}$ and it also decreases $L_{2}$ in general, but for $Re\approx Re_{c}$ , $L_{2}$ could be extremely large compared to $L_{1}$ .

In Miura et al. (Reference Miura, Itano and Sugihara-Seki2014), we have discussed the entry length after which lateral migration has fully developed. Assuming that the lateral force is balanced by a viscous Stokes drag, we derived the conventional expression of the entry length $L_{e}$ as $L_{e}/l^{\prime }\approx 6{\rm\pi}A^{-1}Re^{-1}(d/D)^{-3}$ (Matas et al. Reference Matas, Morris and Guazzelli2004a ; Di Carlo et al. Reference Di Carlo, Irima, Tompkins and Toner2007) when the lateral force is expressed as $F=A(d/D)^{4}{\it\rho}U^{2}D^{2}$ , which was suggested from an asymptotic theory (Ho & Leal Reference Ho and Leal1974; Schonberg & Hinch Reference Schonberg and Hinch1989). Here, $l^{\prime }$ represents a lateral migration distance and $A$ is a constant, depending on the particle position in the cross-section as well as the Reynolds number. Although an alternative expression of $L_{e}/l^{\prime }$ proportional to $(d/D)^{-2}$ instead of $(d/D)^{-3}$ may be obtained for microchannel flows (Di Carlo Reference Di Carlo2009; Choi et al. Reference Choi, Seo and Lee2011), there are no differences between these two expressions in the case of a fixed size ratio $d/D$ as in the present study. Considering the two stage migration properties, we notice that the entry length $L_{e}$ estimated from the above formula may correspond to the length of the first stage, $L_{1}$ . Thus, we here estimate the entry length $L_{e}$ using the magnitude of the lateral force obtained in the present numerical computations and compare them with the values of $L_{1}$ obtained from the trajectories shown in figure 7.

Figures 4(a)–6(a) indicate that $F/{\it\rho}U^{2}d^{2}\approx 4,3$ and $2\times 10^{-3}$ at $Re=144,260$ and 514, respectively, in the region near $(y/(D/2),z/(D/2))=(0.5,0.375)$ . These values lead to $L_{e}/D\approx 50,36$ and 28 at $Re=144,260$ and 514, respectively, if we assume $l^{\prime }\approx D/6$ , i.e. the lateral migration distance is set to be one-sixth of the channel width in the first stage. These estimates show good agreement with the above values of $L_{1}$ obtained for thick lines in figure 7. Thus, $L_{e}$ is indicated to provide the length of the first stage $L_{1}$ . It may be useful to note that the entry length $L_{e}$ , in general, decreases slower than $Re^{-1}$ with increasing $Re$ , since an increase in $Re$ decreases the constant $A$ in the force expression, as shown in the present study.

The main findings of the present study include that the length of the second stage, $L_{2}$ , is much larger than $L_{1}$ or $L_{e}$ . This result suggests that the length of channels required for particle focusing is much larger than $L_{e}$ , which rationalizes the experimental observation of Miura et al. (Reference Miura, Itano and Sugihara-Seki2014); the entry length $L_{e}$ estimated using the above formula was much shorter than that expected from the experimental results of the particle distribution. For microchannel flows at low $Re\;({\leqslant}60)$ , the experimental study of Choi et al. (Reference Choi, Seo and Lee2011) reported that the focusing on the channel face equilibrium position was rather weak in the second stage when the particle distributions were measured at $240D{-}950D$ downstream of the entrance. When the measurements were performed at $640D$ downstream, the fractions of the particles located in the region of $-5^{\circ }\leqslant {\it\theta}\leqslant 5^{\circ }$ around the channel face equilibrium position were less than 30 % for $d/D=0.075$ and about 50 % for $d/D=0.16$ , at $Re=60$ . For smaller $Re$ , the fractions were reduced even when the measurements were performed further downstream. In their study, the length of $L_{e}$ estimated using the above formula described well their experimental observation for the first stage, and $L_{2}$ was assessed to be twice of this length. However, these weak focusing properties in the second stage may imply longer $L_{2}$ than that of the authors’ assessment in their paper. Further studies may be required for the entry length, since this information is one of the most important factors for the design of microchannels for practical applications such as separation, filtration and sorting of particles. We will return to this matter again at the end of the next section.

3.4. Reynolds number dependence of the equilibrium positions

The locations of the channel face and corner equilibrium positions, $y_{e}$ , are plotted against $Re$ in figure 8(a,b), together with the experimental data (Miura et al. Reference Miura, Itano and Sugihara-Seki2014, Reference Miura, Harada, Itano and Sugihara-Seki2015). In figure 8(a), the results of microchannel experiments at lower $Re$ are also plotted for comparable size ratios $(d/D=0.075{-}0.17)$ . The present numerical simulation showed that the channel face equilibrium positions always exist in the range of $Re$ from $20$ to $1000$ and they are stable. As $Re$ increases from small values, micro scale experiments have demonstrated an initial increase in equilibrium position $y_{e}$ (Choi et al. Reference Choi, Seo and Lee2011; Abbas et al. Reference Abbas, Magaud, Gao and Geoffroy2014; Sugihara-Seki Reference Sugihara-Seki2015), while larger millimetre scale experiments have shown a decrease in $y_{e}$ (Miura et al. Reference Miura, Itano and Sugihara-Seki2014, Reference Miura, Harada, Itano and Sugihara-Seki2015). These features are well described by the present numerical simulation.

In figure 8(b), hollow and filled circles denote stable and unstable channel corner equilibrium positions, respectively. The hollow circles in figure 8(b) are in fairly good agreement with the experimental results (Miura et al. Reference Miura, Itano and Sugihara-Seki2014, Reference Miura, Harada, Itano and Sugihara-Seki2015), and exhibit a monotonic increase with increasing $Re$ . This monotonic increase in $y_{e}$ is in agreement with experimental observations of circular tube flows (Segre & Silberberg Reference Segre and Silberberg1962; Matas et al. Reference Matas, Morris and Guazzelli2004a ) and the predictions of perturbation analysis (Asmolov Reference Asmolov1999; Matas et al. Reference Matas, Morris and Guazzelli2009). In circular tube flows, the equilibrium radial position of neutrally buoyant spherical particles, known as the Segre-Silberberg annulus, was reported to move monotonically towards the tube wall with increasing $Re$ up to $Re=1700$ (Segre & Silberberg Reference Segre and Silberberg1962; Matas et al. Reference Matas, Morris and Guazzelli2004a ). Next, we consider the implication of the results shown in figure 8, in relation to the results of previous studies, mainly theoretical and numerical studies on the inertial migration.

Figure 8.

Locations of (a) the channel face equilibrium positions and (b) the channel corner equilibrium positions. The circles represent the results of the present numerical simulations. The experimental results are represented by dots $(d/D=0.11)$ for Miura et al. (Reference Miura, Itano and Sugihara-Seki2014, Reference Miura, Harada, Itano and Sugihara-Seki2015), filled diamonds $(d/D=0.083)$ and filled squares $(d/D=0.17)$ for Sugihara-Seki (Reference Sugihara-Seki2015), and hollow diamonds $(d/D=0.075)$ , hollow triangles $(d/D=0.12)$ and hollow squares $(d/D=0.16)$ for Choi et al. (Reference Choi, Seo and Lee2011), and crosses $(d/D=0.11)$ for Abbas et al. (Reference Abbas, Magaud, Gao and Geoffroy2014). In figure 8(b), hollow and filled circles indicate stable and unstable equilibrium positions, respectively.

There have been many theoretical studies, based on the perturbation method, investigating the origin of the lateral forces on particles suspended in channel flows. Among them, the effect of curvature in Poiseuille velocity profile was first pointed out by Ho & Leal (Reference Ho and Leal1974). From dual perturbation expansions in $Re_{p}$ and $d/D$ with the condition of $Re_{p}\ll (d/D)^{2}$ , they indicated for plane Poiseuille flow that the interaction of the stresslet (symmetric force dipole), generated due to the particle rigidity constraint, and its wall correction with the bulk shear flow produces the inward force, whereas the interaction between the stresslet and the curvature of the bulk velocity profile yields the outward force towards the channel wall. Compared with these effects, the lag velocity of the particle and its slip rotation were found to be less important for neutrally buoyant particles. In accord with this conclusion for small particles at low $Re_{p}$ , the two-dimensional finite element simulations of Feng, Hu & Joseph (Reference Feng, Hu and Joseph1994) indicated that the driving forces of lateral migration for neutrally buoyant particles consist mainly of the wall repulsion and the lift force caused by the velocity profile curvature even for large particle sizes at relatively large Reynolds numbers. Thus, the equilibrium position of a neutrally buoyant spherical particle in channel flow could be determined by a balance between the outward force due to the shear gradient and the inward force due to the wall effect (Ho & Leal Reference Ho and Leal1974; Matas et al. Reference Matas, Morris and Guazzelli2004b ).

For higher $Re$ , Asmolov (Reference Asmolov1999) extended the matched asymptotic approach of Schonberg & Hinch (Reference Schonberg and Hinch1989) up to $Re=1500$ in the condition of $Re_{p}\ll 1$ and $d/D\ll 1$ , and showed that the equilibrium position in plane Poiseuille flow is shifted progressively towards the channel wall as $Re$ increases. The same trend for Poiseuille flow through circular tubes, which was obtained by Matas et al. (Reference Matas, Morris and Guazzelli2009) based on the matched asympotic expansion, explained successfully the experimental observation of the outward shift of the Segre-Silberberg annulus with $Re$ (Segre & Silberberg Reference Segre and Silberberg1962; Matas et al. Reference Matas, Morris and Guazzelli2004a ). This shift of the equilibrium position was understood by taking the above mentioned two competing lift forces into account; the increase in shear-gradient induced lift force with $Re$ is larger than that in wall-effect induced lift force, so that increasing $Re$ shifts the particle equilibrium positions closer to the tube wall (Asmolov Reference Asmolov1999; Matas et al. Reference Matas, Morris and Guazzelli2004b ).

With regard to square channel flows, Hood et al. (Reference Hood, Lee and Roper2015) adopted a perturbation method similar to Ho & Leal (Reference Ho and Leal1974) and showed that the results of Ho & Leal (Reference Ho and Leal1974) hold in a larger parameter space of $Re$ and $d/D$ than the authors assumed. Estimating the stresslet flow field, they concluded that viscous stresses are dominant over the inertial stress even at moderate Reynolds numbers in square channel flows. From the curved velocity profile in the square channel flow, in particular, along the $y$ - ( $z$ -) axis and the diagonal in the channel cross-section, the prediction of the outward shift of the equilibrium position with $Re$ may be also applied for the channel face and corner equilibrium positions if the particles are infinitesimally small. This could account for the present results for the channel face equilibrium positions at low $Re$ and the channel corner equilibrium positions (figure 8). However, the channel face equilibrium positions at high $Re$ exhibited the opposite trend. A possible cause may be the finite size ratio $d/D\;(=0.11)$ and consequently the finite particle Reynolds number ( $Re_{p}\approx 2.5{-}12.3$ ) in the present numerical study and the experiments (Miura et al. Reference Miura, Itano and Sugihara-Seki2014, Reference Miura, Harada, Itano and Sugihara-Seki2015), which violates the assumptions of the matched asymptotic expansion method. According to the study of Hood et al. (Reference Hood, Lee and Roper2015), this result may imply that, due to the finite size effect of the particle, an increase in $Re$ ( ${>}200$ ) increases the lateral force due to the wall effect to a more dominant level than that due to the velocity shear gradient, for a particle placed near the centre of the channel faces. In fact, the present study demonstrates that on the $y$ - (or $z$ -)axis, the lateral forces near the centreline are decreased compared to near the channel wall, as shown in figures 4(a)–6(a).

The importance of the finite particle size for particle migration properties was pointed out in several studies (Feng et al. Reference Feng, Hu and Joseph1994; Di Carlo Reference Di Carlo2009; Matas et al. Reference Matas, Morris and Guazzelli2009; Hood et al. Reference Hood, Lee and Roper2015). In particular, the experimental study of Matas et al. (Reference Matas, Morris and Guazzelli2004a ) for circular tube flows reported the presence of a new equilibrium position closer to the tube centre in addition to the Segre–Silberberg annulus at elevated $Re$ , and the majority of the particles were found to lie on the additional inner annulus for $Re>700$ . Since their matched asymptotic expansion analysis for particles of vanishing size did not provide an additional zero of the lateral force, representing a new equilibrium position, in the range of $Re$ from 1 to 2000, they concluded that the presence of the inner annulus observed experimentally is a result of the finite size of the suspended particle (Matas et al. Reference Matas, Morris and Guazzelli2009).

The shift of the particle positions towards the channel centreline reported in high aspect ratio rectangular microchannel flows (Ciftlik, Ettori & Gijs Reference Ciftlick, Ettori and Gijs2013) may be related to the inward shift of the channel face equilibrium position with high $Re$ ( ${>}200$ ). Ciftlik et al. (Reference Ciftlick, Ettori and Gijs2013) measured fluorescent images of particle streams at downstream positions in the range of $Re$ from 75 to 1500. From their observations in the direction perpendicular to the short channel face, they reported that the positions of the streams are almost constant up to $Re=300$ , and for $300<Re<450$ they migrate towards the channel centreline with increasing $Re$ . Figure S3 of their paper showed that, beyond this value, the focusing positions of the particles approach the channel wall as $Re$ increases up to 1500. The inward shift of the particle focusing positions for $300<Re<450$ is similar to the present result for the channel face equilibrium position, as shown in figure 8(a), although the tendencies in the other ranges of $Re$ are different. These discrepancies may result from differences in experimental conditions, including the aspect ratio of the channel cross-section.

One of the key results of the present study is the large entry length, $L_{2}$ , for particle focusing. As discussed in the previous section, the first stage of particle migration may be described by considering two competing lateral forces: the outward force due to the shear gradient and the inward force due to the wall effect. Thus, the asymptotic expansion analysis in $Re_{p}$ and $d/D$ may explain the particle behaviour even at moderate Reynolds numbers in the first stage (Hood et al. Reference Hood, Lee and Roper2015). However, in the second stage, where the two lateral forces are almost balanced on the heteroclinic orbit, the next order terms of the expansion may be necessary to take into account for describing the migration behaviour precisely. Particle slip velocity as well as slip rotation may play a role in determining the particle lateral motion along the heteroclinic orbit, the contributions of which can be neglected compared to the leading order terms in the first stage. Our results show that at low $Re\;({<}Re_{c})$ , the lateral force acts in the direction from the channel corner equilibrium position to the channel face equilibrium position along all the heteroclinic orbit. This makes the channel face equilibrium position stable and the channel corner equilibrium position unstable. At high $Re\;({>}Re_{c})$ , a saddle point appears on the heteroclinic orbit, and the direction of the lateral force is reversed between the saddle point and the channel corner equilibrium position. This makes the corner equilibrium position stable. It would be important to elucidate the origin of these variations of the lateral force in the second stage of particle migration, since this will enable us to estimate the entry length, $L_{2}$ , for the particle focusing.

Another important issue we need to explore in the near future is the effect of the hydrodynamic interaction between particles, since it may affect the location and the number of equilibrium positions (Humphry et al. Reference Humphry, Kulkami, Weitz, Morris and Stone2010) and be a possible cause for the dispersion of particles observed experimentally (Di Carlo Reference Di Carlo2009; Abbas et al. Reference Abbas, Magaud, Gao and Geoffroy2014). In addition, Abbas et al. (Reference Abbas, Magaud, Gao and Geoffroy2014) found that particles suspended in square channel flows are preferentially located near the channel centre at $Re<1$ , opposed to predictions based on the Stokes flow or extrapolation from inertial migration at higher $Re$ , and this phenomenon was considered to be generated by the shear-induced particle interactions. It may be that the inter-particle interaction affects also the migration behaviour at higher $Re$ . The effect of particle size relative to the channel width is also important, considering the difference in Reynolds number and particle Reynolds number. This effect has been investigated by a regular perturbation analysis of the Stokes equation (Hood et al. Reference Hood, Lee and Roper2015). A recent experimental study in our laboratory for $100\leqslant Re\leqslant 1000$ has demonstrated considerably enhanced focusing properties for larger particles near the equilibrium positions, whereas a clear dependence on particle size has not been observed for the locations of the equilibrium position and the transition behaviours in the range of $0.05\leqslant d/D\leqslant 0.11$ (Miura et al. Reference Miura, Harada, Itano and Sugihara-Seki2015). The significant effect on the particle focusing may reflect the high power dependence of the magnitude of the lateral forces on the size ratio. Further detailed experimental and numerical studies are necessary to elucidate this effect. Another issue to explore is related to the physics near $Re=Re_{c}$ . At $Re\approx Re_{c}$ , we have revealed extremely small lateral forces along the heteroclinic orbit and the presence of a new equilibrium position on the orbit. It is uncertain whether the saddle point and the new equilibrium point appear simultaneously and whether other equilibrium positions are present on the heteroclinic orbit. It would be interesting to carry further investigations on these issues.