3.1. General considerations on macro-transport processes subject to anisotropic diffusion

Note that, due to the convection term by the shear flow, the solution of (2.12) will have the form of a mass cloud moving to infinity at long times. Thus, we hereafter study the solution $P(x^\prime,t)=P(x- \bar {U}t,t)$ in a moving coordinate, where $\bar {U}$ is the drift to be quantified. Then, the transport equation of the probability density function turns into

(3.1)\begin{equation} \frac{\partial P}{\partial t}+\left[ {{Pe}}_f U_f(z)+{{{{Pe}}}_{s}}\sin \theta -{{{\bar{U}}}} \right]\frac{\partial P}{\partial x^\prime}-{{D}_{xx}}(\theta )\frac{{{\partial }^{2}}P}{\partial {{x}^{\prime 2}}}-2{{D}_{xz}}(\theta )\frac{{{\partial }^{2}}P}{\partial x^\prime \partial z}+\mathcal{L}P=0. \end{equation}

Given the slow variation in both the initial distribution and its corresponding solution, it is appropriate to investigate the system's effective dynamics on the extended space and time scales. In this work, we present a small parameter $\epsilon$ as the ratio of the local Péclet number to the global Péclet number. Here, the limit $\epsilon \rightarrow 0$ corresponds to the asymptotically long times during which a single particle samples the entire local spaces ($z$- and $\theta$-spaces) at a longitudinal position before spreading downstream. The local Péclet number can be manifested as the collective influence of translational diffusion along parallel, perpendicular and/or cross directions in the local position space, as well as rotational diffusion in the orientation space. Nevertheless, the combined effect of local aspects is typically overshadowed by the influence of the global (longitudinal) Péclet number. As the convective field tends to an average of zero, the resultant behaviour of $P$ mirrors that characteristic of a pure diffusion process. To address this effect, we resort to the following particular scaling of space and time about $\epsilon$ as:

(3.2)\begin{equation} x^\prime={{\epsilon }^{{-}1}}\xi ,\quad t={{\epsilon }^{{-}2}}\tau. \end{equation}

This scale separation, known as the diffusive scaling (Pavliotis Reference Pavliotis2008, Chapter 13, p. 210), is anticipated for asymptotically long times, to yield the effective drift and diffusivity. The reason for different scales of perturbation in time and space is that this treatment homogenises the local spaces for the derivation of an effective diffusion equation, as could be placed into evidence promptly should we initially approximate an $\epsilon$ scale for both time and space before coming to the present scale separation successively. Now the rescaled probability distribution $P^\epsilon$ as spatio-temporal functions of $\xi$ and $\tau$ satisfies the equation

(3.3)\begin{equation} {{\epsilon }^{2}}\frac{\partial P^\epsilon}{\partial \tau }+\epsilon \left( {{Pe}}_f U_f+{{{{Pe}}}_{s}}\sin \theta -{{{\bar{U}}}} \right)\frac{\partial P^\epsilon}{\partial \xi }-{{\epsilon }^{2}}{{D}_{xx}}\frac{{{\partial }^{2}}P^\epsilon}{\partial {{\xi }^{2}}} -2\epsilon {{D}_{xz}}\frac{{{\partial }^{2}} P^\epsilon}{\partial \xi \partial z}+\mathcal{L} P^\epsilon=0. \end{equation}

Assuming that the probability distribution $P$ possesses an a posteriori unique solution, substitution of an approximate solution of multiple-scale expansions in the limit of $\epsilon \rightarrow 0$

(3.4)\begin{equation} P^\epsilon={{P}_{0}}+\epsilon {{P}_{1}}+{{\epsilon }^{2}}{{P}_{2}}+\textit{O}(\epsilon^3)\end{equation}

into (3.3) gives a hierarchy of equations, belonging to the framework of the classical homogenisation problem.

Let us collect the terms of order 1 as

(3.5)\begin{equation} \textit{O} \left( 1 \right): \mathcal{L}{{P}_{0}}=0. \end{equation}

Since the longitudinal distribution of the initial condition is a function independent of the small-scale variable $\xi$, it is intuitive to conclude that the leading-order behaviour should be a function only of $z$ and $\theta$. Next, we equate to zero the terms of the order of $\epsilon ^1$

(3.6)\begin{equation} \textit{O}( \epsilon ): \mathcal{L}{{P}_{1}}={-}( {{{{Pe}}}_{f}}U_f+{{{{Pe}}}_{s}}\sin \theta -{{{\bar{U}}}})\frac{\partial {{P}_{0}}}{\partial \xi }+2{{D}_{xz}}\frac{{{\partial }^{2}}{{P}_{0}}}{\partial \xi \partial z}. \end{equation}

A separation of variables for the first moment will then yield a cell problem of $b(z,\theta )$, which contains complete local information.

Although the coupling of components of the diffusivity tensor breaks the independence between the global and local spaces, the Onsager-like symmetry relations guarantee the mathematical tractability for the cross-diffusivities, i.e. the transposes of the global, local and coupling diffusion tensors equal themselves. For the order of $\epsilon ^2$, we have

