2.1 Theoretical framework: the method of multiple scales

The time evolution of the probability density $\unicode[STIX]{x1D6FA}(\boldsymbol{x},\boldsymbol{p},t)$ for a dilute swimmer suspension subject to a shearing flow (see figure 1) is given by (Subramanian & Koch Reference Subramanian and Koch2009)

(2.1)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }((\boldsymbol{u}+U_{s}\boldsymbol{p})\unicode[STIX]{x1D6FA})+\unicode[STIX]{x1D735}_{p}\boldsymbol{\cdot }(\dot{\boldsymbol{p}}\unicode[STIX]{x1D6FA})+\frac{1}{\unicode[STIX]{x1D70F}}\left[\unicode[STIX]{x1D6FA}-\int \text{d}\,\boldsymbol{p}^{\prime }\,K(\,\boldsymbol{p}|\boldsymbol{p}^{\prime })\unicode[STIX]{x1D6FA}(\,\boldsymbol{p}^{\prime })\right]-D_{r}\unicode[STIX]{x1D6FB}_{p}^{2}\unicode[STIX]{x1D6FA}=0,\end{eqnarray}$$

where $\boldsymbol{x}$ and $\boldsymbol{p}$ denote the swimmer position and orientation, respectively. The second term on the left-hand side of (2.1) denotes the spatial convection of probability density by both the fluid flow $\boldsymbol{u}$ and swimming velocity $U_{s}\boldsymbol{p}$, and the term involving $\dot{\boldsymbol{p}}$ denotes swimmer rotation by the ambient shear. The terms within brackets model a run-and-tumble process (Berg Reference Berg1993; Subramanian & Koch Reference Subramanian and Koch2009; Nambiar et al. Reference Nambiar, Nott and Subramanian2017), obeying Poisson statistics with a mean run duration $\unicode[STIX]{x1D70F}$; the kernel $K(\,\boldsymbol{p}|\boldsymbol{p}^{\prime })$ characterizes correlation between the pre ($\boldsymbol{p}^{\prime }$)- and post ($\boldsymbol{p}$)-tumble orientations, with $K(\,\boldsymbol{p}|\boldsymbol{p}^{\prime })=1/(4\unicode[STIX]{x03C0})$ for random tumbles. The last term on the left side of (2.1) models the stochastic orientation change due to rotary diffusion, $D_{r}$ being the rotary diffusivity. In writing (2.1), we have neglected any sources of translational diffusion other than that arising from swimming (for instance, a passive Brownian diffusivity). We examine the role of an additional source of translational diffusion in § 3.4.

In the experiments of Rusconi et al. (Reference Rusconi, Guasto and Stocker2014b) and Barry et al. (Reference Barry, Rusconi, Guasto and Stocker2015), 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{\boldsymbol{p}}=B[\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{p}-\boldsymbol{p}(\boldsymbol{E}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p})]+\boldsymbol{W}\boldsymbol{\cdot }\boldsymbol{p}$. Here, $B=(\unicode[STIX]{x1D705}^{2}-1)/(\unicode[STIX]{x1D705}^{2}+1)$ is the Bretherton constant, $E_{ij}=\dot{\unicode[STIX]{x1D6FE}}(z)(\unicode[STIX]{x1D6FF}_{i1}\unicode[STIX]{x1D6FF}_{j3}+\unicode[STIX]{x1D6FF}_{i3}\unicode[STIX]{x1D6FF}_{j1})/2$ and $W_{ij}=\dot{\unicode[STIX]{x1D6FE}}(z)(\unicode[STIX]{x1D6FF}_{i1}\unicode[STIX]{x1D6FF}_{j3}-\unicode[STIX]{x1D6FF}_{i3}\unicode[STIX]{x1D6FF}_{j1})/2$ are the strain rate and vorticity tensors, with $\dot{\unicode[STIX]{x1D6FE}}(z)=\text{d}u/\text{d}z$ 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, $\unicode[STIX]{x1D6FF}_{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 $\unicode[STIX]{x1D705}$ (Leal & Hinch Reference Leal and Hinch1971), and it is via $\unicode[STIX]{x1D705}$ 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. Reference Darnton, Turner, Rojevsky and Berg2007; Das & Lauga Reference Das and Lauga2018), on account of the rapid (counter)-rotation on time scales much shorter than those characterizing shear-induced migration.

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

(2.2)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FA}\boldsymbol{p})+Pe\dot{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D735}_{p}\boldsymbol{\cdot }(\tilde{\dot{\boldsymbol{p}}}\unicode[STIX]{x1D6FA})+\frac{1}{\unicode[STIX]{x1D70F}D_{r}}\left[\unicode[STIX]{x1D6FA}-\int \text{d}\,\boldsymbol{p}^{\prime }\,K(\,\boldsymbol{p}/\boldsymbol{p}^{\prime })\unicode[STIX]{x1D6FA}(\,\boldsymbol{p}^{\prime })\right]-\unicode[STIX]{x1D6FB}_{p}^{2}\unicode[STIX]{x1D6FA}=0,\end{eqnarray}$$

where $\tilde{\dot{\boldsymbol{p}}}=\dot{\boldsymbol{p}}/\dot{\unicode[STIX]{x1D6FE}}(z)$ is a non-dimensional rotation rate. Note that convection by the fluid flow has now been omitted since we only examine the variations of the swimmer concentration in the gradient direction. The concentration variation being orthogonal to the direction of the fluid flow (either imposed or induced, the latter arising from the back coupling due to a swim stress; see Vennamneni et al. (Reference Vennamneni, Garg and Subramanian2018)), this convection term is identically zero. In (2.2), $\unicode[STIX]{x1D716}$, Pe and $\unicode[STIX]{x1D705}$ are dimensionless parameters that primarily determine the different regimes of shear-induced migration. The parameter $\unicode[STIX]{x1D716}=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. Reference Takatori, Yan and Brady2014; Takatori & Brady Reference Takatori and Brady2017), $\unicode[STIX]{x1D716}^{-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 $\unicode[STIX]{x1D70F}D_{r}$ determines the roles of tumbling vis-a-vis rotary diffusion in the swimmer orientation dynamics. Active Brownian particles (ABPs) correspond to $\unicode[STIX]{x1D70F}D_{r}\rightarrow \infty$, and run-and-tumble particles (RTPs) to $\unicode[STIX]{x1D70F}D_{r}\rightarrow 0$ (Tailleur & Cates Reference Tailleur and Cates2009; Saintillan Reference Saintillan2010b); in the latter case, the swimmer Knudsen number is defined as $\unicode[STIX]{x1D716}=U\unicode[STIX]{x1D70F}/H$ instead. Regardless of the particular value of $\unicode[STIX]{x1D70F}D_{r}$, swimmers sample the channel cross-section diffusively for long times, the translational diffusivity being given by $D_{t}=(U_{s}^{2}/(6D_{r}))(\unicode[STIX]{x1D70F}D_{r}/(0.5+\unicode[STIX]{x1D70F}D_{r}))$ (Berg Reference Berg1993; Koch & Subramanian Reference Koch and Subramanian2011; Subramanian & Nott Reference Subramanian and Nott2012) in the general case; for $\unicode[STIX]{x1D70F}D_{r}\rightarrow \infty$, $D_{t}=U_{s}^{2}/(6D_{r})$ as above. In the experiments of Rusconi et al. (Reference Rusconi, Guasto and Stocker2014b), $\unicode[STIX]{x1D716}\sim 0.23$ for the bacterium B. subtilis. The modest value of $\unicode[STIX]{x1D716}$ implies that the time scales characterizing the swimmer orientation dynamics (of $O(\unicode[STIX]{x1D70F})$ or $O(D_{r}^{-1})$) and the diffusive sampling of the channel cross-section (of $O(H^{2}/U_{s}^{2})\unicode[STIX]{x1D70F}$ or $O(H^{2}D_{r}/U_{s}^{2})$) are well separated, their ratio being of $O(\unicode[STIX]{x1D716}^{2})$. Further, since the experimental range of flow rates corresponds to $U_{s}\ll U_{f}(\unicode[STIX]{x1D716}\ll Pe)$, we use the method of multiple scales to analyse the swimmer concentration profiles (Nitsche & Hinch Reference Nitsche and Hinch1997; Subramanian & Brady Reference Subramanian and Brady2004; Kasyap & Koch Reference Kasyap and Koch2012, Reference Kasyap and Koch2014; Vennamneni et al. Reference Vennamneni, Garg and Subramanian2018) for small $\unicode[STIX]{x1D716}$, but with $Pe$ and $\unicode[STIX]{x1D705}$ being arbitrary.

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

(2.3)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}t_{1}}+\unicode[STIX]{x1D716}^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}t_{2}}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FA}\boldsymbol{p})+Pe_{z}\unicode[STIX]{x1D735}_{p}\boldsymbol{\cdot }(\tilde{\dot{\boldsymbol{p}}}\unicode[STIX]{x1D6FA})+\frac{1}{\unicode[STIX]{x1D70F}D_{r}}\left[\unicode[STIX]{x1D6FA}-\int \text{d}\,\boldsymbol{p}^{\prime }\,K(\,\boldsymbol{p}/\boldsymbol{p}^{\prime })\unicode[STIX]{x1D6FA}(\,\boldsymbol{p}^{\prime })\right]\nonumber\\ \displaystyle & & \displaystyle \quad -\,\unicode[STIX]{x1D6FB}_{p}^{2}\unicode[STIX]{x1D6FA}=0,\end{eqnarray}$$

