Migration of confined micro-swimmers subject to anisotropic diffusion

Abstract

Shear-induced migration of elongated micro-swimmers exhibiting anisotropic Brownian diffusion at a population scale is investigated analytically in this work. We analyse the steady motion of confined ellipsoidal micro-swimmers subject to coupled diffusion in a general setting within a continuum homogenisation framework, as an extension of existing studies on macro-transport processes, by allowing for the direct coupling of convection and diffusion in local and global spaces. The analytical solutions are validated successfully by comparison with numerical results from Monte Carlo simulations. Subsequently, we demonstrate from the probability perspective that symmetric actuation does not yield net vertical polarisation in a horizontal flow, unless non-spherical shapes, external fields or direct coupling effects are harnessed to generate steady locomotion. Coupled diffusivities modify remarkably the drift velocity and vertical migration of motile micro-swimmers exposed to fluid shear. The interplay between stochastic swimming and preferential alignment could explain the diverse concentration and orientation distributions, including rheological formations of depletion layers, centreline focusing and surface accumulation. Results of the analytical study shed light on unravelling peculiar self-propulsion strategies and dispersion dynamics in active-matter systems, with implications for various transport problems arising from the fluctuating shape, size and other external or inter-particle interactions of swimmers in confined environments.

1. Introduction

Migration of active matter is basic to our comprehension of diverse physiological environments, exemplified by phenomena like algal blooms (e.g. Durham & Stocker 2012). Specifically, the migration of confined micro-swimmers influences significantly environmental and bio-chemical processes, particularly in the context of dispersion and mixing (e.g. Cates & Tjhung 2018; Gouiller et al. 2021). Nevertheless, a comprehensive theoretical depiction of the fundamental properties of living matter remains unattainable due to its profound complexity. It is plausible, therefore, that overarching principles such as conservation laws and symmetries serve as constraints on the potential collective behaviour exhibited by cells (Marchetti et al. 2013).

The dynamics of micro-swimmers is dominated by viscous forces at low Reynolds numbers, and their self-propulsion can occur only through non-reciprocal fluid manipulation (Purcell 1997). From a kinematic perspective, it is feasible to establish a linear relationship from force and moment matrices towards the velocity and rotation matrices through a grand resistance matrix that encompasses inherently a direct interplay between positional convective translation and orientational rotation (Karrila & Kim 1991). This progress could reveal mechanical and bio-physical properties of living matter at the level of each cell, for locomotion on small scales (Lauga 2011). On the other hand, there is an essential need for a phenomenological Smoluchowski approach on field-based continuum modelling, regarding the non-equilibrium steady state of an active system as emerging from the introduction of non-zero, albeit small, driving forces.

The homogenisation technique has been employed in various applications to investigate the long-term asymptotic dispersion of passive particles (Mei, Auriault & Ng 1996; Pavliotis 2008). Dispersion processes in confined flows, with reactions occurring either in the bulk or on a boundary, were investigated extensively by Ng & Yip (Ng 2006a,b). Wu & Chen (2014) employed the homogenisation method to examine the initial stage and transverse concentration distributions. Despite progress in understanding the dispersion of passive particles, the homogenisation method has not been extended to investigate the dispersion of ellipsoidal micro-swimmers. In this work, an effective strategy involves categorising relaxation processes into fast and slow components. As vigorous applications and extensions, we formulate a homogenised theory that addresses primarily the slow dynamics, with the fast processes considered as sources of noise and damping.

Kumar et al. (2021) examined Taylor dispersion of passive ellipsoids with anisotropic diffusion at high flow Péclet numbers in a Poiseuille flow. Following the progress, asymptotic expressions of mean drift and dispersivity in the limit of large rotary Péclet numbers are derived by Khair (2022). In contrast to passive particles, the rotational dynamics turns out to be crucial, especially for the time-dependent transport of self-propelling micro-swimmers, because a slight deviation of the orientation could lead to markedly distinct distributions over brief time intervals. Anisotropic diffusion of ellipsoidal tracers in suspensions of active particles could display non-Gaussian statistics and dispersive phenomena (Nordanger, Morozov & Stenhammar 2022). Rusconi, Guasto & Stocker (2014) showed that trajectories of bacteria displayed frequent loops in high-shear regions due to the hydrodynamic torque generated by the local shear, corresponding to a shear-induced trapping effect in the high-shear domains. Assuming that the steady distribution mirrors a uniform distribution and translational diffusion is negligible, Bearon, Hazel & Thorn (2011) deduced a steady solution of the probability distribution for spheroidal gyrotactic micro-swimmers in a horizontal Poiseuille flow. Recently, Fung (2023) obtained analytical solutions for steady distributions of motile gyrotactic micro-organisms without translational diffusion in vertical flows. While the spatial and orientational dynamics of passive and active particles could be controlled under various mechanisms through external fields and/or fluctuations (Morris & Brady 1996; Saintillan & Shelley 2013; Hwang & Pedley 2014; Ishikawa & Pedley 2014; Morris 2020; Wang & Cirpka 2021; Wang et al. 2022c; Yang et al. 2022; Zeng, Jiang & Pedley 2022), theory has not yet been developed in a general setting to address the migration of confined micro-swimmers experiencing both anisotropic diffusion and self-propulsion.