(3.7)\begin{equation} \textit{O}( {{\epsilon }^{2}}): \mathcal{L}{{P}_{2}}={-}( {{Pe}}_f U_f+{{{{Pe}}}_{s}}\sin \theta -{{{\bar{U}}}})\frac{\partial {{P}_{1}}}{\partial \xi }+{{D}_{xx}}\frac{{{\partial }^{2}}{{P}_{0}}}{\partial {{\xi }^{2}}} +2{{D}_{xz}}\frac{{{\partial }^{2}}{{P}_{1}}}{\partial \xi \partial z}-\frac{\partial {{P}_{0}}}{\partial \tau }. \end{equation}

These equations will be solved recursively below.

3.2. Perturbations of zeroth order

For the equation of the order of 1, a transient separable solution is

(3.8)\begin{equation} {{P}_{0}}\left( \xi ,z,\theta ,\tau \right)=P_{0}^{\infty }\left( z,\theta \right)c\left( \xi ,\tau \right), \end{equation}

with $\mathcal {L} P_{0}^{\infty } =0$. That is, the historical effect is shown as $c( \xi,\tau )$. At the leading order, the probability density is similar to the classical theory (Brenner & Edwards Reference Brenner and Edwards1993; Jiang & Chen Reference Jiang and Chen2019) subject to isotropic diffusion, albeit with the lateral translational diffusivity opting for its local value.

Specifically, the governing equations of $P_0^\infty$ are

(3.9) \begin{equation} \left. \begin{gathered} \mathcal{L} P_0^{\infty}=0,\\ P_0^{\infty} \left(z,\theta\right)=P_0^{\infty} \left(z,{\rm \pi}-\theta\right), \quad \text{ at } z=0,1, \\ \frac{\partial P_0^{\infty}}{\partial z} \left(z,\theta\right)={-}\frac{\partial P_0^{\infty}}{\partial z} \left(z,{\rm \pi}-\theta\right), \quad \text{ at } z=0,1, \\ \left.P_0^{\infty}\right|_{\theta=0}=\left.P_0^{\infty}\right|_{\theta=2 {\rm \pi}}, \\ \left.\frac{\partial P_0^{\infty}}{\partial \theta}\right|_{\theta=0}=\left.\frac{\partial P_0^{\infty}}{\partial \theta}\right|_{\theta=2 {\rm \pi}}. \end{gathered} \right\} \end{equation}

The asymptotic zeroth moment $P_0^{\infty }$ can be determined explicitly by solving the above equation set. Benchmark solutions of $P_0^{\infty }$ will be shown for some special cases shortly. When it comes to more general cases with particle shape anisotropy and a coupling diffusivity tensor, the zeroth moment can be derived as a series expansion of the eigenfunctions. A Galerkin solution is pursued for $P_0^{\infty }$ as

(3.10)\begin{align} P_0^{\infty}(z,\theta) &=\sum_{i=1}^{\infty} q_i e_i(z,\theta) \nonumber\\ &\equiv \sum_{i=0}^{\infty }\sum_{j=1}^{\infty } \left[ A_{ij} \cos (i{\rm \pi} z)+B_{ij} \cos (i{\rm \pi} z)\cos (j\theta)+ C_{ij} \sin (i{\rm \pi} z)\sin (j\theta)\right], \end{align}

in which $A_{ij}$, $B_{ij}$ and $C_{ij}$ are the coefficients to be obtained by a group of linear equations with the orthogonality of trigonometric functions, and $q_i$ is the abstract coefficient vector. The analytical approximation would be truncated at a finite number, for the computational effectiveness of $P_0^\infty$. Due to the potential non-self-adjoint nature of $\mathcal {L}$, a weak formulation is given in the form of an inner product. By computing in advance the inner product of $e_i$ and $\mathcal {L}e_j$ in a bi-linear form, the problem turns into solving the null space of the local operator. That is, the coefficients are determined by the obtained null vector with an arbitrary constant. The normalisation condition

(3.11)\begin{equation} \langle P_0^\infty \rangle =1, \end{equation}

will ascertain the coefficients to be unique, wherein the averaging operator is defined as

(3.12)\begin{equation} \langle \boldsymbol{\cdot} \rangle \triangleq \int_{0}^{1} \int_{0}^{2{\rm \pi}} \boldsymbol{\cdot} \mathrm{d}\theta\,\mathrm{d}z. \end{equation}

3.3. Perturbations of first order

At $\textit {O}(\epsilon )$, the no-flux boundary condition becomes

(3.13)\begin{equation} \int_0^{2{\rm \pi}} \left[D_{xz}(\theta)\frac{\partial P_0}{\partial \xi} +D_{zz}(\theta)\frac{\partial P_1}{\partial z}-{{Pe}}_s \cos \theta P_1 \right] \,\mathrm{d}\theta=0,\quad \text{at } z=0,1. \end{equation}

With (3.6), (3.13) and corresponding periodic boundary conditions, a centring condition is