where $Pe_{z}=Pe\dot{\unicode[STIX]{x1D6FE}}(z)$ is the local Péclet number. Now, expanding $\unicode[STIX]{x1D6FA}$ as $\unicode[STIX]{x1D6FA}(z,\boldsymbol{p},t;\unicode[STIX]{x1D716})=\unicode[STIX]{x1D6FA}_{0}(z,\boldsymbol{p},t_{1},t_{2})+\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{1}(z,\boldsymbol{p},t_{1},t_{2})+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FA}_{2}(z,\boldsymbol{p},t_{1},t_{2})+\cdots \,$, yields the following set of equations at successive orders in $\unicode[STIX]{x1D716}$

(2.4)$$\begin{eqnarray}\displaystyle & \hspace{-8.39996pt}\displaystyle \mathit{O}(1):Pe_{z}\unicode[STIX]{x1D735}_{p}\boldsymbol{\cdot }(\,\tilde{\dot{\boldsymbol{p}}}G)+\frac{1}{\unicode[STIX]{x1D70F}D_{r}}\left[G-\int \text{d}\,\boldsymbol{p}^{\prime }\,K(\,\boldsymbol{p}/\boldsymbol{p}^{\prime })G(\,\boldsymbol{p}^{\prime })\right]-\unicode[STIX]{x1D6FB}_{p}^{2}G=0, & \displaystyle\end{eqnarray}$$

(2.5)$$\begin{eqnarray}\displaystyle & \hspace{-6.0pt}\displaystyle \mathit{O}(\unicode[STIX]{x1D716}):Pe_{z}\unicode[STIX]{x1D735}_{p}\boldsymbol{\cdot }(\,\tilde{\dot{\boldsymbol{p}}}\unicode[STIX]{x1D6FA}_{1})+\frac{1}{\unicode[STIX]{x1D70F}D_{r}}\!\left[\unicode[STIX]{x1D6FA}_{1}-\int \text{d}\,\boldsymbol{p}^{\prime }\,K(\,\boldsymbol{p}/\boldsymbol{p}^{\prime })\unicode[STIX]{x1D6FA}_{1}(\,\boldsymbol{p}^{\prime })\right]-\unicode[STIX]{x1D6FB}_{p}^{2}\unicode[STIX]{x1D6FA}_{1}=-\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FA}_{0}\boldsymbol{p}),\quad & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.6)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \mathit{O}(\unicode[STIX]{x1D716}^{2}):Pe_{z}\unicode[STIX]{x1D735}_{p}\boldsymbol{\cdot }(\,\tilde{\dot{\boldsymbol{p}}}\unicode[STIX]{x1D6FA}_{2})+\frac{1}{\unicode[STIX]{x1D70F}D_{r}}\left[\unicode[STIX]{x1D6FA}_{2}-\int \text{d}\,\boldsymbol{p}^{\prime }\,K(\,\boldsymbol{p}/\boldsymbol{p}^{\prime })\unicode[STIX]{x1D6FA}_{2}(\,\boldsymbol{p}^{\prime })\right]-\unicode[STIX]{x1D6FB}_{p}^{2}\unicode[STIX]{x1D6FA}_{2}\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{0}}{\unicode[STIX]{x2202}t_{2}}-\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FA}_{1}\boldsymbol{p}).\end{eqnarray}$$

Here, we have written $\unicode[STIX]{x1D6FA}_{0}$ as $\unicode[STIX]{x1D6FA}_{0}=F(z,t_{2})G(\,\boldsymbol{p};Pe_{z})$, with $\int \text{d}\,\boldsymbol{p}\,\unicode[STIX]{x1D6FA}(z,\boldsymbol{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 ($\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t_{1}$) has been omitted in (2.4)–(2.6) . Writing $\unicode[STIX]{x1D6FA}_{1}$ in (2.5) as $\unicode[STIX]{x1D6FA}_{11}\unicode[STIX]{x2202}F/\unicode[STIX]{x2202}z+\unicode[STIX]{x1D6FA}_{12}F$ on account of linearity, and using this form in (2.6) and integrating over orientation space, leads to a drift-diffusion equation for $F$

(2.7)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}F(z,t_{2})}{\unicode[STIX]{x2202}t_{2}}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left(D_{zz}(z;Pe,\unicode[STIX]{x1D705})\frac{\unicode[STIX]{x2202}F(z,t_{2})}{\unicode[STIX]{x2202}z}-V_{z}(z;Pe,\unicode[STIX]{x1D705})F(z,t_{2})\right),\end{eqnarray}$$

where the diffusivity $D_{zz}=-\int \text{d}\,\boldsymbol{p}\cos \unicode[STIX]{x1D703}\unicode[STIX]{x1D6FA}_{11}$ and drift $V_{z}=\int \text{d}\,\boldsymbol{p}\cos \unicode[STIX]{x1D703}\unicode[STIX]{x1D6FA}_{12}$. The diffusivity and drift in (2.7) are scaled with $U_{s}^{2}/D_{r}$ and $U_{s}^{2}/D_{r}H$, respectively. At steady state, for impermeable walls