The interactions between active matter and rigid plates within the confined channel significantly influence a range of biophysical and chemical processes (Shen et al. 2017; Traverso & Michelin 2022; Liao et al. 2023; Zhan, Jiang & Wu 2024). The inclination of active suspensions to accumulate near surfaces is prominent (Rothschild 1963). Hydrodynamic interaction, a potential but non-exclusive mechanism for migration, was addressed by Berke et al. (2008) with regard to Escherichia coli. Alternatively, Li & Tang (2009) proposed a kinematic-based mechanism that elucidates the accumulation phenomenon arising from bacterial collisions. This suggests that hydrodynamic interactions may play only a secondary role in the wall accumulation. The Smoluchowski equation could predict quantitatively the high-shear trapping phenomena (Bearon & Hazel 2015). Nevertheless, in the absence of translational diffusion, non-physical concentration singularities may occur at walls. With periodic boundary conditions in position space, Vennamneni, Nambiar & Subramanian (2020) concentrated on the trapping phenomena within low- and high-shear regions observed experimentally by Rusconi et al. (2014) and Barry et al. (2015), to exclude direct influences of channel plates.

The anisotropic diffusion of ellipsoidal particles is ubiquitous in nature, induced intrinsically by a non-spherical shape (Perrin 1936; Brenner 1967), as well as non-thermal fluctuations of an external force (Thiffeault & Guo 2022). Spatial fluctuations in viscosity give rise to modifications in both translation and rotation (Pedley 2010; Brumley et al. 2015; Chen, Perazzo & Stone 2020; Kamal & Lauga 2023). By abstracting the mechanism underlying speed alterations and reorientation, a field-based model could capture adeptly these dynamic phenomena. Shear-induced migration of active ellipsoids uncovers abundant counter-intuitive phenomena, e.g. low-shear trapping (Rusconi et al. 2014; Vennamneni et al. 2020), upstream rheotaxis (Hill et al. 2007; Kaya & Koser 2009; Omori et al. 2022) and centreline depletion (Rusconi et al. 2014; Ezhilan & Saintillan 2015), etc. Inspired by these recent experimental and theoretical investigations, we attempt to analyse at a population level the shear-induced migration of active ellipsoids subject to coupled diffusivities within a continuum theoretical framework. Intriguingly, how could built-in asymmetries and external driving forces modify the migration of confined micro-swimmers?

The paper is structured as follows. Section 2 presents the Smoluchowski analysis for confined micro-swimmers with anisotropic diffusion, followed in § 3 by asymptotic solutions for steady moments and macroscopic transport coefficients. Effects of anisotropic diffusion, particle shape anisotropy, external fields and self-propulsion on the obtained concentration and orientation distributions are considered in § 4, as well as the transient phenomenological transport coefficients in § 5. Finally, § 6 concludes. Solution procedures of time-dependent distributions can be found in Appendix A, and numerical validations through the Monte Carlo method are shown in Appendix B.

2. Problem formulation

We consider a dilute suspension of confined micro-swimmers in a steady parallel flow subject to anisotropic diffusion. Ellipsoidal micro-swimmers are exemplified with fluctuations in both position and orientation spaces, of which the origins could be athermal and biological. A set of two-dimensional Cartesian coordinates $\boldsymbol {R}=[x, z]^{\rm T}$ with base vectors $({\boldsymbol {e}}_{x}, {\boldsymbol {e}}_{z})$ is introduced, as shown in figure 1. A statistical mechanical description of a particle's motion is embodied in the probability distribution function $P(\boldsymbol {R},\theta, t)$ of finding the swimmer at the right position at time $t$, with $\theta$ standing for the angle of the orientation vector ${\boldsymbol {p}}$ from the positive $z$-axis. The transport problem adopts the form of the Smoluchowski equation in the position and orientation spaces, from the conservation law of mass as

