Population splitting of rodlike swimmers in Couette flow

TL;DR

This paper analyzes how self-propelled rod-like swimmers in a shear flow undergo a population splitting transition, reversing swimming directions and forming two distinct groups, driven by shear rate and characterized by a phase diagram.

Contribution

It introduces the concept of population splitting in active suspensions under shear flow and maps the transition between regimes using a phase diagram based on Peclet numbers.

Findings

Swimmers split into two populations with opposite orientations at high shear rates.

Transition from self-propulsion to shear dominance is discontinuous.

Explains non-monotonic orientation behavior of swimmers near walls.

Abstract

We present a quantitative analysis on the response of a dilute active suspension of self-propelled rods (swimmers) in a planar channel subjected to an imposed shear flow. To best capture the salient features of shear-induced effects, we consider the case of an imposed Couette flow, providing a constant shear rate across the channel. We argue that the steady-state behavior of swimmers can be understood in the light of a population splitting phenomenon, occurring as the shear rate exceeds a certain threshold, initiating the reversal of swimming direction for a finite fraction of swimmers from down- to upstream or vice versa, depending on swimmer position within the channel. Swimmers thus split into two distinct, statistically significant and oppositely swimming majority and minority populations. The onset of population splitting translates into a transition from a self-propulsion-dominated regime to a shear-dominated regime, corresponding to a unimodal-to-bimodal change in the probability distribution function of the swimmer orientation. We present a phase diagram in terms of the swim and flow Peclet numbers showing the separation of these two regimes by a discontinuous transition line. Our results shed further light on the behavior of swimmers in a shear flow and provide an explanation for the previously reported non-monotonic behavior of the mean, near-wall, parallel-to-flow orientation of swimmers with increasing shear strength.

Equations

$$\mathbf{u}(\mathbf{r})=\frac{U_0}{2}\cdot \frac{y+H}{H}\hat{\mathbf{x}}=\dot{\gamma }(y+H)\hat{\mathbf{x}},$$

$$\dot{\gamma }=\frac{U_0}{2H}.$$

$$\frac{\partial \Psi }{\partial t}+{\nabla }_{\mathbf{r}}\cdot ({\mathbf{v}}_{\mathbf{r}}\Psi )+{\nabla }_{\mathbf{p}}\cdot ({\mathbf{v}}_{\mathbf{p}}\Psi )=D_T{\nabla }_{\mathbf{r}}^2\Psi +D_R{\nabla }_{\mathbf{p}}^2\Psi ,$$

$${\mathbf{v}}_{\mathbf{r}}$$

$${\mathbf{v}}_{\mathbf{p}}$$

$$V_s\frac{\partial }{\partial y}(\Psi \sin \theta )-\dot{\gamma }\frac{\partial }{\partial \theta }(\Psi {\sin }^2\theta )=D_T\frac{{\partial }^2\Psi }{\partial y^2}+D_R\frac{{\partial }^2\Psi }{\partial {\theta }^2}.$$

$$\left(D_{T}\frac{\partial}{\partial y}-V_{s}\sin\theta\right)\Psi(y,\theta)\bigg{|}_{y=\pm H}=0.$$

$$2P{e}_s\frac{\partial }{\partial \tilde{y}}(\tilde{\Psi }\sin \theta )-P{e}_f\frac{\partial }{\partial \theta }\left(\tilde{\Psi }{\sin }^2\theta \right)={\xi }^2\frac{{\partial }^2\tilde{\Psi }}{\partial {\tilde{y}}^2}+\frac{{\partial }^2\tilde{\Psi }}{\partial {\theta }^2},$$

$$\left(\xi^{2}\frac{\partial}{\partial\tilde{y}}-2Pe_{s}\sin\theta\right)\tilde{\Psi}(\tilde{y},\theta)\bigg{|}_{\tilde{y}=\pm 1}=0,$$

$$P{e}_s=\frac{1/D_R}{2H/V_s}=\frac{V_s}{2H{D}_R},$$

$$P{e}_f=\frac{1/D_R}{1/\dot{\gamma }}=\frac{U_0}{2H{D}_R},$$

$${\xi }^2=\frac{1/D_R}{H^2/D_T}=\frac{D_T}{H^2{D}_R},$$

$$c(y)={\int }_0^{2\pi }\Psi (y,\theta )\mathrm{d}\theta ,$$

$${\int }_{-1}^1\tilde{c}(\tilde{y})\mathrm{d}\tilde{y}=2.$$

$$⟨\cdots ⟩=\frac{1}{\tilde{c}(\tilde{y})}{\int }_{-\pi /2}^{3\pi /2}(\cdots )\tilde{\Psi }(\tilde{y},\theta )\mathrm{d}\theta ,$$

$$U_4=1-\frac{{\mu }_4}{3{\mu }_2^2},$$

$$\alpha =\frac{\mathrm{d}\ln P{e}_f^{\ast }}{\mathrm{d}\ln P{e}_s}>0,$$

$$⟨{\tilde{p}}_x⟩$$

$$⟨{\tilde{p}}_y⟩$$

Full text

I Introduction

Self-propelled particles, or swimmers, and other active matter systems Ramaswamy (2010); Marchetti et al. (2013); Elgeti et al. (2015); Goldstein (2015); Son et al. (2015); Cates and Tailleur (2015); Rusconi and Stocker (2015); Golestanian et al. (2011); Paxton et al. (2006); Romanczuk et al. (2012); Bechinger et al. (2016) have been at the focus of extensive multidisciplinary research over the last few decades due to their important applications from micro- down to nanoscale Koumakis et al. (2013); Dreyfus et al. (2005); Paxton et al. (2004); Medina-Sánchez et al. (2016); Cheang et al. (2014); Khalil et al. (2014); Vach et al. (2013, 2015) and also their remarkable theoretical and fundamental aspects Purcell (1977); Najafi and Golestanian (2004); Lauga (2011); Zöttl and Stark (2014); Sokolov and Aranson (2012); Golestanian et al. (2007); Elgeti et al. (2015); Marchetti et al. (2013); Ramaswamy (2010); Golestanian et al. (2011); Goldstein (2015); Bechinger et al. (2016); Becker et al. (2003); Najafi and Golestanian (2005); Cates and Tailleur (2015); Son et al. (2015); Berg and Anderson (1973); Berg (2003); Budrene and Berg (1991); Guasto et al. (2012); Rusconi and Stocker (2015); Golestanian et al. (2005); Blair (1995); Romanczuk et al. (2012); Paxton et al. (2006). Swimmers are commonplace in biology and biotechnological applications Ramaswamy (2010); Marchetti et al. (2013); Elgeti et al. (2015); Goldstein (2015); Son et al. (2015); Cates and Tailleur (2015); Rusconi and Stocker (2015); Golestanian et al. (2011); Paxton et al. (2006); Romanczuk et al. (2012); Bechinger et al. (2016); Rusconi and Stocker (2015); Berg and Anderson (1973); Berg (2003); Budrene and Berg (1991); Guasto et al. (2012) with examples ranging from uniflagellated sperm cells Gray (1955); Shack et al. (1974); Woolley (2003); Friedrich et al. (2010); Friedrich and Jülicher (2007), the biflagellate alga Chlamydomonas reinhardtii Ringo (1967); Goldstein (2015); Polin et al. (2009), the colonial alga Volvox carteri Drescher et al. (2009); Brumley et al. (2015) and the multiflagellate bacterium Escherichia coli Berg (2003); Berg and Anderson (1973); Blair (1995) to artificial micro-/nanoswimmers, fabricated using patchy colloids and Janus micro-/nanoparticles Walther and Muller (2013); Valadares et al. (2010). This latter class of swimmers exhibit well-defined, stochastic, self-propelled motion due to asymmetry in their catalytic surface reactions in specific solvent mixtures Jiang et al. (2010); Howse et al. (2007); Volpe et al. (2011); Paxton et al. (2006); Golestanian et al. (2005, 2011); Elgeti et al. (2015); Marchetti et al. (2013); Ramaswamy (2010); Romanczuk et al. (2012); Bechinger et al. (2016); Takagi et al. (2013); Saintillan and Shelley (2007); Hong et al. (2007); their motion closely resembles microorganismal swimming, which involves the hugely more complex mechanisms such as flagellar or ciliary movements.