(3.14) \begin{align} &\left\langle -\left[ P{{e}_{f}}U_f+{{{{Pe}}}_{s}}\sin \theta -{{{\bar{U}}}} \right] \frac{\partial P_{0}}{\partial \xi}+2{{D}_{xz}}\frac{\partial^2 P_{0}}{\partial \xi \partial z} \right\rangle \nonumber\\ &\quad =\langle \mathcal{L} P_1\rangle\nonumber\\ &\quad= \left\langle {{Pe}}_s \cos \theta \frac{\partial P_1}{\partial z} +\frac{\partial \left[ \varOmega(z, \theta) P_1 \right]}{\partial \theta} -2D_{xz}(\theta)\frac{\partial^2 P_0}{\partial \xi \partial z} -D_{zz}(\theta)\frac{\partial^2 P_1}{\partial z^2} -\frac{\partial^2 P_1}{\partial \theta^2} \right\rangle\nonumber\\ &\quad ={-}\int_{0}^{1} \,\mathrm{d}z \int_{0}^{2{\rm \pi}} \frac{\partial}{\partial z} \left[2D_{xz}(\theta)\frac{\partial P_0}{\partial \xi} +D_{zz}(\theta)\frac{\partial P_1}{\partial z}-{{Pe}}_s \cos \theta P_1 \right] \,\mathrm{d}\theta\nonumber\\ &\qquad+\int_{0}^{1} \mathrm{d}z \left[\varOmega(z, \theta) P_1- \frac{\partial P_1}{\partial \theta}\right]_0^{2{\rm \pi}} =0. \end{align}

It follows that the convective term then averages to zero in an appropriate sense, as demonstrated in the next subsection. In other words, the effective behaviour of $P_0$ is that of pure diffusion. Explicitly, with the normalisation condition, we obtain the drift as

(3.15)\begin{equation} {{\bar{U}}}=\left\langle P_{0}^{\infty }\left( {{{{Pe}}}_{f}}U_f+{{{Pe}}_{s}}\sin \theta \right)-2{{D}_{xz}}\frac{\partial P_{0}^{\infty }}{\partial z} \right\rangle. \end{equation}

An a posteriori solution for $P_1$ is available with the centring condition as

(3.16)\begin{equation} {{P}_{1}}\left( \xi ,z,\theta ,\tau \right)={-}b\left( z,\theta \right)\frac{\partial c}{\partial \xi }+P_{0}^{\infty }\left( z,\theta \right)f\left( \xi ,\tau \right), \end{equation}

wherein $b(z,\theta )$ is governed by

(3.17) \begin{equation} \left. \begin{gathered} \mathcal{L}b=P_{0}^{\infty }\left( {{{Pe}}_{f}}U+{{{Pe}}_{s}}\sin \theta -{{{\bar{U}}}} \right)-2{{D}_{xz}}\frac{\partial P_{0}^{\infty }}{\partial z},\\ \left.b\right|_{\theta=\theta_0} =\left.b \right|_{\theta={\rm \pi}-\theta_0}, \quad \text{at } z=0,1,\\ \left. \frac{\partial b}{\partial z}\right|_{\theta=\theta_0} ={-}\left. \frac{\partial b}{\partial z}\right|_{\theta={\rm \pi}-\theta_0},\quad \text{at } z=0,1,\\ \left.b\right|_{\theta=0}=\left.b\right|_{\theta=2 {\rm \pi}}, \\ \left.\frac{\partial b}{\partial \theta}\right|_{\theta=0}=\left.\frac{\partial b}{\partial \theta}\right|_{\theta=2 {\rm \pi}}. \end{gathered} \right\} \end{equation}

The boundary conditions are similar to those of $P_{0}^{\infty }$. Note that $\mathcal {L} [P_{0}( z,\theta )f( \xi,\tau )]=0$, so the arbitrary function $f(\xi,\tau )$ does not indeed enter into the moment equations.

It is obvious from (3.9) that $P_0^\infty$ multiplied by an arbitrary constant constitutes a complementary solution of $b$. Thus, we introduce a normalisation condition as

(3.18)\begin{equation} \langle b_N\rangle =0. \end{equation}

We then devise a decomposition of $b$ as

(3.19)\begin{equation} b(z, \theta)=b_N(z, \theta)+\bar{B} P_0^{\infty}(z, \theta), \end{equation}

where the constant $\bar {B}$ is the integration of $b$ over the cross-section, representing the initial information of probability. Indeed, it is the gradient of $B$ that is unique rather than $B$ itself since it has an additive constant $\bar {B}$, as shown in (3.19). Nevertheless, only $\boldsymbol {\nabla } B$, or alternatively $b_N$ with the normalisation condition (Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003), other than $B$ enters into the computation of dispersivity. Solutions of long-time phenomenological transport coefficients can be derived with $b_N$ solely. Likewise, we will pursue a Galerkin solution of $b_N$ as performed to $P_0^\infty$. The difference lies in that the right-hand side becomes an inhomogeneous source term.

Next of interest to us is utilising $b$, which contains complete local information, to determine the analytical solution of the dispersivity directly.

3.4. Perturbations of second order

With the reflective and periodic boundary conditions, we have