(2.8)$$\begin{eqnarray}F(z)=\aleph \exp \left(\int _{-1}^{z}\,\text{d}z^{\prime }\frac{V_{z}(z^{\prime })}{D_{zz}(z^{\prime })}\right),\end{eqnarray}$$

where

$$\begin{eqnarray}\aleph =\frac{1}{\int _{-1}^{1}\,\text{d}z\exp \left(\displaystyle \int _{-1}^{z}\,\text{d}z^{\prime }{\displaystyle \frac{V_{z}(z^{\prime })}{D_{zz}(z^{\prime })}}\right)}\end{eqnarray}$$

is determined from the normalization condition $\int _{-1}^{1}\,\text{d}z\,F(z)=1$. Determining $D_{zz}$ and $V_{z}$ requires knowledge of the orientation distribution as a function of $Pe$ and $\unicode[STIX]{x1D705}$. For arbitrary $Pe$, the orientation distributions, at different orders in $\unicode[STIX]{x1D716}$, are written as a truncated series in spherical harmonics, with the resulting coefficients being obtained numerically (Nambiar et al. Reference Nambiar, Garg and Subramanian2019a) (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. (Reference Rusconi, Guasto and Stocker2014b), the swimmers are uniformly distributed along $z$ at the initial time, with their orientations following an isotropic distribution. Thereafter, at each time step $\unicode[STIX]{x1D6FF}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 modelled using the method outlined in Ghosh, Samuel & Sinha (Reference Ghosh, Samuel and Sinha2012), whereas run-and-tumble motion is implemented in accordance with Poisson statistics (Berg Reference Berg1993; Krishnamurthy & Subramanian Reference Krishnamurthy and Subramanian2015). The flow-induced rotation at any $z$ is determined using (2.2). Following the update in the swimmer orientations, their positions are updated using a second-order Runge–Kutta 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 18 000 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 $\unicode[STIX]{x1D716}=0.05$, and find this to be sufficient to capture the $\unicode[STIX]{x1D716}\ll 1$ regime that pertains to the multiple scale analysis.

3.1 Swimmer migration towards walls – high-shear trapping

Figure 2.

Swimmer concentration profiles for different Pe: (a) $\unicode[STIX]{x1D705}=\infty$, (b) $\unicode[STIX]{x1D705}=15$. The red-dashed line in (a) denotes the $Pe\rightarrow \infty$ asymptote; the green colour line with circles and magenta colour line with hexagram symbols denote the simulation results for $Pe=10$ and 40, respectively. The inset in (a) highlights the concentration profiles, at the three largest Pe numbers, in the vicinity of the channel centreline. The inset in (b) is the comparison of the numerical profiles for $F^{(2)}(z)=(F-F^{(0)})/Pe^{2}$, with the parabolic depletion predicted by the small-$Pe$ analysis.

In figure 2, we plot the normalized steady state concentration profile [$\bar{F}(z)=F(z,t_{2}\rightarrow \infty )/F(z,0)$], using (2.8) with $D_{zz}$ and $V_{z}$ determined numerically as mentioned in § 2.1, for infinitely slender swimmers ($\unicode[STIX]{x1D705}=\infty$), and for $\unicode[STIX]{x1D705}=15$, over a range of $Pe$. Also shown are the profiles obtained from Langevin simulations, for $\unicode[STIX]{x1D705}=\infty$, for $Pe=10$ and 40, which are in good agreement with the predictions of the multiple scale analysis. All profiles exhibit near-centre 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)+\cdots \,$, 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 2(b), for $\unicode[STIX]{x1D705}=15$ and $Pe\leqslant 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 centreline towards the walls), and for $\unicode[STIX]{x1D705}=\infty$, the profiles appear to asymptote to a singular cusped profile in the limit $Pe\rightarrow \infty$ (see figure 2a). The mechanism underlying near-centre depletion is readily explained. For small $\unicode[STIX]{x1D716}$, 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 centreline (where $Pe_{z}=0$) as shown in figure 3(a,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.

Figure 3.

A schematic highlighting the physical mechanism underlying the wallward migration in the high-shear trapping regime.

Figure 4.

The swimmer diffusivity at the wall $D_{zz}|_{z=-1}$, as a function of $Pe$, for swimmers with $\unicode[STIX]{x1D705}=\infty$ (a), $\unicode[STIX]{x1D705}=15$ (b). The variation of swimmer drift at the wall $V_{z}|_{z=-1}$, against $Pe$, for swimmers with $\unicode[STIX]{x1D705}=\infty$ (c), $\unicode[STIX]{x1D705}=15$ (d).

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 $\unicode[STIX]{x1D705}=\infty$ and 15. For $\unicode[STIX]{x1D705}=\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 $\unicode[STIX]{x1D705}=15$ too, there is an intermediate range of $Pe$ values (${\sim}$ 3 to a little greater than 10) where the aforementioned scaling holds, but there is a deviation at larger $Pe$ values. 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_{sz}^{2}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 Reference Hinch and Leal1972a; Brenner Reference Brenner1974); $t_{c}\sim \langle \unicode[STIX]{x1D703}^{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}=\unicode[STIX]{x0394}z_{+}f_{+}-\unicode[STIX]{x0394}z_{-}f_{-}$, where $\unicode[STIX]{x0394}z_{+}$, $\unicode[STIX]{x0394}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+\unicode[STIX]{x0394}z^{\prime })=U_{s}(Pe_{z}^{-1/3}-Pe_{z+\unicode[STIX]{x0394}z^{\prime }}^{-1/3})$, where $\unicode[STIX]{x0394}z^{\prime }$ is the (dimensionless) gradient projection of the swimmer mean free path. With a swimming velocity of $O(U_{s}Pe_{z}^{-1/3})$ and a decorrelation time of $O(Pe_{z}^{-2/3}D_{r}^{-1})$, one obtains $\unicode[STIX]{x0394}z^{\prime }=U_{s}H^{-1}Pe_{z}^{-1}D_{r}^{-1}$. Using this leads to $V_{z}\sim O(U_{s}^{2}Pe_{z}^{-4/3}H^{-1}D_{r}^{-1}/z)$.

Figure 5.

Comparison of $|V_{z}|/D_{zz}$ against $z$, at different $Pe$ values for swimmers with $\unicode[STIX]{x1D705}=\infty$. The black-solid line represents the large $Pe$ asymptotic profile of $|V_{z}|/D_{zz}$.

Figure 6.

Comparison of normalized concentration profiles: our theory (red line), Langevin simulation (blue line with circles) with the experiments (dotted line in magenta colour) and simulations (black-dashed line) of Rusconi et al. (Reference Rusconi, Guasto and Stocker2014b);$\unicode[STIX]{x1D705}=10$.