Fluid suspensions of active particles are typically confined by rigid boundaries. The near-wall behavior of swimmers has received a lot of attention (see, e.g., Refs. Zhang et al. (2010); Sabass and Seifert (2010); Elgeti and Gompper (2016); Berke et al. (2008); Lauga et al. (2006); Spagnolie and Lauga (2012); Li et al. (2011); Li and Tang (2009); Li and Ardekani (2014); Hernandez-Ortiz et al. (2005); Underhill et al. (2008); Hernandez-Ortiz et al. (2009); Schaar et al. (2015); Hu et al. (2015); Elgeti and Gompper (2009, 2013); Kantsler et al. (2014, 2005); Rusconi et al. (2014); Elgeti and Gompper (2015); Costanzo et al. (2012); Nash et al. (2010); Mathijssen et al. (2016a, b, c); Ezhilan and Saintillan (2015); Tung et al. (2015); Zöttl and Stark (2013); Chilukuri et al. (2015) and the references therein) and turns out to have a profound role in many biological processes, such as biofilm formation O’Toole et al. (2000); Hall-Stoodley et al. (2004) and sperm motility Elgeti et al. (2010); Denissenko et al. (2012); Crowdy (2011); Kantsler et al. (2014), as well as microfluidic setups used for manipulation and separation of active agents, such as bacteria Wu et al. (2009); Qiu et al. (2009); Son et al. (2015).

Very commonly, external flow is present within the confined regions containing active suspensions Sternberg et al. (1999); Pereira et al. (2002). Specifically, near-surface fluid flow is known to cause swimming against the flow, i.e., upstream swimming, which has been observed in experiments Son et al. (2015); Hill et al. (2007); Kaya and Koser (2012, 2009); Rosengarten et al. (1988); Meng et al. (2005); Rusconi et al. (2014); Zöttl and Stark (2013) and the effect attributed to a number of factors, including hydrodynamic surface interactions and/or shear reorientation of swimmers Kaya and Koser (2012, 2009); Ezhilan and Saintillan (2015); Nash et al. (2010); Hill et al. (2007); Tung et al. (2015); Kantsler et al. (2014, 2005); Rusconi et al. (2014); Mathijssen et al. (2016a, b, c); Zöttl and Stark (2013).

In the context of continuum models, Ezhilan and Saintillan Ezhilan and Saintillan (2015) have conducted a thorough study on the response of a confined suspension of active rodlike particles to imposed Poiseuille flow, using a kinetic model in which the active suspension is described by a joint position-orientation probability distribution function (PDF) governed by a non-interacting Smoluchowski equation. Without considering hydrodynamic inter-particle and surface-particle interactions, they have been able to observe many of the important effects arising from confinement and imposed flow coupled with self-propulsion. These include wall accumulation, shear-trapping and upstream swimming (see Ref. Ezhilan and Saintillan (2015) for a comprehensive review of the literature on these effects). It has thus been suggested that hydrodynamic interactions may not be essential to the major phenomena occurring in active suspensions under confinement.

Basing our work on the approach adopted by the above authors, we address the steady-state behavior of a dilute suspension of self-propelled, rodlike, Brownian particles in a planar channel, subjected to an imposed Couette flow. The characteristic feature of the Couette flow is that it presents a constant shear rate across the channel width, modeled here using a linear fluid velocity profile increasing from zero on one of the surfaces (bottom wall) to a finite value on the other surface (top wall). Thus, unlike the pressure-driven Poiseuille flow, where the shear rate varies gradually from zero at the centerline of the channel to maximum values at its boundaries, the Couette model provides a situation where the swimmers experience the same shear-induced torque, independent of their position within the channel. This helps create a more transparent arrangement to investigate the effects of imposed shear on swimmer density and orientation. Three characteristic regions were shown to exist under an imposed Poiseuille flow Ezhilan and Saintillan (2015). These were (a) a centerline depletion region, where self-propulsion and rotational diffusion dominate, (b) an intermediate shear-trapping region, which develops at sufficiently low (high) swim (flow) Péclet numbers due to the interplay between self-propulsion and imposed shear, and (c) swimmer accumulation regions near the walls, where all previous factors play a role. From these three, only the latter two exist when the flow is of the Couette type, as the centerline behavior does not manifest itself in the absence of changes in the shear rate across the channel width.

We shall demonstrate that the steady-state behavior of swimmers in the channel for different strengths of self-propulsion and imposed shear, represented by swim and flow Péclet numbers, respectively, can be described in the light of a population splitting phenomenon. The phenomenon is caused by shear alignment at flows strong enough to flip the swimming direction of a statistically significant (macroscopic) number of active particles from upstream to downstream, or vice versa. In other words, upon increasing the flow Péclet number, a growing fraction of swimmers from a majority population of up- or downstream swimmers will be converting to develop a minority population of swimmers moving in the opposite direction. This behavior can happen near each of the walls, provided that the flow Péclet number is increased beyond a threshold value. It corresponds to a transition from a unimodal to a bimodal PDF of swimmers, and underlies some of the salient features arising from the interplay between shear-induced torque and active self-propulsion, which can be studied and elucidated most clearly in the case of Couette flow. These include the non-monotonic behavior of the mean parallel-to-flow component of the orientation vector of swimmers, showing first an increasing and, then, a decreasing trend toward zero with increasing flow Péclet number. Such a behavior (which has also been found, but remained unexplored, in the case of Poiseuille flow Ezhilan and Saintillan (2015)) portrays the onset of the population splitting phenomenon.

Our results thus indicate that the standard picture for the near-wall behavior of rodlike swimmers in an imposed flow, in which the active particles are perceived to swim either up- or downstream near the walls, may be inaccurate as it is determined based on the mean parallel-to-flow orientation of swimmers. This latter quantity ceases to act as a suitable measure describing the ‘typical’ near-wall behavior of swimmers, when population splitting takes place and gives rise to a minority, but significant and measurable, population of particles swimming in the direction opposite to the majority population; hence, a situation whose characterization requires a knowledge of the higher-order moments of the PDF.