(3.20)\begin{align} \left\langle \mathcal{L} P_2 \right\rangle &=\left\langle -\left[ {{Pe}}_f U_f(z)+{{{{Pe}}}_{s}}\sin \theta -{{{\bar{U}}}} \right]\frac{\partial {{P}_{1}}}{\partial \xi }+{{D}_{xx}}(\theta )\frac{{{\partial }^{2}}{{P}_{0}}}{\partial {{\xi }^{2}}} +2{{D}_{xz}}(\theta )\frac{{{\partial }^{2}}{{P}_{1}}}{\partial \xi \partial z}-\frac{\partial {{P}_{0}}}{\partial \tau } \right\rangle \nonumber\\ &= \left\langle {{Pe}}_s \cos \theta \frac{\partial P_2}{\partial z} +\frac{\partial \left[ \varOmega(z, \theta) P_2 \right]}{\partial \theta} -D_{xz}(\theta)\frac{\partial^2 P_1}{\partial \xi \partial z} -D_{zz}(\theta)\frac{\partial^2 P_2}{\partial z^2} -\frac{\partial^2 P_2}{\partial \theta^2} \right\rangle\nonumber\\ &={-}\int_{0}^{1} \,\mathrm{d}z \int_{0}^{2{\rm \pi}} \frac{\partial}{\partial z} \left[2D_{xz}(\theta)\frac{\partial P_1}{\partial \xi} +D_{zz}(\theta)\frac{\partial P_2}{\partial z}-{{Pe}}_s \cos \theta P_2 \right] \mathrm{d}\theta\nonumber\\ &\quad +\int_{0}^{1} \,\mathrm{d}z \left[\varOmega(z, \theta) P_2- \frac{\partial P_2}{\partial \theta}\right]_0^{2{\rm \pi}} =0. \end{align}

For the equation of the order of $\epsilon ^2$, substitution of solutions of $P_0$ (3.8) and $P_1$ (3.16) into (3.7) gives

(3.21)$$\begin{gather} \mathcal{L}{{P}_{2}}=\left( P{{e}_{f}}U_f+{{{{Pe}}}_{s}}\sin \theta -{{{\bar{U}}}} \right)b\frac{{{\partial }^{2}}c}{\partial {{\xi }^{2}}}-\left( P{{e}_{f}}U_f+{{{{Pe}}}_{s}}\sin \theta -{{{\bar{U}}}} \right) P_{0}^{\infty }\left( z,\theta \right)\frac{\partial f}{\partial \xi } \nonumber\\ +\,{{D}_{xx}}P_{0}^{\infty }\frac{{{\partial }^{2}}c}{\partial {{\xi }^{2}}}-P_{0}^{\infty }\frac{\partial c}{\partial \tau }-2{{D}_{xz}}\frac{\partial b}{\partial z}\frac{{{\partial }^{2}}c}{\partial {{\xi }^{2}}}+2{{D}_{xz}}\frac{\partial P_{0}^{\infty }}{\partial z}\frac{\partial f}{\partial \xi }. \end{gather}$$

With (3.7) and (3.20), we obtain a second centring condition

(3.22)$$\begin{gather} \left\langle \left( {{{{Pe}}}_{f}}U_f+{{{{Pe}}}_{s}}\sin \theta -{{{\bar{U}}}} \right) b\frac{{{\partial }^{2}}c}{\partial {{\xi }^{2}}}-\left( {{{{Pe}}}_{f}}U_f+{{{{Pe}}}_{s}}\sin \theta -{{{\bar{U}}}} \right) P_{0}^{\infty }\frac{\partial f}{\partial \xi } \right. \nonumber\\ \left. \qquad +\,{{D}_{xx}}P_{0}^{\infty }\frac{{{\partial }^{2}}c}{\partial {{\xi }^{2}}}-P_{0}^{\infty }\frac{\partial c}{\partial \tau }-2{{D}_{xz}}\frac{\partial b}{\partial z}\frac{{{\partial }^{2}}c}{\partial {{\xi }^{2}}}+2{{D}_{xz}}\frac{\partial P_{0}^{\infty }}{\partial z}\frac{\partial f}{\partial \xi } \right\rangle =0. \end{gather}$$

It could be seen from (3.16) that $f(\xi,\tau )$ does not contribute to $P_2$ with the help of the first centring condition (3.14). That is, the governing equation of the undetermined function $c(\xi,\tau )$ is

(3.23)\begin{equation} \frac{\partial c}{\partial \tau }-\left\langle \left( {{{{Pe}}}_{f}}U_f+{{{{Pe}}}_{s}}\sin \theta -{{{\bar{U}}}} \right) b+{{D}_{xx}}P_{0}-2{{D}_{xz}}\frac{\partial b}{\partial z} \right\rangle \frac{{{\partial }^{2}}c}{\partial {{\xi }^{2}}}=0. \end{equation}

Note that the normalisation condition of $P_0$ has been utilised in (3.23). As the notation suggests, the migration of micro-swimmers subject to anisotropic diffusion in a confined channel is governed by an effective dispersion equation with the convective term averaged out. It implies evidently that the dispersivity reads

(3.24)\begin{equation} {{\bar{D}}}=\left\langle \left( {{{{Pe}}}_{f}}U_f+{{{{Pe}}}_{s}}\sin \theta -{\bar{U}} \right) b+{{D}_{xx}}P_{0}^{\infty }-2{{D}_{xz}}\frac{\partial b}{\partial z} \right\rangle. \end{equation}

Consequently, the enhanced diffusivity scales with the swim and flow Péclet numbers, and the molecular diffusivity is modified with an addition of the anisotropic effect. For $0 < \epsilon \ll 1$ and long times, the solution $P^\epsilon$ of (3.3) is approximated by $P_0$, a solution of the convection–diffusion equation with slowly varying initial conditions