(2.1)\begin{equation} \frac{\partial P}{\partial t}+ \boldsymbol{\nabla}_{\boldsymbol{R}} \boldsymbol{\cdot} \left[{{Pe}}_f {\boldsymbol{U}_f}(\boldsymbol{R}) P+{{Pe}}_s \boldsymbol{p}P- \boldsymbol{\mathsf{D}} \boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{R}} P\right]+ \frac{\partial}{\partial \theta}\left[\varOmega(z,\theta) P- \frac{\partial P}{\partial \theta}\right]= 0, \end{equation}

where $\boldsymbol {\nabla }_{\boldsymbol {R}}$ is the gradient operator in the position space, ${\boldsymbol {U}_f(\boldsymbol {R})}=[U_f(\boldsymbol {R}),0]^{\rm T}$ is the external velocity. We parameterise all length scales by the width $W^*$ of the channel, flow velocity $U_f^*$ by the mean speed $U_m^*$ over the lateral cross-section, swimming velocity by the swimming speed $V_s^*$ of active ellipsoids, the rate of orientational change $\varOmega ^*$ by the rotational diffusivity $D_\theta ^*$, and dimensional translational diffusivities $(D_{\|}^*$ and $D_{\perp }^*)$ along the parallel and perpendicular axes by $D_\theta ^* W^{* 2}$, with asterisk signs employed to represent the dimensional characteristic quantities. Specifically, the dimensionless parameters are

(2.2)\begin{equation} \left. \begin{gathered} t=t^* D_\theta^*, \quad x=\frac{x^*-U_m^* t^*}{W^*}, \quad z=\frac{z^*}{W^*}, \quad U_f=\frac{U_f^*}{U_m^*}-1, \quad \varOmega=\frac{\varOmega^*}{D_\theta^*},\\ {{Pe}}_s= \frac{V_s^*}{D_\theta^* W^*},\quad {{Pe}}_f=\frac{U_m^*}{D_\theta^* W^*}, \quad D_{\|}=\frac{D_{\|}^{*}}{D_\theta^*W^{*2}},\quad D_{{\perp}}=\frac{D_{{\perp}}^{*}}{D_\theta^*W^{*2}}, \end{gathered} \right\} \end{equation}

wherein $D_\|$ is the dimensionless translational diffusivity parallel to the swimming direction, $D_{\perp }$ is the dimensionless translational diffusivity perpendicular to the swimming direction, ${{Pe}}_s$ is the dimensionless swim Péclet number and ${{Pe}}_f$ is the dimensionless flow Péclet number. Provided inertial motion and particle–particle interaction are absent, the angular velocity could be determined through Jeffery's formulation (Jeffery & Filon 1922; Leal & Hinch 1972; Pedley & Kessler 1992; Lauga 2020), so that the rate of change of orientation is

(2.3)\begin{equation} \varOmega\left({z,\theta}\right)= \frac{{{Pe}}_f}{2}\frac{\mathrm{d}U}{\mathrm{d}z}\left[1+\alpha_0\cos(2\theta)\right]-\lambda\sin\theta, \end{equation}

where $\alpha _0= (AR^2-1)/(AR^2+1)$ is the shape factor characterised by the aspect ratio $AR$, $\lambda =1/2B^*D_\theta ^*$ is the gyrotactic bias parameter and $B^*$ is the reorientation time. Specifically, $\alpha _0=0$ stands for spheroids and $\alpha _0=1$ for infinitely thin rods. The elements of the translational diffusivity tensor are

(2.4)\begin{equation} \boldsymbol{\mathsf{D}}=\left[\begin{array}{cc} D_{x x}(\theta) & D_{x z}(\theta) \\ D_{z x}(\theta) & D_{z z}(\theta) \end{array}\right]=\left[\begin{array}{cc} D_{\|} \sin ^2 \theta+D_{{\perp}} \cos ^2 \theta & \left(D_{\|}-D_{{\perp}}\right) \sin \theta \cos \theta \\ \left(D_{\|}-D_{{\perp}}\right) \sin \theta \cos \theta & D_{\|} \cos ^2 \theta+D_{{\perp}} \sin ^2 \theta \end{array}\right] .\end{equation}