The organization of this paper is as follows: In Section II, we introduce the model used in our study of the shear-induced orientational behavior of rodlike swimmers subject to Couette flow. In Section III, we present the results obtained from our numerical investigations, which we conclude in Section IV.

II Model

II.1 Physical specifications

We consider a dilute suspension of self-propelled, rodlike, Brownian particles (swimmers) confined by two plane-parallel, rigid, no-slip walls forming a long, planar channel of width $2H$, subjected to an imposed Couette flow with constant shear rate, $\dot{\gamma}$, and laminar flow profile

[TABLE]

where $U_{0}$ is the (constant) speed of the moving top plate, ${\mathbf{\hat{x}}}$ the unit vector in the flow (or channel) direction, $y$ the perpendicular coordinate (see Fig. 1), and

[TABLE]

Swimmers are modeled as thin rods (needles) of high aspect ratio giving a Bretherton constant Bretherton (1962) approaching unity. They are described by the orientation vector $\mathbf{p}$, denoting the direction of active swimming at a constant speed $V_{s}$ (in the case of microorganisms, $\mathbf{p}$ may correspond to the tail-head direction as shown schematically in Fig. 1; note, however, that our model is generically applicable also to synthetic self-propelled micro- and nanorods Paxton et al. (2004); Takagi et al. (2013); Saintillan and Shelley (2007); Hong et al. (2007)). We shall be interested primarily in the steady-state properties of the system that exhibit translational symmetry in $x$-direction as well as in the direction perpendicular to the $x-y$ plane and, therefore, depend only on the coordinate $y$ and the orientation angle $\theta$. We neglect hydrodynamic as well as excluded-volume interactions between individual swimmers and between swimmers and the channel walls. This assumption is expected to remain valid in the dilute regime Ezhilan and Saintillan (2015).

II.2 The continuum model

In the continuum description, the joint PDF $\Psi(\mathbf{r},\mathbf{p},t)$ of finding a swimmer at position $\mathbf{r}$ with orientation $\mathbf{p}$ at time $t$ is governed by the Smoluchowski equation Doi and Edwards (1988); Ezhilan and Saintillan (2015)

[TABLE]

where $D_{T}$ and $D_{R}$ are translational and rotational diffusivities, and $\mathbf{v}_{\mathbf{r}}$ and $\mathbf{v}_{\mathbf{p}}$ the deterministic translational and rotational flux velocities, respectively. For rodlike swimmers under the action of the Couette flow, one has

[TABLE]

where $\mathbf{I}$ is the identity tensor. The deterministic translational flux takes contributions from self-propulsion and advection by the flow, while the rotational one takes contributions from the shear rate only. It is the combined action of self-propulsion, advection and shear by the imposed flow, as well as translational and rotational diffusion, that leads to an eventual steady-state distribution function for swimmers within the channel $\Psi=\Psi(y,\theta)$. Using $\mathbf{p}=(p_{x},p_{y})=(\cos\theta,\sin\theta)$, the Smoluchowski equation at steady state reduces to

[TABLE]

Since the walls are rigid and impermeable to fluid flow and active particles, the normal component of the deterministic translational flux must balance the normal component of the stochastic (diffusion) flux at $y=\pm H$, supplementing Eq. (5) with the boundary conditions

[TABLE]

II.3 Non-dimensionalization

Using the channel half-width $H$, and the rotational diffusion timescale $1/D_{R}$, as characteristic scales for length and time, respectively, the non-dimensionalization of Eq. (5) proceeds immediately as

[TABLE]

subject to the boundary conditions

[TABLE]

where $\tilde{y}=y/H$ is the dimensionless coordinate across the channel width and $\tilde{\Psi}=\Psi/c_{0}$ is the dimensionless PDF, with $c_{0}$ being the mean number density of swimmers in the channel. In Eqs. (7) and (8), we have adopted a notation close to that of Ref. Ezhilan and Saintillan (2015) and defined the dimensionless system parameters as follows: (i) the swim Péclet number,

[TABLE]

which gives the ratio of the rotational diffusion timescale $1/D_{R}$, to the swim (across channel) timescale $2H/V_{s}$, (or, equivalently, the ratio of the run length $V_{s}/D_{R}$, to the channel width), (ii) the flow Péclet number,

[TABLE]

which gives the ratio of the rotational diffusion timescale to the timescale for swimmer alignment with the axis of flow $1/\dot{\gamma}$, as the flow would attempt to align the swimmers in/against the direction of flow (i.e., the positive/negative $x$-axis, corresponding to the self-propelled particles swimming downstream/upstream, respectively), and (iii) the dimensionless parameter

[TABLE]

which gives the ratio of the rotational diffusion timescale over the translational diffusion timescale.

The number density of swimmers,

[TABLE]

is also normalized by $c_{0}$, giving the following constraint on the rescaled swimmer density, $\tilde{c}({\tilde{y}})\equiv c(y)/c_{0}$:

[TABLE]

The steady-state properties of swimmers in the channel follow from the solution of Eqs. (7), (8) and (13) in the computational domain $-1\leq\tilde{y}\leq 1$ and $0\leq\theta<2\pi$ (see Appendix A for details). The parameter space of the system is spanned by the three dimensionless parameters $Pe_{f}$, $Pe_{s}$ and $\xi$. The focus of our study will be on the interplay between shear alignment and self-propulsion of swimmers and, as such, we will vary $Pe_{f}$ and $Pe_{s}$ and, unless otherwise noted (Section III.7), we fix the third dimensionless parameter as $\xi=\sqrt{2/3}\simeq 0.82$. This particular value is chosen to match the base parameter values assumed in Ref. Ezhilan and Saintillan (2015) to enable like-for-like (quantitative) comparison between the cases of imposed Poiseuille flow considered in that work and the Couette flow considered here, as will be seen throughout the text (note that the third dimensionless parameter in Ref. Ezhilan and Saintillan (2015) is denoted by $\Lambda$ and its chosen base value $\Lambda=1/6$ corresponds to the mentioned base value for $\xi$, when the relation $\Lambda=\xi^{2}/(4Pe_{s}^{2})$ is used with $Pe_{s}=1$). We should, however, note that these base values imply a channel half-width comparable to particle size, if typical thermal values for the translational and rotational diffusion coefficients are assumed. Therefore, we shall later vary $\xi$ over a wider range, including much smaller values of $\xi$ (much larger channel widths), showing that the key outcomes of our analysis for the shear-induced behavior of swimmers indeed hold over a wide range of dimensionless parameter values (which can themselves be mapped to a broad range of values for the actual system parameters).

III Results and Discussion

III.1 Swimmer density profiles