(3.25)$$\begin{gather} \frac{\partial P_0}{\partial t } +\bar{U} \frac{\partial P_0}{\partial x}= \bar{D} \frac{{{\partial }^{2}}P_0}{\partial {{x}^{2}}}, \end{gather}$$

(3.26)$$\begin{gather}\left. P_0\right|_{t=0} = g(x) P^{(0)}(z,\theta). \end{gather}$$

The important point is that the parabolic equations, cf. (3.9) and (3.17), exhibit independence from the small-scale $\epsilon$. In certain cases, explicit solutions for the two partial differential equations are attainable, as will be shown in § 4.1. Even when this is not the case, the equations lend themselves to rigorous analysis or efficient numerical simulations. Due to the slowly varying coefficients, this is far less computationally expensive than the direct numerical simulation of (2.1). Specifically, (3.15) represents the mean drift velocity of an active Brownian particle, and (3.24) is the dispersivity as a sum of the convective Taylor contribution and the molecular contribution. In summary, the asymptotic phenomenological coefficients are derived precisely with a homogenisation method.

Orientation distribution could be obtained by the averaged zeroth moment over the lateral cross-section as $C^\theta =\int _{0}^{1} P_0^{\infty } (z, \theta )\, \mathrm {d}z$. The vertical flux is the product of normalised vertical polarisation $C^z= \int _{0}^{2{\rm \pi} } P_0^{\infty } (z, \theta ) \,\mathrm {d} \theta$ and vertical concentration $P^z= \int _{0}^{2{\rm \pi} } P_0^{\infty } \cos (\theta )(z, \theta )\, \mathrm {d} \theta / C^z$, with swim Péclet numbers as a coefficient, as

(3.27)\begin{equation} J_z(z)={{Pe}}_s P^z(z) C^z(z)=\int_{0}^{2{\rm \pi}} P_0^{\infty} (z, \theta) {{Pe}}_s\cos \theta\, \mathrm{d} \theta . \end{equation}

Integrating (3.5) from $0$ to $1$ with respect to $\theta$ gives

(3.28)\begin{equation} {{Pe}}_s \cos \theta \frac{\partial C^z}{\partial z} -D_{zz} \frac{\partial^2 C^z}{\partial z^2}=0. \end{equation}

That is, the vertical concentration is governed by the balance between lateral diffusion and convection. Hence, the ensemble-averaged effective vertical migration could be quantified by the ratio of vertical flux to vertical concentration, i.e. normalised vertical polarisation, as

(3.29)\begin{equation} V_m^\infty(z)=\frac{J_z}{C^z}={{Pe}}_s P^z(z). \end{equation}

For further research interest, the transient evolution of these transport coefficients can feature the temporal process in the anisotropic diffusion of active ellipsoids long before the Taylor dispersion regime (Guan et al. Reference Guan, Zeng, Li, Guo, Wu and Wang2021; Debnath et al. Reference Debnath, Jiang, Guan and Chen2022; Wang, Jiang & Chen Reference Wang, Jiang and Chen2022a,  Reference Wang, Jiang and Chen2023). To reflect the transient evolution of $\bar {U}$ and $\bar {D}$, we denote the transient drift velocity and dispersivity as $U_x(t)$ and $D_T(t)$, respectively. Note that we should not take it for granted to write the streamwise drift velocity and dispersivity (Yasuda Reference Yasuda1984; Guan et al. Reference Guan, Zeng, Jiang, Guo, Wang, Wu, Li and Chen2022; Guan & Chen Reference Guan and Chen2024), by reducing the external integral of (3.15) and (3.24) directly. Instead, we should recover a streamwise definition from the time-dependent solutions of moments successively.

4.1. Benchmark solutions without translational diffusion and particle shape anisotropy

First, we resort to a simple benchmark case of spheroidal ($\alpha _0=0$) gyrotactic micro-swimmers without translational diffusion ($D_{\|}=D_\perp =0$). The steady distribution is

(4.2)\begin{equation} \mathcal{L}P_{0}^{\infty }={{Pe}}_s \cos \theta \frac{\partial P_{0}^{\infty }}{\partial z}+\frac{{{Pe}}_f}{2}\frac{\mathrm{d}U}{\mathrm{d}z}\frac{\partial P_{0}^{\infty }}{\partial \theta }-\lambda \frac{\partial }{\partial \theta }\left( P_{0}^{\infty }\sin \theta \right)-\frac{{{\partial }^{2}}P_{0}^{\infty }}{\partial {{\theta }^{2}}}=0. \end{equation}

With an a posteriori condition that the rotational diffusion results in a uniform angular distribution asymptotically, the governing equation of $P_0^\infty$ is independent of $\theta$. In this way, (4.2) reduces to

(4.3)\begin{equation} {{Pe}}_s \cos \theta \frac{\mathrm{d} P_{0}^{\infty }}{\mathrm{d} z}-\lambda \frac{\mathrm{d} }{\mathrm{d} \theta } \left[ P_{0}^{\infty }(z) \sin \theta \right]=0. \end{equation}

Thus, an analytical solution is

(4.4)\begin{equation} P_{0}^{\infty }=\frac{\lambda }{2{\rm \pi} {{Pe}}_s \left[ \exp \left( \lambda /{{{Pe}}_s} \right)-1 \right]}\exp \left( \frac{\lambda z}{{{Pe}}_s} \right). \end{equation}

