2.1. Problem definition and kinetic model

We analyse the dynamics in a dilute suspension of self-propelled slender particles confined between two parallel flat plates and placed in an externally imposed pressure-driven flow as illustrated in figure 1. The channel half-width is denoted by $H$ , and is assumed to be much greater than the characteristic length $L$ of the particles ( $H/L\gg 1$ ), so that the finite size of the particles can be neglected. The external flow follows the parabolic Poiseuille profile

(2.1) $$\begin{eqnarray}\boldsymbol{U}(\boldsymbol{x})=U(z)\,\hat{\boldsymbol{y}}=U_{m}[1-(z/H)^{2}]\hat{\boldsymbol{y}},\end{eqnarray}$$

with maximum velocity $U_{m}$ at the centreline ( $z=0$ ). The shear rate varies linearly with position $z$ across the channel:

(2.2) $$\begin{eqnarray}S(z)=\frac{\text{d}U}{\text{d}z}=-\dot{{\it\gamma}}_{w}\frac{z}{H},\end{eqnarray}$$

where $\dot{{\it\gamma}}_{w}=2U_{m}/H$ is the maximum absolute shear rate attained at the walls ( $z=\pm H$ ).

Figure 1. Problem definition: a dilute suspension of slender active particles with positions $\boldsymbol{x}=(x,y,z)$ and orientations $\boldsymbol{p}=(\sin {\it\theta}\cos {\it\phi},\sin {\it\theta}\sin {\it\phi},\cos {\it\theta})$ is confined between two parallel flat plates ( $z=\pm H$ ) and subject to an imposed pressure-driven parabolic flow.

Following previous models for active suspensions (Saintillan & Shelley Reference Saintillan and Shelley2008a ,Reference Saintillan and Shelley b ), the configuration of the active particles is captured by the probability distribution function ${\it\Psi}(\boldsymbol{x},\boldsymbol{p},t)$ of finding a particle at position $\boldsymbol{x}=(x,y,z)$ with orientation $\boldsymbol{p}=(\sin {\it\theta}\cos {\it\phi},\sin {\it\theta}\sin {\it\phi},\cos {\it\theta})$ at time $t$ , where $\boldsymbol{p}$ also defines the direction of swimming. Conservation of particles is expressed by the Smoluchowski equation (Doi & Edwards Reference Doi and Edwards1986)

(2.3) $$\begin{eqnarray}\frac{\partial {\it\Psi}}{\partial t}+\boldsymbol{{\rm\nabla}}_{x}\boldsymbol{\cdot }(\dot{\boldsymbol{x}}\,{\it\Psi})+\boldsymbol{{\rm\nabla}}_{p}\boldsymbol{\cdot }(\dot{\boldsymbol{p}}\,{\it\Psi})=0,\end{eqnarray}$$

where the translational flux velocity $\dot{\boldsymbol{x}}$ captures self-propulsion with constant velocity $V_{s}$ in the direction of $\boldsymbol{p}$ , advection by the imposed flow and centre-of-mass diffusion with isotropic and constant diffusivity $d_{t}$ :

(2.4) $$\begin{eqnarray}\dot{\boldsymbol{x}}=V_{s}\,\boldsymbol{p}+\boldsymbol{U}(z)-d_{t}\boldsymbol{{\rm\nabla}}_{x}\ln {\it\Psi}.\end{eqnarray}$$

Particle rotations are captured by the angular flux velocity $\dot{\boldsymbol{p}}$ , which includes contributions from the imposed flow via Jeffery’s equation (Jeffery Reference Jeffery1922; Bretherton Reference Bretherton1962), and from rotational diffusion with diffusivity  $d_{r}$ :

(2.5) $$\begin{eqnarray}\dot{\boldsymbol{p}}=S(z)(\hat{\boldsymbol{z}}\boldsymbol{\cdot }\boldsymbol{p})(\unicode[STIX]{x1D644}-\boldsymbol{p}\boldsymbol{p})\boldsymbol{\cdot }\hat{\boldsymbol{y}}-d_{r}\boldsymbol{{\rm\nabla}}_{p}\ln {\it\Psi}.\end{eqnarray}$$

We have assumed that the particles have a high aspect ratio ( $r\rightarrow \infty$ ), leading to a Bretherton constant of ${\it\zeta}=(r^{2}-1)/(r^{2}+1)\approx 1$ , a good approximation for common motile bacteria as well as many self-propelled catalytic micro-rods. Particle–particle hydrodynamic interactions have also been neglected based on the assumption of infinite dilution; such interactions could otherwise be included via an additional disturbance velocity in the expressions for $\dot{\boldsymbol{x}}$ and $\dot{\boldsymbol{p}}$ (Saintillan & Shelley Reference Saintillan and Shelley2008b ). As a result, we expect the distribution of particles to be uniform along the $x$ and $y$ directions, and at steady state the Smoluchowski equation (2.3) for ${\it\Psi}(\boldsymbol{x},\boldsymbol{p},t)={\it\Psi}(z,\boldsymbol{p})$ then simplifies to

(2.6) $$\begin{eqnarray}V_{s}\cos {\it\theta}\frac{\partial {\it\Psi}}{\partial z}-d_{t}\frac{\partial ^{2}{\it\Psi}}{\partial z^{2}}+S(z)\boldsymbol{{\rm\nabla}}_{p}\boldsymbol{\cdot }[\cos {\it\theta}(\unicode[STIX]{x1D644}-\boldsymbol{p}\boldsymbol{p})\boldsymbol{\cdot }\hat{\boldsymbol{y}}\,{\it\Psi}]=d_{r}{\rm\nabla}_{p}^{2}{\it\Psi}.\end{eqnarray}$$

This equation simply expresses the balance of self-propulsion, translational diffusion, particle alignment by the imposed flow and rotational diffusion.

In this work, we treat the translational and rotational diffusivities $d_{t}$ and $d_{r}$ as independent constants, which could result from either Brownian motion or various athermal sources of noise (Drescher et al. Reference Drescher, Dunkel, Cisneros, Ganguly and Goldstein2011; Garcia et al. Reference Garcia, Berti, Peyla and Rafaï2011). The athermal contribution to diffusion may arise due to tumbling or other fluctuations in the swimming actuation of motile micro-organisms, or from fluctuations in the chemical actuation mechanism of catalytic particles. In many active suspensions, such athermal fluctuations are in fact the dominant source of diffusion.

2.2. Boundary conditions

In the continuum limit, the impenetrability of the channel walls is captured by prescribing that the normal component of the translational flux be zero at both walls:

(2.7) $$\begin{eqnarray}\hat{\boldsymbol{z}}\boldsymbol{\cdot }\dot{\boldsymbol{x}}=0\quad \text{at}~z=\pm H.\end{eqnarray}$$

Inserting (2.4) for the translational flux, this leads to a Robin boundary condition for the probability distribution function,

(2.8) $$\begin{eqnarray}d_{t}\frac{\partial {\it\Psi}}{\partial z}=V_{s}\cos {\it\theta}\,{\it\Psi}\quad \text{at}~z=\pm H,\end{eqnarray}$$

expressing the balance of translational diffusion and self-propulsion in the wall-normal direction. Equation (2.8) implies that particles pointing towards a wall ( $\cos {\it\theta}>0$ for the top wall at $z=+H$ ) incur a positive wall-normal gradient ( $\partial {\it\Psi}/\partial z>0$ ), whereas particles pointing away from the wall ( $\cos {\it\theta}<0$ ) incur a negative gradient. This suggests that sorting of orientations should occur and lead to a net polarization towards the walls, accompanied by near-wall accumulation. These effects will indeed be confirmed in § 3. It is important to note that the boundary condition (2.8) requires that the wall-normal swimming flux be balanced by a diffusive flux. In the complete absence of translational diffusion ( $d_{t}=0$ ), the swimming flux can no longer be balanced at the wall: this singular limit, which is ill-posed in our mean-field theory, will not be addressed here. Note also that the balance of the wall-normal fluxes hints at a length scale of $\ell _{a}=d_{t}/V_{s}$ for wall accumulation, as we demonstrate more quantitatively below.

Other types of boundary conditions have been considered in previous works. In particular, several studies have implemented the condition

(2.9) $$\begin{eqnarray}\int _{{\it\Omega}}\hat{\boldsymbol{z}}\boldsymbol{\cdot }\dot{\boldsymbol{x}}\,{\it\Psi}\,\text{d}\boldsymbol{p}=0\quad \text{at}~z=\pm H,\end{eqnarray}$$