Figure 2 shows the numerically calculated swimmer density profiles across the channel for different swim and flow Péclet numbers ($Pe_{s}=0.5$ and 2 in panels a and b, respectively, with $Pe_{f}$ increased from zero to 100). The results confirm the well-established observation that swimmers tend to accumulate at the boundaries, in this case the channel walls at $\tilde{y}=\pm 1$. While active swimming (increased $Pe_{s}$) tends to drive the particles toward the walls, making the profiles more steep near the walls, the imposed shear (increased $Pe_{f}$) tends to align the particles in the horizontal direction, making the profiles more uniform (see, e.g., the nearly uniform density profiles in Fig. 2a, where active self-propulsion is relatively weak). For sufficiently large values of $Pe_{s}$ and $Pe_{f}$ (Fig. 2b), the interplay between wall accumulation and shear alignment effects leads to a plateau-like region around the centerline of the channel and regions with lowered density between the centerline and the walls. The behavior in these latter regions is thus intrinsically different from the centerline depletion found in the case of an imposed Poiseuille flow Ezhilan and Saintillan (2015), where swimmer depletion around the centerline (with vanishing shear rate) is driven by self-propulsion and rotational diffusion. (Note also that, although for the chosen parameter values the curves in Fig. 2 appear to go through nearly the same points close to $\tilde{y}=\pm 0.5$, this feature does not necessarily hold in other ranges of parameters such as larger values of $Pe_{s}$.)

III.2 Mean swimmer orientation profiles

Figure 3 shows how ensemble averages of the parallel-to-flow and perpendicular-to-flow components of the swimmer orientation vector, $\langle\tilde{p}_{x}\rangle=\langle\cos\theta\rangle$ and $\langle\tilde{p}_{y}\rangle=\langle\sin\theta\rangle$, respectively, vary across the channel. Here, $\left\langle\cdots\right\rangle$ denotes the average

[TABLE]

where, with no loss of generality, we have chosen the $[-\pi/2,3\pi/2)$ range for $\theta$, so that the interval can be conveniently divided into a first $[-\pi/2,\pi/2)$ and a second $[\pi/2,3\pi/2)$ half, corresponding to downstream and upstream swimmers, respectively.

When there is no imposed flow ($Pe_{f}=0$), the swimmer orientation is, on average, perpendicular to the direction of channel and toward the walls (i.e., $\langle\tilde{p}_{x}\rangle=0$ and $\langle\tilde{p}_{y}\rangle\neq 0$ as seen in Fig. 3). When Couette flow is imposed, the swimmers take a finite mean parallel-to-flow orientation (see, e.g., the curves for $Pe_{f}=5$ and 10 in Fig. 3). In this case, negative (positive) values of $\langle\tilde{p}_{x}\rangle$ for $\tilde{y}<0$ ($\tilde{y}>0$) are representative of the fact that, on average, swimmers orient in the upstream (downstream) direction in the lower (upper) half of the channel. As noted before, this kind of up- or downstream swimming in imposed flows is already well established (see, e.g., Refs. Son et al. (2015); Hill et al. (2007); Kaya and Koser (2012, 2009); Rosengarten et al. (1988); Meng et al. (2005); Rusconi et al. (2014); Ezhilan and Saintillan (2015)) and can be understood in terms of the shear-induced torques acting on swimmers near the lower or upper wall of the channel (see Section III.3).

It is interesting to note that, once the shear flow is established and gradually strengthened, both components of the mean orientation vector become smaller. The reduction in $\langle\tilde{p}_{y}\rangle$ with increasing shear strength (or flow Péclet number) is a natural consequence of the horizontally directed Couette flow barring swimmer tendency to approach the walls. However, the reduction in $\langle\tilde{p}_{x}\rangle$ with increasing shear strength (see, e.g., the curves for $Pe_{f}>10$ in Fig. 3) is less trivial. A more clear view of the non-monotonic (monotonic) behavior of the mean parallel-to-flow (perpendicular-to-flow) orientation component is given in Fig. 4, where $\langle\tilde{p}_{x}\rangle$ and $\langle\tilde{p}_{y}\rangle$ at the top wall are plotted as a function of the flow Péclet number and for different swim Péclet numbers. As seen, increasing $Pe_{f}$ leads to a monotonic decrease in $\langle\tilde{p}_{y}\rangle$ and a corresponding increase in $\langle\tilde{p}_{x}\rangle$, but only up to a certain threshold (the peaks of the plots in Fig. 4a), beyond which $\langle\tilde{p}_{x}\rangle$ decreases with increasing $Pe_{f}$. Such a behavior has also been found, but remained unexplained, for swimmers in a Poiseuille flow Ezhilan and Saintillan (2015). This behavior in fact portrays the onset of a population splitting phenomenon in swimmer orientation that shall be discussed in the following section (see also Appendix B for more details).

III.3 The population splitting phenomenon

Figure 5 shows the profiles of the PDF, $\tilde{\Psi}(\tilde{y},\theta)$, within the computational domain $-\pi/2\leq\theta<3\pi/2$ and $-1\leq\tilde{y}\leq 1$. As can be seen in the captions, we have chosen $Pe_{s}=\left\{0.25,1,5\right\}$ and $Pe_{f}=\left\{5,20,100\right\}$ as sample parameter values representing (in the same order as they appear in the brackets) active self-propulsion and shear flow of relatively low, medium, and high strength, respectively. As the color bars show, the warmer (cooler) colors represent higher (lower) swimmer density. The symmetry in all profiles arises from the linear structure of Couette flow, resulting in $(\tilde{y},\theta)\rightarrow(-\tilde{y},\pi+\theta)$ symmetry for the solution of the Smoluchowski equation, Eq. (7).

The tendency of swimmers to accumulate at confining boundaries is evident in all frames of Fig. 5: There are always two areas of high swimmer density near the top and bottom wall of the channel. Recalling the definition of $\theta$ discussed after Eq. (14), the first and second halves of the horizontal axes in the profiles of Fig. 5 correspond to down- and upstream swimming, respectively. Hence, the high density region near the top and bottom wall of the channel corresponds to active particles swimming in the downstream and upstream direction, respectively.

It is interesting to note that, in the self-propulsion-dominated regime (represented by typically large values of $Pe_{s}$ such as in Fig. 5c), there is a single population of swimmers on each of the channel walls, swimming in the up- or downstream direction depending on whether they are closer to the bottom or the top wall, respectively (on the top wall in panel c, this single population can be identified as the high-density region around an orientation angle between 0 and $\pi/2$, representing a tilting of swimmer orientation toward the wall). In the shear-dominated regime (represented by typically large values of $Pe_{f}$ such as in Fig. 5g), there are two distinct and well-separated populations of swimmers on each of the walls, swimming in opposite down- or upstream directions (on the top wall in panel g, these populations can be identified as the two high-density regions around the orientation angles 0 and $\pi$, corresponding to down- and upstream swimmers, respectively). Similar features can be seen in other cases (e.g., in panels b and f that appear closer to the self-propulsion-dominated behavior in panel c, and in panels d and h that appear closer to the shear-dominated behavior in panel g). However, the situation becomes more intriguing in the intermediate regime, where neither of the mechanisms are completely dominant (e.g., in panels a, e, i). In these cases, there still appear to be two distinct populations of swimmers on each of the walls but with very different densities. Therefore, the near-wall swimmers are still divided into two distinct populations (with angular separation of around $\pi$), representing a majority and a minority population, which, on the top wall, swim in the down- and upstream directions, respectively. The exact opposite of this situation occurs on the bottom wall. Therefore, the interplay between self-propulsion and imposed shear appears to lead to a discernible change in swimmer distribution from a single population to two separate populations with different up- and downstreaming behavior near the walls. This phenomenon can also occur in other positions across the channel width, depending on the system parameters (compare, e.g., panels c, f and i). We shall refer to this behavior as the population splitting phenomenon, whose characteristic features will be discussed in more detail below.