Figure 1. Ellipsoidal micro-swimmers suspended in a Poiseuille flow subject to anisotropic Brownian diffusion and an external gravity field.

In this context, the present work extends the existing theory concerning macro-transport processes (Jiang & Chen 2021; Wang, Jiang & Chen 2022b; Guan et al. 2023), facilitating the direct interplay of convection and diffusion in the local and global spaces. The transport of active ellipsoids, characterised by anisotropic diffusion, becomes intricate due to the influence of cross-diffusivities.

The no-flux boundary conditions require

(2.5)\begin{equation} \int_0^{2{\rm \pi}} {\boldsymbol{e}}_{z} \boldsymbol{\cdot} \left[ {{Pe}}_f {\boldsymbol{U}_f}(\boldsymbol{R})P+{{Pe}}_s \boldsymbol{p}P- \boldsymbol{\mathsf{D}} \boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{R}} P \right]\, \mathrm{d}\theta=0,\quad \text{at } z=0,1. \end{equation}

Provided that collisions between ellipsoids and channel boundaries are perfectly elastic, we further require in the position space that the incident probability fluxes be compensated by the reflected fluxes as (Bearon et al. 2011; Volpe, Gigan & Volpe 2014; Ezhilan & Saintillan 2015; Jakuszeit, Croze & Bell 2019)

(2.6)$$\begin{gather} P \left({\boldsymbol{R}},\theta, t\right)=P \left({\boldsymbol{R}},{\rm \pi}-\theta, t\right),\quad \text{at } z=0,1, \end{gather}$$

(2.7)$$\begin{gather}\left. \frac{\partial P}{\partial z}\right|_{\theta=\theta_0}={-}\left. \frac{\partial P}{\partial z}\right|_{\theta={\rm \pi}-\theta_0},\quad \text{at } z=0,1. \end{gather}$$

We refer to (2.6–2.7) as reflective boundary conditions in what follows of our derivation. For clarity, no periodic boundary conditions are applied at $z=0$ and $z=1$. It can be verified with simple algebra that the reflective boundary conditions satisfy the no-flux integral condition (2.5). From individual Brownian dynamics, when a micro-swimmer hits the boundary, the swimmer changes symmetrically the swimming direction and bounces off like a billiard. This strong condition might be realised via reasonable control strategies for artificial nano- and micro-robots, yet it is seemingly unrealistic for alive motile micro-organisms. In the orientation space, periodic boundary conditions are imposed as

(2.8)$$\begin{gather} \left.P\right|_{\theta=0} =\left.P\right|_{\theta=2 {\rm \pi}}, \end{gather}$$

(2.9)$$\begin{gather}\left.\frac{\partial P}{\partial \theta}\right|_{\theta=0} =\left.\frac{\partial P}{\partial \theta}\right|_{\theta=2 {\rm \pi}} . \end{gather}$$

The initial probability distribution $P^{(0)}$ is given as

(2.10)\begin{equation} P |_{t = 0} = \frac{1}{2 {\rm \pi}} \delta(x) P^{(0)} (z, \theta) , \end{equation}

where $P^{(0)}$ is a prescribed distribution function.

We introduce a differential operator defined as

(2.11)\begin{equation} \mathcal{L} (\boldsymbol{\cdot})\triangleq {{Pe}}_s \cos \theta \frac{\partial}{\partial z}(\boldsymbol{\cdot}) +\frac{\partial \varOmega(z, \theta)}{\partial \theta}(\boldsymbol{\cdot}) +\varOmega(z, \theta)\frac{\partial}{\partial \theta} (\boldsymbol{\cdot}) -D_{zz}(\theta)\frac{\partial^2}{\partial z^2}(\boldsymbol{\cdot}) -\frac{\partial^2}{\partial \theta^2} (\boldsymbol{\cdot}) , \end{equation}

with reflective boundary conditions in the $z$-space and periodic boundary conditions in the $\theta$-space. Note that the linear operator $\mathcal {L}$ is not self-adjoint if $\varOmega \neq 0$ or ${{Pe}} \neq 0$. The transport equation (2.1) for the probability distribution function could thus be re-written as

(2.12)$$\begin{gather} \frac{\partial P}{\partial t}+\left[{{Pe}}_f U_f(z)+{{Pe}}_s \sin \theta\right] \frac{\partial P}{\partial x}-D_{xx}(\theta)\frac{\partial^2 P}{\partial x^2}-2 D_{xz}(\theta) \frac{\partial^2 P}{\partial x \partial z}+\mathcal{L} P=0, \end{gather}$$

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