where ${\it\Omega}$ denotes the unit sphere of orientation. Equation (2.9) captures the zeroth orientational moment of (2.8) and expresses conservation of active particle concentration. However, this is a weaker condition than (2.8), which can be enforced in a variety of ways and is therefore insufficient as a boundary condition. It is often implemented numerically using a reflection condition on the distribution function. This reflection condition was first used by Bearon, Hazel & Thorn (Reference Bearon, Hazel and Thorn2011) in a two-dimensional model of suspensions of gyrotactic swimmers constrained to a planar domain. Ezhilan, Pahlavan & Saintillan (Reference Ezhilan, Pahlavan and Saintillan2012) also imposed equation (2.9) using the reflection condition for the case of a chemotactic active suspension confined to a thin liquid film, where the primary mechanism for accumulation was chemotaxis as opposed to kinematics. In the absence of external fields, however, the reflection boundary condition allows for a uniform isotropic solution throughout the channel and is therefore unable to capture near-wall accumulation or upstream swimming when a flow is imposed (see appendix A for more details). Kasyap & Koch (Reference Kasyap and Koch2014) also considered chemotactic active suspensions in thin films but used a position/orientation decoupling approximation for the probability distribution function ${\it\Psi}(\boldsymbol{x},\boldsymbol{p},t)$ , allowing them to derive a boundary condition for the number density field expressing the balance of the chemotactic and diffusive fluxes at the boundaries. To our knowledge, the only previously reported use of the boundary condition (2.8) for a confined active suspension was in the work of Elgeti & Gompper (Reference Elgeti and Gompper2013), whose analysis was restricted to equilibrium distributions in the absence of flow and in the limits of narrow channels or weak propulsion.

Finally, it should be kept in mind that the simple boundary condition (2.8) neglects the finite size of the particles and is therefore inaccurate very close to the walls, where steric exclusion prohibits certain particle configurations and should lead to a depletion layer as observed in experiments (Takagi et al. Reference Takagi, Palacci, Braunschweig, Shelley and Zhang2014). The implications of steric exclusion are discussed further in appendix B, where a more detailed boundary condition is derived and enforced on the hypersurface separating allowed from forbidden configurations (Nitsche & Brenner Reference Nitsche and Brenner1990; Schiek & Shaqfeh Reference Schiek and Shaqfeh1995; Krochak, Olson & Martinez Reference Krochak, Olson and Martinez2010). As we show there, the effects of steric exclusion are weak in wide channels ( $H/L\gg 1$ ) such as the ones considered in this work.

2.3. Dimensional analysis and scaling

Dimensional analysis of the governing equations reveals three dimensionless groups:

(2.10a−c ) $$\begin{eqnarray}Pe_{s}=\frac{V_{s}}{2d_{r}H},\quad Pe_{f}=\frac{\dot{{\it\gamma}}_{w}}{d_{r}},\quad {\it\Lambda}=\frac{d_{t}d_{r}}{V_{s}^{2}}.\end{eqnarray}$$

The first parameter $Pe_{s}$ , or swimming Péclet number, can be interpreted as the ratio of the characteristic time scale for a particle to lose memory of its orientation due to rotational diffusion over the time it takes it to swim across the channel width. Equivalently, it is also the ratio of the persistence length of particle trajectories ( $\ell _{p}=V_{s}/d_{r}$ ) over the channel width ( $2H$ ). The second parameter $Pe_{f}$ , or flow Péclet number, compares the same diffusive time scale to the characteristic time for a particle to align under the imposed velocity gradient. The third parameter ${\it\Lambda}$ relates the translational and rotational diffusivities to the swimming speed and is a fixed constant for a given particle type. It can be interpreted as an inverse measure of the strength of propulsion of a swimmer with respect to fluctuations, and the limits of ${\it\Lambda}\rightarrow 0$ and ${\it\Lambda}\rightarrow \infty$ describe the strong and weak propulsion cases, respectively. When ${\it\Lambda}$ is held constant, $Pe_{s}$ also reduces to a measure of the degree of confinement, with $Pe_{s}\rightarrow 0$ and $Pe_{s}\rightarrow \infty$ describing the limits of weak and strong confinement, respectively.

In the following, we non-dimensionalize the governing equations using the characteristic time, length and velocity scales

(2.11a−c ) $$\begin{eqnarray}t_{c}=d_{r}^{-1},\quad \ell _{c}=H,\quad v_{c}=Hd_{r},\end{eqnarray}$$

and also normalize the distribution function ${\it\Psi}$ by the mean number density $n$ defined as

(2.12) $$\begin{eqnarray}n=\frac{1}{2H}\int _{-H}^{H}\int _{{\it\Omega}}{\it\Psi}(z,\boldsymbol{p})\,\text{d}\boldsymbol{p}\,\text{d}z.\end{eqnarray}$$

After non-dimensionalization, the conservation equation (2.6) becomes

(2.13) $$\begin{eqnarray}Pe_{s}\cos {\it\theta}\frac{\partial {\it\Psi}}{\partial z}-2{\it\Lambda}Pe_{s}^{2}\frac{\partial ^{2}{\it\Psi}}{\partial z^{2}}+\frac{Pe_{f}}{2}\,S(z)\boldsymbol{{\rm\nabla}}_{p}\boldsymbol{\cdot }[\cos {\it\theta}(\unicode[STIX]{x1D644}-\boldsymbol{p}\boldsymbol{p})\boldsymbol{\cdot }\hat{\boldsymbol{y}}\,{\it\Psi}]=\frac{1}{2}{\rm\nabla}_{p}^{2}{\it\Psi},\end{eqnarray}$$

where the dimensionless shear rate profile is simply $S(z)=-z$ . The boundary condition (2.8) also becomes

(2.14) $$\begin{eqnarray}\frac{\partial {\it\Psi}}{\partial z}=\frac{1}{2{\it\Lambda}Pe_{s}}\cos {\it\theta}\,{\it\Psi}\quad \text{at}~z=\pm 1.\end{eqnarray}$$

Note that the choice of $H$ for the characteristic length scale is convenient, as it sets the positions of the boundaries to $z=\pm 1$ in the dimensionless system. However, we will see below that alternative length scales are more judiciously chosen in certain limits due to the presence of boundary layers.

2.4. Orientational moment equations

Equation (2.13), together with boundary condition (3.23), cannot be solved analytically in general. While a numerical solution is possible, as we show below, analytical progress can still be made in terms of orientational moments of the distribution function (Saintillan & Shelley Reference Saintillan and Shelley2013). More precisely, we introduce the zeroth, first and second moments of ${\it\Psi}(z,\boldsymbol{p})$ as

(2.15a−c ) $$\begin{eqnarray}c(z)=\langle 1\rangle ,\quad \boldsymbol{m}(z)=\langle \boldsymbol{p}\rangle ,\quad \unicode[STIX]{x1D63F}(z)=\langle \boldsymbol{p}\boldsymbol{p}-\unicode[STIX]{x1D644}/3\rangle ,\end{eqnarray}$$

where the angular brackets $\langle \cdot \rangle$ denote the orientational average

(2.16) $$\begin{eqnarray}\langle h(\boldsymbol{p})\rangle =\int _{{\it\Omega}}h(\boldsymbol{p}){\it\Psi}(z,\boldsymbol{p})\,\text{d}\boldsymbol{p}.\end{eqnarray}$$

The zeroth moment $c(z)$ corresponds to the local concentration of particles. The next two moments are directly related to the polarization vector $\boldsymbol{P}(z)$ and to the nematic order parameter tensor $\unicode[STIX]{x1D64C}(z)$ commonly used in the description of liquid-crystalline systems (Marchetti et al. Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Aditi Simha2013) as

(2.17a,b ) $$\begin{eqnarray}\boldsymbol{m}(z)=c(z)\boldsymbol{P}(z),\quad \unicode[STIX]{x1D63F}(z)=c(z)\unicode[STIX]{x1D64C}(z).\end{eqnarray}$$

Knowledge of these as well as higher moments also allows one to recover the full distribution function as

(2.18) $$\begin{eqnarray}{\it\Psi}(z,\boldsymbol{p})=\frac{1}{4{\rm\pi}}\,c(z)+\frac{3}{4{\rm\pi}}\,\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{m}(z)+\frac{15}{8{\rm\pi}}\,\boldsymbol{p}\boldsymbol{p}\boldsymbol{ : }\unicode[STIX]{x1D63F}(z)+\cdots \,,\end{eqnarray}$$

which can also be interpreted as a spectral expansion of ${\it\Psi}(z,\boldsymbol{p})$ on the basis of spherical harmonics. Near isotropy this expansion converges rapidly, which justifies truncation after a few terms. If only the first three terms corresponding to $c$ , $\boldsymbol{m}$ and $\unicode[STIX]{x1D63F}$ are retained, a closed system of equations can be derived for these variables by taking moments of the conservation equation (2.13) (Baskaran & Marchetti Reference Baskaran and Marchetti2009; Saintillan & Shelley Reference Saintillan and Shelley2013).