III.4 Bimodality of the near-wall swimming

The population splitting phenomenon, and how it emerges and develops upon varying the swim and flow Péclet numbers can be elucidated further by focusing on the PDF of swimmers on the channel walls. As before, we concentrate on the top wall at $\tilde{y}=1$ (all that follows is true for the bottom wall at $\tilde{y}=-1$, if upstream and downstream change places in the wording).

The PDF, $\tilde{\Psi}(\tilde{y}=1,\theta)$, has been plotted in Fig. 6 for four different $Pe_{s}$ and a wide range of $Pe_{f}$. Figure 6a shows the case with non-active rodlike particles ($Pe_{s}=0$). With no driving force for the particles to tend toward the walls, they would exhibit equal chance of aligning in or against the direction of flow (corresponding to the left or right peaks, respectively). This results in a symmetric bimodal distribution (with equal-height peaks) for the particles and vanishing components for their mean orientation vector (see Fig. 4). Note that, in this case, the particles are passive thin rods Doi and Edwards (1988) and, as such, they can only nematically align with the flow, regardless of their positions in the channel.

The situation is different for active swimmers, in which case the near-wall behavior of the particles plays a central role and their population splitting leads to two oppositely swimming majority and minority populations. In the absence of flow ($Pe_{f}=0$), the active particles will, on average, swim straight up toward the (top) channel wall, with no hindering effect from an imposed flow. This results in symmetric unimodal $\tilde{\Psi}$-profiles about the orientation angle $\theta_{\ast}=\pi/2$ (light-gray solid curves in Figs. 6b and c). When the shear is imposed and strengthened, it will attempt to align the swimmers in the direction of flow, leading to shift of this single peak from $\theta_{\ast}=\pi/2$ toward smaller angles (see, e.g., the black dashed curves for $Pe_{f}=5$ in panels b and c). In this regime, self-propulsion is still the dominant mechanism driving the active particles, which form a single statistical population with a unimodal (but not necessarily symmetric) PDF peaked around the mean swimmer orientation, pointing in the downstream direction near the top wall. When $Pe_{f}$ is further increased at a given $Pe_{s}$, we find a second peak, whose height increases and gradually tends toward that of the first peak, as the imposed flow becomes stronger (this is seen more clearly in log-linear scale in the insets of panels b and c). The emergence of this smaller second peak marks the onset of the population splitting phenomenon, as from there onwards, with $Pe_{f}$ increasing further, a growing number of swimmers from the majority population of downstream swimmers will be converting to develop the minority population of upstream swimmers (see the schematic illustration in Fig. 7). The exact opposite of this conversion occurs on the bottom wall. It is thus clear that the PDF of swimmer orientation changes from unimodal to bimodal, exhibiting two peaks that remain well separated (by an angular separation of almost $\pi$) and represent two distinct macroscopic populations of swimmers, as the population splitting takes place and develops further as the flow Péclet number, $Pe_{f}$, is increased. Note also that larger values of $Pe_{f}$ lead generally to narrower PDF peak(s) as stronger shear is more effective in barring particle tilting due to rotational diffusion.

Going back to the results for the parallel-to-flow component of the mean swimmer orientation vector $\langle\tilde{p}_{x}\rangle$ (Section III.2), we found that this quantity exhibits a non-monotonic behavior, eventually, tending to zero across the channel (Fig. 3a) and, in particular, near the channel walls (Fig. 4a) as $Pe_{f}$ is increased. Our foregoing analysis thus clearly shows that describing the near-wall behavior of swimmers in an imposed flow merely based on their mean parallel-to-flow orientation, in which swimmers are perceived to swim either up- or downstream near the walls, may be inaccurate and may even lead to the misleading interpretation that the near-wall, up- or downstream swimming trends are weakened at large shear rates, where $\langle\tilde{p}_{x}\rangle$ tends to zero. This is in fact where the swimmers split into two nearly equal but oppositely swimming populations. Thus, in order to capture the bimodal nature of the swimmer orientation, one needs to consider statistical quantities other than the mean orientation, including higher-order moments of swimmer orientation angle as shall be discussed later.

III.5 The bimodal ratio

The flow strength ($Pe_{f}$) required to generate a certain degree of population splitting would be larger under stronger swimmer self-propulsion (larger $Pe_{s}$). This can be seen from Fig. 8, which shows the bimodal ratio $R_{\Psi}$, defined as the ratio of the minority population of upstream swimmers on the top channel wall (right peak of the profiles in Fig. 6) and the majority population of downstream swimmers (left peak of the profiles in Fig. 6). When there is no active self-propulsion, the two (non-swimming) populations are exactly equal and the ratio is identical to 1, regardless of the flow strength. For finite $Pe_{s}$, it can be seen that the ratio is smaller, for all flow strengths, when self-propulsion is stronger (larger $Pe_{s}$).

Under imposed flow of (theoretically) infinite strength, the two populations will become equal (in this sense like the non-active situation) and in the near-wall $\tilde{\Psi}$-profiles, we would have two Dirac delta functions at $\theta=0$ and $\theta=\pi$. The curves in Fig. 8 are accordingly seen to verge toward $R_{\Psi}=1$ as the flow Péclet number is increased.

It is important to note that, at experimentally feasible values of $Pe_{s}$, like those shown on the graph, the situation is still far from the theoretical limits of zero self-propulsion or infinite flow strength. This has resulted in $R_{\Psi}$ falling well below the limiting value $R_{\Psi}\rightarrow 1$.

III.6 Onset of population splitting: Binder cumulant

As noted in Sections III.3-III.5, the onset and the extent of population splitting and bimodality in swimmer orientation in a Couette flow vary significantly depending on the flow and swim Péclet numbers as well as the vertical position of the swimmers across the channel.

In general, when the PDF takes a bimodal shape, the first-order moments cease to act as suitable quantities to characterize the typical behavior of swimmers. One can instead turn to quantities expressed in terms of the higher-order moments that may reflect the bimodal nature of the PDF. We make use of the notion of the Binder cumulant, commonly used in the context of symmetry-breaking equilibrium phase transitions (as, e.g., occurring within Ising models Binder (1981) with a characteristic uni- to bimodal shape change in the corresponding Boltzmann weight of the mean spin on a lattice). It is defined as

[TABLE]

where $\mu_{2}$ and $\mu_{4}$ are the second- and fourth-order central moments, respectively, of the swimmer orientation angle, $\theta$, at the walls.