The identical $Pe$-scalings of the drift and diffusivity 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 $\unicode[STIX]{x1D705}=\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 centreline leads to $C=0.344$. From figure 2(a), the swimmer concentration profiles are seen to approach the above $Pe$-independent limiting form for $Pe\geqslant 40$. Although, the steady state profile is $Pe$-independent, the time taken to attain steady state is of $O((H^{2}D_{r}/U_{s}^{2})Pe^{4/3})$ and diverges for $Pe\rightarrow \infty$. The limiting concentration profile is singular with a cusp at the centreline ($z=0$). The cusp is, of course, an artefact of the $Pe=\infty$ limit. For any finite $Pe$, the multiple scale analysis breaks down in an asymptotically small interval of $O(U_{s}/D_{r})$ in the vicinity of the centreline (keeping in mind our interpretation of $\unicode[STIX]{x1D716}$ as a swimmer Knudsen number, this layer in the vicinity of the centreline, 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 (Reference Zöttl and Stark2013).

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. (Reference Rusconi, Guasto and Stocker2014b) for $\unicode[STIX]{x1D705}=10$. For the three smallest $Pe$ values (1.25, 2.5 and 5), our theory underpredicts the inhomogeneity compared to the experiments. There is improved agreement at $Pe=10$, and again, at the highest $Pe$ values (25 and 50), the disagreement with the experiments widens. The disagreement at the two highest $Pe$ values, we believe, is because the experiments (Rusconi et al. Reference Rusconi, Guasto and Stocker2014b) 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 Reference Sokolov and Aranson2016; Sokolov et al. Reference Sokolov, Rubio, Brady and Aranson2018), 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_{s}^{2})$. Here, $\bar{u}=(2\bar{\dot{\unicode[STIX]{x1D6FE}}}H/3)$ is the average suspension velocity and $\bar{\dot{\unicode[STIX]{x1D6FE}}}=PeD_{r}$ is 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~\unicode[STIX]{x03BC}\text{m}$, $U_{s}=50~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$, $D_{r}=1~\text{s}^{-1}$ (Rusconi et al. Reference Rusconi, Guasto and Stocker2014b), the entrance length at $Pe=50$ is 1.4 m! For $Pe=10$, the entrance length is approximately 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$ values (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 hydrodynamic 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. Reference Berke, Turner, Berg and Lauga2008; Lauga & Powers Reference Lauga and Powers2009; Chilukuri et al. Reference Chilukuri, Collins and Underhill2014). Wall accumulation in boundary layers with polar orientational order arises even in the absence of hydrodynamic interactions (Elgeti & Gompper Reference Elgeti and Gompper2009, Reference Elgeti and Gompper2016). This mechanism is expected to prevail in relatively confined geometries, and for modest $Pe$, owing to its reliance on polar order (the multiple scale analysis here predicts the polar order to be $O(\unicode[STIX]{x1D716})$ smaller than the nematic order associated with the leading-order orientation distribution $G(\,\boldsymbol{p})$). In the present context, we expect the hydrodynamic wall accumulation to only have a minor influence for the following 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_{r}^{-1})$ 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. (Reference Berke, Turner, Berg and Lauga2008), 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. (Reference Chilukuri, Collins and Underhill2014) 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. (Reference Rusconi, Guasto and Stocker2014b) and Barry et al. (Reference Barry, Rusconi, Guasto and Stocker2015) 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$ values, pointed out above, could be due to multiple reasons. First is the possibility of the assumption $\unicode[STIX]{x1D716}\ll 1$, underlying the multiple scale 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), $\unicode[STIX]{x1D716}\approx 0.23$ in the experiments of Rusconi et al. (Reference Rusconi, Guasto and Stocker2014b). Our simulations, which are not limited by the small $\unicode[STIX]{x1D716}$ assumption, agree well with theory even for $\unicode[STIX]{x1D716}=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 $\unicode[STIX]{x1D716}$ not being small enough is not relevant to the largest $Pe$ values 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 scale analysis. For the high-shear trapping regime considered in this section, the effective Knudsen number is $\unicode[STIX]{x1D716}_{eff}=\unicode[STIX]{x1D716}Pe^{-1/3}$, the exponent reflecting the degree of swimmer alignment in the orientation boundary layer; for low-shear trapping examined next, $\unicode[STIX]{x1D716}_{eff}=\unicode[STIX]{x1D716}Pe^{-1}$ for large $Pe$, being even smaller. 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 RTPs and ABPs. As explained in the introduction, the former involves the swimmer undergoing a sudden large change in its orientation during tumbles (Berg Reference Berg1993, Reference Berg2008), while the latter involves continuous small amplitude fluctuations in their orientation. Our Langevin simulations do suggest a greater inhomogeneity for RTPs than ABPs at the smallest $Pe$ values, 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$ values (the Rusconi et al. (Reference Rusconi, Guasto and Stocker2014b) 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 centreline – low-shear trapping

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

Figure 7.

Swimmer concentration profiles as a function of $Pe$ for different aspect ratios: (a) $\unicode[STIX]{x1D705}=1.25$, (b) 2, (c) 10, (d) 15. In (a), the red-dashed line is the large $Pe$ limiting profile, the solid-black and blue lines with circles represent the concentration profiles obtained from the Langevin simulations for $Pe=4$ and 40, respectively.

Figure 8.

The variation of (a,b) the diffusivity ($D_{zz}$) and (c,d) the absolute drift velocity ($|V_{z}|$) at the channel wall, as a function of $Pe$, for swimmers with different aspect ratios. Here, the $Pe^{4/3}$ and $Pe^{2}$ pre-factors, used to rescale $D_{zz}$ and $V_{z}$, are motivated by the scalings in the high-shear ($Pe^{-4/3}$) and low-shear ($Pe^{-2}$) trapping regimes.

Figure 9.

The variation of (a) $D_{zz}$ and (b) $|V_{z}|$ against $Pe$ for swimmers with $\unicode[STIX]{x1D705}=10$, at different locations across the channel. The black (thin) and blue (thick) lines in (a) and (b), represent the scalings of $D_{zz}$ and $|V_{z}|$ in the high-shear ($Pe^{-4/3}$) and low-shear ($Pe^{-2}$) trapping regimes, respectively. The variation of (c) $D_{zz}$ and (d) $|V_{z}|$ for $\unicode[STIX]{x1D705}=10$ against $0.5Pe_{z}$, for different $Pe$; here, the factor 0.5 in front of $Pe_{z}$ is to ensure that $0.5Pe_{z}$ at the boundary is equal to $Pe$.

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}$, each of them scaled with $Pe^{-4/3}$ and $Pe^{-2}$, against $Pe$ for $\unicode[STIX]{x1D705}=\infty$ and for the $\unicode[STIX]{x1D705}$ values in figure 7. From figure 8(c,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$ values smaller than that corresponding to the zero-crossing (of $V_{z}$), swimmers migrate towards the walls; for larger $Pe$ values, swimmers migrate towards the centreline, leading to a low-shear trapping behaviour. It is important to note that the reversal in migration is actually a function of $Pe_{z}$, occurring for $Pe_{z}$ exceeding a $\unicode[STIX]{x1D705}$-dependent threshold. This is confirmed in figure 9, where the plots of the diffusivity and drift against $Pe$ for $\unicode[STIX]{x1D705}=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 centreline (figure 9a,b). The reversal in drift is also delayed at locations closer to the channel centreline, 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 figure 9(c,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 centreline, there is always a region close to the centreline, 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 centreline regardless of $z$ (this singular profile is subject to the limitations of the multiple scale analysis mentioned earlier in § 3.1). As evident from figure 8, for any finite $\unicode[STIX]{x1D705}$, 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 $\unicode[STIX]{x1D705}$ (10 and 15), there is the emergence of an intermediate scaling regime (evident from figure 8a,c) where $D_{zz},V_{z}\sim Pe^{-4/3}$, already seen in § 3.1. Thus, $\unicode[STIX]{x1D705}=\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 behaviour. Of course, for large enough $\unicode[STIX]{x1D705}$ values, 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 8(b) for $\unicode[STIX]{x1D705}=15$, and in figure 8(d) for both $\unicode[STIX]{x1D705}=10$ and $\unicode[STIX]{x1D705}=15$; owing to the zero-crossing, it takes a larger $Pe$ for the drift to conform to the eventual $O(Pe^{-2})$ scaling behaviour.

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 $\unicode[STIX]{x1D705}$ values of order unity. For large $\unicode[STIX]{x1D705}$, the cross-over $Pe$ is asymptotically large, and marks a transition from an $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 Reference Leal and Hinch1971; Hinch & Leal Reference Hinch and Leal1972a; Brenner Reference Brenner1974), 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{\unicode[STIX]{x1D6FE}}^{-1}\unicode[STIX]{x1D705})$ is comparable to that characterizing rotary diffusion of $O(D_{r}^{-1}Pe^{-2/3})$; this leads to $Pe_{c}\sim O(\unicode[STIX]{x1D705}^{3})$ for $Pe\gg 1$ (this scaling is confirmed in § 3.3.2, where we demarcate high- and low-shear trapping on the $Pe{-}\unicode[STIX]{x1D705}$ 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 $\unicode[STIX]{x0394}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_{r}^{-1})$, lead to a $D_{zz}\sim \langle \unicode[STIX]{x0394}z^{2}\rangle /t_{c}$ which may be written as $O((U_{s}^{2}/D_{r})Pe^{-2})$, confirming the scaling behaviour observed in the low-shear trapping regime; here, $\langle \unicode[STIX]{x0394}z^{2}\rangle$ is the mean square displacement of the swimmer in the gradient direction. Note that, for large $\unicode[STIX]{x1D705}$, swimmers spend an $O(\unicode[STIX]{x1D705}/(D_{r}Pe))$ time remaining nearly aligned with the flow axis at an angle of $O(\unicode[STIX]{x1D705}^{-1})$. This leads to the same estimate for $\unicode[STIX]{x0394}z$ as above, but a faster decorrelation rate of $O(\unicode[STIX]{x1D705}^{2}D_{r})$ (since rotary diffusion only has to induce an angular displacement of $O(\unicode[STIX]{x1D705}^{-1})$), leading to a scaling for $D_{zz}$ of $O((U_{s}^{2}\unicode[STIX]{x1D705}^{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+\unicode[STIX]{x0394}z^{\prime })$. Now, for swimmers with $\unicode[STIX]{x1D705}$ of $O(1)$, using $U_{sz}(z)=U_{s}Pe_{z}^{-1}$, we obtain $V_{z}\sim$$O((U_{s}^{2}H^{-1}D_{r}^{-1}/z)Pe_{z}^{-2})$. The drift and diffusivity scalings above are broadly consistent with those obtained in figure 8(b).

Figure 10.

Variation of swimmer concentration at the wall ($\bar{F}(-1)$) against $Pe$ for swimmers with different aspect ratios.

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(\unicode[STIX]{x1D705})/z$, where $C$ is now a function of $\unicode[STIX]{x1D705}$. One therefore expects a limiting $Pe$-independent profile of the form $(C+1)|z|^{C}$; note that, for any finite $Pe$, the actual profiles break away from this limiting profile in the vicinity of the centreline. The actual profiles exhibit a dip in this region, corresponding to the onset of depletion (wallward migration) at small $Pe_{z}$. Since $C(\unicode[STIX]{x1D705})<0$ for low-shear trapping (owing to the sign reversal of the drift), this limiting profile has a diverging concentration at the centreline. For $\unicode[STIX]{x1D705}=1.25$, $C(\unicode[STIX]{x1D705})\sim -0.238$, and the associated limiting profile compares well to the numerically computed ones in figure 7(a). An interesting consequence of the change in sign of $C(\unicode[STIX]{x1D705})$, leading to low-shear trapping, is that the aforementioned limiting profile fails to be integrable for $C(\unicode[STIX]{x1D705})\leqslant -1$ (this was not 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(\unicode[STIX]{x1D705})$ to equal this threshold value ($-1$) for $\unicode[STIX]{x1D705}\approx 2$, and decrease further for larger $\unicode[STIX]{x1D705}$. The implication of the non-integrability above is that, for $\unicode[STIX]{x1D705}\geqslant 2$, the swimmer concentration profiles must exhibit a ‘centreline collapse’ in the limit $Pe\rightarrow \infty$; in other words, $\lim _{Pe\rightarrow \infty ,\unicode[STIX]{x1D705}\geqslant 2}\bar{F}(z)=\unicode[STIX]{x1D6FF}(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 $\unicode[STIX]{x1D705}=1.25$, the near-wall swimmer concentration approaches a finite limiting value $(\lim _{Pe\rightarrow \infty ,\unicode[STIX]{x1D705}=1.25}\bar{F}(-1)=0.765)$. For $\unicode[STIX]{x1D705}\geqslant 2$, however, the near-wall swimmer concentration decreases monotonically to zero for $Pe\rightarrow \infty$, in accordance with the centreline collapse hypothesized above. The slope characterizing this decrease appears to increase with increasing $\unicode[STIX]{x1D705}$, implying that the collapse increases in intensity for larger $\unicode[STIX]{x1D705}$ (for $\unicode[STIX]{x1D705}=15$, the numerics have not accessed the asymptotic regime yet). The $\unicode[STIX]{x1D705}$-dependence of the collapse does suggest that the drift and diffusivity coefficients must scale differently with $\unicode[STIX]{x1D705}$ for $\unicode[STIX]{x1D705}\gg 1$. Scaling arguments presented above suggested that $D_{zz}\sim O(\unicode[STIX]{x1D705}^{2}Pe^{-2})$ in the low-shear trapping regime, and the large$-\unicode[STIX]{x1D705}$ scaling of the drift must therefore differ from $O(\unicode[STIX]{x1D705}^{2}Pe_{z}^{-2}/z)$ that would emerge from the biased random walk scaling above modified for $\unicode[STIX]{x1D705}\gg 1$. The large-$\unicode[STIX]{x1D705}$ scaling of the drift in this regime will be analysed in detail in a future effort.

Figure 11.

A schematic highlighting the physical mechanism underlying the low-shear trapping regime. The solid-headed and hollow-headed swimmers represent the excess and deficit of swimmer (orientation) probability, respectively.

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 centreline at large $Pe$, beginning from the orientation dynamics of swimmers in a local linear flow description, as shown in figure 11(a). 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 centreline. For any finite $\unicode[STIX]{x1D705}$, and sufficiently large $Pe$ ($Pe\gg 1$ for $\unicode[STIX]{x1D705}$ of $O(1)$; $Pe\gg \unicode[STIX]{x1D705}^{3}$ for $\unicode[STIX]{x1D705}\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 (Reference Leal and Hinch1971), the distribution across the orbits is still controlled by asymptotically weak rotary diffusion, although this does not 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 11(a). 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 11(b) 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 sizes 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. These sources and sinks balance the divergence of the flux of the probability density; this flux is along Jeffery orbits for large $Pe$, and is the Jeffery angular velocity times the probability density. The swimming-induced source causes the flux to increase through quadrants II and III, attaining a maximum along the positive gradient axis; the sink in quadrants I and IV causes a decrease thereafter, with the flux attaining a minimum along the negative gradient axis. The flux being proportional to the probability density, leads to an excess of orientations in quadrants I and II relative to quadrants III and IV. This excess probability implies a net drift into the low-shear regions. The migration towards the low-shear regions continues till the swimmers reach sufficiently close to the channel centre where $Pe_{z}$ eventually falls below $\unicode[STIX]{x1D705}^{3}$. From this location onward, and up until the centreline, 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 centreline and the wall, with this location shifting towards the centreline with increasing $Pe$.

3.3 Depletion index and the $Pe{-}\unicode[STIX]{x1D705}$ migration portrait

Figure 12.

Schematic illustrating the depletion index ($DI$) calculation using steady state swimmer concentration profiles for $\unicode[STIX]{x1D705}=15$, for three different values of $Pe$ (100, 400 and 1000). (a) Area ($A_{1}$) under the initial uniform distribution (black colour dash-dotted line) from the wall ($z=-1$) up to its intersection with the steady state swimmer concentration profile ($z_{in}$). (b) Area ($A_{2}$) under the steady state profile from $z=-1$ to $z_{in}$. (c) The differential area $A_{1}-A_{2}$. The symbols $A_{1}$, $A_{2}$ (the respective shaded areas) and $z_{in}$ corresponding to $Pe=100,400$ and 1000, are represented in blue, red and black colour, respectively. The depletion index is expressed as: $DI=[2(A_{1}-A_{2})/A]\times 100$, where $A$ is the total area under any (initial uniform distribution or steady state concentration) profile and the factor of 2 accounts for the symmetry about the centreline.

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. Figure 12 defines the $DI$ using swimmer concentration profiles for $\unicode[STIX]{x1D705}=15$ (in a half-channel width) for three different $Pe$. $DI<0$ (${>}$ 0) corresponds to high-shear (low-shear) trapping. Next, in § 3.3.2, we organize the migration behaviour on a $Pe{-}\unicode[STIX]{x1D705}$ plane (for small $\unicode[STIX]{x1D716}$), thereby delineating the regimes of high-shear and low-shear trapping across a numerically determined phase boundary.

3.3.1 The depletion index

Figure 13.

Plot of $\unicode[STIX]{x0394}DI$ ($100-DI$) against $1/Pe$ for swimmers of different aspect ratios. The subscript ‘$cc$’ in the aspect ratio ($\unicode[STIX]{x1D705}_{cc}$) denotes ‘centreline collapse’. The inset represents the comparison of $DI_{Exp}$ (defined separately; see text), as predicted by the analysis here, and the experimental results of Rusconi et al. (Reference Rusconi, Guasto and Stocker2014b).

To analyse swimmer migration at large $Pe$, $\unicode[STIX]{x0394}DI=100-DI$ is plotted against $1/Pe$ for swimmers of different aspect ratios in figure 13. The plot discriminates between high ($\unicode[STIX]{x0394}DI>100$) and low-shear trapping behaviour ($\unicode[STIX]{x0394}DI<100$), and, in addition, helps differentiate between centreline collapse ($\unicode[STIX]{x0394}DI\rightarrow 0$ for $1/Pe\rightarrow 0$ for $\unicode[STIX]{x1D705}\geqslant 2$) and non-collapse ($\unicode[STIX]{x0394}DI$ approaches a finite value as for $1/Pe\rightarrow 0$ for $\unicode[STIX]{x1D705}<2$) behaviour within the low-shear trapping regime. Figure 13 clearly illustrates the singular nature of infinitely slender swimmers. For $\unicode[STIX]{x1D705}=\infty$, $\unicode[STIX]{x0394}DI$ plateaus at a finite value (${>}100$) for $1/Pe\rightarrow 0$. In contrast, for any finite $\unicode[STIX]{x1D705}\geqslant 2,\unicode[STIX]{x0394}DI$ approaches zero in the same limit, and the rate of approach, rather counter-intuitively, increases with increasing $\unicode[STIX]{x1D705}$. For large $\unicode[STIX]{x1D705}$ values, this accelerated approach towards the origin begins after an intermediate extended plateau with $\unicode[STIX]{x0394}DI>100$ (the same plateau value as for $\unicode[STIX]{x1D705}=\infty$). A kink in the $\unicode[STIX]{x0394}DI$-curves results as they cross the homogeneous state given by $DI=0$ (for instance, see the $\unicode[STIX]{x0394}DI$-curve for $\unicode[STIX]{x1D705}=15$ in figure 13). This kink is an artefact of our definition of $DI$ – although suited for characterizing centreline collapse, the point $z_{in}$ used in the definition, to calculate $A_{1}$ and $A_{2}$ above changes discontinuously as the wall concentration crosses the homogeneous value (see the locations of $z_{in}$ for $Pe=100$ and 400 in figure 12c). The latter happens during the transition from a near-centre depletion to a near-centre excess profile with increasing $Pe$, and, as a result, the profiles on the homogeneous line ($DI=0$) in figure 13 are near-centre excess profiles with wall concentration equal to the homogeneous value. This kink does not alter the general character of the $DI$-curves, however. Also note that Rusconi et al. (Reference Rusconi, Guasto and Stocker2014b), likely motivated by experimental considerations, had defined a 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 13 is again unchanged by this altered definition.

It is worth reiterating the behaviour of the depletion index, $\unicode[STIX]{x0394}DI$, as a function of the swimmer aspect ratio, starting from infinitely slender swimmers. Regardless of $\unicode[STIX]{x1D705}$, $\unicode[STIX]{x0394}DI$ starts from 100 for $Pe\ll 1$. For $\unicode[STIX]{x1D705}=\infty$, $\unicode[STIX]{x0394}DI$ increases with increasing $Pe$, approaching a constant ($\unicode[STIX]{x0394}DI=111.1$) for $Pe\rightarrow \infty$. For large but finite $\unicode[STIX]{x1D705}$, $\unicode[STIX]{x0394}DI$ again increases with increasing $Pe$, asymptoting to the aforementioned plateau value in the range $1\ll Pe\ll \unicode[STIX]{x1D705}^{3}$, and thereafter, decreases rapidly with further increase in $Pe$. This behaviour continues with decreasing $\unicode[STIX]{x1D705}$, 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 $\unicode[STIX]{x1D705}$. For the ABPs examined here, this trend continues until a critical aspect ratio of approximately 2, when $\unicode[STIX]{x0394}DI$ approaches zero only at a logarithmic rate for large $Pe$; that is, for $\unicode[STIX]{x1D705}\approx 2$, $\unicode[STIX]{x0394}DI$ is $O(1/\ln Pe)$ for $Pe\gg 1$. For $\unicode[STIX]{x1D705}<2$, $\unicode[STIX]{x0394}DI$ approaches a finite positive plateau (less than 100) for $Pe\rightarrow \infty$, with this plateau value approaching 100 as $\unicode[STIX]{x1D705}$ approaches unity (spherical swimmers). It is worth noting that the infinite-aspect-ratio assumption is used routinely in the rheological context, when analysing swimmer suspensions, and generally yields good agreement between experiments and theory (Saintillan Reference Saintillan2010a,Reference Saintillanb; López et al. Reference López, Gachelin, Douarche, Auradou and Clément2015; Nambiar et al. Reference Nambiar, Nott and Subramanian2017, Reference Nambiar, Phanikanth, Nott and Subramanian2019b). In the present context, however, even for moderately large aspect ratios, $\unicode[STIX]{x1D705}\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. (Reference Rusconi, Guasto and Stocker2014b); 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 13 for $\unicode[STIX]{x1D705}=10$. While the disagreement appears to be 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 Reference Sokolov and Aranson2016; Sokolov et al. Reference Sokolov, Rubio, Brady and Aranson2018) 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 centreline collapse as $Pe\rightarrow \infty$).

3.3.2 The $Pe{-}\unicode[STIX]{x1D705}$ migration portrait

Figure 14 organizes the high-shear and low-shear trapping regimes on the $Pe{-}\unicode[STIX]{x1D705}$ 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 \unicode[STIX]{x1D705}^{3}$; the pre-factor 0.07 given in figure 14 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 $\unicode[STIX]{x1D705}\rightarrow 1$ can transform from a near-centre depletion to a near-centre 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 $\unicode[STIX]{x1D6FA}=F/4\unicode[STIX]{x03C0}+B\unicode[STIX]{x1D6FA}_{1}+\cdots \,$, one finds $Pe_{c}\sim 4.5$ (see appendix C). Therefore, as $\unicode[STIX]{x1D705}\rightarrow 1$, ABPs 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 centreline, already seen in the numerical profiles. This transition for near-spherical swimmers has again been validated using numerical predictions based on the multiple scale analysis, and Langevin simulations.

Figure 14.

Shear-induced migration behaviour on the $Pe{-}\unicode[STIX]{x1D705}$ plane. The red and blue circles, respectively, represent the high-shear trapping and low-shear trapping regimes (based on the profiles obtained from the multiple scale analysis, with partial verification from Langevin simulations), and the grey dashed line is an approximate boundary separating the two. The solid black lines denote the near-sphere and large-aspect-ratio asymptotes for this boundary. The diamonds and triangles respectively, corresponding to the different trapping regimes of the algae Heterosigma and Dunaliella observed by Barry et al. (Reference Barry, Rusconi, Guasto and Stocker2015).

Also shown in figure 14 are the results from the experiments of Barry et al. (Reference Barry, Rusconi, Guasto and Stocker2015) corresponding to both high-shear and low-shear trapping of the phytoplankton algae Heterosigma and Dunaliella. The effective swimmer aspect ratios ($\unicode[STIX]{x1D705}$) in the experiments were determined from a best fit of the theoretical orientation distribution, obtained by solving the Fokker–Planck equation with $\unicode[STIX]{x1D705}$ as a fitting parameter, with the experimentally determined one. Note that the $\unicode[STIX]{x1D705}$ values 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. (Reference Barry, Rusconi, Guasto and Stocker2015), for $Pe\leqslant 5$, Heterosigma exhibits high-shear trapping. However, Dunaliella, exhibits high-shear trapping until $Pe\approx 5$, and transitions to low-shear trapping at $Pe\approx 17$. As seen in figure 14, our predictions are consistent with these observations; in particular, we too predict a transition from high-shear to low-shear trapping, with increasing $Pe$, for Dunaliella. An important point to note here is that the $\unicode[STIX]{x1D716}$ values for the algal microswimmers above are close to unity (for Heterosigma, $\unicode[STIX]{x1D716}\sim 1.12$, and for Dunaliella ${\sim}0.717$), calling into question the validity of the method of multiple scale analysis. However, as mentioned earlier, the domain of validity of the analysis is controlled by $\unicode[STIX]{x1D716}_{eff}$, and not $\unicode[STIX]{x1D716}$. For large $Pe$ values, $\unicode[STIX]{x1D716}_{eff}\sim \unicode[STIX]{x1D716}Pe^{-1}$, and remains small even when $\unicode[STIX]{x1D716}$ is of order unity. This then allows for the trapping regimes of the aforementioned algae to be identified using the multiple time scale analysis. The consistency between the experiments of Barry et al. (Reference Barry, Rusconi, Guasto and Stocker2015) and our theory serves as an important validation of the analysis and arguments presented here.

3.4 Migration pattern of swimmers with passive translational diffusivity

The analysis described thus far shows 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 Reference Leal and Hinch1971; Hinch & Leal Reference Hinch and Leal1972a,Reference Hinch and Lealb; Leal & Hinch Reference Leal and Hinch1972, Reference Leal and Hinch1973; Hinch & Leal Reference Hinch and Leal1973), 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, the orientation dynamics turn out to be much more crucial for active particles, being responsible for the onset of a low-shear trapping regime! Indeed, Nitsche & Hinch (Reference Nitsche and Hinch1997) have previously used the method of multiple scales to examine shear-induced migration in a dilute suspension of passive anisotropic particles, and only found a high-shear trapping regime. The latter arises on account of the anisotropy of the mobility, and thence, the Brownian diffusivity; the lower diffusivity (mobility) in the transverse direction, of flow-aligned particles close to the walls, leads to accumulation. This 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 fibre suspensions have failed to observe the predicted high-shear trapping (Nitsche & Hinch Reference Nitsche and Hinch1997) at lower volume fractions (Strednak et al. Reference Strednak, Shaikh, Butler and Guazzelli2018). Importantly, since the 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 (for sufficiently large $Pe$, $\unicode[STIX]{x0394}DI$ decays as $O(Pe^{-1})$ for any finite $\unicode[STIX]{x1D705}$, and as $O(Pe^{-1/3})$ for $\unicode[STIX]{x1D705}=\infty$). This is, of course, in sharp contrast to the active case examined here.

In light of the contrasting nature of the active and passive migration behaviour, it is worth examining (briefly) the effect of an independent source of translational diffusion (thermal or otherwise) on shear-induced migration. This corresponds to an additional term of the form (Kim & Karrila Reference Kim and Karrila2005)

(3.1)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x1D701}D_{r}}{U_{s}^{2}}k_{B}T\left[M_{\Vert }\boldsymbol{p}\boldsymbol{p}+M_{\bot }(\boldsymbol{I}-\boldsymbol{p}\boldsymbol{p})\right]\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D701}\bar{D}\left[\frac{\left(M_{\Vert }^{\prime }+2M_{\bot }^{\prime }\right)}{3}\boldsymbol{I}+\left(M_{\Vert }^{\prime }-M_{\bot }^{\prime }\right)\left(\,\boldsymbol{p}\boldsymbol{p}-\frac{\boldsymbol{I}}{3}\right)\right],\end{eqnarray}$$