In the problem of interest to us here, symmetries dictate that the only non-zero components of $\boldsymbol{m}$ and $\unicode[STIX]{x1D63F}$ are $m_{z}$ and $D_{zz}=-2D_{xx}=-2D_{yy}$ in the absence of flow. When a flow is applied in the $y$ direction, $m_{y}$ and $D_{yz}=D_{zy}$ are also expected to become non-zero, and $D_{yy}$ need no longer be equal to $D_{xx}$ . The governing equations for these variables can be obtained as

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle Pe_{s}\frac{\text{d}m_{z}}{\text{d}z}-2{\it\Lambda}Pe_{s}^{2}\frac{\text{d}^{2}c}{\text{d}z^{2}}=0, & \displaystyle\end{eqnarray}$$

(2.20) $$\begin{eqnarray}\displaystyle & \displaystyle Pe_{s}\frac{\text{d}D_{zz}}{\text{d}z}-2{\it\Lambda}Pe_{s}^{2}\frac{\text{d}^{2}m_{z}}{\text{d}z^{2}}+\left(\frac{1}{6{\it\Lambda}}+1\right)m_{z}=-\frac{1}{10}Pe_{f}S(z)m_{y}, & \displaystyle\end{eqnarray}$$

(2.21) $$\begin{eqnarray}\displaystyle & \displaystyle Pe_{s}\frac{\text{d}D_{yz}}{\text{d}z}-2{\it\Lambda}Pe_{s}^{2}\frac{\text{d}^{2}m_{y}}{\text{d}z^{2}}+m_{y}=\frac{2}{5}Pe_{f}S(z)m_{z}, & \displaystyle\end{eqnarray}$$

(2.22) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{4}{15}Pe_{s}\frac{\text{d}m_{z}}{\text{d}z}-2{\it\Lambda}Pe_{s}^{2}\frac{\text{d}^{2}D_{zz}}{\text{d}z^{2}}+3D_{zz}=\frac{4}{7}Pe_{f}S(z)D_{yz}, & \displaystyle\end{eqnarray}$$

(2.23) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{2}{15}Pe_{s}\frac{\text{d}m_{z}}{\text{d}z}-2{\it\Lambda}Pe_{s}^{2}\frac{\text{d}^{2}D_{yy}}{\text{d}z^{2}}+3D_{yy}=-\frac{3}{7}Pe_{f}S(z)D_{yz}, & \displaystyle\end{eqnarray}$$

(2.24) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{5}Pe_{s}\frac{\text{d}m_{y}}{\text{d}z}-2{\it\Lambda}Pe_{s}^{2}\frac{\text{d}^{2}D_{yz}}{\text{d}z^{2}}+3D_{yz}=Pe_{f}S(z)\left(\frac{1}{10}c+\frac{5}{14}D_{zz}-\frac{2}{7}D_{yy}\right). & \displaystyle\end{eqnarray}$$

$D_{xx}$

$D_{yy}$

$D_{zz}$

$\unicode[STIX]{x1D63F}$

$Pe_{f}=0$

Boundary conditions for these variables are also readily obtained by taking moments of (3.23), yielding

(2.25) $$\begin{eqnarray}\frac{\text{d}c}{\text{d}z}=\frac{1}{2{\it\Lambda}Pe_{s}}m_{z},\end{eqnarray}$$

(2.26a,b ) $$\begin{eqnarray}\frac{\text{d}m_{z}}{\text{d}z}=\frac{1}{2{\it\Lambda}Pe_{s}}\left(D_{zz}+\frac{1}{3}c\right),\quad \frac{\text{d}m_{y}}{\text{d}z}=\frac{1}{2{\it\Lambda}Pe_{s}}D_{yz},\end{eqnarray}$$

(2.27a−c ) $$\begin{eqnarray}\frac{\text{d}D_{zz}}{\text{d}z}=\frac{2}{15{\it\Lambda}Pe_{s}}m_{z},\quad \frac{\text{d}D_{yy}}{\text{d}z}=-\frac{1}{15{\it\Lambda}Pe_{s}}m_{z},\quad \frac{\text{d}D_{yz}}{\text{d}z}=\frac{1}{10{\it\Lambda}Pe_{s}}m_{y},\end{eqnarray}$$

all to be enforced at $z=\pm 1$ . For symmetry reasons, we expect $c$ , $m_{y}$ , $D_{yy}$ and $D_{zz}$ to be even functions of $z$ , whereas $m_{z}$ and $D_{yz}$ are expected to be odd functions.

Integrating equation (2.19) and making use of the boundary condition (2.25) easily shows that (2.19) can be replaced by

(2.28) $$\begin{eqnarray}m_{z}-2{\it\Lambda}Pe_{s}\frac{\text{d}c}{\text{d}z}=0\end{eqnarray}$$

at every point in the channel, underlining the direct relation between transverse polarization and concentration gradients. We also note that the normalization condition (2.12) on the distribution function translates into an integral condition on the concentration field expressing conservation of the total particle number:

(2.29) $$\begin{eqnarray}\int _{-1}^{1}c(z)\,\text{d}z=2.\end{eqnarray}$$

As we discuss next, solution of the system (2.19)–(2.24) subject to the boundary conditions (2.25)–(2.27) and to the integral constraint (2.29) is possible under certain assumptions, and provides results that are in excellent quantitative agreement with the full numerical solution of the Smoluchowski equation (2.13) over a wide range of values of the Péclet numbers.

3.1. Theory based on moment equations

In the absence of flow, the moment equations derived in § 2.4 only involve $c$ , $m_{z}$ and $D_{zz}$ , and simplify to

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle m_{z}-2{\it\Lambda}Pe_{s}\frac{\text{d}c}{\text{d}z}=0, & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle Pe_{s}\frac{\text{d}D_{zz}}{\text{d}z}-2{\it\Lambda}Pe_{s}^{2}\frac{\text{d}^{2}m_{z}}{\text{d}z^{2}}+\left(\frac{1}{6{\it\Lambda}}+1\right)m_{z}=0, & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{4}{15}Pe_{s}\frac{\text{d}m_{z}}{\text{d}z}-2{\it\Lambda}Pe_{s}^{2}\frac{\text{d}^{2}D_{zz}}{\text{d}z^{2}}+3D_{zz}=0, & \displaystyle\end{eqnarray}$$

(3.7a,b ) $$\begin{eqnarray}\frac{\text{d}m_{z}}{\text{d}z}=\frac{1}{2{\it\Lambda}Pe_{s}}\left(D_{zz}+\frac{1}{3}c\right),\quad \frac{\text{d}D_{zz}}{\text{d}z}=\frac{2}{15{\it\Lambda}Pe_{s}}m_{z}\quad \text{at}~z=\pm 1.\end{eqnarray}$$

Using this set of equations, we first proceed to derive a relation between the values of the concentration and wall-normal polarization at the boundaries. First, we integrate equation (3.6) across the channel width and use the second boundary condition in (3.7) to arrive at

(3.8) $$\begin{eqnarray}\int _{-1}^{1}D_{zz}(z)\,\text{d}z=0.\end{eqnarray}$$

Now, combining (3.4) and (3.5), integrating from $z$ to 1 and making use of the first boundary condition gives

(3.9) $$\begin{eqnarray}D_{zz}-2{\it\Lambda}Pe_{s}\frac{\text{d}m_{z}}{\text{d}z}+\frac{6{\it\Lambda}+1}{3}c=2{\it\Lambda}\,c(1).\end{eqnarray}$$

This relation can be integrated once more across the channel width. Using condition (3.8) together with the parity properties of $c$ and $m_{z}$ , this simplifies to

(3.10) $$\begin{eqnarray}c(\pm 1)=\left(1+\frac{1}{6{\it\Lambda}}\right)\mp Pe_{s}m_{z}(\pm 1),\end{eqnarray}$$

providing a simple relation between concentration and polarization at the walls. Inserting this relation into the first condition in (3.7) yields a new set of boundary conditions that does not involve the concentration:

(3.11a,b ) $$\begin{eqnarray}\frac{\text{d}m_{z}}{\text{d}z}=\frac{1}{2{\it\Lambda}Pe_{s}}\left(D_{zz}+\frac{6{\it\Lambda}+1}{18{\it\Lambda}}\right)\mp \frac{1}{6{\it\Lambda}}m_{z},\quad \frac{\text{d}D_{zz}}{\text{d}z}=\frac{2}{15{\it\Lambda}Pe_{s}}m_{z}\quad \text{at}~z=\pm 1.\end{eqnarray}$$

Equations (3.5) and (3.6), together with these boundary conditions, form a coupled system of second-order linear ordinary differential equations for $m_{z}$ and $D_{zz}$ that can be solved analytically. Once these variables are known, the concentration profile is easily obtained from the polarization by integration of (3.4) along with condition (3.10).

Figure 2. Equilibrium distributions in the absence of flow and for various swimming Péclet numbers $Pe_{s}$ (with ${\it\Lambda}=1/6$ ), obtained by numerical solution of (2.13) using finite volumes: (a) concentration  $c$ , (b) wall-normal polarization  $m_{z}$ , and (c) wall-normal nematic order parameter  $D_{zz}$ .