Further coupling of the external field (e.g. force fields) between the global (streamwise) and local spaces could be introduced in the meantime, which does not break the generality of the manipulation discussed hereafter by absorbing the external field into the velocity terms. It is supposed the coupled diffusivities are symmetrical,

(2.14)\begin{equation} D_{xz}(\theta)=D_{zx}(\theta), \end{equation}

as embodied in (2.4). When these coupled diffusing terms are set to zero, the resulting transport problem reduces to the special case in the existing dispersion analysis of active particles (Jiang & Chen 2019). As a substantial extension of the previous dispersion theory, the present case is more general by including the direct coupling between global and local transport processes.

3. Asymptotic solutions of macro-transport processes with anisotropic diffusion

Here, we proceed to evaluate the anisotropic diffusion of micro-swimmers by a homogenisation method. Since only the effects of local translational diffusivity in the lateral direction enter into the local operator, any coupling fluctuations of micro-swimmers can be neglected safely, to the leading order. That is, the relevant diffusivity component at the leading order only arises in the vertical direction, while the cross-diffusivities act in the sequence of higher orders. Some formal transformations are considered below.

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 2008, 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 1993; Jiang & Chen 2019) 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 2002; Manela & Frankel 2003), 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. 2021; Debnath et al. 2022; Wang, Jiang & Chen 2022a,  2023). 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 1984; Guan et al. 2022; Guan & Chen 2024), 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. Migration of confined micro-swimmers

In this section, we will investigate the shear-induced migration of confined micro-swimmers subject to anisotropic diffusion in a Poiseuille flow

(4.1)\begin{equation} U_f(z)={-}6z^2+6z-1. \end{equation}

The external flow field exerts two effects on the transport processes: (a) streamwise convection leading to spatial non-uniformity; (b) shear-induced rotation causing orientational redistribution. This introduces inherently multi-scale effects for the evolution of the shear-induced migration of micro-swimmers.

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. 2011; Vennamneni et al. 2020). 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 (2015) with extreme wall accumulation for a Robin boundary condition and by Ezhilan & Saintillan (2015) 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. (2020) 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. 2014). 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$.

5. Transient solutions and evolution of phenomenological transport coefficients

To comprehend the transient impacts of anisotropic diffusion, we investigate the temporal evolution of typical phenomenological transport coefficients. Non-spherical shape anisotropy and external gravitational forces are employed to induce vertical locomotion at a population level. Specifically, coupled diffusivities are revealed to influence remarkably the transient drift velocity $U_x(t)$ and dispersivity $D_T(t)$ at high orders.

Unlike spherical particles in a Poiseuille flow, elongated micro-swimmers display a persistent non-zero drift for asymptotically long times, as demonstrated in figure 8(a). Initially, the micro-swimmers are displaced on the centreline where flow convection is most pronounced, resulting in a notably positive drift relative to the mean flow speed. A sharp decline in transient drift is expected as micro-swimmers diffuse stochastically and swim actively, causing disparate sampling of individual streamlines with unequal probabilities over time. Consequently, this preferential sampling leads to a modest temporal increase in drift during the transitional phase. For the asymptotic dispersion period, the drift converges asymptotically to a constant value, as predicted by (3.15). In the scenario involving elongated micro-swimmers with robust self-propulsion, the upstream rheotaxis surpasses distinctly the flow convection for asymptotically long times, a trend that is augmented with higher cross-diffusivity $D_{xz}$. Utilising (3.15), it becomes apparent that the inclusion of a non-zero $D_{\|}$, corresponding to the isotropic condition in the present case, indeed magnifies the overall cross-diffusivity. A monotonic increase in the effective diffusivity, known as dispersivity, is observed due to the combined effects of swimming dispersion and shear dispersion. Specifically, anisotropic diffusion modifies dispersivity by means of streamwise diffusivity and cross-diffusivities. Figure 8(b) displays a notably increased diffusivity, by orders of magnitude, relative to the overall translational diffusivity. For the isotropic case, cross-diffusivities play no role in determining dispersivity; hence, effective diffusion is reinforced as overall streamwise diffusivity ($D_{xx}$) increases. With the predominance of shear dispersion subject to anisotropic diffusion, an increase in $D_{\perp }$ leads to a reduction in dispersivity, which arises from the adverse influence of associated coupling terms.

