Shear induced migration of microswimmers in pressure-driven channel flow

TL;DR

This paper models how microswimmers migrate in pressure-driven channel flow, revealing regimes of wall or centerline accumulation depending on shear rate and swimmer shape, and validates findings with simulations and experiments.

Contribution

It introduces a multiple scales method to derive a drift-diffusion model for swimmer concentration, explaining shear-induced migration behaviors and recent experimental results.

Findings

Swimmers with high aspect ratio migrate towards walls at moderate shear.

At high shear, swimmers tend to concentrate at the channel center.

Low-aspect-ratio swimmers show a Pe-independent profile at high shear.

Abstract

We study the shear induced migration of microswimmers (primarily, active Brownian particles or ABP's) in a plane Poiseuille flow. For wide channels characterized by $U_b/HD_r \ll 1$, the separation between time scales characterizing the swimmer orientation dynamics (of O($D^{-1}_r$)) and those that characterize migration across the channel (of O($H^{2}D_r/U^{2}_b$)), allows for use of the method of multiple scales to derive a drift-diffusion equation for the swimmer concentration profile; here, $U_b$ is the swimming speed, $H$ is the channel half-width, and $D_r$ is the swimmer rotary diffusivity. The steady state concentration profile is a function of the P\'eclet number, $Pe = U_{f}/(D_r H)$ ($U_f$ being the channel centerline velocity), and the swimmer aspect ratio $\kappa$. Swimmers with $ \kappa \gg 1$ (with $ \kappa \sim$ O(1)), in the regime $1 \ll \textit{Pe} \ll \kappa^3$ ($Pe\sim$ O(1)), migrate towards the channel walls, corresponding to a high-shear trapping behavior. For $Pe \gg \kappa^3 $ ($Pe \gg $ 1 for $\kappa \sim$ O(1)), however, swimmers migrate towards the centerline, corresponding to a low-shear trapping behavior. Interestingly, within the low-shear trapping regime, swimmers with $\kappa < 2$ asymptote to a $Pe$-independent concentration profile for large $Pe$, while those with $\kappa \geq 2$ exhibit a `centerline-collapse' for $Pe \to \infty$. The prediction of low-shear-trapping, validated by Langevin simulations, is the first explanation of recent experimental observations [Barry $\textit{et al}$. (2015)]. We organize the high-shear and low-shear trapping regimes on a $Pe-\kappa$ plane, thereby highlighting the singular behavior of infinite-aspect-ratio swimmers.

Equations

$$\frac{\partial \Omega }{\partial t}+U_s{\nabla }_{\mathbf{x}}\cdot (\Omega \mathbf{p})+{\nabla }_p\cdot (\dot{\mathbf{p}}\Omega )+\frac{1}{\tau }\left[\Omega -\int d{\mathbf{p}}^{\prime }K(\mathbf{p}|{\mathbf{p}}^{\prime })\Omega ({\mathbf{p}}^{\prime })\right]-D_r{\nabla }_p^2\Omega =0,$$

$$\frac{\partial \Omega }{\partial t}+ϵ{\nabla }_{\mathbf{x}}\cdot (\Omega \mathbf{p})+Pe\dot{\gamma }{\nabla }_p\cdot (\widetilde{\dot{\mathbf{p}}}\Omega )+\frac{1}{\tau D_r}\left[\Omega -\int d{\mathbf{p}}^{\prime }K(\mathbf{p}/{\mathbf{p}}^{\prime })\Omega ({\mathbf{p}}^{\prime })\right]-{\nabla }_p^2\Omega =0,$$

$$\frac{\partial \Omega }{\partial t_1}+{ϵ}^2\frac{\partial \Omega }{\partial t_2}+ϵ{\nabla }_{\mathbf{x}}\cdot (\Omega \mathbf{p})+P{e}_z{\nabla }_p\cdot (\widetilde{\dot{\mathbf{p}}}\Omega )+\frac{1}{\tau D_r}\left[\Omega -\int d{\mathbf{p}}^{\prime }K(\mathbf{p}/{\mathbf{p}}^{\prime })\Omega ({\mathbf{p}}^{\prime })\right]-{\nabla }_p^2\Omega =0,$$

$$\text{O}(1)$$

$$\text{O}(ϵ)$$

$$\text{O}({ϵ}^2)$$

$$\frac{\partial F(z,t_2)}{\partial t_2}=\frac{\partial }{\partial z}\left(D_{zz}(z;Pe,\kappa )\frac{\partial F(z,t_2)}{\partial z}-V_z(z;Pe,\kappa )F(z,t_2)\right),$$

$$F(z)=\aleph \exp \left({\int }_{-1}^{z}d{z}^{\prime }\frac{V_z(z^{\prime })}{D_{zz}(z^{\prime })}\right),$$

$$\aleph =\frac{1}{{\int }_{-1}^{1}dz\exp \left({\int }_{-1}^{z}d{z}^{\prime }\frac{V_z(z^{\prime })}{D_{zz}(z^{\prime })}\right)}$$

$${\nabla }_p\cdot (\widetilde{\dot{\mathbf{p}}}G)$$

$${\mathcal{L}}_y$$

$${\mathcal{L}}_z$$

$${\nabla }_p\cdot (\widetilde{\dot{\mathbf{p}}}G)=-3\kappa \sqrt{\frac{2\pi }{15}}(Y_2^{-1}-Y_2^1)G$$

$$G(\mathbf{p})=\sum _{l=0}^{\infty }\sum _{m=-l}^l{a}_{l,m}{Y}_l^m(\mathbf{p})$$

$${\nabla }_p\cdot (\widetilde{\dot{\mathbf{p}}}G)=\sum _{l=0}^{\infty }\sum _{m=-l}^l{a}_{l,m}{\nabla }_p\cdot (\widetilde{\dot{\mathbf{p}}}{Y}_l^m)=\sum _{l=0}^{\infty }\sum _{m=-l}^l{a}_{l,m}[-3B\sqrt{\frac{2\pi }{15}}(Y_2^{-1}-Y_2^1)Y_l^m$$

$$+\frac{i}{2}{\mathcal{L}}_y|Y_l^m>\left\{1+B\left(-\frac{2}{3}\sqrt{\pi }{Y}_0^0+\frac{8}{3}\sqrt{\frac{\pi }{5}}{Y}_2^0\right)\right\}+B\sqrt{\frac{2\pi }{15}}\left(Y_2^{-1}+Y_2^1\right){\mathcal{L}}_z|Y_l^m>],$$

$$\sum _{l=0}^{\infty }\sum _{m=-l}^l{a}_{l,m}$$

$$+BP{e}_z\sqrt{\frac{2\pi }{15}}\left(Y_2^{-1}+Y_2^1\right){\mathcal{L}}_z|Y_l^m>+l(l+1)Y_l^m$$

$$+\frac{1}{\tau D_r}\left(Y_l^m-\sum _{s=0}^{\infty }\sum _{t=-s}^s\int d{\mathbf{p}}^{\prime }{A}_s\frac{4\pi (-1{)}^t{Y}_s^t(\mathbf{p})Y_s^{-t}({\mathbf{p}}^{\prime })Y_l^m({\mathbf{p}}^{\prime })}{(2s+1)}\right)]=0.$$

$$A_n=\frac{\int d(\mathbf{p}\cdot {\mathbf{p}}^{\prime })\beta \exp (\beta (\mathbf{p}\cdot {\mathbf{p}}^{\prime }))/(4\pi \sinh \beta )P_n(\mathbf{p}\cdot {\mathbf{p}}^{\prime })}{\int d(\mathbf{p}\cdot {\mathbf{p}}^{\prime })P_n(\mathbf{p}\cdot {\mathbf{p}}^{\prime })P_n(\mathbf{p}\cdot {\mathbf{p}}^{\prime })}\text{for}n\ge 1.$$

$$i{\mathcal{L}}_y|Y_l^m>$$

$${\mathcal{L}}_z|Y_l^m>$$

$$\sum _{l=0}^{\infty }\sum _{m=-l}^l{a}_{l,m}$$

$$\left\{1+B\left(-\frac{2}{3}\sqrt{\pi }{Y}_0^0+\frac{8}{3}\sqrt{\frac{\pi }{5}}{Y}_2^0\right)\right\}+P{e}_{z}B\sqrt{\frac{2\pi }{15}}\left(Y_2^{-1}+Y_2^1\right)m{Y}_l^m$$

$$+\frac{Y_l^m}{\tau D_r}\left(1-\frac{A_l}{4\pi }\right)+l(l+1)Y_l^m]=0,$$

$$\sum _{l=0}^{\infty }\sum _{m=-l}^l{a}_{l,m}$$

$$+P{e}_z\left(\frac{B}{3}\right)\sqrt{\frac{4\pi }{5}}\left[k_{lm}^{\prime }<r,s|Y_2^0|l,m+1>-k_{lm}^{\prime \prime }<r,s|Y_2^0|l,m-1>\right]$$

$$+P{e}_{z}B\sqrt{\frac{2\pi }{15}}\left[(m-3)<r,s|Y_2^{-1}|l,m>+(m+3)<r,s|Y_2^1|l,m>\right]$$

$$\displaystyle+\delta_{lr}\delta_{ms}\left[\frac{1}{\tau D_{r}}\left(1-\frac{4\pi A_{l}}{(2l+1)}\right)+l(l+1)\right]\Bigr{\}}=0.$$

$$<r,s|Y_a^b|l,m>=\int d\mathbf{p}{Y}_r^{s^{\ast }}(\mathbf{p})Y_a^b(\mathbf{p})Y_l^m(\mathbf{p}),$$

$$\int d\mathbf{p}\Omega =\int d\mathbf{p}({\Omega }_0+ϵ{\Omega }_1+{ϵ}^2{\Omega }_2+\cdots )=F,$$

$$a_{0,0}=1/\sqrt{4\pi }.$$

$$\sum _{l=1}^{\infty }\sum _{m=-l}^l$$

$$+P{e}_z\left(\frac{B}{3}\right)\sqrt{\frac{4\pi }{5}}\left[a_{lm}^{\prime }<r,s|Y_2^0|l,m+1>-a_{lm}^{\prime \prime }<r,s|Y_2^0|l,m-1>\right]$$

$$+P{e}_{z}B\sqrt{\frac{2\pi }{15}}\left[(m-3)<r,s|Y_2^{-1}|l,m>+(m+3)<r,s|Y_2^1|l,m>\right]$$

$$\displaystyle+\delta_{lr}\delta_{ms}\left[\frac{1}{\tau D_{r}}\left(1-\frac{A_{l}4\pi}{(2l+1)}\right)+l(l+1)\right]\Bigr{\}}$$

$$=P{e}_{z}B\sqrt{\frac{3}{10}}\left[<r,s|Y_2^{-1}|0,0>-<r,s|Y_2^1|0,0>\right].$$

$$P{e}_z{\nabla }_p\cdot (\widetilde{\dot{\mathbf{p}}}{\Omega }_1)+\frac{1}{\tau D_r}\left[{\Omega }_1-\int d{\mathbf{p}}^{{}^{\prime }}K(\mathbf{p}/{\mathbf{p}}^{{}^{\prime }}){\Omega }_1({\mathbf{p}}^{{}^{\prime }})\right]-{\nabla }_p^2{\Omega }_1=-\cos \theta G\frac{\partial F}{\partial z}-\cos \theta F\frac{\partial G}{\partial z}.$$

Full text

1 Introduction