Solving these equations yields complicated expressions for $c$ , $m_{z}$ and $D_{zz}$ that are omitted here for brevity. The profiles, which are illustrated in figure 2 and will be discussed in more detail below, reveal one important finding: while a significant wall-normal polarization exists in the near-wall region, nematic alignment is relatively weak throughout the channel for ${\it\Lambda}\gtrsim 0.1$ . This suggests seeking a yet simpler solution that neglects nematic order altogether. If the moment expansion (2.18) is truncated after two terms, the equations for $c$ and $m_{z}$ simplify to

(3.12a,b ) $$\begin{eqnarray}m_{z}-2{\it\Lambda}Pe_{s}\frac{\text{d}c}{\text{d}z}=0,\quad -2{\it\Lambda}Pe_{s}^{2}\frac{\text{d}^{2}m_{z}}{\text{d}z^{2}}+\left(\frac{1}{6{\it\Lambda}}+1\right)m_{z}=0,\end{eqnarray}$$

subject to the conditions

(3.13a,b ) $$\begin{eqnarray}\frac{\text{d}m_{z}}{\text{d}z}=\frac{c}{6{\it\Lambda}Pe_{s}}\quad \text{at}~z=\pm 1\quad \text{and}\quad \int _{-1}^{1}c(z)\,\text{d}z=2.\end{eqnarray}$$

Solving these equations is straightforward and provides elegant expressions for the concentration and polarization profiles:

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle c(z)=\frac{B[6{\it\Lambda}\cosh B+\cosh Bz]}{6{\it\Lambda}B\cosh B+\sinh B}, & \displaystyle\end{eqnarray}$$

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle m_{z}(z)=\frac{6{\it\Lambda}Pe_{s}B^{2}\sinh Bz}{3(6{\it\Lambda}B\cosh B+\sinh B)}, & \displaystyle\end{eqnarray}$$

(3.16) $$\begin{eqnarray}B^{-1}={\it\Lambda}Pe_{s}\sqrt{\frac{12}{1+6{\it\Lambda}}}\end{eqnarray}$$

defines the dimensionless decay length of the excess concentration at the walls. In dimensional terms, this decay length is given by $B^{-1}H=\ell _{a}\sqrt{3/(1+6{\it\Lambda})}$ where $\ell _{a}=d_{t}/V_{s}$ . In the limit of strong propulsion ( ${\it\Lambda}\ll 1$ ), it simplifies to $\sqrt{3}\,\ell _{a}$ . In the limit of weak propulsion ( ${\it\Lambda}\gg 1$ ), it becomes $\ell _{d}/\sqrt{2}$ where $\ell _{d}=\sqrt{d_{t}/d_{r}}$ is a purely diffusive length scale. For Brownian particles, $\ell _{d}$ is typically of the order of the particle size  $L$ , though this may not be the case for active particles subject to athermal sources of noise. Next, we focus more precisely on these two limits by rescaling the governing equations with the appropriate scales identified here.

3.2. Strong propulsion limit: ${\it\Lambda}\rightarrow 0$

In the limit of small ${\it\Lambda}$ , the above discussion suggests rescaling the Smoluchowski equation using the accumulation length scale $\ell _{a}$ , yielding

(3.17) $$\begin{eqnarray}\cos {\it\theta}\frac{\partial {\it\Psi}}{\partial z}-\frac{\partial ^{2}{\it\Psi}}{\partial z^{2}}={\it\Lambda}{\rm\nabla}_{p}^{2}{\it\Psi},\end{eqnarray}$$

subject to the boundary condition

(3.18) $$\begin{eqnarray}\frac{\partial {\it\Psi}}{\partial z}=\cos {\it\theta}\,{\it\Psi}\quad \text{at}~z=\pm H^{\ast }.\end{eqnarray}$$

Here, $H^{\ast }=(2{\it\Lambda}Pe_{s})^{-1}$ is the channel half-height rescaled by the accumulation length scale $\ell _{a}$ . We recall that $d_{t}\rightarrow 0$ is an ill-posed limit in our theory and the analysis presented here really corresponds to the limit of ${\it\Lambda}\rightarrow 0$ with finite ${\it\Lambda}Pe_{s}$ (finite $H^{\ast }$ ) or, in terms of length scales, to the limit of $\ell _{p}\rightarrow \infty$ with finite $\ell _{a}$ . In dimensional terms, this corresponds to the limit of $d_{r}\rightarrow 0$ with finite $d_{t}$ . The gap-averaged isotropy constraint is now expressed as

(3.19) $$\begin{eqnarray}\int _{-H^{\ast }}^{H^{\ast }}{\it\Psi}\,\text{d}z=\frac{H^{\ast }}{2{\rm\pi}}.\end{eqnarray}$$

The leading-order solution corresponding to ${\it\Lambda}=0$ , which was previously obtained by Elgeti & Gompper (Reference Elgeti and Gompper2013), is written

(3.20) $$\begin{eqnarray}{\it\Psi}^{(0)}(z,{\it\theta})=\frac{H^{\ast }\cos {\it\theta}}{4{\rm\pi}\sinh (H^{\ast }\cos {\it\theta})}\exp (z\cos {\it\theta}),\end{eqnarray}$$

and it is easily seen that it satisfies zero wall-normal flux pointwise throughout the channel. In particular, it shows that wall accumulation is possible even in the absence of rotational diffusion and is simply a result of a coupling between self-propulsion, translational diffusion and confinement. This solution can then be corrected to order $O({\it\Lambda})$ by solving the first-order inhomogeneous equation

(3.21) $$\begin{eqnarray}\cos {\it\theta}\,{\it\Psi}^{(1)}(z,{\it\theta})-\frac{\partial {\it\Psi}^{(1)}}{\partial z}={\rm\nabla}_{p}^{2}\int _{-H^{\ast }}^{z}{\it\Psi}^{(0)}(z,{\it\theta})\,\text{d}z,\end{eqnarray}$$

subject to boundary condition (3.18). An exact analytical solution to this equation can again be obtained but is cumbersome and omitted here for brevity.

3.3. Weak propulsion limit: ${\it\Lambda}\rightarrow \infty$

In the limit of large ${\it\Lambda}$ , the Smoluchowski equation is rescaled using the diffusive length scale $\ell _{d}$ as

(3.22) $$\begin{eqnarray}\frac{1}{\sqrt{{\it\Lambda}}}\cos {\it\theta}\frac{\partial {\it\Psi}}{\partial z}-\frac{\partial ^{2}{\it\Psi}}{\partial z^{2}}={\rm\nabla}_{p}^{2}{\it\Psi},\end{eqnarray}$$

subject to

(3.23) $$\begin{eqnarray}\frac{\partial {\it\Psi}}{\partial z}=\frac{1}{\sqrt{{\it\Lambda}}}\cos {\it\theta}\,{\it\Psi}\quad \text{at}~z=\pm H^{\dagger },\end{eqnarray}$$

where $H^{\dagger }=(2\sqrt{{\it\Lambda}}Pe_{s})^{-1}$ . The leading-order solution in the limit of ${\it\Lambda}\rightarrow \infty$ is uniform and isotropic and corresponds to the case of a passive particle: ${\it\Psi}^{(0)}(z,{\it\theta})=1/4{\rm\pi}$ . It can be corrected asymptotically using a regular perturbation expansion in powers of ${\it\Lambda}^{-1/2}$ :

(3.24) $$\begin{eqnarray}{\it\Psi}(z,{\it\theta})={\it\Psi}^{(0)}(z,{\it\theta})+{\it\Lambda}^{-1/2}\,{\it\Psi}^{(1)}(z,{\it\theta})+{\it\Lambda}^{-1}\,{\it\Psi}^{(2)}(z,{\it\theta})+\cdots \,.\end{eqnarray}$$

Recursively solving for higher-order terms yields

(3.25) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\Psi}^{(1)}(z,{\it\theta})=\frac{3}{4{\rm\pi}\sqrt{2}\cosh (\sqrt{2}H^{\dagger })}\sinh (\sqrt{2}z)\cos {\it\theta},\qquad & \displaystyle\end{eqnarray}$$

(3.26) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\Psi}^{(2)}(z,{\it\theta})=\left[-\frac{1}{15}\frac{\cosh \sqrt{2}z}{\cosh (\sqrt{2}H^{\dagger })}+\frac{\tanh (\sqrt{2}H^{\dagger })}{5\sqrt{3}}\frac{\cosh (\sqrt{6}z)}{\sinh (\sqrt{6}H^{\dagger })}\right]\left(\cos ^{2}{\it\theta}-\frac{1}{3}\right),\qquad & \displaystyle\end{eqnarray}$$

${\it\Psi}^{(1)}$