For a non-periodic random variable running over the whole real axis, the central moments can be calculated using standard textbook definitions; hence, the ratio $\mu_{4}/\mu_{2}^{2}$ will coincide with the standard kurtosis $\kappa$, and the Binder cumulant will be expressed in terms of the excess kurtosis $\kappa_{e}=\kappa-3$ as $U_{4}=-\kappa_{e}/3$. In this case, the Binder cumulant vanishes ($U_{4}=0$) for a unimodal Gaussian PDF. It takes non-zero values when the PDF is bimodal and, in the limiting case of a distribution consisting of two delta-function peaks, it gives $U_{4}=2/3$.

In the present case, however, the angular random variable $\theta$ is defined over a circle and the PDF is periodic in $\theta$, i.e., $\tilde{\Psi}(\tilde{y},\theta)=\tilde{\Psi}(\tilde{y},\theta+2\pi)$. Naively applying the standard definitions to calculate the central moments over a full period leads to results that generally depend on the choice of the period. This problem can be circumvented systematically by resorting to directional statistics Mardia and Jupp (2000). However, for the sake of simplicity, we take a heuristic and more intuitive route to construct an interval-independent Binder cumulant as follows. For an arbitrary choice of angular interval $[\theta_{1},\theta_{1}+2\pi)$ on the top or bottom wall, we determine the majority peak location $\theta_{\mathrm{max}}$ inside the interval. We evaluate the mean orientation angle, $\overline{\theta}$, by taking a symmetric interval equal to a full period around $\theta_{\mathrm{max}}$. This step is performed because the location of the peak varies and shifts toward smaller (larger) angles on the top (bottom) wall upon increasing the flow Péclet number. The result for $\overline{\theta}$ is then used to evaluate the central moments $\mu_{n}=\overline{(\theta-\overline{\theta})^{n}}$ by taking a symmetric interval equal to a full period around $\overline{\theta}$. The resulting central moments and the Binder cumulant will be independent of $\theta_{1}$. Using this procedure, the Binder cumulant can be shown to have the following limiting values: $U_{4}=2/5$ for $Pe_{f}=Pe_{s}=0$, where the distribution function is uniform, $U_{4}\rightarrow 0$ for $Pe_{f}=0$ and $Pe_{s}\gg 1$, where the distribution takes a sharp unimodal shape, $U_{4}\rightarrow 1/3$ for $Pe_{s}=0$ and $Pe_{f}\gg 1$, where the distribution takes a bimodal shape with two sharp peaks at an angular separation equal to half the period (note that $U_{4}$ in this latter case is half the limiting value it takes for a non-periodic random variable mentioned above).

The computed Binder cumulant is shown in Fig. 9 as a function of the flow Péclet number for different values of the swim Péclet number. In the absence of self-propulsion ($Pe_{s}=0$, green symbols), $U_{4}$ decreases monotonically with $Pe_{f}$ from a value close to 2/5 (uniform PDF) on the one end to 1/3 in the other (two sharp peaks), in accordance with the expected limiting values. For a finite $Pe_{s}$, the morphological variations of the PDF from a unimodal to a bimodal shape is clearly indicated by a drop in $U_{4}$ at intermediate $Pe_{f}$, while $U_{4}$ stays close to its limiting values on the two ends (zero for unimodal PDF at small $Pe_{f}$ and 1/3 for two sharp peaks at large $Pe_{f}$). The larger the swim Péclet number, the larger is the drop in the Binder cumulant at its minimum. The latter occurs at a swim Péclet number nearly equal (within our margin of error) to the value representing the onset of population splitting, indicating that this onset can also be identified by a large negative Binder cumulant.

We note that excess kurtosis, $\kappa_{e}$ (which, as noted previously, is proportional to the Binder cumulant, $U_{4}=-\kappa_{e}/3$, for non-periodic random variables) has long been considered as a measure for unimodality versus bimodality in statistics literature, where its usefulness and limitations have also been discussed extensively Darlington (1970); Moors (1986); Balanda and Macgillivray (1988). In the case of two-point distributions (equivalent to the two-state model introduced in Appendix B), the Binder cumulant can be shown to take large negative values as the distribution becomes less symmetric (i.e., as the probability density at one of the points or, equivalently, the bimodal ratio, $R_{\Psi}$, is decreased), and, eventually, tend to negative infinity as the onset of bimodality is approached ($R_{\Psi}\rightarrow 0$) Darlington (1970). Although such a behavior is consistent with the trend seen in Fig. 9, further attempt to analyze the behavior of the shown Binder cumulant requires an account of the cyclic nature of the orientation angle, which goes beyond the scope of the present work.

III.7 Phase diagram

In Fig. 10, we present a phase diagram in terms of the inverse flow Péclet number $1/Pe_{f}$, the swim Péclet number $Pe_{s}$, and the third dimensionless parameter of the system $\xi$ (see Eq. (11)). The plots in this figure are shown in log-log scale and symbols represent the inverse values of the flow Péclet number at the onset of population splitting, $Pe_{f}^{\ast}$ (computed either directly by studying the $\tilde{\Psi}$-profiles or by examining the Binder cumulant plots, with closely matching results), as a function of $Pe_{s}$ for fixed $\xi=0.041,0.41,0.82$ in the main set, and as a function of $\xi$ for fixed $Pe_{s}=0.1$ and 1 in the inset. Therefore, the area below (above) each curve represents the bimodal (unimodal) regime. As noted in Section II.3, the base value of $\xi=0.82$ was used to enable a quantitive comparison between our results and those in the case of an imposed Poiseuille flow Ezhilan and Saintillan (2015), as discussed throughout the previous sections. Figure 10 presents our results for the population-splitting behavior of swimmers over a much wider range of values for $\xi$. The curves that indicate the boundary between the uni- and bimodal regimes appear to take approximately linear forms in the log-log representation of Fig. 10. This can suggest possible scaling-like behavior for $Pe_{f}^{\ast}$, especially in the regime $Pe_{s}<1$; in the regime of parameters explored here, the exponent,

[TABLE]

appears to be weakly dependent on $\xi$ and given roughly as $\alpha\simeq 1$ for $Pe_{s}<1$.

We also note that the lines separating the unimodal and bimodal regimes in the phase diagram can be identified as lines of discontinuous transitions in analogy with the standard terminology in the context of equilibrium phase transitions Kardar (2007). To further elucidate this point, we consider the behavior of the loci $\theta_{\ast}$ of the peaks of the PDF (see Fig. 6), corresponding to the most probable orientations of swimmers near the walls. We also consider both positive and negative shear rates and, thus, allow $Pe_{f}$ to take both positive and negative values, with the negative values corresponding to a situation where the direction of the fluid flow in the channel is reversed.