Biologically active microswimmers, a subset of the low-Reynolds-number dwellers, have motivated research spanning over several decades [Hancock (1953); Gray & Hancock (1955); Brennen & Winet (1977); Pedley & Kessler (1992); Koch & Subramanian (2011); Guasto et al. (2012); Marchetti et al. (2013)]. Owing to their intrinsic activity, suspensions of such microswimmers exhibit a host of interesting phenomena such as enhanced tracer diffusion [Wu & Libchaber (2000); Underhill et al. (2008); Leptos et al. (2009); Krishnamurthy & Subramanian (2015); Patteson et al. (2016); Stenhammar et al. (2017)], long-ranged orientational order [Saintillan & Shelley (2007, 2008); Underhill & Graham (2011); Stenhammar et al. (2017); Nambiar et al. (2019a)], negative viscosities [López et al. (2015); Bechtel & Khair (2017); Nambiar et al. (2017, 2019b); Takatori et al. (2014a, b); Takatori & Brady (2015)], collective motion/bacterial turbulence [Dombrowski et al. (2004); Underhill et al. (2008); Subramanian & Koch (2009); Subramanian et al. (2011); Wensink et al. (2012); Marchetti et al. (2013)], among others. The underlying mechanisms driving most of these phenomena may be understood from studying quiescent swimmer suspensions, or those subjected to homogeneous shear flows. However, nearly all biologically relevant problems involve motile microorganisms in complex flow environments [Guasto et al. (2012); Rusconi et al. (2014a)]. Complex flows, with or without boundaries, are known to alter the nutrient landscape [Taylor & Stocker (2012)] influence swimmer rheotaxis [Stocker et al. (2006); Fu et al. (2012)], result in shear induced migration [Rusconi et al. (2014b); Bearon & Hazel (2015); Barry et al. (2015)], in turn leading to shear-banding instabilities [Guo et al. (2018); Laxminarsimharao et al. (2018)], and possibly, influencing biofilm formation [Rusconi et al. (2010, 2014a); Kim et al. (2014)]. The study in this paper examines the orientation dynamics and transport of microorganisms in one of the simplest inhomogeneous shearing flows - the pressure-driven flow between a pair of parallel plates.

Microorganisms, bacteria such as E. coli and B. subtilis, algae such as C. reinhardtii [Berg (2008); Elgeti et al. (2015); Guasto et al. (2012)], and even artificial swimmers, exhibit an intrinsic stochasticity in their swimming motion owing to a tendency to reorient as they swim. The intrinsic stochasticity may be biologically motivated [Berg (2008); Elgeti et al. (2015); Guasto et al. (2012)], or a consequence of reaction with the solvent [Ebbens et al. (2014); Moran & Posner (2017)] or driven by an external field [Chaturvedi et al. (2010); Fischer & Ghosh (2011); Peyer et al. (2013)]. The reorientations occur either as discrete finite-amplitude events leading to run-and-tumble dynamics, or as continuous small-amplitude events leading to active Brownian dynamics. Regardless of its origin or detailed nature, the interplay between the stochastic dynamics of the swimmer and an imposed inhomogeneous shearing flow has been shown to result in migration patterns [Chilukuri et al. (2014); Rusconi et al. (2014b); Chilukuri et al. (2015); Barry et al. (2015); Ezhilan & Saintillan (2015); Bearon & Hazel (2015); Sokolov & Aranson (2016); Sokolov et al. (2018)] that stand in sharp contrast to those known for suspensions of passive particles [Gadala-Maria & Acrivos (1980); Leighton & Acrivos (1987); Koh et al. (1994); Nott & Brady (1994); Nitsche & Hinch (1997); Strednak et al. (2018)]. Experiments and simulations on suspensions of passive particles in pressure-driven channel flow reveal, irrespective of particle geometry (rigid spheres and/or slender fibers), migration towards the channel center [Koh et al. (1994); Strednak et al. (2018); Nott & Brady (1994)]. Microswimmers on the other hand, depending on swimmer geometry and flow characteristics, exhibit migration both towards the wall (high-shear trapping) and channel center (low-shear trapping) [Rusconi et al. (2014b); Barry et al. (2015); Bearon & Hazel (2015)]. Recent experiments in a pressure-driven microchannel flow have found slender bacteria to migrate towards the walls [Rusconi et al. (2014b)], whereas relatively round algae have been found to migrate towards the wall or the centerline depending on the particular species and the prevailing shear rate [Barry et al. (2015)]. Earlier computations by Ezhilan & Saintillan (2015), for infinitely slender swimmers have only confirmed the existence of high-shear trapping over a range of $Pe$. Here, $Pe=U_{f}/(D_{r}H)$ is the rotary Péclet number, the shear rate ($U_{f}/H$) measured in units of the inverse rotary diffusivity ($D^{-1}_{r}$); $U_{f}$ and $H$ are, respectively, the centerline velocity and channel half-width, and $D_{r}$ the swimmer diffusivity. An independent computational investigation by Bearon & Hazel (2015) revealed regimes of both high and low-shear trapping, but this was in contrast to the authors’ analysis which, similar to Ezhilan & Saintillan (2015), predicted high-shear trapping regardless of swimmer aspect ratio ($\kappa$); this apparent discrepancy remains unexplained. There have been additional simulations [Chilukuri et al. (2014, 2015)] examining the spatial distribution, and the dispersion in the flow direction, of microswimmers (modeled as dumbbells), in pressure driven channel flow. The channel in these cases was, however, quite narrow, and as a result, confinement effects, together with wall-mediated hydrodynamic interactions were dominant. In contrast, the experiments mentioned above used channels much wider than the single swimmer dimension.

From the above summary, it is clear that a comprehensive understanding of the interplay between swimmer shape and shear induced migration in a channel flow is lacking. In particular, there has been no attempt to interpret and explain the observations of Rusconi et al. (2014b) and Barry et al. (2015) based on a common theoretical frame work. In this work, we use both theory and Langevin simulations, to comprehensively characterize shear induced migration of microswimmers in the $Pe-\kappa$ plane, rationalizing the aforementioned observations in the process. Two distinct classes of swimmer concentration profiles are identified in this parameter plane, corresponding to high-shear and low-shear trapping behavior. In §2, the theoretical framework and the simulation protocol adopted are described. In §2.1, beginning with the equation governing the swimmer probability density in position-orientation space, and using the separation of the fast (O$(D_{r}^{-1})$) reorientation time scale and the much slower (O$(H^{2}D_{r}/U_{s}^{2})$) time scale for migration across the channel ($U_{s}$ here being the swimming speed), we derive a drift-diffusion equation for the cross-stream swimmer concentration profile with the aid of the method of multiple scales. This is followed by the description of the simulation scheme in §2.2; the scheme obtains the concentration profile by numerically integrating the Langevin equations of motion for the individual swimmers with periodic boundary conditions. The steady state swimmer concentration profiles obtained from the multiple scales analysis are presented in §3; high-shear trapping is discussed in §3.1 and low-shear trapping in §3.2. In each case comparisons are drawn with both experiments [Rusconi et al. (2014b); Barry et al. (2015)], and the Langevin simulations described in §2.2. In §3.3, the swimmer depletion index, a measure of the spatial inhomogeneity of the swimmer concentration, is plotted as a function of Pe for swimmers of different aspect ratios, highlighting both the singular behavior of infinitely slender swimmers, and the centerline-collapse that occurs for swimmers with $\kappa$ finite but greater than (approximately) 2, for $Pe\rightarrow\infty$, in the low-shear trapping regime. Next, the shear induced migration behavior is organized on the $Pe-\kappa$ plane, demarcating the low and high-shear trapping regimes. Finally, in §4, we present concluding remarks, that include scaling arguments for the threshold governing the transition from active to passive shear induced migration patterns, and directions for future work.

2 Theoretical framework and the Langevin simulation protocol

In §2.1, we use the method of multiple scales to derive the swimmer concentration profiles as a function of the transverse coordinate in plane Poiseuille flow. Next, in §2.2, we describe the scheme adopted to simulate a discrete system of swimmers in the same flow.

2.1 Theoretical framework: the method of multiple scales

The time evolution of the probability density $\Omega(\mathbf{x},\mathbf{p},t)$ for a dilute swimmer suspension subject to a shearing flow (see figure 1), is given by [Subramanian & Koch (2009)]:

[TABLE]

where $\mathbf{x}$ and $\mathbf{p}$ denote the swimmer position and orientation. The second term on the left hand side of (1) denotes spatial convection of the probability density owing to swimming with speed $U_{s}$, and the term involving $\dot{\mathbf{p}}$ denotes swimmer rotation by the ambient shear. The terms within brackets model a run-and-tumble process [Berg (1993); Subramanian & Koch (2009); Nambiar et al. (2017)], obeying Poission statistics with a mean run duration $\tau$; the kernel $K(\mathbf{p}|\mathbf{p}^{\prime})$ characterizes correlation between the pre ($\mathbf{p}^{\prime}$)- and post ($\mathbf{p}$)-tumble orientations, with $K(\mathbf{p}|\mathbf{p}^{\prime})=1/(4\pi)$ for random tumbles. The last term on the left side of (1) models the stochastic orientation change due to rotary diffusion, $D_{r}$ being the rotary diffusivity.

In the experiments of Rusconi et al. (2014b) and Barry et al. (2015), the ratio of the channel width ($H$) to the swimmer size ($L$) is approximately $40$. Thus, it is reasonable to assume that the swimmer senses a local linear flow, with its rotation governed by the Jeffery equation; therefore, $\dot{\mathbf{p}}=B\left[\mathbf{E}\cdot\mathbf{p}-\mathbf{p}(\mathbf{E}:\mathbf{p}\mathbf{p})\right]+\mathbf{W}\cdot\mathbf{p}$. Here, $B=(\kappa^{2}-1)/(\kappa^{2}+1)$ is the Bretherton constant, $E_{ij}=\dot{\gamma}(z)(\delta_{i1}\delta_{j3}+\delta_{i3}\delta_{j1})/2$ and $W_{ij}=\dot{\gamma}(z)(\delta_{i1}\delta_{j3}-\delta_{i3}\delta_{j1})/2$ are the strain rate and vorticity tensors, with $\dot{\gamma}(z)=du/dz$ being the shear rate at $z$, and $u(z)=U_{f}(1-(z/H)^{2})$ for the plane Poiseuille flow under consideration. In the above, $\delta_{ij}$ is the Kronecker delta, and the coordinate system used appears in figure 1. In using Jeffery’s equation above, we model the swimmers as equivalent spheroids with an aspect ratio $\kappa$ [Leal & Hinch (1971)], and it is via $\kappa$ that the swimmer geometry enters the analysis, determining the nature of the shear induced migration. This implicitly assumes an axisymmetric cross-section, which is a reasonable assumption, at least for bacteria [Darnton et al. (2007); Das & Lauga (2018)], on account of the rapid (counter)-rotation on time scales much shorter than those characterizing shear induced migration.

Non-dimensionalizing (1) using $D_{r}^{-1}$, $U_{f}$, and $H$ as scales for time, velocity and length, respectively, yields:

[TABLE]