${\it\Psi}^{(2)}$

3.4. Numerical results and discussion

Figure 2 shows the full numerical solution for the concentration $c$ , wall-normal polarization $m_{z}$ and nematic order parameter $D_{zz}$ obtained by finite-volume solution of the Smoluchowski equation (2.13) as described in appendix B. Here, we fix the value of ${\it\Lambda}$ and focus on the effect of $Pe_{s}$ , which is an inverse measure of confinement. The concentration profiles shown in figure 2(a) exhibit significant accumulation of particles near the boundaries, especially at low values of $Pe_{s}$ . As anticipated, this accumulation is accompanied by polarization towards the boundaries as a direct consequence of the boundary condition (2.25), as well as by a weak nematic alignment. As $Pe_{s}$ increases, the spatial heterogeneity and anisotropy near the walls progressively extend through the entire channel as the two boundary layers thicken and eventually merge. Further increase in the swimming Péclet number leads to a flattening of the profiles, which is especially significant when $Pe_{s}>1$ . This flattening is a direct consequence of the scaling of translational diffusion with $Pe_{s}^{2}$ in (2.13), causing it to overwhelm self-propulsion, which scales with $Pe_{s}$ . The influence of ${\it\Lambda}$ is illustrated in figure 3, where it is seen to be similar to that of $Pe_{s}$ : increasing ${\it\Lambda}$ leads to a thickening of the boundary layers and flattening of the concentration profiles, again due to the scaling of translational diffusion with ${\it\Lambda}$ in (2.13).

Figure 3. Equilibrium distributions in the absence of flow and for various values of ${\it\Lambda}$ (with $Pe_{s}=0.25$ ), obtained by numerical solution of (2.13) using finite volumes: (a) concentration  $c$ , (b) wall-normal polarization  $m_{z}$ , and (c) wall-normal nematic order parameter  $D_{zz}$ . Solutions based on moment equations are nearly identical, as illustrated in figure 4.

Figure 4. The relative r.m.s. error for the concentration between the finite-volume solution and the two-moment analytical solution (3.14) for different values of ${\it\Lambda}$ . Solutions based on moment equations are nearly identical to the finite-volume solution for sufficiently large values of ${\it\Lambda}$ .

The finite-volume numerical solution of the full conservation equation (2.13) is in excellent quantitative agreement with the two- and three-moment approximations derived previously, which are not shown in figure 2 as they are nearly indistinguishable over the entire channel width as long as ${\it\Lambda}\gtrsim 0.1$ . The root-mean-square (r.m.s.) error between the two-moment solution of (3.4) and the finite-volume solution is indeed plotted in figure 4, where it remains below $10^{-3}$ for all values of $Pe_{s}$ considered here when ${\it\Lambda}\gtrsim 0.1$ . This finding may seem quite surprising considering the strong approximation made when truncating expansion (2.18) after only two terms, and strongly validates the use of approximate moment equations such as (2.19)–(2.24) when modelling active suspensions, at least in the absence of flow. For very small values of ${\it\Lambda}$ , however, nematic alignment at the walls becomes significant, as seen in figure 3(c), so that the nematic tensor can no longer be neglected and the two-moment solution loses its accuracy; in this case, the alternative expressions derived in the small ${\it\Lambda}$ limit in § 3.2 can be used instead.

Figure 5. Wall accumulation in the absence of flow as a function of $Pe_{s}$ (at ${\it\Lambda}=1/6$ ): (a) concentration $c(\pm 1)$ at the walls; (b) boundary layer thickness ${\it\delta}$ , defined as the distance from the wall where $c(1-{\it\delta})=1$ ; and (c) fraction ${\it\delta}^{\ast }$ of particles inside the boundary layer, defined as the integral of $c(z)$ over the boundary layer thickness. The solid line shows the theoretical prediction based on the two-moment solution (3.14), and symbols show full numerical results using finite volumes.

The influence of $Pe_{s}$ on wall accumulation is analysed more quantitatively in figure 5, showing the values of the wall concentration $c(\pm 1)$ , the boundary layer thickness ${\it\delta}$ defined as the distance from the wall where $c(1-{\it\delta})=1$ , and the fraction ${\it\delta}^{\ast }$ of particles inside the boundary layer defined as

(3.27) $$\begin{eqnarray}{\it\delta}^{\ast }=\int _{1-{\it\delta}}^{1}c(z)\,\text{d}z.\end{eqnarray}$$

Analytical expressions for these quantities can be derived from the two-moment solution (3.14). In particular, the boundary layer thickness is obtained as

(3.28) $$\begin{eqnarray}{\it\delta}(Pe_{s})=1-\frac{1}{B}\log \left\{\frac{\sinh B}{B}\pm \left[\left(\frac{\sinh B}{B}\right)^{2}-1\right]^{1/2}\right\},\end{eqnarray}$$

which has the two limits

(3.29a,b ) $$\begin{eqnarray}\lim _{Pe_{s}\rightarrow 0}{\it\delta}(Pe_{s})=0\quad \text{and}\quad \lim _{Pe_{s}\rightarrow \infty }{\it\delta}(Pe_{s})=1-\frac{1}{\sqrt{3}}.\end{eqnarray}$$

Similarly, the fraction of particles inside the boundary layer is given by

(3.30) $$\begin{eqnarray}{\it\delta}^{\ast }(Pe_{s})=1-\frac{6{\it\Lambda}B(1-{\it\delta})\cosh B+\sinh [B(1-{\it\delta})]}{6{\it\Lambda}B\cosh B+\sinh B},\end{eqnarray}$$

and has the same limits as ${\it\delta}(Pe_{s})$ when $Pe_{s}\rightarrow 0$ and  $\infty$ .

As shown in figure 5(a), the wall concentration reaches its maximum in the limit of $Pe_{s}\rightarrow 0$ , and steadily decreases towards 1 as $Pe_{s}$ increases due to the smoothing effect of translational diffusion. This is accompanied by an increase in the boundary layer thickness ${\it\delta}$ , which asymptotes at high values of $Pe_{s}$ . The fraction ${\it\delta}^{\ast }$ of particles near the walls shows a similar trend, but interestingly also exhibits a weak maximum for $Pe_{s}\approx 1.135$ when wall accumulation due to self-propulsion and translational diffusion are of similar magnitudes; at this value of $Pe_{s}$ , ${\it\delta}^{\ast }\approx 0.46$ corresponding to nearly half the particles being trapped near the walls. As previously observed in figure 4, excellent agreement is obtained between the two-moment approximation and the numerical solution of the full governing equations.

4.1. Weak-flow limit: regular asymptotic expansion

We now proceed to analyse the effects of an external pressure-driven flow, first focusing on the case of a weak flow for which $Pe_{f}\ll 1$ . Since the parameter ${\it\Lambda}$ is fixed for a given type of swimmers, we keep it constant in the rest of the paper and focus on the effects of $Pe_{s}$ and $Pe_{f}$ . The form of the governing equations suggests seeking an approximate solution as a regular expansion of the moments of the distribution function in powers of $Pe_{f}$ . The leading-order $O(Pe_{f}^{0})$ solution corresponding to the absence of flow was previously calculated in § 3. It is henceforth denoted by $c^{(0)}$ , $\boldsymbol{m}^{(0)}$ , $\unicode[STIX]{x1D63F}^{(0)}$ , and we recall that $m_{y}^{(0)}=D_{yz}^{(0)}=0$ . Inspection of the moment equations (2.19)–(2.24) reveals that the interaction of the applied shear profile $S(z)$ with this leading-order solution perturbs $m_{y}$ and $D_{yz}$ at order $O(Pe_{f})$ . On the other hand, $c$ , $m_{z}$ , $D_{zz}$ and $D_{yy}$ are only perturbed by the flow at order $O(Pe_{f}^{2})$ due to its interaction with $m_{y}$ and $D_{yz}$ . Based on these observations, we expand the solution as

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle c(z)=c^{(0)}(z)+Pe_{f}^{2}\,c^{(2)}(z)+O(Pe_{f}^{3}), & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle m_{z}(z)=m_{z}^{(0)}(z)+Pe_{f}^{2}m_{z}^{(2)}(z)+O(Pe_{f}^{3}), & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle D_{zz}(z)=D_{zz}^{(0)}(z)+Pe_{f}^{2}D_{zz}^{(2)}(z)+O(Pe_{f}^{3}), & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle D_{yy}(z)=D_{yy}^{(0)}(z)+Pe_{f}^{2}D_{yy}^{(2)}(z)+O(Pe_{f}^{3}), & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle m_{y}(z)=Pe_{f}m_{y}^{(1)}(z)+O(Pe_{f}^{3}), & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle D_{yz}(z)=Pe_{f}D_{yz}^{(1)}(z)+O(Pe_{f}^{3}). & \displaystyle\end{eqnarray}$$

$m_{y}$