In order to show a wide range of negative and positive values for $Pe_{f}$ in a single frame, we plot $\theta_{\ast}$ on the top wall as a function of the redefined variable $\Pi_{f}={\mathrm{sgn}}(Pe_{f})\left(\log_{10}|Pe_{f}|^{-1}+C\right)$ with $C=4$ in Fig. 11. This change of (abscissa) variable enables a plot on a logarithmic scale with the branches, representing negative (positive) values of $Pe_{f}$, appearing separately on the left (right) half of the frame. Note that, for this latter purpose, any positive value of $C$ can be used, provided that $\log_{10}|Pe_{f}|^{-1}+C$ remains positive over the range plotted for $Pe_{f}$; in the figure, we vary $|Pe_{f}|$ from $4\times 10^{-3}$ up to $4\times 10^{3}$ and we have thus conveniently set $C=4$. Also, in the mentioned representation, large values of $|Pe_{f}|$ are mapped to small values for $|\Pi_{f}|$, and vice versa. The narrow region around $\Pi_{f}=0$ corresponds to excessively high shear rates $|Pe_{f}|\gg 1$ that cannot be treated numerically. In spite of this, and given that the states with $Pe_{f}=\pm\infty$ would be indistinguishable, one can formally extrapolate the curves shown in Fig. 11 down to $\Pi_{f}=0$ from both sides, in which case the right and left branches will be joined. It should be pointed out that the current approach based on low-Reynolds-number hydrodynamics Purcell (1977); Happel and Brenner (1983) breaks down at very large $Pe_{f}$, and the region around $\Pi_{f}=0$ is expected to remain outside the regime of applicability of our results.

Figure 11 is obtained for fixed $Pe_{s}=1$ and $\xi=0.82$ but it clearly shows the characteristic features discussed in the preceding sections: $\theta_{\ast}$ is single-valued and tends to $\pi/2$ when $|\Pi_{f}|\rightarrow\infty$ or, equivalently, when $|Pe_{f}|\rightarrow 0$; $\theta_{\ast}$ remains single-valued in the unimodal regime above (below) a threshold value for $|\Pi_{f}|$ ($|Pe_{f}|$), and becomes multi-valued with two separate branches in the bimodal regime below (above) that threshold. The threshold value of $|\Pi_{f}|$, corresponding to the transition point separating the uni- and bimodal regimes, is indicated by a dashed vertical line in Fig. 11 (occurring at $Pe_{f}\simeq 4.1$ for the given parameter values). The discontinuous nature of the transition is thus evident from this graph.

The lower and upper branches for positive shear rate, corresponding to $\Pi_{f}>0$ (right half of the plot), give the most probable swimmer orientation near the (top) wall. The lower (upper) branch corresponds to the majority (minority) population and can be interpreted as the globally (locally) most probable state. The situation is the exact opposite for negative shear rate, corresponding to $\Pi_{f}<0$ (left half of the plot).

The results in Fig. 11 may thus appear as indicating the possibility of a hysteresis-like behavior for the most probable orientation of swimmers with respect to the flow Péclet number in a channel, where the direction of the Couette flow can be reversed. This kind of behavior might indeed be possible in reality, that is, swimmers in a majority population, oriented in the direction of a sufficiently large (positive) shear rate ($\theta_{\ast}\simeq 0$), may resist flipping at once to form a majority population in the opposite direction ($\theta_{\ast}\simeq\pi$) in accordance with the reversed flow. It is important to note, however, that interpretation of our results in Fig. 11 in terms of a hysteresis-like effect and, in general, any analogies with equilibrium phase transitions, should be done with caution (see Section IV).

IV Conclusion and discussion

We present a thorough study of the effects of an externally imposed Couette flow on the steady-state behavior of a dilute suspension of active (self-propelled) rodlike swimmers in a planar channel with rigid no-slip walls. Using a continuum description for non-interacting active rods Ezhilan and Saintillan (2015), we present a detailed investigation of the interplay between the imposed shear and the intrinsic self-propulsion of the swimmers in the channel.

In the self-propulsion-dominated regime, the swimmers form a single population, strongly accumulated near the channel walls and swimming in the down- and upstream directions near the top and bottom walls, respectively, in accordance with the directions intuitively implied from the linear structure of the Couette flow. Upon increasing the shear rate, the swimmers become more strongly aligned with the flow, resulting in larger magnitudes for the mean parallel-to-flow component of the orientation vector of the down- (up-) streaming swimmers near the top (bottom) wall. This behavior, however, holds only up to a threshold shear rate (or flow Péclet number, $Pe_{f}$), beyond which the shear-induced torque acting on the particles becomes strong enough to reverse the swimming direction of a finite fraction of swimmers from the direction dictated by the flow; this fraction of swimmers thus forms a minority, but statistically significant and macroscopic, population of particles swimming in the opposite direction to the majority population. Near the top wall, the conversion is from the downstream (majority) to upstream (minority), and vice versa near the bottom wall. This behavior, which we refer to as the population splitting phenomenon, corresponds to a morphological uni- to bimodal change in the joint (position-orientation) probability distribution function (PDF) of the swimmers.

We thus argue that the standard picture for the near-wall behavior of swimmers in a shear flow, in which the swimmers are perceived, based on their mean parallel-to-flow orientation, to swim up- or downstream may be inaccurate and misleading, as it implies overall up- or downstream swimming near the walls, missing the fact that beyond the onset of population splitting, swimmers form two discernible, oppositely swimming, statistical populations. Consequently, while one might expect stronger imposed flow (larger $Pe_{f}$) to lead to larger magnitudes for the parallel-to-flow component of the swimmer orientation, this latter quantity increases only up to the threshold value of $Pe_{f}$ and takes the reverse path of decreasing in magnitude as the population splitting takes place and develops further with increasing $Pe_{f}$. The inadequacy of the mean swimmer orientation to describe the ‘typical’ behavior of swimmers is nothing but a direct consequence of the bimodal nature of the swimmer orientation, necessitating the knowledge of higher-order moments of the PDF of swimmer orientation. We show that a conveniently defined Binder cumulant (used commonly in the context of equilibrium phase transitions Binder (1981)) can identify the uni- to bimodal transition of the PDF (or, the population splitting of swimmers) by exhibiting a global minimum of large negative value at the flow Péclet number corresponding to the onset of this transition. It is also worthwhile to note that in the limits, where the imposed shear (self-propulsion) effects are completely dominant, the PDF tends to the one found in the non-active (no-shear) case with two equal peaks (one symmetric peak); hence, assuming these two limiting cases per se, a uni- to bimodal transition, going through a population-splitting threshold for swimmer orientation, could be envisaged. Our analysis places this intuitive picture (which, however and to the best of our knowledge, remained unexplored) in a systematic framework and elucidates its nature.

We map out a phase diagram in terms of the flow and swim Péclet numbers and show that the uni- and bimodal regimes are separated by a discontinuous transition line. Our results also suggest the possibility of a hysteresis-like behavior for the most probable orientation of swimmers with the flow Péclet number in a channel, where the direction of the flow can be reversed. It is important to note, however, that a hysteresis-like behavior would be possible in situations where the time-scale of swimmer response is much longer than that of flow reversal (in equilibrium phase transitions, such a mechanism is provided by the inter-particle couplings that lead to free energy barriers between stable and/or meta-stable states of the system). Addressing this issue thus requires a more sophisticated modeling approach to incorporate several key factors such as (time-dependent) swimmer response to variations in the external flow, interactions between swimmers and interactions between swimmers and the walls, which have been neglected in our current analysis. Hence, interpretation of our results in terms of a hysteresis-like effect and making any analogies with equilibrium phase transitions should be done with caution.