where $\tilde{\dot{\mathbf{p}}}=\dot{\mathbf{p}}/\dot{\gamma}(z)$ is a non-dimensional rotation rate. In (2), $\epsilon$, Pe and $\kappa$ are dimensionless parameters that primarily determine the different regimes of shear induced migration. The parameter $\epsilon=U_{s}/D_{r}H$, defined as the ratio of the swimmer mean free path $U_{s}/D_{r}$ to $H$, may be regarded as a swimmer Knudsen number. Note that in earlier efforts [Takatori et al. (2014b); Takatori & Brady (2017)], $\epsilon^{-1}$, which may be written as $U_{s}H/D_{t}$, with $D_{t}\sim U_{s}^{2}/(6D_{r})$ being the translational diffusivity, has been interpreted as the swim Péclet number. The parameter $\tau D_{r}$ determines the roles of tumbling vis-a-vis rotary diffusion in the swimmer orientation dynamics. Active Brownian particles (ABP’s) correspond to $\tau D_{r}\rightarrow\infty$, and run-and-tumble particles (RTP’s) to $\tau D_{r}\rightarrow 0$ [Tailleur & Cates (2009); Saintillan (2010b)]; in the latter case, the swimmer Knudsen number is defined as $\epsilon=U\tau/H$ instead. Regardless of the particular value of $\tau D_{r}$, swimmers sample the channel cross-section diffusively for long times, the translational diffusivity being given by $D_{t}=(U^{2}_{s}/(6D_{r}))(\tau D_{r}/(0.5+\tau D_{r}))$ [Berg (1993); Koch & Subramanian (2011); Subramanian & Nott (2012)] in the general case; for $\tau D_{r}\to\infty$, $D_{t}=U^{2}_{s}/(6D_{r})$ as above. In the experiments of Rusconi et al. (2014b), $\epsilon\sim 0.23$ for the bacterium B. subtilis, and in Barry et al. (2015), $\epsilon\sim 0.29$ and 0.15 for the algal species Heterosigma and Dunaliella, respectively. For these small values of $\epsilon$, the time scales characterizing the swimmer orientation dynamics (of O($\tau$) or O($D^{-1}_{r}$)) and the diffusive sampling of the channel cross-section (of O($H^{2}/U^{2}_{s}\tau$) or O($H^{2}D_{r}/U^{2}_{s}$)) are well separated, their ratio being of O($\epsilon^{2}$). Further, since the experimental range of flow rates correspond to $U_{s}\ll U_{f}(\epsilon\ll Pe)$, we use the method of multiple scales to analyze the swimmer concentration profiles [Subramanian & Brady (2004); Nitsche & Hinch (1997); Kasyap & Koch (2014, 2012); Laxminarsimharao et al. (2018)] for small $\epsilon$, but with $Pe$ and $\kappa$ being arbitrary.

As a first step, we split the time derivative in (2) in terms of the fast intrinsic time scale $t_{1}=t$ and the slow diffusive time scale $t_{2}=\epsilon^{2}t$, resulting in $\epsilon$:

[TABLE]

where $Pe_{z}=Pe\dot{\gamma}(z)$ is the local Péclet number. Now, expanding $\Omega$ as $\Omega(z,\mathbf{p},t;\epsilon)=\Omega_{0}(z,\mathbf{p},t_{1},t_{2})+\epsilon\Omega_{1}(z,\mathbf{p},t_{1},t_{2})+\epsilon^{2}\Omega_{2}(z,\mathbf{p},t_{1},t_{2})+\ldots$, yields the following set of equations at successive orders in $\epsilon$:

[TABLE]

Here, we have written $\Omega_{0}$ as: $\Omega_{0}=F(z,t_{2})G(\mathbf{p};Pe_{z})$, with $\int d\mathbf{p}\Omega(z,\mathbf{p},t_{2})=F(z,t_{2})$ denoting the $z$-dependent (normalized) swimmer concentration profile. The aforementioned ansatz assumes a quasi-steady response ($t_{1}\rightarrow\infty$), with the swimmer orientation distribution having relaxed to a $Pe_{z}$-dependent local equilibrium; accordingly, the fast time derivative ($\partial/\partial t_{1}$) has been omitted in (4-2.1) . Writing $\Omega_{1}$ in (5) as $\Omega_{11}\partial F/\partial z+\Omega_{12}F$ on account of linearity, and using this form in (2.1) and integrating over orientation space, leads to a drift-diffusion equation for $F$:

[TABLE]

where the diffusivity $D_{zz}=-\int d\mathbf{p}\cos\theta\Omega_{11}$ and drift $V_{z}=\int d\mathbf{p}\cos\theta\Omega_{12}$. At steady state, for impermeable walls:

[TABLE]

where

[TABLE]

is determined from the normalization condition $\int_{-1}^{1}dzF(z)=1$. Determining $D_{zz}$ and $V_{z}$ requires knowledge of the orientation distribution as a function of $Pe$ and $\kappa$. For arbitrary $Pe$, the orientation distributions, at different orders in $\epsilon$, are written as a truncated series in spherical harmonics, with the resulting coefficients being obtained numerically [Nambiar et al. (2019a)] (see appendix A), while for small $Pe$, the orientation distribution is determined analytically via a perturbation expansion (see appendix B).

2.2 Langevin simulation scheme

The objective is to simulate a population of $N$ swimmers subject to a parabolic flow profile, with the purpose of determining the cross-stream swimmer concentration profile. Thus, there are three degrees of freedom per swimmer: two orientation degrees of freedom and one spatial (the swimmer $z$ coordinate). In other words, while we keep track of the swimmer orientation in three dimensions, only the cross-stream location $z$ of each swimmer is recorded. Following the experimental protocol of Rusconi et al. (2014b), the swimmers are uniformly distributed along $z$ at the initial time, with their orientations following an isotropic distribution. Thereafter, at each time step $\delta t$, the swimmer orientation changes on account of either rotary diffusion or run-and-tumble motion, and due to the rotation by the flow. Here, the rotary diffusion process is modeled using the method outlined in Ghosh et al. (2012), whereas, run-and-tumble motion is implemented in accordance with Poisson statistics [Berg (1993); Krishnamurthy & Subramanian (2015)]. The flow-induced rotation at any $z$ is determined using (2). Following the update in the swimmer orientations, their positions are updated using an $RK2^{nd}$ order scheme, while also accounting for periodic boundary conditions in the $z$-direction. To avoid artefacts arising from a jump in shear rate which the swimmer encounters while swimming across the cross-stream (periodic) boundaries, the flow profile is chosen to be a double Poiseuille flow, which implies that the spatial extent of the simulation domain is $z\in[0,4H]$. Nevertheless, we have verified that both single and double Poiseuille flow profiles yield the same steady state concentration profiles. Every result reported here is obtained from simulating a set of ten independent initial conditions, with a given run having 18000 swimmers; every run has been averaged over each of the four half-widths ($H$) of the double Poiseuille flow profile. For most runs, we fix $\epsilon=0.05$, and find this to be sufficient to capture the $\epsilon\ll 1$ regime that pertains to the multiple scales analysis.

3 Results and Discussion

Herein, the results for the steady state cross-stream swimmer concentration profiles are presented primarily for ABP’s. While RTP’s exhibit analogous behavior for small to O(1) $Pe$’s, the location of the boundary separating the high-shear and low-shear trapping regimes differs. The resulting implications, and the associated sensitivity of this boundary to the precise swimming mechanism, will be reported in a later communication. We first discuss swimmer migration towards walls leading to high-shear trapping, highlighting the ranges of $\kappa$ and $Pe$ that correspond to this behavior (§3.1), before moving on to examining migration towards the centerline leading to low-shear trapping (§3.2).

3.1 Swimmer migration towards walls - High-shear trapping

In figure 2, we plot the normalized steady state concentration profile [$\bar{F}(z)=F(z,t_{2}\to\infty)/F(z,0)$], using (8) with $D_{zz}$ and $V_{z}$ determined numerically as mentioned in §2.1, for infinitely slender swimmers ($\kappa=\infty$), and for $\kappa=15$, over a range of $Pe$. Also shown are the profiles obtained from Langevin simulations, for $\kappa=\infty$, for $Pe=10$ and 40, which are in good agreement with the predictions of the multiple scales analysis. All profiles exhibit near-center depletion caused by wallward migration of swimmers. For small $Pe$, a power series expansion of the swimmer probability density in $Pe$, of the form $F(z)=F^{(0)}+Pe^{2}F^{(2)}(z)+\ldots$, yields a parabolic depletion profile at O($Pe^{2}$), the quadratic scaling being consistent with invariance to flow reversal (see appendix B). As shown in the inset of figure 2b, for $\kappa=15$ and $Pe\leq 0.5$, the numerical profile for $F^{(2)}(z)=(F-F^{(0)})/Pe^{2}$ agrees well with the small-$Pe$ analytical solution. With increasing $Pe$, there is a qualitative change in the shape of the concentration profiles (the profile curvature changing sign as one moves from the centerline towards the walls), and for $\kappa=\infty$, the profiles appear to asymptote to a singular cusped profile in the limit $Pe\rightarrow\infty$ (see figure 2a). The mechanism underlying near-center depletion is readily explained. For small $\epsilon$, the swimmer orientation distribution rapidly attains a $Pe$-dependent equilibrium at each $z$, with the peak of this distribution aligning with the flow direction with increasing $Pe$. Since the local Péclet number, $Pe_{z}$, increases away from the centerline (where $Pe_{z}=0$) as shown in figures 3a and b, swimmers are, on average, more flow-aligned near the walls. The resulting reduction in the gradient component of the swimming velocity leads to a net wallward migration, eventually leading to an inhomogeneous steady state where this migration is balanced by an opposing diffusive flux.

To further investigate the nature of the wallward migration, we plot, in figure 4, $D_{zz}$ and $V_{z}$ at the channel wall as a function of Pe for $\kappa=\infty$ and 15. For $\kappa=\infty$, both $D_{zz}$ and $V_{z}$, for any non-zero $z$ scale as $Pe^{-4/3}$ for $Pe\gg 1$; being determined by $Pe_{z}$, the onset of this asymptotic scaling regime is postponed to a progressively larger $Pe$ with decreasing $z$ (not shown). For $\kappa=15$ too, there is an intermediate range of $Pe$’s ($\sim$ 3 to a little greater than 10) where the aforementioned scaling holds, but there is a deviation at larger $Pe$’s. As will be seen in §3.2, this deviation is a signature of the impending low-shear trapping regime. The observed scalings may be rationalized by noting that $D_{zz}\sim U^{2}_{z}t_{c}$, where $U_{sz}\sim$ O$(U_{s}Pe^{-1/3})$ is the swimming speed, projected along the gradient direction, of a swimmer in the large-$Pe$ orientation boundary layer, with an angular extent of O($Pe^{-1/3}$), that forms around the flow direction [Hinch & Leal (1972a); Brenner (1974)]; $t_{c}\sim\langle\theta^{2}\rangle/D_{r}\sim Pe^{-2/3}/D_{r}$ is the time over which rotary diffusion causes the swimmer orientation to fluctuate across the flow axis, leading to a de-correlation in its gradient-directed motion. The scaling for $V_{z}$ may be obtained by from arguments based on a gradient-directed biased random walk. One may write $V_{z}=\Delta z_{+}f_{+}-\Delta z_{-}f_{-}$, where $\Delta z_{+}$, $\Delta z_{-}$ represent the random walk jumps with frequencies $f_{+}$, $f_{-}$ along the positive and negative gradient directions, respectively. Expressing the random walk jumps as the ratio of swimming velocity in that direction to the associated frequency, one has $V_{z}=U_{sz}(z)-U_{sz}(z+\Delta z^{\prime})=U_{s}(Pe^{-1/3}_{z}-Pe^{-1/3}_{z+\Delta z^{\prime}})$, where $\Delta z^{\prime}$ is the (dimensionless) gradient-projection of the swimmer mean free path. With a swimming velocity of O($U_{s}Pe^{-1/3}_{z}$) and a decorrelation time of O($Pe^{-2/3}_{z}D^{-1}_{r}$), one obtains $\Delta z^{\prime}=U_{s}H^{-1}Pe^{-1}_{z}D^{-1}_{r}$. Using this leads to $V_{z}\sim$ O$(U^{2}_{s}Pe^{-4/3}_{z}H^{-1}D^{-1}_{r}/z)$.