Figure 8. Transient distributions of phenomenological transport coefficients. The swimmers are initially displaced at $z=0.5$ and uniformly oriented. Common parameters are ${{Pe}}_s=1$, ${{Pe}}_f=10$, $\alpha _0=0.9$ and $\lambda =0$.

Finally, our focus shifts to gyrotactic micro-swimmers exposed to an external gravitational field, as shown in figure 9. For simplicity, we maintain the perpendicular diffusivity at zero and manipulate the parallel diffusivity to elucidate the influence of anisotropic diffusion. As evidenced previously in figure 2, gyrotactic micro-swimmers accumulate near the wall in the dispersion regime and exhibit a pronounced preference for upstream migration. Consequently, we anticipate that gyrotactic micro-swimmers will experience downstream convection with a positive drift as initially discharged at the centreline, later exhibit transient rises due to stochastic swimming and eventually converge to a constant negative value due to upstream rheotaxis. Upstream migration becomes more pronounced with increasing gyrotaxis, as shown in figure 9(a). Enhanced diffusion, exhibiting a monotonic time-dependent behaviour, is observed in comparison with the original translational diffusion, as depicted in figure 9(b). Notably, in the exceptional scenario of pronounced gyrotaxis ($\lambda =5$) coupled with low translational diffusivities, a non-monotonic dependence of transient dispersivity on time and a significant reduction in asymptotic dispersivity are observed. This phenomenon arises because nearly all micro-swimmers accumulate near the wall regions where flow convection is minimal, resulting in limited opportunities for bulk streamline sampling through Brownian diffusion.

Figure 9. Transient distributions of phenomenological transport coefficients under gravity. The swimmers are initially displaced at $z=0.5$ and uniformly oriented. Common parameters are ${{Pe}}_s=1$, ${{Pe}}_f=10$, $\alpha _0=0.9$, $D_{\perp }=0$ and $\lambda =0.5$ (unless otherwise declared).

6. Conclusion

This work establishes a population-level theoretical framework aimed at modelling macro-transport processes of elongated micro-swimmers involving anisotropic diffusion. In the absence of shape anisotropy and stochastic diffusion, we present a detailed examination of a benchmark solution of the zeroth moment for spheroidal gyrotactic micro-swimmers.

Utilising the analytical solutions, we demonstrate that symmetric actuation in a horizontal flow does not lead to net vertical locomotion of micro-swimmers from the probability perspective, unless non-spherical shapes, external fields or direct coupling effects are utilised to induce vertical polarisation. The inherent asymmetries and external fields can exert substantial influence on the observed concentration and orientation profiles of self-propelling micro-swimmers, as well as phenomenological transport coefficients.

Multiple patterns of shear-induced migration are elucidated with varying anisotropic diffusivities. In particular, highly elongated particles, moving slowly through a confined channel, tend to align themselves in a manner that amplifies the adverse effects of the strain rate, thereby adapting to balance the pronounced fluid rotation. It is intriguing to observe a distinct and persistent peak of the concentration distribution with moderate gyrotaxis, surpassing the nearby average values on both channel walls. Diverse collective behaviours encompassing upstream migration and shear alignment, as well as the formation of depletion layers, centreline focusing and surface accumulation, have been observed for various swim and flow Péclet numbers.

Notably, as the flow strength increases, distinctive phenomena of shear trapping are characterised by a rich array of symmetric depletion layers around the centreline. This implies that shear-induced migration exhibits a high degree of independence from the intricacies of the propulsion system and reorientation mechanisms. We also observe and explain the remarkable phenomenon whereby a reduction in the perpendicular diffusivity increases the effective diffusivity caused by the shear dispersion mechanism. A substantial decrease in steady dispersivity has been reported, deviating from the expectations of classical Taylor dispersion theory.

The theoretical analysis offers valuable insights into understanding distinct self-propulsion strategies and the transport dynamics within active-matter systems. As a first step, the current study has the potential to tackle a multitude of transport problems stemming from external or inter-particle interactions of spatially varying shapes, sizes and others for swimmers in constrained environments. For example, micro-swimmers can generate stress in the fluid, potentially causing collective motion at higher concentrations (Nordanger et al. 2022). Jeffery orbits, typically associated with Newtonian fluids, undergo notable changes due to viscoelastic stresses induced by the activity (Choudhary, Nambiar & Stark 2023). The conceptualisation of the rheological dynamics in dense active suspensions is absent from the present theoretical framework. Future work could encompass extensions to incorporate high-concentration interactions, non-Newtonian effects and realistic boundary conditions.