in the dimensionless equation governing the phase-space probability density (equation (2.2) in § 2); here, $M_{\Vert }$ and $M_{\bot }$ denote the $\unicode[STIX]{x1D705}$-dependent translational mobilities along and transverse to the swimmer axis, with $M_{\Vert }^{\prime }=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}LM_{\Vert }$, $M_{\bot }^{\prime }=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}LM_{\bot }$; the averaged mobility, $(M_{\Vert }+2M_{\bot })/3$ corresponds to the passive translational diffusivity for an isotropic orientation distribution ($Pe=0$). The parameter $\bar{D}=(k_{B}TD_{r}/6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}LU_{s}^{2})$ represents the ratio of the active to the passive diffusivities, with $k_{B}$ being the Boltzmann constant, $T$ the temperature and $\unicode[STIX]{x1D702}$ the suspension viscosity (Kim & Karrila Reference Kim and Karrila2005). The additional parameter $\unicode[STIX]{x1D701}$ accounts for the possible athermal origin of the translational diffusion. Thus, $\unicode[STIX]{x1D701}\geqslant 1$, with $\unicode[STIX]{x1D701}=1$ corresponding to the translational diffusion having a thermal origin. Application of the method of multiple scale analysis again leads to the drift-diffusion equation (2.7), but with additional passive contributions to both the diffusivity and drift given by