Note that the benchmark solution (4.4) satisfies the reflective boundary conditions. We are aware that an identical formalism of this benchmark has also been proposed previously (Bearon et al. Reference Bearon, Hazel and Thorn2011; Vennamneni et al. Reference Vennamneni, Nambiar and Subramanian2020). For distributions of non-gyrotactic micro-swimmers ($\lambda =0$), (4.4) asymptotes to uniformity as

(4.5)\begin{equation} P_{0}^{\infty }=\frac{1} {2{\rm \pi}}. \end{equation}

Subsequently, the steady drift could be derived as

(4.6)\begin{equation} {{{U}_s^\infty}}= P_{0}^{\infty }\left( {{Pe}}_f U+{{{Pe}}_{s}}\sin \theta \right) ={\frac{\lambda {{{{Pe}}}_{f}}U_f(z){{\exp}\left({\lambda z/{{{Pe}}_s}}\right)}}{{{{{Pe}}}_{s}}\left( {\exp \left({\lambda /{{Pe}}_s}\right)}-1 \right)}}. \end{equation}

With the vertical migration velocity (3.29), it is interesting to see whether gyrotactic micro-swimmers tend to swim vertically within a horizontal flow in the presence of reorientation. The steady vertical polarisation equals zero asymptotically as

(4.7)\begin{equation} V_m^\infty=0. \end{equation}

This may seem simple, but is in fact far reaching in that symmetric actuation leads to no ensemble-averaged vertical movement at a population level. Conversely, imposed symmetry breaking would bring some exceptional deviations from the prediction. That is, non-spherical shapes, external body forces and direct coupling effects in a noisy environment can all be exploited to generate effective locomotion of micro-swimmers.

4.2. Spatial patterns: effects of anisotropic diffusion at leading order

Gyrotaxis, which arises from the combined effects of gravity and shear, can induce an overall upward tendency of the micro-swimmers. Figure 2 illustrates the influence of anisotropic diffusion on steady distributions of the spatial and orientational probability. The exerted driving force is as weak as $\lambda = 0.5$, and the parallel translational diffusivity $D_{\|}$ is tuned with $D_{\perp }$ fixed as zero as a hallmark of anisotropy on the vertical diffusivity $D_{zz}$, since the streamwise diffusivity and cross-diffusivity do not factor into $P_0^\infty$. For strong translational diffusion, the concentration distribution approaches (yet will never reach) uniformity. Concurrently, the orientational distribution exhibits a preference for upward migration since gravity is sufficiently strong. In the limit of weak translational diffusion, the orientational distribution becomes fairly uniform while the concentration distribution shows peculiar accumulation at walls. With moderate diffusivity ($10^{-5} \le D_{\|}\le 10^{-3}$), the concentration distribution exhibits wall accumulation due to gyrotaxis, as shown in figure 2(a). The peak of $C^\theta$ occurs around $3{\rm \pi} /2$ in figure 2(b), signifying a prevalence of upstream swimming. The amalgamation of these observations suggests that micro-swimmers rotate to swim against the flow near walls where the flow is relatively weak.

Figure 2.

Steady distributions of concentration and orientation subject to gyrotaxis. Common parameters are ${{Pe}}_s=0.1$, ${{Pe}}_f=10$, $\alpha _0=0$, $\lambda =0.5$ and $D_{\perp }=0$.

Figure 3 illustrates an internal symmetry breaking originating from particle shape anisotropy, which in turn leads to the generation of steady non-zero vertical locomotion from the probability perspective. Using (2.3) with $\lambda =0$, the angular velocity of micro-swimmers is governed by the local flow vorticity and the shape-induced strain rate. For highly elongated particles ($\alpha _0=0.9$) swimming slowly (${{Pe}}_s=0.1$) through the channel, it becomes evident that micro-swimmers tend to align themselves in a manner that amplifies effectively the impact of the strain rate, thereby sampling a position preferentially to compete against the pronounced fluid rotation. Specifically, the contribution of shape-induced rotation is maximally negative when the orientation angle $\theta$ assumes values of either ${\rm \pi} /2$ or $3{\rm \pi} /2$, in accordance with the orientation depicted in figure 3(b).

Figure 3.

Steady distributions of concentration and orientation subject to a shape-induced strain rate. Common parameters are ${{Pe}}_s=0.1$, ${{Pe}}_f=10$, $\alpha _0=1$ and $\lambda =0$.

Various patterns of concentration and orientation distributions occur with different anisotropic diffusivities. We direct our attention to the influence of anisotropic diffusion in a comparably strong flow field (${{Pe}}_f=10$). In figure 3(a), for a large perpendicular diffusivity, swimmers tend to form a symmetric cup-like concentration focused around the centreline. Conversely, pairs of depletion layers emerge, displaced symmetrically about $z=0.5$, with their maxima diminishing in amplitude while converging towards the wall as $D_{\perp }$ increases. A subtle transition from downstream to upstream swimming becomes noticeable in figure 3(b) from a probabilistic standpoint in response to stronger translational diffusion. In contrast to the gyrotactic cases discussed previously, the pronounced external flow reduces significantly the wall accumulation for elongated micro-swimmers. This effect arises from the rapid rotation imposed by the strong external flow on the particles near the wall, causing them to swim inward and compelling the elongated particles within the bulk to preferentially align themselves with the flow direction. The cup-like focusing pattern results from a delicate balance between shear alignment and diffusion. Thus a slight alteration in the diffusivity tensor can readily transition it into a bi-modal trapping pattern.