The identical drift and diffusivity scalings above imply that the ratio $V_{z}/D_{zz}$ is independent of $Pe$ for large $Pe$, and hence, so is the concentration profile. We find $V_{z}/D_{zz}\propto 1/z$, consistent with the scaling arguments above, with a constant of proportionality ($C$) of order unity. This leads to a large-$Pe$ functional form of the concentration profile given by $\bar{F}=(1+C)|z|^{C}$. In figure 5, for $\kappa=\infty$, we compare the variation of $|V_{z}|/D_{zz}$ against $z$, for different $Pe$. A fit to the large-$Pe$ variation away from the centerline leads to $C=0.344$. From figure 2a, the swimmer concentration profiles are seen to approach the above $Pe$-independent limiting form for $Pe\geq 40$. Although, the steady state profile is $Pe$-independent, the time taken to attain steady state is of O($(H^{2}D_{r}/U^{2}_{s})Pe^{4/3}$) and diverges for $Pe\rightarrow\infty$. The limiting concentration profile is singular with a cusp at the centerline ($z=0$). The cusp is, of course, an artefact of the $Pe=\infty$ limit. For any finite $Pe$, the multiple scales analysis breaks down in an asymptotically small interval of O($U_{s}/D_{r}$) in the vicinity of the centerline (keeping in mind our interpretation of $\epsilon$ as a swimmer Knudsen number, this layer in the vicinity of the centerline, where ballistic swimmer trajectories play an important role, may be termed a Knudsen layer in traditional kinetic theory parlance). Note that both the swimmer orientation distribution and concentration in this Knudsen layer may be influenced by the bound trajectories observed for deterministic swimmers observed in Zöttl & Stark (2013).