(3.2)$$\begin{eqnarray}\displaystyle & \displaystyle D_{zz}^{p}=\unicode[STIX]{x1D701}\bar{D}\int \text{d}\,\boldsymbol{p}\left[\frac{(M_{\Vert }^{\prime }+2M_{\bot }^{\prime })}{3}+\frac{(M_{\Vert }^{\prime }-M_{\bot }^{\prime })}{3}(3\cos ^{2}\unicode[STIX]{x1D703}-1)\right]G(\,\boldsymbol{p}), & \displaystyle\end{eqnarray}$$

(3.3)$$\begin{eqnarray}\displaystyle & \displaystyle V_{z}^{p}=\unicode[STIX]{x1D701}\bar{D}\int \text{d}\,\boldsymbol{p}\left[\frac{(M_{\Vert }^{\prime }+2M_{\bot }^{\prime })}{3}+\frac{(M_{\Vert }^{\prime }-M_{\bot }^{\prime })}{3}(3\cos ^{2}\unicode[STIX]{x1D703}-1)\right]\unicode[STIX]{x2202}G(\,\boldsymbol{p})/\unicode[STIX]{x2202}z. & \displaystyle\end{eqnarray}$$

Choosing $D_{r}=1~\text{s}^{-1}$, $L\approx 9~\unicode[STIX]{x03BC}\text{m}$, $U_{s}\approx 50~\unicode[STIX]{x03BC}\text{m}$, $T\approx 300~\text{K}$ and $\unicode[STIX]{x1D702}\approx 10^{-3}~\text{Ns}~\text{m}^{-2}$, which are typical values for microorganisms in an aqueous medium, yields $\bar{D}\approx 9.765\times 10^{-6}$ (Rusconi et al. Reference Rusconi, Guasto and Stocker2014b; López et al. Reference López, Gachelin, Douarche, Auradou and Clément2015). In figure 15, we compare the variation of $\unicode[STIX]{x0394}DI$ against $1/Pe$, for the $\unicode[STIX]{x1D705}$ values shown in figure 13 (which appears again as figure 15(a) with $\unicode[STIX]{x1D701}=0$), with the cases $\unicode[STIX]{x1D701}=1$ and 10. While the swimmer migration remains similar to the active case (figure 15a) for low to moderate $Pe$, the passive translational diffusivity given by (3.2) becomes dominant at large $Pe$ (since the active diffusivity approaches zero for $Pe\rightarrow \infty$), causing $\unicode[STIX]{x0394}DI$ to approach 100 (the homogeneous state) in the limit $Pe\rightarrow \infty$. While this trend, on the whole, is similar to that reported by Ezhilan & Saintillan (Reference Ezhilan and Saintillan2015), the nature of the approach to the homogeneous state is a function of $\unicode[STIX]{x1D705}$. The approach to homogeneity may be from either the high-shear or low-shear trapping regime depending on $\unicode[STIX]{x1D705}$. With increasing $\unicode[STIX]{x1D701}$, the approach from the high-shear trapped regime extends to a progressively larger range of aspect ratios, with a concomitant collapse of the depletion index curves towards the homogeneous state.