At finite translational diffusivities, we analyse the effect of strong isotropic diffusion with various ellipsoidal shapes and gyrotaxis strengths. In figure 4, the funnel-shaped concentration distribution is featured for non-gyrotactic elongated micro-swimmers. Recall that the perfectly spherical active particles would sample a uniform distribution in the $z$-space and exhibit no net vertical polarisation. The more the ellipsoids are elongated, the more enhanced the non-uniformity becomes. Both ends of the funnel-shaped distributions in figure 4(a) are anticipated as a result of the polarisation. Since the polarisation distribution is anti-symmetric about $z=0.5$, we consider only the upper half-plane. Due to shear-induced migration, the ellipsoidal micro-swimmers near the centreline tend to migrate upwards, i.e. towards high-shear regions. When micro-swimmers approach the upper wall, the reflective condition turns the swimming direction opposite, as indicated in figure 4(b). This observation for the reflective boundary conditions deviates from the high-shear trapping, as reported by Bearon & Hazel (Reference Bearon and Hazel2015) with extreme wall accumulation for a Robin boundary condition and by Ezhilan & Saintillan (Reference Ezhilan and Saintillan2015) with singular cusped profiles for a periodic double-Poiseuille boundary condition. The nomenclature of high-shear trapping with regard to this typical kind of steady concentration distribution, once deemed suitable for delineating this phenomenon, now appears less appropriate. More precisely, this constitutes a distinctive transport process moulded by the interplay between shape- and shear-induced migration, under the reflective boundary conditions. For spherical gyrotactic particles subject to strong isotropic diffusion, the deviation of the concentration distributions from uniformity in figure 4(c) and upswimming behaviour in figure 4(d) are remarkable and significantly enhanced with the gyrotaxis strengths. It is interesting to investigate the balance between gyrotactic upswimming and self-propulsion. With small $\lambda$, we have shown in figure 2 that the orientational distribution peaks in the upward direction. The swimmers migrate upwards rapidly near the centreline while exhibiting almost no net vertical locomotion near the wall. As gyrotaxis grows, a steady concentration layer of micro-swimmers forms with appreciable thickness due to the wall reflection and shear trapping. The thickness of this accumulation layer declines rapidly when gyrotaxis or diffusion dominates self-propulsion.

Figure 4.

Steady distributions of concentration and vertical migration velocity with different particle shape anisotropies and gyrotaxis: (a) vertical concentration distribution and (b) polarisation distribution for non-gyrotactic micro-swimmers ($\lambda =0$); (c) vertical concentration distribution and (d) polarisation distribution for spherical micro-swimmers ($\alpha _0=0$). Common parameters are ${{Pe}}_s=0.1$, ${{Pe}}_f=5$, $D_{\|}=10^{-2}$ and $D_{\perp }=10^{-2}$.

Using (3.5), the change of orientation originates from three parts: (a) relative strength of the lateral translational diffusivity to the rotational diffusivity, (b) angular velocity under the external effect of flow orientation and gravity field and (c) reflection when bouncing at walls. As illustrated in figure 5, the micro-swimmers exhibit three distinct equilibrium orientations. A possible first state could be peaking at approximately $3{\rm \pi} /2$, indicative of horizontal swimming against the prevailing flow; the second displays a pronounced tendency to sample inclined angles relative to the upright direction; the third orientation represents a transitional mode positioned between the above two states. When strong rotational diffusion dominates, the equilibrium balance is dominated by gyrotaxis and self-propulsion. Steady concentration layers emerge as a consequence of gyrotactic trapping when gyrotaxis exhibits relative weakness. Conversely, for sufficiently robust gravitational forces, micro-swimmers accumulate invariably towards the wall, disregarding the wall reflection, and align themselves with the upstream direction, as illustrated in figure 5(a,b). Note that the transition from low to high gyrotaxis signifies a distinct concentration distribution, in which a persistent characteristic peak becomes conspicuously pronounced, surpassing the nearby average values on both sides of the wall. This observation aligns with the upswimming behaviour depicted in figure 5(c), wherein the trapping layers correspond to the negligible vertical migration velocity in the wall regions, and wall accumulations correspond to the approach of the peak of polarisation distribution towards $z=1$.

Figure 5.

Steady distributions with different translational diffusivities and gyrotaxis. Common parameters are ${{Pe}}_s=0.1$, ${{Pe}}_f=10$ and $\alpha _0=0$.