Our analysis is based on a few simplifying assumptions such as neglect of hydrodynamic interactions, either between the active particles themselves, or between the particles and the channel walls. In dense swimmer suspensions, such inter-particle interactions, including also steric interactions, would have to be included (see, e.g., Refs. Wioland et al. (2013); Lushi et al. (2014); Saintillan and Shelley (2008a, b, 2013); Ezhilan and Saintillan (2015); Li and Ardekani (2014); Hernandez-Ortiz et al. (2005); Underhill et al. (2008); Hernandez-Ortiz et al. (2009); Zöttl and Stark (2013)). Wall hydrodynamic interactions have been suggested to have a prominent role in the well-established migration of active particles toward walls (see, e.g., Refs. Berke et al. (2008); Lauga et al. (2006); Spagnolie and Lauga (2012); Li and Ardekani (2014); Mathijssen et al. (2016a, b, c); Zöttl and Stark (2013)), while there are several studies that suggest otherwise (see, e.g., Refs. Elgeti and Gompper (2015); Costanzo et al. (2012); Ezhilan and Saintillan (2015); Li et al. (2011); Li and Tang (2009) and references therein). This work is thus an example among many simulation works that demonstrate wall accumulation of swimmers, while wall hydrodynamic interactions are neglected.

Another interesting problem is the role of particle shape and its finite aspect ratio. Here, we have modeled the swimmers as thin rods (needles) and, as such, our analysis ignores hydrodynamic effects, such as oscillatory motion of particle orientation in Jeffery orbits Jeffery (1922), arising as a result of finite particle aspect ratio. This and other factors such as swimmer chirality are among the directions we plan to pursue in our future research.

Acknowledgements.

M.K., J.A. and A.F. acknowledge the School of Physics, IPM, Tehran, for hospitality and support during the completion of this work. A.N. acknowledges Iran Science Elites Federation for partial support. We thank M.R. Shabanniya for useful discussions.

Appendix A Numerical methods

The steady-state properties of swimmers in the channel follow from the solution of Eq. (7) that can be evaluated numerically in the rescaled domain $-1\leq\tilde{y}\leq 1$ and $0\leq\theta<2\pi$ using the boundary conditions (8), the constraint (13) and the fact that $\tilde{\Psi}(\tilde{y},\theta=0)=\tilde{\Psi}(\tilde{y},\theta=2\pi)$. All numerical simulations in this work have been performed using COMSOL Multiphysics v5.2. In all cases reported here, the relative numerical error was determined by direct validation of the numerical solutions to be less than $1\%$. We used structured (square-shaped) mesh elements and cubic discretization in all simulations. A mesh-independency study was performed for the test case of swimmer density on the top wall ($\tilde{y}=1$). We observed the results converging, with an error margin below $1\%$, with quadratic discretization. For our case (cubic), the convergence rate was, expectedly, even better.

Appendix B Mean near-wall swimmer orientation

Here, we show how the non-monotonic (monotonic) behavior of the mean parallel-to-flow (perpendicular-to-flow) orientation component $\langle\tilde{p}_{x}\rangle$ ($\langle\tilde{p}_{y}\rangle$) as a function of $Pe_{f}$ (Fig. 4) can be understood on a semi-quantitative basis using the population spitting phenomenon. We focus on the top channel wall at $\tilde{y}=1$ and fisrt consider the unimodal regime of small to intermediate $Pe_{f}$ at a fixed $Pe_{s}>0$. When there is no imposed shear ($Pe_{f}=0$), the near-wall PDF $\tilde{\Psi}(\tilde{y}=1,\theta)$ of swimmer orientation is unimodal and symmetric, being centered at $\theta_{\ast}\simeq\pi/2$. It thus follows that $\langle\tilde{p}_{x}\rangle=\langle\cos\theta\rangle=0$, while $\langle\tilde{p}_{y}\rangle=\langle\sin\theta\rangle$ is finite and depends on the exact value of $Pe_{s}$ (see Figs. 3 and 4). As noted in Section III.4, upon increasing $Pe_{f}$ and before the population splitting sets in, the unimodal distribution for swimmer orientation becomes increasingly more asymmetric and its peak location shifts to smaller angles, i.e., $\theta_{\ast}$ tends to zero. Hence, $\langle\tilde{p}_{x}\rangle$ increases in magnitude, while $\langle\tilde{p}_{y}\rangle$ decreases.

We now turn to the bimodal regime, where swimmer orientation is sharply peaked around the orientation angles $\theta_{\ast}=\theta_{u}$ (for upstream minority) and $\theta_{d}$ (for downstream majority) with an angular separation of $\theta_{u}-\theta_{d}\simeq\pi$ (in fact, $\theta_{u}\simeq\pi$ and $\theta_{d}\simeq 0$). A rough estimate for $\langle\tilde{p}_{x}\rangle$ and $\langle\tilde{p}_{y}\rangle$ can be obtained by adopting an idealized two-state model, in which the swimmer orientation angle on the top wall can only take one of the two values $\{\theta_{u},\theta_{d}\}$, representing the minority and majority populations, respectively, with no angular dispersion around these peak values. Hence, we find

[TABLE]

where $R_{\Psi}$ is the bimodal ratio at the top wall (Section III.5). Thus, within the mentioned two-state model, $\langle\tilde{p}_{y}\rangle$ vanishes as $\theta_{d}$ tends to zero, which, together with what we found above for the unimodal regime, suggests that $\langle\tilde{p}_{y}\rangle$ should monotonically decrease to zero as $Pe_{f}$ is increased in accord with our numerical results (Fig. 4b). The situation is found to be different for $\langle\tilde{p}_{x}\rangle$; since $\theta_{d}$ is close to zero in the bimodal regime, $\cos\theta_{d}\simeq 1$ and $\langle\tilde{p}_{x}\rangle$ can be further approximated as $\langle\tilde{p}_{x}\rangle\simeq(1-R_{\Psi})/(1+R_{\Psi})$, which finds its maximum near the onset of population splitting ($R_{\Psi}=0$) and decreases to zero as the minority population is strengthened ($R_{\Psi}\rightarrow 1$).

Our foregoing analysis clearly demonstrates that the increasing (decreasing) regime for $\langle\tilde{p}_{x}\rangle$ as a function of $Pe_{f}$, as shown in Fig. 4a, is associated with the unimodal (bimodal) regime, respectively. Therefore, the non-monotonic behavior of $\langle\tilde{p}_{x}\rangle$ signifies the underlying population splitting phenomenon. This same mechanism can explain the non-monotonic near-wall behavior of $\langle\tilde{p}_{x}\rangle$ found in the case of an imposed Poiseuille flow Ezhilan and Saintillan (2015), as the latter represents a linear flow profile near the walls. The foregoing argument can be used for swimmers on the bottom wall by replacing minority with majority, upstream with downstream, and $R_{\Psi}$ with $1/R_{\Psi}$.