Figure 15.

Variation of $\unicode[STIX]{x0394}DI~(100-DI)$ against $1/Pe$ for swimmers of different aspect ratios. The parameter $\unicode[STIX]{x1D701}$ accounts for the athermal origin of the translational diffusivity (a) $0$, (b) $1$ and (c) $10$; here, $\unicode[STIX]{x1D701}=0$ represents the absence of passive translational diffusivity ($\unicode[STIX]{x0394}DI$ of figure 13) and $\unicode[STIX]{x1D701}=1$ corresponds to the translational diffusivity having a thermal origin. The black dash-dotted line corresponds to the homogeneous state ($DI=0$). The inset in (c) represents the magnified view of $\unicode[STIX]{x0394}DI$ versus $1/Pe$, highlighting the approach towards the homogeneous state at large $Pe$.

Figure 16.

Variation of $\unicode[STIX]{x0394}DI~(100-DI)$ against $1/Pe$ for swimmers of different aspect ratios, based on the mean passive translational diffusivity as $Pe\rightarrow \infty$.

Figure 15(a–c) shows that the inclusion of passive translational diffusion can have a non-trivial effect on activity-induced shear migration. In particular, a large enough $\unicode[STIX]{x1D701}$ value can impede the appearance of the low-shear trapping regime. To illustrate this, we consider the choice made in Ezhilan & Saintillan (Reference Ezhilan and Saintillan2015); the authors used a dimensionless parameter, $\unicode[STIX]{x1D6EC}=D_{t}D_{r}/U_{s}^{2}$, to incorporate the effects of a scalar translational diffusivity. Most of their results correspond to $\unicode[STIX]{x1D6EC}=1/6$, which makes the translational diffusivity equal to the swimming-induced diffusivity ($U_{s}^{2}/6D_{r}$) in the absence of a flow. The translational diffusion in the present analysis is a tensor owing to which there no direct analogue to the parameter $\unicode[STIX]{x1D6EC}$. Nevertheless, one may choose a suitably averaged translational diffusion to replace $D_{t}$ in $\unicode[STIX]{x1D6EC}$ above. Since the effects of translational diffusion are most important at large $Pe$, such an average ($\langle D_{t}\rangle$) may be based on the infinite$-Pe$ (but $\unicode[STIX]{x1D705}$-dependent) orientation distribution derived by Leal & Hinch (Reference Leal and Hinch1971). Thus, replacing $D_{t}$ by $\langle D_{t}\rangle$ in $\unicode[STIX]{x1D6EC}$ above, and choosing $\unicode[STIX]{x1D6EC}=1/6$, leads to the depletion index curves shown in figure 16; note that the fixed value of $\unicode[STIX]{x1D6EC}$ implies that $\unicode[STIX]{x1D701}$ is now a function of $\unicode[STIX]{x1D705}$, and thus, each of the depletion curves in figure 16 corresponds to a different $\unicode[STIX]{x1D701}$. Although not shown, in contrast to the near-centre depletion profiles seen in figure 2, the large-$Pe$ concentration profiles along these curves begin to approach uniformity in the bulk, with a narrow depletion layer around the centreline. The choice of the orientational average in $\langle D_{t}\rangle$ does not alter the above conclusion; we have verified that use of the isotropic diffusivity $(D_{\Vert }+2D_{\bot })/3$, appropriate for $Pe=0$, leads to depletion index curves similar to those shown in figure 16 (the isotropic average and the $Pe=\infty$ Leal–Hinch average serve to bound the arbitrary-$Pe$ translational diffusivity).