$D_{yz}$

$O(Pe_{f})$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle Pe_{s}\frac{\text{d}D_{yz}^{(1)}}{\text{d}z}-2{\it\Lambda}Pe_{s}^{2}\frac{\text{d}^{2}m_{y}^{(1)}}{\text{d}z^{2}}+m_{y}^{(1)}=\frac{2}{5}S(z)m_{z}^{(0)}, & \displaystyle\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{Pe_{s}}{5}\frac{\text{d}m_{y}^{(1)}}{\text{d}z}-2{\it\Lambda}Pe_{s}^{2}\frac{\text{d}^{2}D_{yz}^{(1)}}{\text{d}z^{2}}+3D_{yz}^{(1)}=S(z)\left(\frac{1}{10}c^{(0)}+\frac{1}{2}D_{zz}^{(0)}\right), & \displaystyle\end{eqnarray}$$

(4.9a,b ) $$\begin{eqnarray}\frac{\text{d}m_{y}^{(1)}}{\text{d}z}=\frac{1}{2{\it\Lambda}Pe_{s}}D_{yz}^{(1)},\quad \frac{\text{d}D_{yz}^{(1)}}{\text{d}z}=\frac{1}{10{\it\Lambda}Pe_{s}}m_{y}^{(1)}\quad \text{at}~z=\pm 1.\end{eqnarray}$$

Note that the forcing terms on the right-hand sides of equations (4.7) and (4.8) are known and capture the interaction of the local shear rate $S(z)$ with the equilibrium distributions in the absence of flow.

Figure 6. Effect of a weak applied flow: leading-order $O(Pe_{f})$ corrections of (a) streamwise polarization $m_{y}$ and (b) shear nematic alignment $D_{yz}$ for different values of the swimming Péclet number, obtained by numerical solution of (4.7)–(4.9).

A numerical solution of equations (4.7)–(4.9) is plotted in figure 6 for different values of $Pe_{s}$ . At low values of the swimming Péclet number, figure 6(a) shows an upstream polarization ( $m_{y}<0$ ) near the boundaries, and a downstream polarization ( $m_{y}>0$ ) near the centre of the channel. The upstream polarization, which has previously been observed in both experiments and simulations and is at the origin of the well-known phenomenon of upstream swimming, is a simple and direct consequence of the shear rotation of the particles near the wall, which tend to point towards the walls in the absence of flow, as explained in § 3. This interaction is encapsulated in the right-hand side in (4.7). The downstream polarization near the centreline is a more subtle effect arising from self-propulsion through the first term on the left-hand side of (4.7). As $Pe_{s}$ increases and the boundary layers thicken, upstream swimming becomes weaker near the boundaries due to the weaker wall-normal polarization there; however, $m_{y}$ is also observed to become negative across the entire channel owing to the thickening of the polarized boundary layers into the bulk of the channel, as previously shown in figure 2(b).

The mean streamwise swimming velocity $\overline{V}_{y}$ of the active particles with respect to the imposed flow can be defined in terms of the polarization as

(4.10) $$\begin{eqnarray}\overline{V}_{y}=\frac{1}{2}\int _{-1}^{1}Pe_{s}\,m_{y}(z)\,\text{d}z=\frac{Pe_{s}Pe_{f}}{2}\int _{-1}^{1}m_{y}^{(1)}(z)\,\text{d}z=Pe_{s}Pe_{f}\,\overline{m}_{y}^{(1)}.\end{eqnarray}$$

An expression for $\overline{m}_{y}^{(1)}$ can be derived based on the moment equations. We first take the cross-sectional average of (4.7) and use the first boundary condition to obtain

(4.11) $$\begin{eqnarray}\overline{m}_{y}^{(1)}=-\frac{1}{5}\int _{-1}^{1}z\,m_{z}^{(0)}(z)\,\text{d}z.\end{eqnarray}$$

Since $m_{z}^{(0)}$ is an odd function of $z$ with $m_{z}^{(0)}(z)\geqslant 0$ for $z\geqslant 0$ , the integrand on the right-hand side is always positive across the channel, and therefore the mean upstream polarization is negative: $\overline{m}_{y}^{(1)}<0$ . This also implies that $\overline{V}_{y}<0$ , i.e. there is a net upstream flux of particles against the mean flow for all values of ${\it\Lambda}$ and $Pe_{s}$ in the weak-flow limit. Using (3.5) for $m_{z}^{(0)}(z)$ , we can rewrite the right-hand side as

(4.12) $$\begin{eqnarray}\overline{m}_{y}^{(1)}=-\frac{1}{5\displaystyle \left(\frac{1}{6{\it\Lambda}}+1\right)}\left[2{\it\Lambda}Pe_{s}^{2}\int _{-1}^{1}z\,\frac{\text{d}^{2}m_{z}^{(0)}}{\text{d}z^{2}}\,\text{d}z-Pe_{s}\int _{-1}^{1}z\,\frac{\text{d}D_{zz}^{(0)}}{\text{d}z}\,\text{d}z\right].\end{eqnarray}$$

After integration by parts and application of the boundary condition on $m_{z}^{(0)}(z)$ together with (3.8), this simplifies to

(4.13) $$\begin{eqnarray}\overline{m}_{y}^{(1)}=-\frac{2Pe_{s}}{15\left(\displaystyle \frac{1}{6{\it\Lambda}}+1\right)}[c^{(0)}(1)-6{\it\Lambda}Pe_{s}\,m_{z}^{(0)}(1)].\end{eqnarray}$$

Recalling that $c^{(0)}(1)$ and $m_{z}^{(0)}(1)$ are related via (3.10), we obtain two expressions for the mean streamwise swimming velocity in terms of either the concentration or the wall-normal polarization at the top wall in the absence of flow:

(4.14) $$\begin{eqnarray}\overline{V}_{y}=-\frac{4{\it\Lambda}}{5}Pe_{s}^{2}Pe_{f}[c^{(0)}(1)-1]=-\frac{2}{15}Pe_{s}^{2}Pe_{f}[1-6{\it\Lambda}Pe_{s}m_{z}^{(0)}(1)].\end{eqnarray}$$

Since the concentration at the wall in the absence of flow always exceeds the mean when $Pe_{s}>0$ , equation (4.14) again confirms that $\overline{V}_{y}<0$ . If we further make use of the simplified two-moment analytical solution (3.14) for the concentration profile, we arrive at a simple expression for the mean upstream velocity in terms of the swimming and flow Péclet numbers:

(4.15) $$\begin{eqnarray}\overline{V}_{y}=-\frac{4{\it\Lambda}}{5}Pe_{s}^{2}Pe_{f}\left[\frac{B\cosh B-\sinh B}{6{\it\Lambda}B\cosh B+\sinh B}\right].\end{eqnarray}$$

This simple analytical prediction for $\overline{V}_{y}$ will be tested against numerical simulations at arbitrary $Pe_{f}$ in § 4.2, where it will be shown to provide an excellent estimate for the swimming flux up to $Pe_{f}\approx 2$ .

Figure 7. Equilibrium concentration profiles (at  ${\it\Lambda}=1/6$ ) for (a)  $Pe_{s}=0.25$ (strong wall accumulation) and (b)  $Pe_{s}=1.0$ (weak accumulation) and for various values of the flow Péclet number $Pe_{f}$ , obtained by finite-volume solution of the governing equation (2.13).

The effects of the external flow on nematic alignment are also illustrated in figure 6(d), where $D_{yz}$ is found to vary almost linearly across the channel width and has the same sign as the external shear rate profile $S(z)$ . The right-hand side in (4.8) provides a simple explanation for these findings, where we see that shear nematic alignment results primarily from the interaction of the flow with the concentration profile and with the wall-normal nematic alignment. As $Pe_{s}$ increases, shear nematic alignment decreases due to the decrease in $c$ and $D_{zz}$ inside the boundary layers as seen in figure 2(a,c), and to self-propulsion through the first term on the left-hand side of (4.8).

4.2. Strong-flow limit: scaling analysis

As we shall see in § 4.3 and figure 7, the regime of high flow Péclet number is also quite interesting, as it can result in a depletion near the channel centreline surrounded by regions where particles become trapped. The thickness of this depletion region will be found to decrease with increasing flow strength, suggesting the presence of another boundary layer near $z=0$ in the limit of $Pe_{f}\gg 1$ . Insight into this regime can be gained by analysing the behaviour of the governing equation (2.13) for $Pe_{f}\gg 1$ and $Pe_{s}\ll 1$ . If the swimming Péclet number is low, the wall boundary layers are very thin and have negligible impact on the dynamics in the bulk of the channel. Inspection of (2.13) suggests that, in the outer region away from both the channel walls and the centreline, the dominant balance is between shear alignment and rotational diffusion:

(4.16) $$\begin{eqnarray}\frac{Pe_{f}}{2}\,S(z)\boldsymbol{{\rm\nabla}}_{p}\boldsymbol{\cdot }[\cos {\it\theta}(\unicode[STIX]{x1D644}-\boldsymbol{p}\boldsymbol{p})\boldsymbol{\cdot }\hat{\boldsymbol{y}}\,{\it\Psi}]\approx \frac{1}{2}{\rm\nabla}_{p}^{2}{\it\Psi}.\end{eqnarray}$$

In this region, the concentration is expected to be nearly uniform, and the particle orientation distribution is primarily nematic as a result of the competition between the local shear rate and rotational diffusion (as would occur in a passive rod suspension). This corresponds to the shear-trapping region where cross-streamline migration is very weak owing to the strong alignment with the flow.

However, as we move closer and closer to the centreline, the local shear rate decreases, causing a concomitant decrease in shear alignment and increase in cross-streamline migration due to self-propulsion. This transition corresponds to the edge of the central boundary layer from which particles are depleted, and the position ${\it\delta}_{D}$ of this transition region (or half-thickness of the depletion layer) can be estimated by balancing the magnitudes of the terms describing self-propulsion and shear alignment in (2.13), i.e.

(4.17) $$\begin{eqnarray}\frac{Pe_{s}}{{\it\delta}_{D}}\sim \frac{Pe_{f}}{2}{\it\delta}_{D},\end{eqnarray}$$

from which we find

(4.18) $$\begin{eqnarray}{\it\delta}_{D}\approx C\sqrt{{\it\chi}},\end{eqnarray}$$

where the prefactor $C$ is a numerical constant and where we have defined

(4.19) $$\begin{eqnarray}{\it\chi}=\frac{Pe_{s}}{Pe_{f}}=\frac{V_{s}}{2\dot{{\it\gamma}}_{w}H}.\end{eqnarray}$$

The dimensionless group ${\it\chi}$ can be interpreted as the ratio of the time scale $\dot{{\it\gamma}}_{w}^{-1}$ it takes a particle to align with the flow over the characteristic time scale $2H/V_{s}$ it takes it to swim across the channel width. If ${\it\chi}$ is small, particles align with the flow much faster than they can cross the channel, leading to significant shear trapping. On the other hand, if ${\it\chi}$ is large, particles cross the channel much faster than they align with the flow and shear trapping does not occur. As we show in appendix C, this scaling for ${\it\delta}_{D}$ can indeed also be derived by considering the individual trajectories of deterministic swimmers released from the centreline, which can be shown to become trapped at a distance of the order of ${\it\delta}_{D}$ . It will also be shown to agree quite well with numerical results in § 4.3, where we will find that ${\it\delta}_{D}\approx 2.404\sqrt{{\it\chi}}$ provides an excellent estimate for the thickness of the depletion layer when $Pe_{s}\lesssim 0.25$ and $Pe_{f}\gtrsim 50$ .

To gain further understanding of the effect of shear rate on the intensity of depletion, we rescale lengths by ${\it\delta}_{D}$ inside the central boundary layer to rewrite the governing equation (2.13) as

(4.20) $$\begin{eqnarray}\frac{{\it\Gamma}}{C}\cos {\it\theta}\frac{\partial {\it\Psi}}{\partial z}-2{\it\Lambda}\frac{{\it\Gamma}^{2}}{C^{2}}\frac{\partial ^{2}{\it\Psi}}{\partial z^{2}}-\frac{C{\it\Gamma}}{2}z\,\boldsymbol{{\rm\nabla}}_{p}\boldsymbol{\cdot }[\cos {\it\theta}(\unicode[STIX]{x1D644}-\boldsymbol{p}\boldsymbol{p})\boldsymbol{\cdot }\hat{\boldsymbol{y}}\,{\it\Psi}]=\frac{1}{2}{\rm\nabla}_{p}^{2}{\it\Psi},\end{eqnarray}$$

where the dimensionless group ${\it\Gamma}=\sqrt{Pe_{s}Pe_{f}}$ emerges as the most significant parameter governing the profile of the depletion layer. Unsurprisingly, we find that self-propulsion and shear rotation have the same magnitude upon rescaling. In this region, self-propulsion, which scales with ${\it\Gamma}$ , has the effect of enhancing depletion by driving particles away from the centreline; this competes against translational diffusion, scaling with ${\it\Gamma}^{2}$ , which has the effect of smoothing concentration gradients and thus hampers depletion. This suggests the following dependence of the concentration profile on $Pe_{f}$ . As flow strength is increased from small values, the depletion layer forms and continually narrows according to (4.18) for ${\it\delta}_{D}$ . As long as ${\it\Gamma}<1$ , self-propulsion dominates translational diffusion and increasing $Pe_{f}$ (and therefore ${\it\Gamma}$ ) enhances depletion. This trend reverses when ${\it\Gamma}\sim O(1)$ , when translational diffusion starts to overcome self-propulsion, leading to a subsequent decrease in the strength of depletion for ${\it\Gamma}>1$ . This qualitative explanation for the non-monotonic dependence of the strength of depletion upon ${\it\Gamma}$ (and hence upon the mean shear rate of the imposed Poiseuille flow) is consistent with the experimental observations of Rusconi et al. (Reference Rusconi, Guasto and Stocker2014), and is also borne out by numerical solutions of the governing equations, as we describe next.

4.3. Arbitrary flow strengths: finite-volume calculations and discussion

We now test and extend the key predictions from the weak-flow asymptotics and strong-flow scaling analysis from the preceding sections by performing finite-volume numerical simulations of the governing equation (2.13) for arbitrary values of $Pe_{s}$ and $Pe_{f}$ using the algorithm of appendix C. Typical concentration profiles are illustrated in figure 7 for various values of $Pe_{f}$ , and for the two values of $Pe_{s}=0.25$ and $1.0$ corresponding to cases where wall accumulation in the absence of flow is strong and weak, respectively. In both cases, the leading effect of the external flow on $c$ is to decrease wall accumulation. This trend is easily understood as a result of the alignment of the particles with the flow, which reduces wall-normal polarization and thereby hinders accumulation. This decrease in accumulation also results in a net increase in the concentration in the central parts of the channel and in the flattening of the profiles in the strong-flow limit. When $Pe_{s}$ is small, as in figure 7(a), a depletion layer is also observed to form near the channel centreline and to progressively narrow with increasing $Pe_{f}$ , in agreement with the theoretical predictions of § 4.2. At high values of $Pe_{f}$ , the three distinct regions identified in § 4.2 (wall accumulation, shear trapping and centreline depletion) in fact become clearly visible. However, if the swimming Péclet number is increased to $Pe_{s}=1.0$ , as in figure 7(b), the thickening of the wall boundary layers suppresses shear trapping and depletion at the centreline, leading to a nearly uniform concentration profile in the strong-flow limit.

Figure 8. Equilibrium streamwise and wall-normal polarization profiles (at  ${\it\Lambda}=1/6$ ) for (a,c) $Pe_{s}=0.25$ and (b,d) $Pe_{s}=1.0$ and for various values of the flow Péclet number $Pe_{f}$ , obtained by finite-volume solution of the governing equation (2.13). The streamwise polarization $m_{y}$ is shown in (a,b), and the wall-normal polarization $m_{z}$ in (c,d).

Corresponding profiles for the wall-normal and streamwise polarization are also shown in figure 8. As expected, rotation of the particles by the flow causes a decrease in the wall-normal polarization, and also results in a non-zero streamwise polarization $m_{y}$ , as previously discussed in § 4.1. This streamwise polarization is especially strong in the near-wall region, where $m_{y}$ is negative, indicating upstream swimming. It is significantly weaker near the centre of the channel, where it is found to be positive for $Pe_{s}=0.25$ but remains negative across the entire channel when $Pe_{s}=1.0$ owing to the overlap of the two wall boundary layers.

Figure 9. Effect of swimming and flow Péclet numbers on: (a) wall concentration $c(\pm 1)$ , (b) streamwise polarization $m_{y}(\pm 1)$ at the channel walls, and (c) streamwise polarization $m_{y}(0)$ at the channel centreline.

These trends are made more quantitative in figure 9, showing the dependence of $c(\pm 1)$ , $m_{y}(\pm 1)$ and $m_{y}(0)$ on the swimming and flow Péclet numbers. As previously discussed, the wall concentration is seen to decrease with increasing flow strength irrespective of the value of $Pe_{s}$ , and asymptotically tends to 1 in the strong-flow limit as the concentration profiles flatten. Figure 9(b) shows that the streamwise polarization at the walls is always negative, which implies that the active particles always swim upstream near the boundaries. Interestingly, we find that there is maximum upstream swimming at $Pe_{f}\approx 10$ , and the upstream motion is reduced at higher values of the flow Péclet number. The streamwise polarization at the channel centreline shows complex trends, as shown in figure 9(c). As predicted by the weak-flow asymptotic analysis of § 4.1, $m_{y}(0)$ is found to be positive for low values of $Pe_{s}$ and negative for high values of $Pe_{s}$ . Its absolute value increases with flow strength in both cases up to $Pe_{f}\approx 10$ , beyond which further increasing flow strength reduces the polarization. The decrease in both $m_{y}(\pm 1)$ and $m_{y}(0)$ at high $Pe_{f}$ is a likely consequence of the dominant effect of the shear alignment term in (2.13), which promotes nematic rather than polar order.