The magnitude of anisotropic diffusion impacts the migration of micro-swimmers only through the alteration of $D_{zz}$. That is, $D_{\|}$ is dominant for the upswimming mode of orientation, while $D_{\perp }$ becomes advantageous for the upstream swimming mode. For large gyrotaxis ($\lambda =10$), micro-swimmers accumulate remarkably near the wall while displaying a pronounced preference for alignment against the flow direction. In this case, the existence of parallel diffusivity works indeed as an amplifier of the overall vertical translational diffusivity, thereby smoothing the sharp accumulation near the wall. Under moderate gyrotaxis conditions ($\lambda =5$), the balance mechanism becomes more intricate with the potential gyrotactic trapping or wall accumulation. Here, translational diffusion is governed primarily by $D_{\|}$, and anisotropic diffusion emerges into effect when $D_{\|} \neq D_{\perp }$. With an increase in the magnitude of the overall $D_{zz}$, the concentration distribution becomes more flattened, transitioning from wall accumulation to gyrotactic trapping if gyrotaxis is fixed. During the approach to uniformity, a transitional characteristic peak is observed within the range of moderate values for $D_{zz}$, e.g. the red line for $10^{-3}< D_{zz}<10^{-2}$ in figure 5(a).

4.3. Diversity in concentration distributions: effects of flow convection and self-propulsion

Figure 6 shows the effect of flow strength on steady distributions of concentration and orientation with weak self-propulsion and isotropic translational diffusion. Diversity of steady concentration distributions is revealed as a result of different shear strengths in figure 6(a). When the flow speed is comparable to the swim speed, micro-swimmers exhibit a flattened nearly uniform distribution whereas tiny accumulations at the wall and slight focusing near the centreline are observed. As shown in figure 6(b), the orientation distributions with different flow Péclet numbers are qualitatively consistent, aligning preferentially with the flow direction either along the upstream or downstream direction. With weak ${{Pe}}_f$, this non-uniformity of orientation is impaired remarkably. Interestingly, the upstream tendency is sampled a bit more from the probability perspective, consistent with the mechanism of upstream rheotaxis for wall accumulations similar to figure 5.

Figure 6.

Steady distributions of concentration and orientation as functions of flow Péclet numbers. Common parameters are ${{Pe}}_s=0.1$, $D_{\|}=10^{-4}$, $D_{\perp }=10^{-4}$, $\lambda =0$ and $\alpha _0=0.9$.

When the flow strength grows, shear trapping appears as a peculiar phenomenon featuring symmetric depletion layers about the centreline. For periodic boundary conditions, Vennamneni et al. (Reference Vennamneni, Nambiar and Subramanian2020) have discussed intensively and rationalised the trapping phenomenon. In the present work, due to the imposed boundary conditions, trapping would be more complicated and we define three potential patterns as centreline focusing, wall accumulation and shear trapping. Note that the bi-modal distribution of depletion layers is a hallmark of shear trapping, trapped either in the low- or high-shear regions. Of interest for us to elucidate are distributions and mechanisms in this new system with effective boundary conditions. As demonstrated in figure 3(a), the shear-induced migration of elongated micro-swimmers could form a cup-shaped focus at the centre for strong translational diffusion. In contrast, funnel-shaped shear trapping with reduced wall accumulation is found for low translational diffusivities and great flow strengths. As illustrated in figure 6(a), the trapped depth converges towards the centreline with increasing flow Péclet number. That is, the strong flow rotates rapidly the particles with the largest shear rates at $z=0 \textrm { and } 1$, thereby depleting the near-wall micro-swimmers and rushing them inward. On the other hand, the micro-swimmers on the centreline migrate with an inclined angle to the flow direction, only stable when swimming away from the centreline since they would otherwise experience an anti-symmetric flow vorticity. Combining these two mechanisms, the bi-modal trapping distribution is formed. Note that the orientation still peaks around ${\rm \pi} /2$ and $3{\rm \pi} /2$, exhibiting upstream rheotaxis and shear alignment in figure 6(b), whereas the width of the Gaussian-like distribution decreases with ${{Pe}}_f$.

Steady distributions of concentration and orientation for moderate values of ${{Pe}}_f$ as functions of swim Péclet numbers are presented in figure 7. With small ${{Pe}}_s$, the micro-swimmers become trapped near the central axis and walls, resulting in a symmetrical depletion layer around $z=0.5$. Within the near-wall domains, a balance is achieved between stochastic swimming and wall reflection at a fixed depth, thereby trapping particles preferentially aligned with the flow. In the central region, the mechanism of shear-induced alignment contributes to various sampling modes of micro-swimmers. Note that the trapping depth converges towards the boundary plates as self-propulsion intensifies. As ${{Pe}}_s$ increases further, the trapping transitions into a pair of uni-modal distributions symmetric about $z=0.5$. When self-propulsion is as strong as or even surpasses flow convection, noteworthy wall accumulations and subtle centreline focusing phenomena are observed. With increasing swim Péclet numbers, the characteristic alignment of ellipsoidal micro-swimmers is compressed considerably, resulting in a nearly uniform distribution in both the position and orientation spaces. That is, the strength of shear-induced trapping is influenced significantly by the swim Péclet numbers, with its highest manifestation occurring at intermediate levels of activity (Rusconi et al. Reference Rusconi, Guasto and Stocker2014). However, the active swimming could not alter qualitatively the steady distributions, yet eliminate quantitatively the deviations between the near-wall regions and the central area.

Figure 7.

Steady distributions of concentration and orientation as functions of the swim Péclet numbers. Common parameters are ${{Pe}}_f=10$, $D_{\|}=10^{-4}$, $D_{\perp }=10^{-4}$, $\lambda =0$ and $\alpha _0=0.9$.