In figure 6, we compare our theoretical and simulation results (which are in mutual agreement regardless of $Pe$) with those reported in Rusconi et al. (2014b) for $\kappa=10$. For the three smallest $Pe$’s (1.25, 2.5 and 5), our theory under predicts the inhomogeneity compared to the experiments. There is improved agreement at $Pe=10$, and again, at the highest $Pe$’s (25 and 50), the disagreement with the experiments widens. The disagreement at the two highest $Pe$’s, we believe, is because the experiments [Rusconi et al. (2014b)] sampled under-developed concentration profiles. The channel length available for development of the concentration profile corresponds to a single arm of the serpentine tube used in the experiments, and is approximately 5 cm (note that, in light of the later experiments by Aronson and co-workers [Sokolov & Aranson (2016); Sokolov et al. (2018), in a curvilinear geometry, it is likely that the bends of the serpentine tube significantly disrupt the concentration profile developed in the straight sections). At large $Pe$, the entrance length may be estimated based on the time scale for attainment of steady state mentioned above, and is O($0.06\bar{u}H^{2}D_{r}Pe^{4/3}/U^{2}_{s}$). Here, $\bar{u}=(2\bar{\dot{\gamma}}H/3)$ is the average suspension velocity and $\bar{\dot{\gamma}}=PeD_{r}$ is the the mean shear rate; the numerical pre-factor (0.06) in the above estimate is determined from the Langevin simulations. Based on this estimate, and using $H=212.5\mu$m, $U_{s}=50\mu$m $s^{-1}$, $D_{r}=1s^{-1}$ [Rusconi et al. (2014b)], the entrance length at $Pe=50$ is 1.4m! For $Pe=10$, the entrance length is about 3 cm, which is smaller than the aforementioned length of a straight section. This would suggest that the measured profile at $Pe=10$ is indeed fully developed. The fact that our simulation results at intermediate times closely match shapes of the experimental profiles at two largest $Pe$’s (not shown here), reinforces the above assertion of the latter profiles not being fully developed ones. A second reason for deviation from theory could be the role of wall interactions, which are not included in the present analysis, the focus here being on bulk mechanisms for migration. Theory and simulations based on image induced interactions of pusher-type swimmers, and experiments on smooth swimming bacterial strains have shown swimmer accumulation at the boundaries in the absence of an imposed flow [Berke et al. (2008); Lauga & Powers (2009); Chilukuri et al. (2014)]. However, in the present context we expect this effect to only have a minor influence for two reasons. First, the time scale $t_{w}\sim$ O($L/U_{s}$), over which the wall interaction is significant, is comparable to the time scale for orientation decorrelation of O$(D^{-1}_{r})$ for wild-type swimmers under consideration, and the latter acts to disrupt any hydrodynamically induced accumulation for small $Pe$; in contrast, for the smooth swimmers used in Berke et al. (2008), the decorrelation time scale given by O($D_{r-Brownian}^{-1}$), was greater by at least two orders of magnitude. At large $Pe$, rotary diffusion induced decorrelation is weak, and one then has a sustained image-induced interaction between nearly flow-aligned swimmers. Wall interactions, even in this case, are expected to only induce a concentration enhancement within a region of O($L$) from the walls. Using Brownian dynamics simulations for a confined channel geometry, with hydrodynamic interactions between the swimmers and wall included, Chilukuri et al. (2014) showed that wall accumulation is indeed localized to the immediate vicinity of boundaries, consistent with the above arguments, and that this accumulation is further suppressed when the swimmers are subject to a pressure driven channel flow. Given that the observations of Rusconi et al. (2014b) and Barry et al. (2015) were in the central portion of the channel, one expects the swimmer concentration to remain largely unaffected by wall interactions.

The disparity at the smallest $Pe$’s, pointed out above, could be due to multiple reasons. First is the possibility of the assumption $\epsilon\ll 1$, underlying the multiple scales analysis, breaking down. Based on the microchannel geometry and the known mean free path for B. subtilis (estimated from the decorrelation time that is in turn inferred from the measured translational diffusivities), $\epsilon\approx 0.23$ in the experiments of Rusconi et al. (2014b). Our simulations, which are not limited by the small $\epsilon$ assumption, agree well with theory even for $\epsilon=0.2$, ruling out the lack of smallness of the swimmer Knudsen number as a cause of theory-experiment incongruity. Note that the issue of $\epsilon$ not being small enough is not relevant to the largest $Pe$’s since the stronger flow alignment of the swimmers with the unaltered decorrelation time implies that the projection of the mean free path in the cross-stream direction shrinks with increasing $Pe$, increasing the range of validity of the multiple scales analysis. For the high-shear trapping regime considered in this section, the effective Knudsen number is $\epsilon_{eff}=\epsilon Pe^{-1/3}$, the exponent reflecting the degree of swimmer alignment in the orientation boundary layer; for low-shear trapping examined next, $\epsilon_{eff}=\epsilon Pe^{-1}$ for large $Pe$. Next, it is also possible that the particular choice of the swimming mechanism determines the detailed concentration profile; for instance, there might be a significant difference between the degree of inhomogeneity in the swimmer concentration between RTP’s and ABP’s. As explained in the introduction, the former involves the swimmer undergoing a sudden large change in its orientation during tumbles [Berg (1993, 2008)], while the latter involves continuous small amplitude fluctuations in their orientation. Our Langevin simulations do suggest a greater inhomogeneity for RTP’s than ABP’s at the smallest $Pe$’s, with this difference diminishing in magnitude with increasing $Pe$. Thus, the difference between the swimming mechanism of B. subtilis, relative to the pure ABP-dynamics assumed in the theoretical framework here, may cause the discrepancy between experiment and theory at the smallest $Pe$’s (the Rusconi et al. (2014b) experiments only measured the translation diffusivities, and do not therefore have detailed information with regard to the relative importance of tumbling vis-a-vis rotary diffusion). The sensitivity of the shear induced migration pattern to the swimming mechanism does deserve a more detailed study, and this will be reported in a future effort.

3.2 Swimmer migration towards the channel centerline - Low-shear trapping

Figure 7 shows the $Pe$-dependent swimmer concentration profiles for a range of swimmer aspect ratios, including the case $\kappa=15$ examined earlier (now over a larger range in $Pe$ extending upto 2500). For small to intermediate $Pe$, one observes high-shear trapping regardless of $\kappa$, consistent with the results of the small-$Pe$ analysis, and as seen before, the corresponding concentration profiles exhibit a single minimum at the centerline. At large $Pe$, however, there is a qualitative change with the concentration profiles now exhibiting a pair of maxima symmetrically disposed about the centerline; the maxima increase in amplitude, while converging towards the centerline with increasing $Pe$ (consistent with the earlier finite element computations of Bearon & Hazel (2015)). The aforementioned transition occurs across a $\kappa$-dependent threshold $Pe$. The Langevin simulations reinforce the existence of the transition in the concentration profiles; as shown in figure 7a, the profile at $Pe=$4 exhibits a depletion of swimmers near the centerline, and the one at $Pe=$ 40 exhibits a near-center excess.

To understand the transition from high-shear to low-shear trapping with increasing $Pe$, in figure 8, we plot $D_{zz}|_{z=-1}$ and $|V_{z}|_{z=-1}$, appropriately scaled, against $Pe$ for $\kappa=\infty$ and for the $\kappa$’s in figure 7. From figures 8c and d, $V_{z}$ is seen to undergo a change in sign with increasing $Pe$ (the associated zero-crossing appears as a sharp dip in the log-log plots), which is responsible for the reversal in migration. For $Pe$’s smaller than that corresponding to the zero-crossing (of $V_{z}$), swimmers migrate towards the walls; for larger $Pe$’s, swimmers migrate towards the centerline, leading to a low-shear trapping behavior. It is important to note that the reversal in migration is actually a function of $Pe_{z}$, occurring for $Pe_{z}$ exceeding a $\kappa$-dependent threshold. This is confirmed in figure 9, where the plots of the diffusivity and drift against $Pe$ for $\kappa=10$, at different transverse locations ($z$), confirm the onset of the asymptotic O($Pe^{-2}$) scaling regime at progressively larger $Pe$ for $z$ approaching the centerline (figures 9a and b). The reversal in drift is also delayed at locations closer to the channel centerline, with the earliest reversal occurring at the walls (the one shown in figure 9b). To further confirm the role of $Pe_{z}$ (rather than $Pe$), in figures 9c and d, $D_{zz}$ and $V_{z}$ are plotted as a function of $Pe_{z}$. This collapses the different $D_{zz}$ curves into a single one, and also collapses the zero-crossings of the different $V_{z}$ curves onto a single location corresponding to a critical $Pe_{z}$. Since $Pe_{z}$ always becomes arbitrarily small close enough to the centerline, there is always a region close to the centerline, regardless of $Pe$, where swimmers continue to exhibit a wallward migration. This leads to the non-monotonicity, with a pair of maxima bracketing a central dip in the concentration profiles shown in figure 7; thus, the central portions of all the profiles in figure 7, corresponding to low-shear trapping, still resemble the depletion profiles seen earlier is §3.1. It is only the limiting excess profile, for $Pe=\infty$, where swimmers migrate towards the centerline regardless of $z$ (this singular profile is subject to the limitations of the multiple scales analysis mentioned earlier in §3.1). As evident from figure 8, for any finite $\kappa$, both $D_{zz}$ and $V_{z}$ exhibit an O($Pe^{-2}$) scaling in the low-shear trapping regime that emerges for sufficiently large $Pe$, although for the larger values of $\kappa$ (10 and 15), there is the emergence of an intermediate scaling regime (evident from figures 8a and c) where $D_{zz},V_{z}\sim Pe^{-4/3}$, already seen in §3.1. Thus, $\kappa=\infty$ is the only case where the O($Pe^{-4/3}$) scaling regime persists in the limit $Pe\rightarrow\infty$, as indicated by the black-dashed lines in figure 8. Infinitely slender swimmers therefore represent a singular limit with regard to shear induced migration behavior. Of course, for large enough $\kappa$’s, the $Pe$ corresponding to the transition to low-shear trapping may be large enough as to be numerically inaccessible. This is seen to be the case in figure 8b for $\kappa$=15, and in figure 8d for both $\kappa=10$ and $\kappa=15$; owing to the zero-crossing, it takes a larger $Pe$ for the drift to conform to the eventual O($Pe^{-2}$) scaling behavior.

From the results presented thus far, it is seen that the transition from high to low-shear trapping occurs at a threshold $Pe$ of order unity for $\kappa$’s of order unity. For large $\kappa$, the cross-over $Pe$ is asymptotically large, and marks a transition from O($Pe^{-4/3}$) to an O($Pe^{-2}$) scaling regime for the drift and diffusivity coefficients (see figure 9). From classical work on the orientation dynamics of passive anisotropic Brownian particles in a simple shear flow [ Leal & Hinch (1971); Hinch & Leal (1972a); Brenner (1974)], the former scaling regime implies the localization of the swimmer orientation distribution in an O($Pe^{-1/3}$) boundary layer around the flow direction where there is a balance between the rotary diffusion and flow-induced rotation. In contrast, when flow-induced rotation along Jeffery orbits occurs on a time scale much shorter than rotary diffusion, the latter acts as a regular perturbation; there is no orientational boundary layer, and one expects a scaling regime involving an integral powers of $Pe^{-1}$. Thus, the aforementioned cross-over $Pe$ must correspond to the point where the flow-induced (Jeffery) rotation, on a time scale of O($\dot{\gamma}^{-1}\kappa$) is comparable to that characterizing rotary diffusion of O($D^{-1}_{r}Pe^{-2/3}$); this leads to $Pe_{c}\sim$O($\kappa^{3}$) for $Pe\gg 1$ (this scaling is confirmed in §3.3.2, where we demarcate high and low-shear trapping on the $Pe-\kappa$ plane). Now, at leading order, one has purely flow induced Jeffery rotations of the swimmer, with an orbital time period $t_{jeff}\sim 1/(D_{r}Pe)$. The Jeffery rotations coupled with swimming lead to a swimmer trajectory, that along the cross-stream direction, has a bounded but oscillatory character with an amplitude $\Delta z$ of O($U_{s}/(D_{r}Pe)$). Rotary diffusion disrupts the exact periodicity of the leading order trajectory, causing an occasional slip of O($U_{s}/(D_{r}Pe)$). These random displacements in the $z$-direction, occurring over a time scale of O($D^{-1}_{r}$), leads to a $D_{zz}\sim\langle\Delta z^{2}\rangle/t_{c}$ which may be written as O($(U^{2}_{s}/D_{r})Pe^{-2}$), confirming the scaling behavior observed in the low-shear trapping regime; here, $\langle\Delta z^{2}\rangle$ is the mean square displacement of the swimmer in the gradient direction. Note that, for large $\kappa$, swimmers spend an O$(\kappa/(D_{r}Pe))$ time remaining nearly aligned with the flow axis at an angle of O($\kappa^{-1}$). This leads to the same estimate for $\Delta z$ as above, but a faster decorrelation rate of O($\kappa^{2}D_{r}$) (since rotary diffusion only has to induce an angular displacement of O($\kappa^{-1}$)), leading to a scaling for $D_{zz}$ of O($(U^{2}_{s}\kappa^{2}/D_{r})Pe^{-2}$). Based on the biased random walk argument as mentioned in §3.1, we can express $V_{z}=U_{sz(z)}-U_{sz(z+\Delta z^{\prime})}$. Now, for swimmers with $\kappa$ of O(1), using $U_{sz(z)}=U_{s}Pe^{-1}_{z}$, we obtain $V_{z}\sim$ O($(U^{2}_{s}H^{-1}D^{-1}_{r}/z)Pe^{-2}_{z}$). The drift and diffusivity scalings above are broadly consistent with those obtained in figure 8b.

For large $Pe$, similar to the high-shear trapping case, $V_{z}$ and $D_{zz}$ exhibit the same O($Pe^{-2}$) scaling, their ratio being independent of $Pe$ and of the form $C(\kappa)/z$, where $C$ is now a function of $\kappa$. One therefore expects a limiting $Pe$-independent profile of the form $(C+1)|z|^{C}$; note that, for any finite $Pe$, this profile breaks down in the vicinity of the centerline, where the actual finite-$Pe$ profiles exhibit a dip, corresponding to the onset of depletion (wallward migration) at small $Pe_{z}$. Since $C(\kappa)<0$ for low-shear trapping (owing to the sign-reversal of the drift), this limiting profile has a diverging concentration at the centerline. For $\kappa=1.25$, $C(\kappa)\sim-0.238$, and the associated limiting profile compares well to the numerically computed ones in figure 7a. An interesting consequence of the change in sign of $C(\kappa)$, leading to low-shear trapping, is that the aforementioned limiting profile fails to be integrable for $C(\kappa)\leq-1$ (this wasn’t an issue in the high-shear trapping case, since the swimmer concentration in all profiles, including the one corresponding to $Pe=\infty$, is bounded below by zero). We find $C(\kappa)$ to equal this threshold value (-1) for $\kappa\approx 2$, and decrease further for larger $\kappa$. The implication of the non-integrability above is that, for $\kappa\geq 2$, the swimmer concentration profiles must exhibit a ‘centerline collapse’ in the limit $Pe\rightarrow\infty$; in other words, $\displaystyle{\lim_{Pe\rightarrow\infty,\kappa\geq 2}}\bar{F}(z)=\delta(z)$. Figure 10 shows the variation of swimmer concentration at the wall $(\bar{F}(-1))$ against $Pe$ for the aspect ratios considered in figure 7. For $\kappa=$ 1.25, the near-wall swimmer concentration approaches a finite limiting value $\left(\displaystyle{\lim_{Pe\rightarrow\infty,\kappa=1.25}\bar{F}(-1)=0.765}\right)$. For $\kappa\geq 2$, however, the near-wall swimmer concentration decreases monotonically to zero for $Pe\rightarrow\infty$, in accordance with the centerline-collapse hypothesized above. The slope characterizing this decrease appears to increase with increasing $\kappa$, implying that the collapse increases in intensity for larger $\kappa$ (for $\kappa=$ 15, the numerics have not accessed the asymptotic regime yet). The $\kappa$-dependence of the collapse does suggest that the drift and diffusivity coefficients must scale differently with $\kappa$ for $\kappa\gg 1$. Scaling arguments presented above suggested that $D_{zz}\sim$O$(\kappa^{2}Pe^{-2})$ in the low-shear trapping regime, and the large$-\kappa$ scaling of the drift must therefore differ from O($\kappa^{2}Pe^{2}_{z}/z$) that would emerge from the biased-random-walk scaling above modified for $\kappa\gg 1$. The large-$\kappa$ scaling of the drift in this regime will be analyzed in detail in a future effort.

Unlike the rather intuitive picture underlying high-shear trapping (see figure 3 in §3.1), the physical mechanism underlying low-shear trapping is more subtle. Here, we present arguments that support migration of swimmers towards the centerline at large $Pe$, beginning from the orientation dynamics of swimmers in a local linear flow description, as shown in figure 11a. For the parabolic flow under consideration, the local linear approximation yields a simple shear flow with a shear rate that decreases from a maximum at the wall to zero at the centerline. For any finite $\kappa$, and sufficiently large $Pe$ ($Pe\gg 1$ for $\kappa$ of O(1); $Pe\gg\kappa^{3}$ for $\kappa\gg 1$), swimmers at any non-zero $z$ execute Jeffery rotations at leading order, with a time period commensurate with the $z$-dependent shear rate of the simple shear flow above. The resulting distribution of azimuthal angles within any Jeffery orbit is symmetrical about the flow direction (as first shown by Leal & Hinch (1971), the distribution across the orbits is still controlled by asymptotically weak rotary diffusion, although this doesn’t affect the argument for the direction of drift given below). At the next order, O($1/Pe$), rotary diffusion breaks the aforementioned symmetry, leading to an excess of orientations in the extensional quadrants, and an equal and opposite deficit in the compressional ones (this antisymmetry is a consequence of Brownian motion being a weak but regular perturbation; in the singular case, relevant to the high-shear trapping regime, almost the entire swimmer probability is confined to the extensional quadrants). The aforementioned distributions at O(1) and O($1/Pe$) are sketched in figure 11a. Owing to the decrease of $Pe_{z}$, the amplitudes of the excess and deficit above increase as one moves away from the wall. The lower figure in 11b contains sketches of the local simple shear flows at three different $z$-locations, with the local excess and deficit in orientations in accordance with the arguments above (the swimmer lengths shown denote the amplitude of the excess or deficit). This picture may now be used to argue, in a quadrant-wise manner, for the swimming-induced changes the orientation distribution. At any $z$, in quadrant I, lesser excess swims in from the stronger shear flow below, and a greater excess exits into the weaker shear flow above; in quadrant III, the exact opposite occurs. Similarly, in quadrant II (IV), a greater (lesser) deficit exits into the weaker (stronger) shear flow above (below), with a lesser (greater) deficit swimming in from the stronger (weaker) shear flow below (above), leading to a net excess (deficit) of swimmers. As a result, there is a swimming-induced source of probability in quadrants II and III, and a sink in quadrants I and IV. In turn, this implies that the azimuthal flux of probability will, on a clockwise traversal, increase through quadrants II and III, and then decrease through IV and I. Accounting for the depleted orientation probability in the compressional quadrants (II and IV), this flux-profile leads to an excess of swimmers in quadrants I and II, implying a net drift towards low-shear regions. This migration towards the low-shear regions continues till the swimmers reach sufficiently close to the channel center where $Pe_{z}$ eventually falls below $\kappa^{3}$. From this location onward, and up until the centerline, the swimmers revert to a wallward migration. Therefore, the steady state concentration profile is such that it exhibits a maximum at a location intermediate between the centerline and the wall, with this location shifting towards the centerline with increasing $Pe$.

3.3 Depletion index and the $Pe-\kappa$ migration portrait

Having discussed the two kinds of migration patterns mentioned above, in §3.3.1, we quantify the same in terms of the depletion index $DI$, a non-dimensional measure of the inhomogeneity of the swimmer concentration field: $DI<0~{}(>0)$ corresponds to high-shear (low-shear) trapping. Next, in §3.3.2, we organize the migration behavior on a $Pe-\kappa$ plane (for small $\epsilon$), thereby delineating the regimes of high-shear and low-shear trapping across a numerically determined phase boundary.

3.3.1 The depletion index (DI)

We define $DI=2(A_{1}-A_{2})/A$, where $A_{1}$ represents the area under the initial uniform distribution from the wall up to its intersection with the steady state swimmer concentration profile ($z_{in}$), $A_{2}$ denotes the analogously defined area under the steady state profile from the channel wall to $z_{in}$ and $A$ is the total area under either profile. For depletion profiles, conservation of swimmer number does imply that equality of areas between the wall and $z_{in}$, and between $z_{in}$ and centerline. To analyze swimmer migration at large $Pe$, $\Delta DI=100-DI$ is plotted against $1/Pe$ for swimmers of different aspect ratios in figure 12. The plot discriminates between high ($\Delta DI>100$) and low-shear trapping behavior ($\Delta DI<100$), and in addition, helps differentiate between centerline-collapse ($\Delta DI\rightarrow 0$ for $1/Pe\rightarrow 0$ for $\kappa\geq 2$) and non-collapse ($\Delta DI$ approaches a finite value as for $1/Pe\rightarrow 0$ for $\kappa<2$) behavior within the low-shear trapping regime. Figure 12 clearly illustrates the singular nature of infinitely slender swimmers. For $\kappa=\infty$, $\Delta DI$ plateaus at a finite value ($>100$) for $1/Pe\rightarrow 0$. In contrast, for any finite $\kappa\geq 2,\Delta DI$ approaches zero in the same limit, and the rate of approach, rather counter intuitively, increases with increasing $\kappa$. For large $\kappa$’s, this accelerated approach towards the origin begins after an intermediate extended plateau with $\Delta DI>100$ (the same plateau value as for $\kappa=\infty$). Note that Rusconi et al. (2014b), likely motivated by experimental considerations, had defined depletion index ($DI_{Exp}$) based on the area between the actual and uniform profile in the central half width of the channel. The general character of figure 12 is, however, unchanged regardless of the definition.

It is worth reiterating the behavior of the depletion index, $\Delta DI$, as a function of the swimmer aspect ratio starting from infinitely slender swimmers. Regardless of $\kappa$, $\Delta DI$ starts from 100 for $Pe\ll 1$. For $\kappa=\infty$, $\Delta DI$ increases with increasing $Pe$, approaching a constant ($\Delta DI=111.1$) for $Pe\rightarrow\infty$. For large but finite $\kappa$, $\Delta DI$ again increases with increasing $Pe$, asymptoting to the aforementioned plateau value in the range $1\ll Pe\ll\kappa^{3}$, and thereafter, decreases rapidly with further increase in $Pe$. This behavior continues with decreasing $\kappa$, with the $Pe$-interval corresponding to the intermediate plateau above progressively shrinking in extent, and the eventual rate of approach to zero at large $Pe$ also decreasing with decreasing $\kappa$. For the ABP’s examined here, this trend continues until a critical aspect ratio of approximately 2, when $\Delta DI$ approaches zero only at a logarithmic rate for large $Pe$; that is, for $\kappa=2$, $\Delta DI$ is O($1/\ln{Pe}$) for $Pe\gg 1$. For $\kappa<2$, $\Delta DI$ approaches a finite positive plateau (less than 100) for $Pe\rightarrow\infty$, with this plateau value approaching 100 as $\kappa$ approaches unity (spherical swimmers). It is worth noting that the infinite-aspect-ratio assumption is used routinely in the rheological context, when analyzing swimmers suspensions, and generally yields good agreement in between experiments and theory [Saintillan (2010a, b); López et al. (2015); Nambiar et al. (2017, 2019b)]. In the present context, however, even for moderately large aspect ratios, $\kappa\sim 4-6$ (relevant in the biological context), at large enough $Pe$, the swimmers can, in fact, end up being trapped in the low-shear regions. Thus, the effective aspect ratio parameter appears to be of considerable importance in the context of shear-induced migration.

Although our theoretical predictions and simulations are in mutual agreement, there is nevertheless a disagreement with the experiments of Rusconi et al. (2014b); aspects of this disagreement were already elaborated on in §3.1. In order to illustrate the discrepancy, we have plotted $DI_{Exp}$, against $Pe$, as an inset in figure 12 for $\kappa=10$. While the disagreement appears a matter of detail, in the sense that the $DI_{Exp}$ in the experiments peaks at a smaller $Pe$ ($\approx$ 10), and starts decreasing earlier, it is not! The interpretation in the original experiments was that the $DI_{Exp}$ would asymptote to zero for larger $Pe$. As discussed earlier, this decay of the $DI_{Exp}$ to zero, implying the approach of the swimmer concentration profile towards homogeneity, is a consequence of the experimental profiles not being fully developed on account of the limited residence time (ironically enough, the apparent peaking of the $DI_{Exp}$ at a finite shear rate in these experiments was motivation for subsequent experiments by Aronson and co-workers [Sokolov & Aranson (2016); Sokolov et al. (2018)], in a curvilinear geometry, where the inhomogeneity in swimmer concentration was argued to persist at infinite $Pe$, albeit for entirely different reasons). As argued above, the $DI_{Exp}$ does not approach zero, for any aspect ratio, in the limit of large $Pe$. Thus, the apparent decay of the theoretical $DI_{Exp}$ is a signature of the transition to a low-shear trapping regime, and the theoretical $DI_{Exp}$ will eventually asymptote to -1, for the aspect ratio of 10 considered (which allows for a centerline collapse as $Pe\rightarrow\infty$).

3.3.2 The $Pe-\kappa$ migration portrait

Figure 13 organizes the high-shear and low-shear trapping regimes on the $Pe-\kappa$ plane, with an approximate boundary separating the two. We differentiate the high-shear trapping from the low-shear trapping, based on the slope of the swimmer concentration profile at the channel wall. Across the phase boundary, there is a change in sign of this slope, arising from the reversal in the drift; see figure 8 in §3.2 (note that this criterion suffices in the present context where the focus is on shear induced migration in the bulk; inclusion of wall interactions, would require a modification of this criterion. This is not an issue for comparison with the experiments, however, since the swimmer concentrations were monitored well away from the boundaries). The large aspect ratio asymptote for the aforementioned boundary between the two trapping regimes has already been determined based on the scaling arguments given in §3.2, and is of the form $Pe_{c}\propto\kappa^{3}$; the pre-factor 0.07 given in figure 13 is determined from a best fit for the curve separating the numerically determined depletion and excess profiles. Now, due to the absence of any preferential alignment under shear, spherical-swimmer concentration profiles are homogeneous regardless of $Pe$, in the dilute non-interacting regime under consideration. Nevertheless, the small but finite inhomogeneity that arises in the limit $\kappa\rightarrow 1$ can transform from a near-center depletion to a near-center excess at a finite $Pe$. This is indeed the case, and based on an expansion of the swimmer probability density in the small Bretherton constant $B$, as $\Omega=F/4\pi+B\Omega_{1}+\cdots$, one finds $Pe_{c}\sim 4.5$ (see appendix C). Therefore, as $\kappa\to 1$, ABP’s exhibit high and low-shear trapping profiles across $Pe_{c}=4.5$, with the analytical forms for the latter profiles also exhibiting the characteristic pair of maxima on either side of the centerline, already seen in the numerical profiles. This transition for near-spherical swimmers has again been validated using numerical predictions based on the multiple scales analysis, and Langevin simulations.

Also shown in figure 13 are the results from the experiments of Barry et al. (2015) corresponding to both high-shear and low-shear trapping of the phytoplankton algae Heterosigma and Dunaliella. The effective swimmer aspect ratios ($\kappa$) in the experiments were determined from a best fit of the theoretical orientation distribution, obtained by solving the Fokker-Planck equation with $\kappa$ as a fitting parameter, with the experimentally determined one. Note that the $\kappa$’s determined in the manner above were functions of $Pe$, implying that the swimmer shape changed on account of the shear induced deformation. Based on the concentration profiles and depletion index plots ($DI_{Exp}$) reported in Barry et al. (2015), for $Pe\leq 5$, Heterosigma exhibit high-shear trapping, whereas, Dunaliella, show high-shear trapping till $Pe\sim 5$ but transition to low-shear trapping at at $Pe\sim 12.5$. Similar to the observation of Barry et al. (2015), for the algae Dunaliella, we too observed the transition from the high-shear to low-shear trapping with increase in $Pe$. The experimental data points are thus consistent with the trapping regimes identified in figure 13, and this consistency serves as an important validation of the analysis and arguments presented here.

4 Conclusions

Motivated by recent microfluidic experiments [Rusconi et al. (2014b); Barry et al. (2015)], in this work we have used both theory and simulations to systematically study shear induced migration in a dilute suspension of microswimmers subject to a plane Poiseuille flow in a wide channel. On the theoretical front, we have used the method of multiple time scales to derive a drift-diffusion equation governing the steady state cross-stream swimmer concentration profile. Based on the swimmer geometry, as characterized by its equivalent aspect ratio $\kappa$, and a rotary Péclet number $Pe=U_{f}/(D_{r}H)$ that characterizes the relative efficiency of stochastic reorientation and flow-induced rotation, we delineate high-shear and low-shear trapping regimes. Our theoretical predictions are reinforced by the results of Langevin simulations, and establish for the first time, the existence of a low-shear trapping regime at sufficient large $Pe$ for any finite $\kappa$. The transition from high-shear to low-shear trapping appears consistent with the experimental data available from Barry et al. (2015); the biflagellated algal species, Chlamydomonas and Dunaliella, used in the experiments exhibit a low-shear trapping regime at high shear rates, with the latter species going through a high-shear trapping regime at lower shear rates. The approach of microswimmers to any surface is the first step towards colonization and biofilm formation. The existence of both high- and low-shear trapping regimes points to the possibility of manipulating fluid shear to possibly reinforce or retard biofilm formation; thus, the high and low-shear trapping regimes in figure 13 may perhaps be equivalently interpreted as regimes where hydrodynamics favors and opposes biofilm formation. A richer range of concentration patterns and dynamical regimes, relative to those identified earlier in Rusconi et al. (2014b) and Bearon & Hazel (2015), is expected in the presence of an ambient chemical gradient that causes chemotactic swimmers to drift towards a boundary (as a possible precursor for biofilm formation), while a high-$Pe$ shear drives migration in the opposite direction. This particular aspect will be pursued in the near future.

The organization of shear induced migration of active swimmers as a $Pe-\kappa$ phase portrait has important implications. The fact that swimmers on either side of the phase boundary in figure 13 migrate in opposite directions allows, in principle, for a microfluidic shape sorter ($\kappa$ being a surrogate for the swimmer shape). A shape-based sorting might be feasible in the low-shear trapping regime alone by exploiting the centerline collapse (for $Pe\rightarrow\infty$) that occurs for $\kappa\geq 2$, (application of the multiple scales analysis shows that results for the plane Poiseuille flow derived herein, including the threshold aspect ratio for centerline collapse, directly carry over to the pipe Poiseuille case). A second interesting implication of the centerline collapse is for the dispersion of microswimmers in the flow direction (this has been examined, in the limit of confined channels, by Chilukuri et al. (2015)). It is known from classical Taylor-dispersion theory that, in the limit where the microswimmers uniformly sample the channel cross-section, the mean square displacement in the flow direction scales as O($Pe_{t}^{2}D_{zz}t$), and therefore, increases with increasing $Pe_{t}$; here, $Pe_{t}=U_{f}H/D_{zz}=Pe^{3}/\epsilon^{2}$, is the Péclet number based on the swimmer translational diffusivity. For large-aspect-ratio microswimmers in the low-shear trapping regime identified here, the extent of the channel cross-section sampled decreases sharply with increasing $Pe_{t}$ (or increasing $Pe$ with $\epsilon$ fixed), pointing to a possible non-trivial scaling of the flow-induced dispersion in the limit of large $Pe_{t}$.

An interesting consequence of the detailed analysis here is the remarkable sensitivity of the shear induced migration to subtle changes in the orientation dynamics. These differences in orientation dynamics were originally identified in the context of passive anisotropic particles, in a series of studies by Hinch and Leal [Leal & Hinch (1971); Hinch & Leal (1972a); Leal & Hinch (1972); Hinch & Leal (1972b); Leal & Hinch (1973); Hinch & Leal (1973)], who were motivated by the need to understand the subtle role of weak Brownian motion in strong shearing flows, and its implications for suspension rheology. In fact, they turn out to much more crucial for active particles, being responsible for the onset of a low-shear trapping regime! Indeed, Nitsche & Hinch (1997) have previously used the method of multiple scales to examine shear induced migration in a dilute suspension of passive anisotropic particles. Such particles do exhibit a high-shear trapping on account of the anisotropy of the mobility, and thence, the Brownian diffusivity. The trapping arises because of the lower transverse diffusivity(mobility) of flow-aligned particles close to the walls, but it is a relatively weak effect owing to the mobility anisotropy (defined as the ratio of the longitudinal to the transverse mobility coefficients) saturating at 2 even for infinitely slender particles; experiments on fiber suspensions have failed to observe the predicted high-shear trapping [Nitsche & Hinch (1997)] at lower volume fractions [Strednak et al. (2018)]. Since the Brownian mobility, and therefore, the diffusivity at any non-zero $z$, asymptotes to a finite value for $Pe\rightarrow\infty$, the concentration profile must become spatially homogeneous in this limit; in sharp contrast to the active case examined here.

Assuming an independent source of translational diffusion (thermal or otherwise), with components $D^{t}_{\parallel}$ and $D^{t}_{\perp}$ along and transverse to the flow direction, in addition to the translational diffusion arising from swimming activity considered here, one may use scaling arguments to predict the transition from the active shear induced migration analyzed here to the migration behavior of passive anisotropic particles analyzed by Nitsche & Hinch (1997). This transition is most relevant to large aspect ratio swimmers, since swimmers with moderate $\kappa$ would have transitioned to their asymptotic low-shear trapping regime (collapse or non-collapse) at fairly moderate values of $Pe$. To estimate the threshold for the active-passive transition, we balance $D^{t}_{\perp}$ with the high-$Pe$ asymptote for the translational diffusivity in the high-shear trapping regime, given by $\sim 0.22Pe^{-4/3}U^{2}_{s}/D_{r}$ (see figure 9); we consider the diffusivity at the wall, which gives an upper bound for the threshold. The balance yields $Pe_{c}\sim[0.22U_{s}^{2}/(D_{r}\,D^{t}_{\perp})]^{3/4}$. Equating this threshold to $\kappa^{3}$, one obtains a threshold aspect ratio $\kappa_{c}\sim[0.22U_{s}^{2}/(D_{r}\,D^{t}_{\perp})]^{1/4}$ (as must be the case, the same expressions for $Pe_{c}$ and $\kappa_{c}$ result on using $D_{zz}\sim Pe^{-2}\kappa^{2}U^{2}_{s}/D_{r}$ corresponding to the low-shear trapping regime). The significance of $\kappa_{c}$ above is that the concentration profiles for swimmers with $\kappa>\kappa_{c}$ would transition to homogeneity (for $Pe\rightarrow\infty$) directly from the high-shear trapping regime, while the concentration profiles for swimmers with $\kappa<\kappa_{c}$ will asymptote to homogeneity from the low-shear trapping regime. When $D^{t}_{\perp}$ has a thermal origin, using the familiar Stokes-Einstein expression with the translational mobility for a slender body, one has $D^{t}_{\perp}=k_{b}T\ln\kappa/(4\pi\eta L)$ [Dhont (1996)] with $k_{b}$ being the Boltzmann constant, $T$ being the temperature, and $\eta$ being the suspension viscosity. This yields $\kappa_{c}\sim$O(10) for $L\approx 9\mu$m, $U_{s}\approx 50\mu$m, $\kappa=10$, $T\approx 300$ and $\eta\approx 10^{-3}$ for the aqueous medium used to suspend the microorganisms [Rusconi et al. (2014b); López et al. (2015)]. Thus, the B. subtilis in Rusconi et al. (2014b) experiments are expected to transition to a passive-migration behavior at a $Pe$ that is about two orders of magnitude higher than the highest $Pe$ sampled in the experiments. The larger algae used in Barry et al. (2015) should yield much larger $\kappa_{c}$ values.

Finally, it must be noted that figure 13 depicts migration behavior specifically for ABP’s. One expects a similar phase portrait for other swimmer types (RTP’s, for instance), although the aforementioned sensitivity to the orientation dynamics obviously points to different $Pe-\kappa$ scalings for the phase boundary (and a different threshold aspect ratio for a possible centerline collapse). Such an analysis, and its implications towards flow-based separation of microswimmers, with different swimming characteristics, will be reported separately.

Acknowledgements

L.N.Rao would like to thank Science and Engineering Research Board, India (Grant No. PDF/2017/002050) for the financial support.

Appendix A Numerical scheme for the swimmer probability densities and concentration profile.

In the following, we describe the numerical scheme adopted to evaluate the swimmer orientation probability densities $G(\mathbf{p})$ and $\Omega_{1}(\mathbf{p})$, at O(1) and O($\epsilon$), respectively, and thence the steady state swimmer concentration profile of (7), mentioned in §2.1. To solve (4), for $G(\mathbf{p})$, we use a spherical coordinate system with its polar axis aligned with the $z$-axis (see figure 1). As a result, the components of the swimmer orientation vector are given by $p_{1}=\sin\theta\cos\phi$, $p_{2}=\sin\theta\sin\phi$, $p_{3}=\cos\theta$. The flow-induced rotation term in (4), $\nabla_{p}\cdot(\tilde{\dot{\mathbf{p}}}G)$, may now be written as:

[TABLE]

where the operators [Messiah (1962); Doi & Edwards (1978)]

[TABLE]

Expressing (LABEL:EQ:GExpansion) in terms of spherical harmonics, we have:

[TABLE]

Now, expanding

[TABLE]

and substituting in (11), one may write:

[TABLE]

where $Y_{l}^{m}(\mathbf{p})=\sqrt{(2l+1)/(4\pi)(l-m)!/(l+m)!}P_{n}^{m}(\cos\theta)\exp{(im\phi)}$, represent the spherical harmonics and $P_{n}^{m}$ are the associated Legendre functions [Abramowitz & Stegun (1965)] and $\mathcal{L}_{i}|Y_{l}^{m}>=\mathcal{L}_{i}(Y_{l}^{m})$.

Next, we substitute the expansion for $G(\mathbf{p})$ above in the remaining terms of (4), with the tumbling kernel being given by $\displaystyle K(\mathbf{p}/\mathbf{p}^{\prime})=\beta\exp{(\beta(\mathbf{p}\cdot\mathbf{p}^{\prime}))}/(4\pi\sinh\beta)=\sum_{n=0}^{\infty}A_{n}P_{n}(\mathbf{p}\cdot\mathbf{p}^{\prime})$ [Subramanian & Koch (2009); Nambiar et al. (2017)]; here, $\beta$ is the correlation parameter characterizing the tumbles, and $P_{n}$ represent the Legendre polynomial of the $n$ th degree, with argument being the cosine of the angles between the pre- and post-tumble operations. Further, employing the addition theorem of spherical harmonics and that $\nabla^{2}_{p}Y_{l}^{m}=-l(l+1)Y_{l}^{m}$ [Abramowitz & Stegun (1965)], (4) may finally be written as the following summation over the spherical harmonics:

[TABLE]

In the above equation, $A_{0}=1/4\pi$ (from the conservation of swimmers during tumble events: $\int d\mathbf{p}K(\mathbf{p}/\mathbf{p}^{\prime})=\int d\mathbf{p}^{\prime}K(\mathbf{p}/\mathbf{p}^{\prime})=1$) and

[TABLE]

Using the following identities in (15) [Doi & Edwards (1978); Messiah (1962)]

[TABLE]

and after simplification, we get:

[TABLE]

where $k_{lm}^{\prime}=\sqrt{l(l+1)-m(m+1)}$ and $k_{lm}^{\prime\prime}=\sqrt{l(l+1)-m(m-1)}$.

Now, multiplying the above equation with the conjugate spherical harmonic $Y_{r}^{s^{*}}(\mathbf{p})$, then integrating it over the unit sphere, and using the orthogonality of spherical harmonics, we obtain the following infinite sequence of linear equations for the $a_{l,m}$’s:

[TABLE]

In the above equation, the Clebsch-Gordan coefficients

[TABLE]

will be evaluated using the Racah formula [Arfken & Weber (1999); Doi & Edwards (1978); Messiah (1962)].

Using $\Omega_{0}=G(\mathbf{p})F(z)$ and (12) in the following normalization condition (see §2.1)

[TABLE]

we get

[TABLE]

Substituting (20), in (18), we have:

[TABLE]

An appropriately truncated version of (21) is solved numerically to obtain the $a_{l,m}$’s, subject to convergence of the orientation distribution.

Now substituting $\Omega_{0}=G(\mathbf{p};Pe_{z})F(z,t_{2})$ in the right hand side of (5), after simplification, we obtain:

[TABLE]

To evaluate $\partial G/\partial z$ in the equation above, we first differentiate (4) with respect to $z$:

[TABLE]

Substituting for $G$ on the right hand side, using (12), and using the expansion

[TABLE]

in (23), we get:

[TABLE]

Differentiating (19) with respect to $z$, using (24) and using the orthogonality of spherical harmonics, we get

[TABLE]

Now, following a procedure analogous to that leading to (21), the sequence of linear equations governing the $b_{l,m}$’s in (25) may be written as:

[TABLE]

where, $k_{l^{\prime}m^{\prime}}^{\prime}=\sqrt{l^{\prime}(l^{\prime}+1)-m^{\prime}(m^{\prime}+1)}$ and $k_{l^{\prime}m^{\prime}}^{\prime\prime}=\sqrt{l^{\prime}(l^{\prime}+1)-m^{\prime}(m^{\prime}-1)}$. One can solve the system of linear equations (27) to determine $dG/dz$ as a truncated version of (24).

To solve for $\Omega_{1}$ in (5), on account of linearity, we first write:

[TABLE]

and equating the coefficients of $\partial F/\partial z$ and $F$ on both sides, we get the following equations governing $\Omega_{11}$ and $\Omega_{12}$:

[TABLE]

Expanding

[TABLE]

[TABLE]

using (28) in (19) and using the orthogonality of spherical harmonics, we get

[TABLE]

From (12), (24), (19), (27), (31), (32), and following the numerical scheme adopted for solving (23), we obtain the unknown coefficients $c_{l,m}$ and $d_{l,m}$.

Using (28), (31), (32) in (2.1), integrating over the orientation degrees of freedom and using $\int d\mathbf{p}\Omega_{0}(\mathbf{p})=F$ (from (19)), we get the drift-diffusion equation describing the evolution of swimmer concentration profile

[TABLE]

in terms of the slow time variable $t_{2}$, and with the drift and diffusivity coefficients being defined by:

[TABLE]

and

[TABLE]

The drift-diffusion equation (34) appears as (7) in §2.1. The truncation in the above sequences of equations is such as to ensure converged values for the $D_{zz}$ and $V_{z}$. The truncated system has $l=L_{max}$. For $\kappa=1.25,~{}2$, and 4, reported in the manuscript, we chose $L_{max}=40$; for the other aspect ratios, we chose $L_{max}=40$ and 80 for $Pe<200$ and $Pe\geq 200$, respectively.

We calculate the steady state swimmer concentration profile satisfying (34), while imposing the zero-flux condition at boundaries, which leads to

[TABLE]

at $z=\pm 1$.

The solution of the above first order ordinary differential equation can be expressed as:

[TABLE]

which appears as (8) in §2.1. In (38), the normalization constant $\aleph$ can be calculated while imposing the condition$\int_{-1}^{1}dzF(z)=1$, which leads to:

[TABLE]

We evaluate the integrals in (38) and (39) using Simpson’s rule.

Appendix B Small-$Pe$ expansion for the swimmer concentration profile

Herein, we obtain closed form analytical expression for the swimmer concentration, to O($Pe^{2}$). The orientation probability density at each order (in $\epsilon$) in the multiple scales analysis (see §2.1) is further expanded in powers of $Pe$ to O($Pe^{2}$). Thus, the probability densities of O(1) ($G$ in (4)), and at O($\epsilon$) ($\Omega_{1}$ in (5)) are expanded to O($Pe^{2}$), with the imposition of the normalization condition: $\int d\textbf{p}\Omega=F(z)$. Here and in §C, we choose a coordinate system different from the one used in §2.1 (figure 1), such that the polar ($\theta$) and azimuthal ($\phi$) angles are measured from the negative vorticity and flow axes, respectively. We first consider the probability density at leading order in the multiple scales analysis, and expand it about the isotropic orientation distribution, for $Pe\ll 1$, as:

[TABLE]

with $G_{0}=1/4\pi$.

At successive orders in $Pe$, we obtain:

[TABLE]

We solve (41) for $G_{1}$, by using the modified Green’s function $\mathcal{G}^{I}_{M}(\mathbf{p}/\mathbf{p{{}^{\prime}}})$ [Subramanian & Koch (2009)] which satisfies:

[TABLE]

where $\delta(\mathbf{p})$ denotes the Dirac delta function in orientation space, with the forcing in (43) being a localized source at an orientation $\mathbf{p}^{\prime}$, together with a compensating uniformly distributed sink over the remainder of the unit sphere. The modified Green’s function result given in Subramanian & Koch (2009), derived for the rotary diffusion of swimmers can be easily generalized to the case of run-and-tumble plus rotary diffusion, so as to arrive at the following expression for $\mathcal{G}^{I}_{M}(\mathbf{p}/\mathbf{p^{\prime}})$ as:

[TABLE]

where $D_{n}=\left[\frac{1}{\tau D_{r}}(1-\frac{4\pi A_{n}}{2n+1})+n(n+1)\right].$

Using the expression for $\mathcal{G}^{I}_{M}(\mathbf{p}/\mathbf{p^{\prime}})$ in (44), the orientation probability density $G_{1}$, at O($Pe$), may be written as:

[TABLE]

Similarly, using the modified Green’s function, the solution of (42) may also be expressed as:

[TABLE]

Again, using (44), (45) in the above equation, expressing $\nabla_{p}\cdot(G_{1}\mathbf{\tilde{\dot{p}}})$ in terms of spherical harmonics, and using the orthogonality of the latter, we get:

[TABLE]

From (40), (45) and (47), we obtain the small $Pe$ of the leading order orientation probability density, $G(\mathbf{p})$, to O($Pe^{2}$) as:

[TABLE]

Now, considering the equation at O($\epsilon$), given by (5), and expanding the unknown probability density as:

[TABLE]

leads to the following three equations at successive orders in $Pe$:

[TABLE]

[TABLE]

The solutions of (50) and (B) may again be expressed as a convolution involving the modified Greens function. Using the modified Green’s function (44) and $G_{0}=1/4\pi$ in (50), we have:

[TABLE]

Using (45) for $G_{0}$, (53) for $\Omega_{10}$, and the modified Green’s function (44), the unknown probability density in (B) can be expressed as:

[TABLE]

Similarly, using (47) for $G_{2}$, (54) for $\Omega_{11}$ and the modified Green’s function (44), the unknown probability density in (52) can be expressed as:

[TABLE]

In the above equation, we only mention the coefficients of $(Y_{1}^{-1}+Y_{1}^{1})$, since only these contribute to the drift, diffusivity, and hence the swimmer concentration profile (see (57) below).

Considering (2.1), integrating over the orientation degrees of freedom and using the normalization condition $\int d\mathbf{p}\Omega_{0}=F$ (see (19)), we get:

[TABLE]

Expressing $p_{3}$ in terms of spherical harmonics, we have:

[TABLE]

From orthogonality, only the coefficients of $(Y_{1}^{-1}+Y_{1}^{1})$, in $\Omega_{10},\Omega_{11}$ and $\Omega_{12}$, survive. Substituting for $\Omega_{10},\Omega_{11}$, $\Omega_{12}$ from (53), (54), (55), the shear rate profile for the parabolic flow, $\dot{\gamma}=-2z$, and evaluating the orientation-space integrals, we get:

[TABLE]

where

[TABLE]

In (58), the odd orders in $Pe$ don’t contribute, as must be the case, so the concentration profile is invariant to flow reversal. We therefore expand the swimmer concentration as $F=F^{(0)}+Pe^{2}F^{(2)}+\cdots$. Substituting this expansion in (58), imposing the zero-flux condition at the boundaries ($z=\pm 1$), we get:

[TABLE]

At successive orders in $Pe$, we obtain:

[TABLE]

Using (60) in (61), the leading order inhomogeneity (${F}^{(2)}$) is independent of the correction to the diffusivity and the O($Pe^{2}$) drift alone dictates the inhomogeneity. It is nevertheless worth noting that $D_{zz}^{(2)}$ exhibits a profile consistent with physical arguments - the fact that it is negative, and larger closer to the walls, where the transverse diffusivity is suppressed due to increasing flow alignment of the swimmer. Although not relevant to the steady state, $D_{zz}^{(2)}$ will influence the transient evolution from a generic initial condition.

Solving (60) and (61) for $F^{(0)}$ and $F^{(2)}$, while imposing the normalization condition of swimmer concentration $\int_{-1}^{1}dzF(z)=\int_{-1}^{1}dz(F^{(0)}+Pe^{2}F^{(2)}+\cdots)=1$, we get:

[TABLE]

where

[TABLE]

with $B$ as the Bretherton constant and $D_{i}$ being given by (44). The comparison of the numerical and the analytical profiles of $F^{(2)}(z)$ for swimmers with $\kappa=15$ ($B=0.991$) and for different $Pe$ is given in figure 2a (inset) of §3.1.

Appendix C Concentration profile for the near-spherical swimmers

In the following, we obtain an analytical expression for the swimmer concentration profile, for an arbitrary $Pe$, but for near-spherical swimmers which correspond to an asymptotically small Bretherton constant ($B\ll 1$). As we did in §B, we start from the small-$\epsilon$ expansion of probability density ($\Omega$) in the multiple scales analysis, and expand the probability densities at each order in $\epsilon$, $G$ in (4) at O(1), and $\Omega_{1}$ in (5) at O($\epsilon$), to O($B$), while also imposing the normalization condition: $\int d\textbf{p}\Omega=F(z)$. To begin with, consider the governing equation (4) for $G$:

[TABLE]

where $\mathbf{\dot{\tilde{p}}^{I}}=\mathbf{W}\cdot\mathbf{p}/\dot{\gamma}$, $\mathbf{\dot{\tilde{p}}^{II}}=B[\mathbf{E}\cdot\mathbf{p}-\mathbf{p}(\mathbf{E}:\mathbf{p}\mathbf{p})]/\dot{\gamma}$, represent the rotation rates of the swimmer at O(1) and O($B$), respectively, due to the local ambient vorticity and rate of strain. In (64), for near-spherical swimmers, we expand $G$ in the parameter $B$, about the spatially homogeneous isotropically oriented base state

[TABLE]

At leading order, one obtains the isotropically oriented state of spherical swimmers that spin at a uniform rate with the ambient vorticity. Thus, $G_{0}=1/4\pi$. At O($B$), we have the following governing equation for $G_{1}$:

[TABLE]

To solve the above equation for $G_{1}$, we again seek a modified Green’s function $\mathcal{G}^{II}_{M}(\mathbf{p}/\mathbf{p}^{\prime})$, that satisfies

[TABLE]

Here, $\mathcal{G}^{II}_{M}(\mathbf{p}/\mathbf{p}^{\prime})$ differs from $\mathcal{G}^{I}_{M}(\mathbf{p}/\mathbf{p}^{\prime})$, given by (43), in that it includes the effect of vorticity-induced rotation in addition to the orientation decorrelation due to rotary diffusion and run-and-tumble dynamics. Writing down the modified Green’s function as a series in spherical harmonics

[TABLE]

substituting in (67), and following the procedure mentioned in appendix B, we obtain

[TABLE]

where $R_{n,m}=\left[\left(1-4\pi A_{n}/({2n+1})\right)/\tau D_{r}+n(n+1)-imPe_{z}/2\right].$ Note that the azimuthal degeneracy is broken by the presence of the $m$-dependent flow contribution. The effect of this flow contribution is to decrease the translational swimmer diffusivity at large $Pe$. For, large $Pe$, the spherical swimmer rotates with a period of O($D^{-1}_{r}Pe^{-1}$). The rotation and swimming lead to a mean free path of O($U_{s}Pe^{-1}D^{-1}_{r}$). Taken together with the decorrelation time of O($D_{r}^{-1}$), this implies a scaling of $(U_{s}^{2}/D_{r})Pe^{-2}$ for the gradient component of the diffusivity. Thus, the rapid rotation induced by the ambient vorticity leads to an asymptotically small translational diffusivity at large $Pe$.

Using the convolution integral involving the modified Green’s function given by (68), one may write the O($B$) swimmer probability density $G_{1}$, governed by (66), in the following final form:

[TABLE]

Next, considering the swimmer probability density at O($\epsilon$) in the multiple scales analysis, and expanding in $B$:

[TABLE]

we obtain the following equations at successive orders in $B$:

[TABLE]

Substituting $G_{0}=1/4\pi$ and using (67), the solution of (72) can be expressed as:

[TABLE]

Similarly, from (74), (70) and (67), we obtain the solution of (C) as:

[TABLE]

As in §B, only the terms involving $(Y_{1}^{-1}+Y_{1}^{1})$ contribute to the diffusivity ($D_{zz}$) and drift ($V_{z}$). So for simplicity, we do not mention the coefficients of $Y_{3}^{\pm 1},Y_{3}^{\pm 3}$ in (75). Integrating (2.1) over the orientation degrees of freedom and using $\int d\mathbf{p}G=1$ (from (19)), we get:

[TABLE]

Using (71), (74), (75), in the above equation and integrating, gives:

[TABLE]

where

[TABLE]

In the above equation, $D^{(0)}_{zz}$ scales as O($1/Pe^{2}$) for large $Pe$, consistent with the scaling arguments above. From (81), we see that the O($B$) drift at the wall, for near-spherical swimmers, changes sign when $Pe_{z}^{2}=84$. This may be rewritten as $\dot{\gamma}|_{z=-1}^{2}Pe^{2}=84$, whence, one obtains the threshold $Pe$ for drift reversal as: $Pe_{c}=\sqrt{21}\sim 4.5$. This critical $Pe$ is mentioned in §3.3.2, and appears in the $Pe-\kappa$ phase portrait in figure 13.

Now, to obtain the steady state swimmer concentration profile we expand

[TABLE]

substituting in the normalization condition $\int_{-1}^{1}dzF(z)=1,$ yields $F^{(0)}=1/2$ and

[TABLE]

Using $F^{(0)}=1/2$, substituting (82) in (77), and equating the coefficients of $B$ on both sides, we get:

[TABLE]

From (78), (81) and (83), the solution of (84), is:

[TABLE]

Thus from (82), (85) and using $F^{(0)}=1/2$, we obtain the following steady state concentration profile for near-spherical swimmers:

[TABLE]

In figure 14, we compare the numerical and analytical results of $F^{(1)}(z)$ for near-spherical swimmers at $Pe=40$. Expectedly, the numerical profiles approach the analytical one for $\kappa\rightarrow 1$. More importantly, the spatial structure of the swimmer concentration profile is similar to that obtained for finite $B$- a pair of concentration maxima on either side of the channel centerline.