Figure 10. (a) Magnitude of the average upstream swimming velocity $|\overline{V}_{y}|$ as a function of $Pe_{f}$ for different values of $Pe_{s}$ (at  ${\it\Lambda}=1/6$ ), and (b) dependence of $|\overline{V}_{y}|/Pe_{f}$ on $Pe_{s}$ for different values of $Pe_{f}$ . Symbols show finite-volume numerical simulations, and dotted lines show the theoretical prediction of (4.14).

The dependence of the average streamwise swimming velocity $\overline{V}_{y}$ defined in (4.10) on both Péclet numbers is shown in figure 10, where numerical results are compared with the weak-flow theoretical prediction of (4.14). Consistent with figure 9(b) for the streamwise polarization at the walls, we find that $\overline{V}_{y}<0$ , and that $|\overline{V}_{y}|$ first increases nearly linearly with $Pe_{f}$ in agreement with the predictions of § 4.1. This increase persists up to $Pe_{f}\approx 10$ , beyond which $|\overline{V}_{y}|$ starts decreasing again. Excellent quantitative agreement is found with (4.14) for $Pe_{f}\lesssim 2.0$ . This is confirmed in figure 10(b), showing the dependence of $|\overline{V}_{y}|/Pe_{f}$ on swimming Péclet number: the upstream velocity is found to increase with $Pe_{s}$ , primarily as a result of the corresponding increase in swimming speed of individual particles, and a collapse of all the curves onto the theoretical prediction of (4.14) is observed when $Pe_{f}\lesssim 2.0$ .

As seen in figure 7(a), shear trapping and centreline depletion are observed in the central portion of the channel at high flow Péclet number if $Pe_{s}$ is sufficiently low. This is illustrated more clearly in figure 11, where concentration and wall-normal polarization profiles are shown in the central portion of the channel for various values of the flow Péclet number and for $Pe_{s}=0.125$ . This value was chosen to match the experiments of Rusconi et al. (Reference Rusconi, Guasto and Stocker2014), where the following parameters were reported: $V_{s}=50~{\rm\mu}\text{m}$ , $d_{r}=1~\text{s}^{-1}$ and $2H=400~{\rm\mu}\text{m}$ . As seen in figure 11(a), increasing $Pe_{f}$ from zero first results in a decrease in the concentration at the centreline, corresponding to the formation of the depletion layer. As the concentration decreases, the width of the depletion layer is also found to decrease. This trend continues up to $Pe_{f}\approx 20$ , above which the concentration at the centreline starts increasing again, even though the depletion layer keeps narrowing. These trends are in very good agreement with the experiments of Rusconi et al. (Reference Rusconi, Guasto and Stocker2014), who also reported a non-monotonic dependence of the strength of depletion on shear rate; in fact, the profiles shown in figure 11 are very similar to the experimental profiles at equivalent values of  $Pe_{f}$ . The trends in the concentration profile are easily understood based on figure 11(b) for the wall-normal polarization, which reflects the net swimming velocity across the channel and provides insight into cross-streamline migration. Indeed, the polarization profiles exhibit peaks on both sides of the depletion layer, corresponding to a strong migration away from the centre. These peaks increase in magnitude and also shift towards the centreline as flow strength increases and the depletion layer narrows. Beyond those peaks, $m_{z}$ quickly decays to zero where the concentration profiles plateau in accordance with (2.28) and shear trapping of the particles takes place.

Figure 11. (a) Concentration profiles in the central portion of the channel for $Pe_{s}=0.125$ and various values of the flow Péclet number $Pe_{f}$ , obtained by finite-volume solution of (2.13). (b) Corresponding profiles of the wall-normal polarization  $m_{z}$ .

Figure 12. (a) Depletion layer thickness ${\it\delta}_{D}$ , defined as the distance from the centreline where the wall-normal polarization reaches its maximum, as a function of $\sqrt{{\it\chi}}=\sqrt{Pe_{s}/Pe_{f}}$ . (b) Depletion index $A_{D}$ defined in (4.21) as a function of ${\it\Gamma}=\sqrt{Pe_{s}Pe_{f}}$ .

These trends are tested more quantitatively against the strong-flow scaling analysis of § 4.2 in figure 12. We first define the thickness ${\it\delta}_{D}$ of the depletion layer as the distance from the centreline where $m_{z}$ reaches its maximum, when such a maximum exists. Based on the analysis of § 4.2, we expect ${\it\delta}_{D}$ to scale linearly with $\sqrt{{\it\chi}}=\sqrt{Pe_{s}/Pe_{f}}$ in strong flows, and this is indeed confirmed in figure 12(a). We find that ${\it\delta}_{D}$ can only be defined when $\sqrt{{\it\chi}}\lesssim 0.16$ or $Pe_{f}\gtrsim 40Pe_{s}$ , which corresponds to the shear-trapping regime. Best agreement with the scaling prediction is obtained in the low- $Pe_{s}$ and high- $Pe_{f}$ limit, and a linear least-squares fit to the data for $Pe_{s}\leqslant 0.25$ and $Pe_{f}\geqslant 50$ shows that ${\it\delta}_{D}\approx 2.404\sqrt{{\it\chi}}$ . As $Pe_{s}$ increases, the numerical results depart from this prediction, primarily due to the thickening of the wall boundary layers, which causes them to interact with the parts of the channel where shear trapping and depletion occur. We further quantify the shape of the depletion layer by introducing a depletion index $A_{D}$ measuring the amount of particles depleted from the centre due to trapping in high-shear regions:

(4.21) $$\begin{eqnarray}A_{D}=\int _{0}^{{\it\delta}_{D}}c(z)\,\text{d}z-{\it\delta}_{D}c({\it\delta}_{D}).\end{eqnarray}$$

As we argued in § 4.2 based on (4.20), the shape of the depletion layer is expected to depend upon ${\it\Gamma}=\sqrt{Pe_{s}Pe_{f}}$ , and indeed the numerical data for the depletion index for various values of $Pe_{s}$ and $Pe_{f}$ are found to collapse onto a master curve when plotted versus ${\it\Gamma}$ in figure 12(b). In agreement with the trends observed in figure 11(a), the depletion index shows a non-monotonic dependence on ${\it\Gamma}$ , with maximum depletion occurring for ${\it\Gamma}\approx 2$ .

Figure 13. Schematic summary of the dynamics in the limits of $Pe_{s}\ll 1$ and $Pe_{f}\gg 1$ . The channel can be roughly divided into three regions: (A) near the walls, particles accumulate in a boundary layer of thickness ${\it\delta}\sim {\it\Lambda}Pe_{s}$ ; (B) away from the walls and centreline, strong nematic alignment by the flow leads to shear trapping and a nearly uniform concentration profile; (C) near the centreline, particle propulsion leads to a depletion layer of thickness ${\it\delta}_{D}\sim {\it\Gamma}$ . Only the left half of the channel $z\in [-1,0]$ is shown; the corresponding other half can be obtained by symmetry and by noting that $m_{z}$ is an even function of $z$ , whereas $m_{y}$ and $D_{yz}$ are both odd functions.

The dynamics in the limits of $Pe_{s}\ll 1$ and $Pe_{f}\gg 1$ are summarized schematically in figure 13, where the channel can be roughly divided into three distinct regions. Region (A), with thickness ${\it\delta}\sim {\it\Lambda}Pe_{s}$ , abuts the channel wall and is characterized by wall accumulation and a net polarization towards the wall. These effects occur even in the absence of flow, and are in fact mitigated by the flow, which tends to decrease the wall concentration and rotate particles to induce upstream polarization. Away from both the wall and the channel centreline is region (B), where the concentration profile is nearly uniform and shear trapping occurs. Here, polarization is weak but there is a strong nematic alignment of the particles due to the applied shear. The local shear rate decreases in magnitude as we approach the centreline and enter region (C), which has a characteristic thickness of ${\it\delta}_{D}\sim \sqrt{Pe_{s}/Pe_{f}}$ . In this last region, particles are depleted due to a net polarization towards the walls, which drives migration away from the centre but is counterbalanced by translational diffusion. Increasing $Pe_{s}$ causes both regions (A) and (C) to widen, up to a point where they merge and the three regions can no longer be distinguished. Increasing $Pe_{f}$ , on the other hand, tends to weaken wall accumulation but does not change the thickness of region (A), while it also causes the narrowing of region (C).