1. Introduction
The individual and collective behaviour of motile micro-organisms in flow environments covers abundant interesting phenomena with various applications, from hydrodynamic focusing by gyrotaxis as a novel method of cell concentration or separation (Kessler Reference Kessler1985a), to gyrotactic trapping of motile algae responsible for thin phytoplankton layers in marine turbulence (Durham, Kessler & Stocker Reference Durham, Kessler and Stocker2009). Unlike passive particles, self-propelled micro-organisms swim through a coupled process of translation and rotation, respectively in the position and orientation space, which could shape macroscopic features of the environmental landscape. The gyrotactic method, an emerging technique of concentrating or separating cells, combines the activity of cells and distinguishes them based on swimming behaviour and morphology (Kessler Reference Kessler1985b). For micro-algae bioreactors, algal suspensions are commonly pumped in an upright pipe to enhance the mixing of nutrients and prevent exposure to light (Bees & Croze Reference Bees and Croze2010). Through vertical pipe flows, surprising gyrotaxis-relevant structures have been observed experimentally and explored theoretically (Kessler Reference Kessler1984, Reference Kessler1985a,Reference Kesslerb; Pedley & Kessler Reference Pedley and Kessler1987).
In upright pipe flows, Kessler (Reference Kessler1984) proposed a simplified theory of gyrotaxis accounting for the hydrodynamic focusing of motile algae corresponding with experimental evidence (Kessler Reference Kessler1985a,Reference Kesslerb). Stable equilibrium distributions of swimming orientation are investigated by Pedley & Kessler (Reference Pedley and Kessler1987) with a vertical Poiseuille pipe flow. Bees & Croze (Reference Bees and Croze2010) and Hill & Bees (Reference Hill and Bees2002) developed the fundamental generalised Taylor dispersion (GTD) framework of Frankel & Brenner (Reference Frankel and Brenner1989, Reference Frankel and Brenner1993) to study the dispersion of gyrotactic micro-organisms in unbounded shear flows. Bearon, Hazel & Thorn (Reference Bearon, Hazel and Thorn2011) widened the applicability of the GTD model in a vertical pipe flow with two steps of mean operations, and obtained the drift and diffusivity in position space, referred to as the ‘two-step GTD model’ subsequently. As it is mathematically challenging and computationally expensive to obtain solutions of the Smoluchowski equation directly, probability density functions are approximated in terms of orientational moments, e.g. with a closure approximation that totally neglects nematic order (Baskaran & Marchetti Reference Baskaran and Marchetti2009; Saintillan & Shelley Reference Saintillan and Shelley2013), based on the observation that significant wall-normal polarisation is usually accompanied by relatively weak nematic alignment (Ezhilan & Saintillan Reference Ezhilan and Saintillan2015; Theillard & Saintillan Reference Theillard and Saintillan2019). Differing from the well-known continuum models like the Pedley–Kessler and two-step GTD models, Jiang & Chen (Reference Jiang and Chen2019, Reference Jiang and Chen2020) took the streamwise coordinate to represent the unique global space, and obtained steady solutions without approximations made to the fundamental Smoluchowski equation. Recent research progress on the transport of active particles in confined flows has led to a renewed interest in pre-asymptotic active dispersion (Schweitzer Reference Schweitzer2003; Jiang & Chen Reference Jiang and Chen2021). With time-dependent solutions directly from the Smoluchowski equation (Jiang & Chen Reference Jiang and Chen2021), the temporal evolution of dispersion of active particles through a vertical pipe is investigated in this work.
Although extensive efforts have been conducted for the transient diffusion of one- and two-dimensional motions of active particles in unbounded position space, there has been little quantitative analysis of the dispersion of self-propelled micro-organisms with translational motion and stochastic rotation in three dimensions. For example, Howse et al. (Reference Howse, Jones, Ryan, Gough, Vafabakhsh and Golestanian2007) derived a two-dimensional projected expression of mean square displacement (MSD) under the combined action of rotation and translation. The asymptotic analyses upon the theoretical solution (Howse et al. Reference Howse, Jones, Ryan, Gough, Vafabakhsh and Golestanian2007) will assign an initially diffusive regime (i.e. MSD linear in time), an intermediate ballistic regime (i.e. MSD in quadratic relation with time) and an asymptotic dispersion regime (i.e. diffusivity linear in time). These basic regimes of transient diffusion of active particles are collaborated perfectly with numerical simulations (Volpe, Gigan & Volpe Reference Volpe, Gigan and Volpe2014; Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016) and observed in experiments (Zheng et al. Reference Zheng, ten Hagen, Kaiser, Wu, Cui, Silber-Li and Löwen2013). For active particles with fluctuating swimming speed in two dimensions, Peruani & Morelli (Reference Peruani and Morelli2007) expressed the dispersion explicitly relevant to the characteristic time of speed fluctuation and orientation. In simple shear flows, Brownian motions of self-propelled particles are identified by rich superdiffusive scaling laws, with rotational diffusion unconfined or confined in one and two spatial dimensions along a linear channel (ten Hagen, van Teeffelen & Löwen Reference ten Hagen, van Teeffelen and Löwen2011a; Ghosh et al. Reference Ghosh, Misko, Marchesoni and Nori2013). Particular attention is drawn to the temporal scaling laws of MSD (ten Hagen, Wittkowski & Löwen Reference ten Hagen, Wittkowski and Löwen2011b). Furthermore, Sandoval et al. (Reference Sandoval, Marath, Subramanian and Lauga2014) derived explicitly a long-time expression of the full MSD tensor for active particles with the account of run-and-tumble dynamics. In the presence of external hydrodynamics and irregular geometries, series solutions can be given with the homotopy analysis method (Liao Reference Liao2004; Apaza & Sandoval Reference Apaza and Sandoval2020). Directly from the moment equations, we will derive theoretically an explicit benchmark solution of dispersion for active Brownian particles by decoupling from the shear dispersion. The general scaling laws in the absence of flow reveal various dispersion regimes, from the initial diffusive regime to the pre-asymptotic superdiffusive regime and eventually the dispersive regime. Richer dispersion regimes are present when the flow strength, flow direction, gyrotaxis and particle anisotropy participate jointly in the balance.
For a solute dispersing through a tube, Taylor (Reference Taylor1953) demonstrated that the centre of the injected material moves at the mean flow velocity and eventually follows a Gaussian distribution with an enhanced dispersion coefficient, honoured as the Taylor dispersion mechanism. The diffusive–advective coupling of orientable Brownian particles rotated and elongated in unbounded Poiseuille flow is explained clearly by Foister & Ven (Reference Foister and Ven1980), while the interaction between Taylor dispersion and shear velocity field was delivered later by Frankel & Brenner (Reference Frankel and Brenner1993). For self-propelled Brownian particles, richer physical mechanisms are doomed under the competitive effects of biased swimming, shear dispersion and Brownian-rotational diffusion. Rusconi, Guasto & Stocker (Reference Rusconi, Guasto and Stocker2014) investigated the shear-induced trapping under random orientation and cell alignment via microfluidic experiments. Bearon & Hazel (Reference Bearon and Hazel2015) described the trapping of bacteria in high-shear regions with a continuum framework in a channel flow in good agreement with experimental results. Interestingly, non-spherical active particles under shear-induced transport could also be trapped in low-shear regions (Barry et al. Reference Barry, Rusconi, Guasto and Stocker2015; Vennamneni, Nambiar & Subramanian Reference Vennamneni, Nambiar and Subramanian2020). Upstream swimming and non-monotonic dispersion of active Brownian particles are recovered by Peng & Brady (Reference Peng and Brady2020), systematically investigating the effect of translational diffusion, activity, flow strength and reorientation time on active dispersion. Nevertheless, owing to the lack of intensive study on pre-asymptotic active dispersion in pipe flows, little systematic research into the mutual role of swimming-induced dispersion and shear-induced dispersion is drawn. This indicates a need to understand the various perceptions of the predominance of self-propelled swimming that exist among abundant living organisms. In this work, we will first resort to a benchmark case of pre-asymptotic active dispersion in extreme shear where Brownian rotation or shear dispersion dominates alone, prior to general discussions under competitive effects. Sensitivity analysis upon active swimming, external flow field and particle anisotropy are elucidated.
The active particles considered here are self-propelled swimmers with orientable Brownian diffusion, exposed to an external vertical Poiseuille flow field. For bottom-heavy micro-organisms, e.g. Chlamydomonas nivalis, Heterosigma akashiwo and Dunaliella salina, joint efforts of gravitational and viscous torques could result in directed locomotion (Jeffery Reference Jeffery1922; Fenchel & Finlay Reference Fenchel and Finlay1984, Reference Fenchel and Finlay1986), referred to as gyrotaxis. Interesting phenomena are observed in pipes, for example, gyrotactic focusing on the axis in a downwelling pipe flow and wall accumulation at the periphery in an upwelling pipe flow (Kessler Reference Kessler1985a,Reference Kesslerb, Reference Kessler1986; Pedley & Kessler Reference Pedley and Kessler1992). In high-shear regions, transient thin phytoplankton is observed in situ and simulated numerically with an underlying gyrotactic trapping mechanism first proposed by Durham et al. (Reference Durham, Kessler and Stocker2009). This can be realised when strong viscous torques persistently override the gravitational torques so that gyrotaxis-driven upswimming is interrupted for transient periods from hours to weeks (Durham & Stocker Reference Durham and Stocker2012). Recently, Fung, Bearon & Hwang (Reference Fung, Bearon and Hwang2022) uncovered the interactions of gyrotactic active particles between orientational diffusion and shear dispersion by reducing the Smoluchowski equation, extending the work of Bearon & Hazel (Reference Bearon and Hazel2015) and Vennamneni et al. (Reference Vennamneni, Nambiar and Subramanian2020) with inclusion of biased motility. From the fundamental Smoluchowski equation, an additional settling speed of gyrotactic swimmers is found efficient to maintain the presence of thin layers (Wang, Jiang & Chen Reference Wang, Jiang and Chen2022). However, much uncertainty still exists about the temporal evolution of hydrodynamic accumulation. This work attempts to explore the potential transitions from one concentrative phenomenon to the other during pre-asymptotic dispersion. Given that the trajectory of an individual gyrotactic micro-organism depends on hydrodynamic and gravitational torques, the rate of change of orientation 
$\boldsymbol {p}$ is governed generally by (Leal & Hinch Reference Leal and Hinch1972; Pedley & Kessler Reference Pedley and Kessler1992; Bearon Reference Bearon2022)
\begin{equation} \dot{\boldsymbol{p}}=\frac{1}{2 B^\ast}\,[\boldsymbol{k}-(\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{p}) \boldsymbol{p}]+\frac{1}{2}\,\omega^\ast{\wedge} \boldsymbol{p}+ \alpha_0 \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{E}^\ast \boldsymbol{\cdot} (\boldsymbol{I}-\boldsymbol{p}\boldsymbol{p}), \end{equation}
where 
$B^\ast$ is the gyrotactic orientation time quantifying the balance between hydrodynamic torque against viscous resistance, 
$\boldsymbol {k}$ is the unit vector pointing vertically upwards, 
$\omega ^\ast$ is the vorticity, 
$\alpha _0$ is the shape factor, with 
$\alpha _0 = 0$ representing spheres while 
$\alpha _0 = 1$ denotes infinitely elongated thin rods, and 
$\boldsymbol {E}^\ast$ is the rate-of-strain tensor. Populations of swimming micro-organisms that are initially oriented randomly and discharged uniformly in the ambient flow could react to external fields and even change their shapes periodically (Elgeti & Gompper Reference Elgeti and Gompper2015; Walker et al. Reference Walker, Ishimoto, Moreau, Gaffney and Dalwadi2022). Based on the formulated dynamics of trajectory, the swimming micro-organisms could be guided to form a radial accumulation through a vertical pipe. The main question of interest to us is: how does this hydrodynamic focusing come into being? More interestingly, are there any newly observed regimes before approaching this final form?
This work is organised as follows. The transport formulation directly from the Smoluchowski equation is given in § 2. From the fundamental moment equations for the simple case without a background flow, we first derive a theoretical solution of MSD in three dimensions in § 3, and then report the temporal scaling laws validated against the benchmark results. Section 4 presents the pre-asymptotic dispersion of spherical active particles in downwelling and upwelling pipe flows under the competitive effects of shear dispersion and Brownian-rotational dispersion, where transient effective trapping is revealed as a newly observed transitional regime before hydrodynamic focusing. The effect of particle anisotropy on vortical rotation and shear alignment of gyrotactic active particles is investigated in § 5. Finally, § 6 concludes.
2. Formulation and solution
Consider a dilute suspension of self-propelled swimmers through a vertical pipe with radius 
$a^{\ast }$, as shown in figure 1. With gyrotaxis-driven upswimming and convection along the axial direction, the active particles are also subject to translational swimming and Brownian rotation in the projection of the radial direction. Fundamentally different from the traditional gyrotactic trapping in widely investigated horizontal flows, the transient gyrotaxis-relevant interaction through a vertical pipe shows a multi-dimensional effect responsible for richer individual behaviour and collective phenomena. In the position space, the velocity profile under polar coordinates 
$(r, \psi, z)$ can be expressed with base vectors 
$(\boldsymbol {e}_r, \boldsymbol {e}_{\psi }, \boldsymbol {e}_z)$ as 
$\boldsymbol {U} = U (r)\, \boldsymbol {e}_z,$ where 
$U$ is the flow speed.
Figure 1. Schematic of active particles in a downwelling flow through a vertical pipe. (a) Vertical pipe of radius 
$a^*$ with pointwise active particles initially oriented randomly and discharged uniformly at the origin point 
$O$ of the pipe. (b) Enlarged view of the point source of active particles. The swimming direction can be adjusted to be aligned with certain directions or distributed uniformly. (c) Descriptions of 
$\boldsymbol {p}$ in the cylindrical and spherical coordinate systems, where 
$\boldsymbol {p}$ is the swimming orientation with the polar angle 
$\theta$ and azimuthal angle 
$\phi$.
A set of spherical coordinates 
$(\rho, \theta, \phi )$ with base vectors 
$(\boldsymbol {e}_{\rho }, \boldsymbol {e}_{\theta }, \boldsymbol {e}_{\phi })$ is introduced, where 
$\theta$ is the polar angle between 
$\boldsymbol {p}$ and 
$\boldsymbol {e}_{\psi }$, and 
$\phi$ is the azimuthal angle between 
$\boldsymbol {k}$ and the projection of 
$\boldsymbol {p}$ onto the 
$z\unicode{x2013} r$ plane, as shown in figure 1(c). In the orientation space, the swimming direction vector can be expressed as
\begin{equation} \boldsymbol{p}=\boldsymbol{e}_{\rho} = p_r (\theta, \phi)\,\boldsymbol{e}_r + p_{\psi} (\theta, \phi)\, \boldsymbol{e}_{\psi} + p_z (\theta, \phi)\, \boldsymbol{e}_z, \end{equation}in which
\begin{equation} p_r =- \sin \theta \sin \phi, \quad p_{\psi} = \cos \theta, \quad p_z =- \sin \theta \cos \phi. \end{equation}With dimensionless variables and parameters defined as
\begin{equation} \left.\begin{gathered} t = t^{{\ast}} D^{{\ast}}_{{r}}, \quad r = \frac{r^{{\ast}}}{a^{{\ast}}},\quad z = \frac{z^{{\ast}}}{a^{{\ast}}} - {Pe}_{{f}}\,t, \quad U = \frac{U^{{\ast}}}{U^{{\ast}}_{{m}}} - 1, \\ {Pe}_{{s}} =\frac{V_{{s}}^{{\ast}}}{D^{{\ast}}_{{r}} a^{{\ast}}}, \quad {Pe}_{{f}} =\frac{U^{{\ast}}_{{m}}}{D^{{\ast}}_{{r}} a^{{\ast}}}, \quad D_{{t}} = \frac{D^{{\ast}}_{{t}}}{D^{{\ast}}_{{r}} a^{{\ast} 2}}, \quad \lambda = \frac{1}{2 B^{{\ast}} D^{{\ast}}_{{r}}}, \end{gathered}\right\} \end{equation}the governing equation in the moving coordinate system at the speed of mean flow velocity reads
\begin{equation} \frac{\partial P}{\partial t} + \boldsymbol{\nabla}_R \boldsymbol{\cdot} \left[\left({Pe}_{{f}}\,U\boldsymbol{e}_z +{Pe}_{{s}}\,\boldsymbol{p} \right) P - D_{{t}}\,\boldsymbol{\nabla}_R P \right] + \boldsymbol{\nabla}_p \boldsymbol{\cdot} \left(\dot{\boldsymbol{p}_{{r}}} P - \boldsymbol{\nabla}_p P \right) = 0, \end{equation}
where 
$\boldsymbol{\nabla} _R = (\partial / \partial r) \boldsymbol {e}_r + (1/r) (\partial / \partial \psi ) \boldsymbol {e}_{\psi } + (\partial /\partial z) \boldsymbol {e}_z$ , 
$\boldsymbol{\nabla} _p = (\partial /\partial \theta ) \boldsymbol {e}_{\theta } + (1/ \sin \theta )(\partial /\partial \phi ) $ 
$\boldsymbol {e}_{\phi }$, 
${Pe} _{{s}}$ is the swimming Péclet number normalised by the rotational diffusivity, 
${Pe} _{{f}}$ is the flow Péclet number, 
$D_{{t}}$ is the ratio of the translational diffusivity to rotational diffusivity, and 
$\lambda$ is the bias parameter.
For perfectly elastic collisions between particles and solid boundaries (Bearon et al. Reference Bearon, Hazel and Thorn2011; Volpe et al. Reference Volpe, Gigan and Volpe2014; Ezhilan & Saintillan Reference Ezhilan and Saintillan2015; Jakuszeit, Croze & Bell Reference Jakuszeit, Croze and Bell2019), reflective boundary conditions at the pipe walls read
$$\begin{gather} P (1, \psi, z, \theta, \phi, t)= P (1, \psi, z, \theta, - \phi, t), \end{gather}$$
$$\begin{gather}\frac{\partial P}{\partial r} (1, \psi, z, \theta, \phi, t) =- \frac{\partial P}{\partial r} (1, \psi, z, \theta, - \phi, t), \end{gather}$$which guarantee a zero probability flux normal to walls:
\begin{align} \int^{\rm \pi}_0 \sin \theta \,\mathrm{d} \theta \int^{2 {\rm \pi}}_0 \,\mathrm{d} \phi (\boldsymbol{e}_r \boldsymbol{\cdot} \boldsymbol{J}_R) = \int^{\rm \pi}_0 \sin \theta \,\mathrm{d} \theta \int^{\rm \pi}_{- {\rm \pi}} \mathrm{d} \phi \left( - {Pe}_{{s}} \sin \theta \sin \phi\,P - D_{{t}} \frac{\partial P}{\partial r} \right) =0. \end{align}The swimming flux is balanced not only by reflections (2.5) considered in Bearon et al. (Reference Bearon, Hazel and Thorn2011), Ezhilan & Saintillan (Reference Ezhilan and Saintillan2015) and Jiang & Chen (Reference Jiang and Chen2019), but also by the radial diffusive flux in (2.6) (Jiang & Chen Reference Jiang and Chen2020).
On the axis (
$r = 0$), it is reasonable to require
\begin{equation} P |_{r = 0} \neq \infty . \end{equation}
For 
$\psi$, the periodic boundary conditions are imposed as
\begin{equation} P |_{\psi = 0} = P |_{\psi = 2 {\rm \pi}},\quad \left. \frac{\partial P}{\partial \psi} \right|_{\psi = 0} = \left.\frac{\partial P}{\partial \psi} \right|_{\psi = 2 {\rm \pi}}. \end{equation}
Likewise, a finite probability condition for 
$\theta$ reads
\begin{equation} P |_{\theta = {\rm \pi}/{2}} \neq \infty, \quad P |_{\theta =- ({\rm \pi}/{2})} \neq \infty, \end{equation}
and for 
$\phi$, the periodic boundary conditions are imposed as
\begin{equation} P |_{\phi = {\rm \pi}} = P |_{\phi =- {\rm \pi}},\quad \left. \frac{\partial P}{\partial \phi} \right|_{\phi = {\rm \pi}} = \left.\frac{\partial P}{\partial \phi} \right|_{\phi =- {\rm \pi}} . \end{equation}
An initial probability distribution 
$P^{(0)}$ is given as
\begin{equation} P |_{t = 0} = P^{(0)} ( r,\psi,z,\theta, \phi) . \end{equation}The macroscopic transport characteristics can be derived rigorously from the moments of the probability distribution of particles. The local moments of probability density function (p.d.f.) are defined in the classical form (Aris Reference Aris1956; Brenner & Edwards Reference Brenner and Edwards1993) as
\begin{equation} P_n (r, \theta, \phi, t) \triangleq \int^{\infty}_{- \infty} z^n\,P (r, z, \theta, \phi, t)\,\mathrm{d} z,\quad n = 0, 1, \ldots. \end{equation}
The zeroth-order moment 
$P_0$ is the marginal p.d.f. of a function of 
$r$, 
$\theta$ and 
$\phi$, and therefore can be viewed as the number density distribution of active particles in the corresponding orientation space (Ezhilan & Saintillan Reference Ezhilan and Saintillan2015; Jiang & Chen Reference Jiang and Chen2019).
The global moments 
$M_n$, i.e. the moments of the mean probability density distribution 
$\bar {P}$, are given as
\begin{equation} M_n (t) \triangleq \int^{\infty}_{- \infty} z^n\,\bar{P}\,\mathrm{d} z = 2 {\rm \pi}\, \bar{P}_n(t), \quad n = 0, 1, \ldots, \end{equation}and
\begin{equation} \bar{P}(z,t) \triangleq \int ^{1}_0 r\,\mathrm{d} r \int^{\rm \pi}_0 \sin \theta \,\mathrm{d}\theta \int^{2 {\rm \pi}}_0 \mathrm{d} \phi\, P (r, z, \theta, \phi, t). \end{equation}
We use a bar to denote integration over the radial and orientation space. Obviously, 
$M_0=1$. In the position space, quantitative results of number density distribution are reported by introducing an orientation-space-mean p.d.f.
\begin{equation} \langle P_0 \rangle_O (r,t) \triangleq \int_{0}^{\rm \pi} \sin \theta \,\mathrm{d} \theta \int_{0}^{2{\rm \pi}} \mathrm{d} \phi\, P_0 (r, \theta, \phi, t), \end{equation}
representing an accurate description of the zeroth-order moment in the 
$r$–
$\psi$ plane (here averaged over 
$\psi$).
We can calculate the basic macroscopic transport characteristics with the aid of global moments. The drift 
$U_d$ is related to the first-order moments and dispersivity 
$D_T$ to the second-order moments, as
$$\begin{gather} U_d (t) \triangleq \frac{\mathrm{d} \mu_x}{\mathrm{d} t} = \frac{\mathrm{d} M_1}{\mathrm{d} t}, \end{gather}$$
$$\begin{gather}D_T (t) \triangleq \frac{1}{2}\,\frac{\mathrm{d} \sigma^2}{\mathrm{d} t} = \frac{1}{2}\,\frac{\mathrm{d} M_2}{\mathrm{d} t} - M_1\,\frac{\mathrm{d} M_1}{\mathrm{d} t}, \end{gather}$$
where 
$\mu _x$ and 
$\sigma ^2$ are the mean displacement and MSD, respectively:
\begin{equation} \mu_x \triangleq \frac{M_1}{M_0} = M_1 \end{equation}and
\begin{equation} \sigma^2 \triangleq \frac{M_2}{M_0} - \frac{M^2_1}{M^2_0} = M_2 - M_1^2 . \end{equation}The asymptotically long-time values corresponding to the macroscopic transport coefficients have been explored recently by Jiang & Chen (Reference Jiang and Chen2020), extending the classical dispersion problem devised by Taylor (Reference Taylor1953).
We also derive the skewness of the p.d.f. to depict the asymmetries of the probability distribution, which is of appreciable effects during the initial and anomalous diffusion stage. The skewness 
$\gamma _1$ is expressed as
\begin{equation} \gamma_1 \triangleq \frac{\kappa_3}{\sigma^3}, \end{equation}
with the third-order cumulant 
$\kappa _3$ of the p.d.f. defined as
\begin{equation} \kappa_3 \triangleq \frac{M_3}{M_0} - 3\,\frac{M_2 M_1}{M_0^2} + 2\,\frac{M_1^3}{M_0^3} = M_3 -3 M_2 M_1 +2 M_1^3. \end{equation}For further research interest, the temporal evolution of the above macroscopic transport coefficients can feature the longitudinal transport process in the pre-asymptotic anomalous diffusion stage before the dispersion regime. Employing the bi-orthogonal expansion method (Strand, Kim & Karrila Reference Strand, Kim and Karrila1987; Brezinski Reference Brezinski1991; Nambiar et al. Reference Nambiar, Phanikanth, Nott and Subramanian2019; Jiang & Chen Reference Jiang and Chen2021), we can derive the local moments by series expansions of surface spherical harmonics and Bessel functions with the aid of the Galerkin method to deal with the associated eigenvalue problem, as detailed in Appendix A.
Besides, the normalised polarisation vector 
$\boldsymbol {P}(r,t)$ and nematic order parameter 
$\boldsymbol{\mathsf{Q}}(r,t)$ commonly referred to in liquid-crystalline systems could be defined in terms of orientational moments of the p.d.f. (Ezhilan & Saintillan Reference Ezhilan and Saintillan2015; Yan & Brady Reference Yan and Brady2015; Theillard & Saintillan Reference Theillard and Saintillan2019; Peng & Brady Reference Peng and Brady2020) as
$$\begin{gather} \boldsymbol{P}(r,t) \triangleq \int_{0}^{\rm \pi} \sin \theta \,\mathrm{d} \theta \int_{0}^{2{\rm \pi}} \mathrm{d} \phi\, \boldsymbol{p}(\theta, \phi)\, P_0(r, \theta, \phi, t) / \langle P_0 \rangle_O (r,t), \end{gather}$$
$$\begin{gather}\boldsymbol{\mathsf{Q}}(r,t) \triangleq \int_{0}^{\rm \pi} \sin \theta \,\mathrm{d} \theta \int_{0}^{2{\rm \pi}} \mathrm{d} \phi\, \left(\boldsymbol{pp}-\frac{\boldsymbol{I}}{3}\right) P_0(r, \theta,\phi, t) / \langle P_0 \rangle_O (r,t), \end{gather}$$
which would help to clarify the configuration of the swimmer in various flow conditions and particle anisotropies. Due to symmetries in the present problem, the components of interest are the normalised wall-normal polarisation 
$P^{r}$, streamwise polarisation 
$P^{z}$, and shear nematic alignment parameter 
$Q^{rz}$. We will also focus on orientational moments
$$\begin{gather} \boldsymbol{m}(r,t) \triangleq \langle P_0 \rangle_O (r,t)\,\boldsymbol{P}(r,t), \end{gather}$$
$$\begin{gather}\boldsymbol{\mathsf{D}}(r,t) \triangleq \langle P_0 \rangle_O (r,t)\,\boldsymbol{\mathsf{Q}}(r,t), \end{gather}$$for cases with strong polarisation to avoid numerical dispersion, where similar superscripts are adopted to denote the components.
Based on the above formulation, pre-asymptotic dispersion of active particles with random swimming, gyrotactic torque, external flow convection, and rotational and translational diffusion can be modelled from a single-particle Smoluchowski equation. Through a vertical pipe, the idealised pointwise swimmers are initially oriented randomly and discharged uniformly at the origin point of pipe. Their swimming orientation is instantaneously variable subject to the internal energy stored in the biological body, the external stimulation of flow field, and the reflective effect of the boundary wall. Simultaneous efforts of gravitational torques and hydrodynamic torques could form directed locomotion of gyrotactic swimmers, known as gyrotaxis. Kessler (Reference Kessler1984, Reference Kessler1985b, Reference Kessler1986) and Pedley & Kessler (Reference Pedley and Kessler1992) for the first time proposed the well-known hydrodynamic focusing phenomenon in pipe flows, including gyrotactic focusing on the axis in the downwelling regions and wall accumulation at the periphery in the upwelling regions. In recent years, Durham et al. (Reference Durham, Kessler and Stocker2009) and Durham & Stocker (Reference Durham and Stocker2012) proposed that one of the mechanisms behind the red tides of harmful algae could be due to gyrotactic trapping in horizontal shear flows, which will be realised when strong viscous torques persistently override the gravitational torques so that gyrotaxis-driven upswimming is interrupted. Projecting this case into a vertical pipe wherein the upswimming characteristics of gyrotactic swimmers and hydrodynamic focusing are in the orthogonal directions, it would be interesting to investigate what gyrotactic trapping in horizontal flows could be like with such an upright configuration, which is commonly adopted for micro-algae bioreactors for cultivation and separation.
Through a vertical pipe, the swimmers are revealed in this work to be vortically rotated under shear–gyrotaxis interaction to form a pre-asymptotic concentration peak in the radial direction near the wall before the approach to directed locomotion of swimming-induced hydrodynamic focusing. Here, we explain the physics behind the transition from transient concentration peaks in the near-wall region towards hydrodynamic focusing in upwelling and downwelling flows with the illustration in figure 2. For point gyrotactic swimmers discharged on the axis, the velocity gradient of intrinsic Poiseuille flow increases from zero, which means that the particles tend to move outwards in upwelling flows, with a symmetric mechanism of gyrotactic focusing in downwelling flows (Kessler Reference Kessler1985a). The linear increment of velocity gradient of the Poiseuille pipe flow allows the free swimmers trapped with rapid rotations in the high-shear regions away from the pipe axis and the wall. Owing to the gyrotaxis-driven upswimming, on average the swimmers perform stably outward locomotion persistently. This transient shear–gyrotaxis interaction of pointwise active particles is responsible for richer phenomena (e.g. gyrotactic focusing, wall accumulation, pre-asymptotic concentration peaks, and potential transitions from one to the other) observed in upright pipe flows with different shear strengths through a vertical pipe.
Figure 2. Illustration of transient shear–gyrotaxis interactions of active particles through a vertical pipe with the reflective boundary conditions. In an upwelling Poiseuille pipe flow at the left hand, the gyrotactic swimmer pointing inwards to the axis (
$-{\rm \pi} <\phi <0$) is unstable since the joint efforts of the gravitational torque 
$T_G$ and viscous torque 
$T_\nu$ rotate the particles to move outwards and upwards. With the temporal evolution of gyrotaxis, the self-propelled particles would be trapped in the middle region away from the pipe wall before their eventual accumulation at the periphery of the pipe. In a downwelling Poiseuille pipe flow at the right hand, gyrotactic swimmers are instead stable with inward locomotion (
$0<\phi <{\rm \pi}$) in specific shear regions. The competition between the pre-asymptotic near-wall peak and hydrodynamic focusing results in the rich phenomena observed in vertical pipe flows with different shear strengths. Transient effective trapping of pointwise active particles through a vertical pipe is revealed to be quite different from that in a horizontal flow as commonly observed.
Because transient symmetry breaking and other deviations from normality are known to be rather persistent, it is of basic interest to explore the following questions.
(i) What is the origin of hydrodynamic focusing?
(ii) How do active swimming, external flow and particle anisotropy influence this evolution?
These questions could be answered only if the time-dependent solutions of the Smoluchowski equation through a vertical pipe are obtained. Although the Pedley–Kessler and two-step GTD models successfully gave some remarkable physical features, inconsistencies and restrictions from experimental observations and numerical simulations have been reported broadly (Bearon et al. Reference Bearon, Hazel and Thorn2011; Croze et al. Reference Croze, Sardina, Ahmed, Bees and Brandt2013; Croze, Bearon & Bees Reference Croze, Bearon and Bees2017; Fung, Bearon & Hwang Reference Fung, Bearon and Hwang2020; Jiang & Chen Reference Jiang and Chen2020). The time-dependent solution directly from the Smoluchowski equation enables a rational treatment of a broader range of parameter space, without the need to take the Péclet number and variance of shear stress small as practised in previous analytical methods. Attention will first be given to the pre-asymptotic spreading of spherical gyrotactic swimmers in the absence of external flow, to reveal the general scaling laws of dispersion with the predominance of only self-propulsion and diffusion. Hydrodynamic force and torque will later participate in the competition where the direction of the external flow could enhance or hinder the gyrotactic rotation, giving birth to richer phenomena. The origin of hydrodynamic focusing would be analysed by elucidating the predominance of swimming in the absence of background flow, the competitive effects of shear dispersion and Brownian rotation, and influences of particle anisotropy, in sequence. In the paradigm of dispersion framework, we will present the number density distribution, polarisation, nematic order parameter, drift velocity, dispersivity and skewness of non-spherical particles as functions of activity, particle anisotropy and external flow strength.
3. Predominance of self-propelled swimming: benchmark solution without background flows
The case with the predominance of self-propelled swimming is of interest in its own right. This section is devised to depict the effect of swimming without any ambient flow. It is taken to decouple the shear-induced dispersion and deduce a benchmark solution through a superposition of the active Brownian particle dynamics and translational contribution.
3.1. Number density, polarisation and nematic order parameters
The effect of swimming-induced dispersion on the local distribution of passive and self-propelled Brownian particles is illustrated in figure 3 with different swimming abilities. For passive particles (
${Pe} _s=0$), the local distribution witnesses a peak on the axis due to the initial displacement at 
$r=0$. Its approach to uniformity is much slower than that of active particles since the translational diffusion is usually negligible compared to Brownian rotational diffusion. Richer physical structures are anticipated for active particles with different swimming Péclet numbers. Although weak swimming locomotion makes no qualitative difference to the local distribution in contrast with passive particles, the rotational diffusion of swimming orientation results in a more effective approach to transverse uniformity under the swimming-induced dispersion process. The spherical active particles with stronger locomotion tend to concentrate radially near the pipe centre at relatively smaller times. The maximum point of local distribution is closely related to 
${Pe} _s$, as shown in figures 3(a,b). The particles will eventually reach an equilibrium state, and the swimmer with a larger swimming speed disperses faster under the simultaneous effects of swimming-induced dispersion and translational diffusion, as shown in figure 3(c). There is no doubt that all the curves will eventually coincide with each other as a uniform line within the limit of the dispersion regime.
Figure 3. The predominance of swimming is revealed during the approach to dispersion, in which the active swimming dispersion mechanism cooperates with the classical shear dispersion mechanism. Orientation-space-mean p.d.f.s 
$\langle P_0 \rangle _O$ of spherical active particles is shown with respect to 
$r$ without a background flow for different swimming Péclet numbers at (a) 
$t=0.5$, (b) 
$t=1$, and (c) 
$t=2$. All swimmers are initially randomly oriented and discharged with uniform distribution at the origin of the pipe. Common parameters: 
${Pe} _f=0$, 
$D_t=0.01$, 
$\alpha _0=0$ and 
$\lambda =0$.
We further explore the effect of activity upon distributions of number density, polar order and nematic order parameters to give a more mechanistic understanding of the dispersion phenomena without a background flow. Swimming dispersion is enhanced with the activity, which results in a uniform distribution for asymptotically large times. Unique concentration peaks in the middle regions are observed for swimmers with comparably strong activity at 
$t=0.5$, as shown in figure 4(a). Wall-normal components of the polarisation vector 
$m^r=-\int _0^{2{\rm \pi} } \mathrm {d}\phi \int _0^{{\rm \pi} } \mathrm {d} \theta \, P_0 \sin ^2\theta \sin \phi$ could well depict the mean swimming direction of active particles along the radial direction. Additionally, the wall-normal nematic order parameter 
$D^{rr}=\int _0^{2{\rm \pi} } \mathrm {d}\phi \int _0^{{\rm \pi} } \mathrm {d} \theta \, P_0 \sin \theta$ 
$(\sin ^2\theta \sin ^2\phi -\frac {1}{3})$helps to analyse the radial alignment of active particles irrespective of swimming orientation. For spherical active particles with medium motility, the concentration peak is accompanied by polarisation in the central region and nematic alignment. This is not a coincidence because the concentrative phenomena in the middle region give rise to a large diffusion flux, which is balanced by the wall-normal swimming flux. As a result of Brownian diffusion, the swimmers preferentially turn outwards radially to smooth out the non-uniformity. In contrast, active particles with strong activity reach the equilibrium state in advance for this benchmark case.
Figure 4. Distributions in the absence of background flow at 
$t=0.5$ as functions of swimming Péclet number: (a) number density, (b) wall-normal polarisation vector, and (c) wall-normal nematic order parameter. Pointwise active particles are initially oriented randomly and displaced uniformly at the origin of the pipe. Common parameters: 
${Pe} _f=0$, 
$D_t=0.01$, 
$\alpha _0=0$ and 
$\lambda =0$.
3.2. Dispersion
With the understanding of mean displacement being zero due to the symmetry of random walk activity, the MSD in the axial direction 
$\sigma ^2$ can characterise the dispersion of active particles. Using the well-known Ornstein–Uhlenbeck formulation (Uhlenbeck & Ornstein Reference Uhlenbeck and Ornstein1930), quantitative descriptions are made precisely by calculating the MSD in two dimensions (Howse et al. Reference Howse, Jones, Ryan, Gough, Vafabakhsh and Golestanian2007; Ghosh et al. Reference Ghosh, Misko, Marchesoni and Nori2013; Volpe et al. Reference Volpe, Gigan and Volpe2014; Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016). For passive Brownian particles, the two-dimensional MSD is 
$4 D_t t$ for asymptotically large times. For self-propelled Brownian particles, the normalised two-dimensional theoretical MSD reads (Howse et al. Reference Howse, Jones, Ryan, Gough, Vafabakhsh and Golestanian2007; Peruani & Morelli Reference Peruani and Morelli2007; Volpe et al. Reference Volpe, Gigan and Volpe2014)
\begin{equation} {MSD}=(4D_t+{Pe}_s^2)t+\frac{{Pe}_s^2}{2}\left[\exp{(-2t)}-1\right]. \end{equation}
This is formally equivalent to the Ornstein–Uhlenbeck formula with respect to the MSD of Brownian particles with inertia (Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016), accounting for the ballistic (
${MSD} \propto t^2$) and diffusive (
${MSD} \propto t$) regimes (Uhlenbeck & Ornstein Reference Uhlenbeck and Ornstein1930). By solving analytically the Langevin equation, ten Hagen et al. (Reference ten Hagen, van Teeffelen and Löwen2011a) depicted the overdamped Brownian motion of self-propelled particles with the translational and rotational diffusion confined to one or two dimensions.
We can initiate these calculations with the concentration moments to three-dimensional active particles to derive a benchmark value of dispersion 
$\sigma ^2=\langle x^2 \rangle$ without background flows in vertical pipes:
\begin{equation} \sigma^2=\left(2 D_t+ \frac{{Pe}_s^2}{3}\right)t+ \frac{{Pe}_s^2}{6} \left[\exp{(-2t)}-1\right].\end{equation}
Differentiating (3.2) with respect to 
$t$ according to (2.18), we deduce a benchmark solution of dispersivity
\begin{equation} D_T= D_t+ \frac{{Pe}_s^2}{6}-\frac{{Pe}_s^2}{6} \exp{(-2t)}.\end{equation} As shown in figure 5, pre-asymptotic dispersion of spherical active particles through a pipe without a background flow is illustrated. The numerical computation of (2.20) and (2.18) collaborates perfectly with the benchmark solution (3.2) and (3.3), and computation is also validated against Lagrangian particle simulation in Appendix B. The active Brownian particle model adopted in the present work assumes a constant swimming speed, with the direction of velocity oriented smoothly via rotational Brownian diffusion, which is applicable for self-phoretic colloids or smooth swimmers such as E. coli mutants (Cates & Tailleur Reference Cates and Tailleur2015). For passive particles, the dispersion basically follows a diffusive scaling with respect to 
$t$. In contrast, the dispersion of self-propelled active particles experiences an exponential decay, showing both ballistic and diffusive temporal scaling laws. Within the short time scales of 
$t \ll 1$, (3.2) reduces to 
$\sigma ^2 \sim 2 D_t t$ when diffusion dominates over self-propulsion. This means that the motion of active particles in the initial regime can indeed be featured by an explicit diffusive law, which is deduced implicitly in figure 5. In the case of weak self-propulsion, this diffusive scaling has also been reported with experimental observations (Zheng et al. Reference Zheng, ten Hagen, Kaiser, Wu, Cui, Silber-Li and Löwen2013). Noting that the dimensional reorientation time exerted by the Brownian rotational diffusion on the active particles is characterised by 
$t_R^{\ast } \equiv 1/D_r^*$, 
$t \sim 1$ means that the dimensional time 
$t^{\ast }$ is of the same order as 
$t_R^{\ast }$. It is more intuitive to examine (3.2) in its dimensional form as 
$\sigma ^2=2 D_t^{\ast } t^{\ast }+V_s^{\ast 2}t^{\ast }t_R^{\ast }/(3a^{\ast 2})-V_s^{\ast 2}t_R^{\ast 2}/(6a^{\ast 2})+O(t^{\ast }) \sim 2 D_t^{\ast } t^{\ast }+ D_r^{\ast 2} t^{\ast 2}/6 \times V_s^{\ast 2}/$ 
$(a^{\ast 2} D_r^{\ast 2})$, where the second term on the right-hand side could be rewritten as the product of the non-dimensional swimming Péclet number 
${Pe} _{{s}} = {V_{{s}}^{\ast }}/{D^{\ast }_{{r}} a^{\ast }}$ and time 
$t = t^{\ast } D^{\ast }_{{r}}$ with a ballistic law. That is, for relatively larger times of characteristic reorientation temporal scales (
$t \sim 1$), we have deduced 
$\sigma ^2 \sim 2 D_t t+ ({{Pe} _s^2}/{6})t^2$. This superdiffusive law is shown clearly in figure 5(a) with tiny translational diffusivity and non-zero swimming Péclet numbers during the intermediate regimes. In the absence of flow, we demonstrate that non-diffusive scaling is produced solely with the random swimming of active particles. For asymptotically large times (
$t\gg 1$), the dispersion follows an enhanced diffusive scaling 
$\sigma ^2 \sim (2 D_t+ {{Pe} _s^2}/{3})t$ and reaches the so-called dispersion regime. The corresponding enhanced diffusivity, which effectively increases in magnitude compared to the translational diffusivity, is the generalised Taylor dispersivity, as plotted in figure 5(b). To sum up, the general scaling laws reveal various dispersion regimes in the absence of ambient flow, at least from the initial diffusive regime, to the pre-asymptotic superdiffusive regime and eventually the dispersive regime. Richer dispersion regimes are doomed to present when the flow strength, flow direction, gyrotaxis and particle anisotropy participate in the balance.
Figure 5. General scaling laws without background flows: (a) MSD and (b) dispersivity. Ballistic and diffusive scaling laws are examined for the temporal evolution of the effective diffusivity, dedicating the initial regime, pre-asymptotic regime and dispersion regime. Numerically computed (lines) and benchmark (symbols) values of MSD and 
$D_T$ for active particles are illustrated with swimming Péclet numbers 
${Pe} _s=0$ (circles), 
${Pe} _s=0.1$ (triangles), 
${Pe} _s=0.5$ (diamonds), 
${Pe} _s=1$ (squares) and 
${Pe} _s=2$ (inverted triangles). Computation of (2.20) and (2.18) collaborates perfectly with the benchmark solution (3.2) and (3.3), respectively. All swimmers are initially oriented randomly, and discharged with uniform distribution at the origin of the pipe. Common parameters: 
${Pe} _f=0$, 
$D_t=10^{-4}$, 
$\alpha _0=0$ and 
$\lambda =0$.
4. Competition of shear dispersion and Brownian rotation: unique effective trapping
In this section, we further consider the pre-asymptotic dispersion of spherical active particles in a Poiseuille pipe flow under the competitive effects of shear dispersion and Brownian-rotational dispersion. We will depict the respective asymptotic cases in weak shear (
${Pe} _f/{Pe} _s\ll {}1$), strong shear (
${Pe} _f/{Pe} _s\gg {}1$) and general shear flows (arbitrary 
${Pe} _f/{Pe} _s$) subsequently.
The steady-state local distribution of active particles in pipe flows is elucidated by Jiang & Chen (Reference Jiang and Chen2020), who explain the cell accumulation phenomena with the underlying mechanisms of hydrodynamic focusing in general downwelling and upwelling pipe flows. This section pays special attention to the pre-asymptotic dispersion of active particles in vertical Poiseuille pipe flows with extreme shear, and then discusses the case of arbitrary shear flows without loss of generality.
4.1. Active dispersion in weak shear: Brownian rotation dominates
For 
${Pe} _f\ll 1$, the case of pre-asymptotic dispersion in weak shear where the Brownian rotational diffusion is predominant over the effects of shear is discussed in the position space and orientation space.
In the position space, the radial distribution of the orientation-space-mean p.d.f.s 
$\langle P_0 \rangle _O$ of non-gyrotactic and gyrotactic active particles differ quantitatively exposed to weak shear as shown in figure 6. A peak of concentration in the middle region away from the wall is observed in figure 6(a), with a smooth transition towards gyrotactic focusing on the axis of the pipe within the dispersion regime. All the swimmers are discharged as a central point source, while they tend to focus on the axis during the pre-asymptotic status owing to gravitational torques and vortical rotation in downwelling flows, which means that the radial distribution of 
$\langle P_0 \rangle _O$ forms a peak at 
$r=0$. For large times, the peak values of 
$\langle P_0 \rangle _O$ of gyrotactic active particles are comparably greater than non-gyrotactic ones with corresponding parameter configurations, because gyrotaxis drives the particles upwards remarkably, and prevents the distribution from being completely uniform. The self-propelled particles in downwelling flows eventually reach an equilibrium state, with gyrotactic swimmers focusing at the pipe centre, and non-gyrotactic ones approaching uniformity at asymptotically large times, as illustrated in figures 6(a) and 6(b), respectively. Gyrotactic focusing in upwelling flows is plotted in figure 6(c), where particles finally concentrate at the pipe wall in contrast to the uniform line of non-gyrotactic active particles in figure 6(b).
Figure 6. Orientation-space-mean p.d.f.s 
$\langle P_0 \rangle _O$ of spherical active particles with respect to 
$r$ in weak shear as functions of time with different flow directions: (a) gyrotactic active particles (
$\lambda =1$) in a downwelling Poiseuille flow; (b) non-gyrotactic active particles (
$\lambda =0$) in a downwelling Poiseuille flow; and (c) gyrotactic active particles (
$\lambda =1$) in an upwelling Poiseuille flow. The approach to hydrodynamic focusing of gyrotactic micro-organisms and uniformity of non-gyrotactic active particles is of particular interest. Peculiar concentration peaks in the middle region of upright pressure-driven flows are observed for only gyrotactic active flows. All swimmers are initially oriented randomly and discharged with uniform distribution at the origin of the pipe. Common parameters: 
${Pe} _f=0.1$, 
${Pe} _s=1$, 
$D_t=0$ and 
$\alpha _0=0$.
4.2. Active dispersion in strong shear: shear dispersion dominates
In the limit of 
${Pe} _f \gg 1$, the case of pre-asymptotic dispersion in strong shear where the effect of shear is predominant over the Brownian-rotational diffusion is discussed.
In the position space, the orientation-space-mean p.d.f. of active particles peaks at the centre of the pipe within the dispersion regime as shown in figures 7(a,b). Under the effect of gyrotaxis, the peak values of active particles are remarkably larger for the greater axial velocity on the axis of intrinsic Poiseuille flow. The distribution of 
$\langle P_0 \rangle _O$ of non-gyrotactic swimmers will eventually approach uniform lines (cf. Jiang & Chen Reference Jiang and Chen2020, figure 3), while gyrotactic swimmers in the downwelling flows eventually sample a biased peak distribution at 
$r=0$. The competitive effect of upwelling convection and gyrotaxis-driven upswimming results in much richer physical structures. A peak of concentration is observed once again in the middle region away from the pipe wall during the pre-asymptotic stage, with a smooth transition towards the wall accumulation, as shown in figure 7. Discharged uniformly in the orientation and position space at 
$r=0$, the active particles tend to move outwards and form a crest in the middle of the radial position. In this case, the time scale to reach an equilibrium state is comparably longer owing to the combined swimming locomotion of shear dispersion and gyrotaxis.
Figure 7. Orientation-space-mean p.d.f.s 
$\langle P_0 \rangle _O$ of spherical active particles with respect to 
$r$ in strong shear as functions of time with different flow directions: (a) gyrotactic active particles (
$\lambda =1$) in a downwelling Poiseuille flow; (b) non-gyrotactic active particles (
$\lambda =0$) in a downwelling Poiseuille flow; and (c) gyrotactic active particles (
$\lambda =1$) in an upwelling Poiseuille flow. The flow strength and direction influence directly the concentration distribution of active particles, namely gyrotactic focusing on the axis in downwelling flows and wall accumulation at the periphery for upwelling cases. Interestingly, a natural transition from effective trapping towards hydrodynamic focusing is revealed inherently within such configurations. All swimmers are initially oriented randomly and discharged with uniform distribution at the origin of the pipe. Common parameters: 
${Pe} _f=10$, 
${Pe} _s=1$, 
$D_t=0$ and 
$\alpha _0=0$.
Shear trapping of pointwise active particles through a vertical pipe is revealed to be quite different from that in a horizontal flow as commonly observed, as shown in figures 2 and 6(a). Note that a unique collective phenomenon as a pre-asymptotic transitional mode before hydrodynamic focusing, named transient effective trapping, is also observed during the temporal evolution in the strong upwelling flow at relatively large times, e.g. figure 7(c) at 
$t=8.0$. Interestingly, the competition between shear trapping and wall accumulation results in a fluctuated trapping layer. This transient ‘vertical trapping layer’ appears when the gravitational torques persistently override the viscous torques so that the upswimming of gyrotactic self-propelled particles is disrupted by the shear rotation and trapped at a distance to the wall. The trapping of gyrotactic swimmers in high-shear regions is later substituted by hydrodynamic focusing. This implies that gyrotactic trapping is only a transient phenomenon without the continuous supply of energy e.g. from the external flow field.
We would like to dig further into the orientational characteristics of the evolutional origin of hydrodynamic focusing. For micro-organisms like phytoplankton, the magnitude of swimming velocity is usually small compared with the flow speed. However, even low activity still triggers strong concentrative behaviour and surprising spatio-temporal patterns, as shown in figure 8. We specify three typical moments of time to illustrate the dynamics of three regimes: (a) initial regime at 
$t=3$; (b) pre-asymptotic regime at 
$t=8$; and (c) dispersion regime at 
$t=50$. The number density profile shown in figure 8(a) reveals the temporal evolution of gyrotactic active particles with a smooth transition, from initial convective diffusion, transient effective trapping towards the final wall accumulation. Wall-normal polarisation 
$P^{r}=-\int _0^{2{\rm \pi} } \mathrm {d}\phi \int _0^{{\rm \pi} } \mathrm {d} \theta \, P_0 \sin ^2\theta \sin \phi / \langle P_0\rangle _O$ and streamwise polarisation 
$P^{z}=-\int _0^{2{\rm \pi} } \mathrm {d}\phi \int _0^{{\rm \pi} } \mathrm {d} \theta \, P_0 \sin ^2\theta \cos \phi / \langle P_0\rangle _O$ could well capture the mean direction of swimming, shown in figures 8(b,c). As anticipated, the unique concentration peak is accompanied by strong polarisation in high-shear regions away from the central axis and pipe wall. Note that during this pre-asymptotic regime characterised by shear trapping, strong upstream rheotaxis also occurs in contrast to uniform downstream polarisation for asymptotic dispersion regimes. The shear nematic alignment parameter 
$Q^{rz}=\int _0^{2{\rm \pi} } \mathrm {d}\phi \int _0^{{\rm \pi} } \mathrm {d} \theta \, P_0 \sin ^2 \theta \sin \phi \cos \phi / \langle P_0\rangle _O$ helps to analyse the shear alignment of active particles irrespective of swimming orientation, shown in figure 8(d). The active swimming is revealed to move outwards stably near the central axis because of gyrotactic balance, as illustrated in figure 2, whereas the trapping layer witnesses an unstable configuration behaving as tumbling. This indicates an evolutional origin of hydrodynamic focusing, with the newly observed regime named transient effective trapping revealed as a pre-asymptotic transitional mode.
Figure 8. Evolutional origin of hydrodynamic focusing in an upwelling Poiseuille flow: (a) number density; (b) wall-normal polarisation vector; (c) wall-normal nematic order parameter; and (d) shear nematic alignment parameter. Pointwise active particles are initially oriented randomly and displaced uniformly at the origin of the pipe. Three typical time scales specified in different lines are 
$t=3$ (initial regime), 
$t=8$ (effective trapping in the pre-asymptotic regime), and 
$t=50$ (dispersion regime). Common parameters: 
${Pe} _f=10$, 
${Pe} _s=1$, 
$D_t=0$, 
$\alpha _0=0$ and 
$\lambda =1$.
4.3. Double dispersion mechanisms: arbitrary Péclet numbers
On the basis of the above discussions, we come to a more general case of arbitrary Péclet numbers, where the competitive physics of swimming dispersion and shear dispersion are in effect simultaneously. In this subsection, we investigate how the unique effective trapping is induced by gyrotaxis and active dispersion. Double dispersion mechanisms are characterised by 
${Pe} _f$ and 
${Pe} _s$ fixed, respectively.
Macroscopic transport coefficients, including the drift 
$U_d$, dispersivity 
$D_T$ and skewness 
$\gamma _1$, are also calculated for arbitrary flow Péclet numbers in figure 9. The drift denotes the first-order time derivative term of the moment 
$M_1$ under a moving coordinate at the speed of mean flow velocity. Different from the persistent zero drift for particles dispersing without a background flow or the transient non-zero drift for non-gyrotactic active particles in parallel flows (Guan et al. Reference Guan, Zeng, Li, Guo, Wu and Wang2021, Reference Guan, Zeng, Jiang, Guo, Wang, Wu, Li and Chen2022; Jiang & Chen Reference Jiang and Chen2021), we find persistent non-zero values of 
$U_d$ with gyrotactic swimmers. The non-monotonic variation of 
$U_d$ results from the combined action of convection, Brownian rotation, and gyrotaxis-driven upswimming. The ups and downs in figure 9(a) show the transitions of interacting competitive mechanisms. With the increment of flow Péclet numbers, the value of the drift 
$U_d$ within the dispersion regime increases. For strong shear flows (e.g. 
${Pe} _f=5$ and 10), the convection dominates the upswimming of gyrotactic active particles to result in a downward drift. Overall, the non-zero drift represents that the mass centre of the cloud of gyrotactic swimmers moves beyond the mean flow speed on average from the pre-asymptotic dispersion regime.
Figure 9. Macroscopic transport characteristics as functions of flow Péclet number under the effect of downwelling shear flows: (a) drift, (b) dispersivity, and (c) skewness. By fixing 
${Pe} _s$ as unity, the flow strength significantly enhances the Taylor dispersion of active particles, with the combined action of active dispersion. Persistent asymmetries are reported during the pre-asymptotic regime. Interestingly, transitions from positive to negative skewness are illustrated in weak shear for the interchanges of various balance mechanisms. All gyrotactic swimmers are initially oriented randomly and discharged with uniform distribution at the origin of the pipe. Common parameters: 
$D_t=0$, 
$\alpha _0=0$ and 
$\lambda =1$.
The dispersivity is a key measurement for active dispersion, consisting of translational diffusion and swimming-induced diffusion according to the generalised Taylor dispersion theory. Of particular interest to us is the temporal evolution of gyrotactic active particles exposed to pipe flows with different 
${Pe} _f$ values, and we thus remove the influence of translational diffusion for clarity. Unlike the ballistic and diffusive time scalings specified by the benchmark solution (3.3) in figure 5(b), richer variations of dispersivity are observed with the addition of flows as shown in figure 9(b). In weak shear, the gyrotaxis-driven upswimming is predominant over the advection of the downwelling Poiseuille flow as a conflicting mechanism. During the initial stage, the downwelling flow flushes the particles on the axis with a maximum velocity, while it is too weak to change the upswimming characteristics of gyrotactic swimmers on average. Therefore, the dispersivity for 
${Pe} _f$ experiences an almost monotonic increase with time. The balance of this dominance will be broken with a slight increment in flow Péclet number (which is comparable to the swimming Péclet number), i.e. the dispersivity can experience a sudden drop due to the flushing of opposite flow velocity. During the pre-asymptotic period, the upswimming of active particles near the axis unfolds this upstream swimming against the flow. The swimming-induced diffusion eventually leads to an effectively enhanced dispersivity with the continuous temporal effect of gyrotaxis. Giant active Taylor dispersion has been reported with hydrodynamic gradients (Dehkharghani et al. Reference Dehkharghani, Waisbord, Dunkel and Guasto2019). When a strong downwelling flow is applied to the gyrotactic swimmers, much more complicated variations of dispersivity are observed. The time scale to reach an equilibrium state is remarkably delayed.
The skewness of the p.d.f. is plotted in figure 9(c). During the initial stage, the convection flushes the particles on the axis at twice the speed of the mean flow velocity. Thus negative skewness of the local p.d.f. of active particles is presented for downstream convection, with the same reason as for passive particles (Aris Reference Aris1956; Aminian et al. Reference Aminian, Bernardi, Camassa, Harris and McLaughlin2016). When the gyrotaxis-driven upswimming overrides the convection of relatively weak downwelling flows, the axial p.d.f. reveals a tendency of the right-skewed distribution shown as the transient positive values in figure 9(c). However, the gyrotactic focusing on the axis, along with the effect of convection, leads to the rapid transition to negative skewness. For asymptotically large times, the skewness approaches zero for all flow Péclet numbers. That is, the mean concentration distribution of active particles approaches normality asymptotically in the streamwise direction like passive particles (Chatwin Reference Chatwin1970). Also, the time scale to reach an equilibrium state is sensitive to the flow intensity and swimming speed.
5. Particle anisotropy matters: implications on cultivation and separation of motile cells
Shape profoundly impacts the swimming behaviour, concentrative phenomena and equilibrium status on the transport of gyrotactic active particles. Pre-asymptotic active dispersion of swimmers with a specified aspect ratio under the combined effect of biased upswimming, vortical and Brownian rotation, and shear alignment are elucidated. For ellipsoidal active particles (
$\alpha _0>0$), shear-induced alignment is present due to noticeable strain effects of the ambient fluid (Ezhilan & Saintillan Reference Ezhilan and Saintillan2015). Preferential alignment and spatial heterogeneity of elongated bacteria have been reported through experimental investigations (Rusconi et al. Reference Rusconi, Guasto and Stocker2014; Barry et al. Reference Barry, Rusconi, Guasto and Stocker2015). Streamwise alignment (Pedley & Kessler Reference Pedley and Kessler1992) is remarkable for infinitely thin rods (
$\alpha _0=1$). In this section, particular attention is paid to the high-shear regions where transitions from pre-asymptotic near-wall peaks of concentration to wall accumulation are observed. The competing physics of the transient shear–gyrotaxis interaction gives rise to richer individual behaviour observed in pipe flows, e.g. shear alignment, vortical rotation and gyrotaxis-driven upswimming of ellipsoidal active particles.
5.1. Number density, polarisation and nematic order parameters
The orientation-space-mean number density distribution 
$\langle P_0 \rangle _O$ of ellipsoidal active particles with respect to 
$r$ in an upwelling flow is shown in figure 10 with different aspect ratios (
$\alpha _0=0$, 0.5, 0.75 and 
$1$). During the pre-asymptotic stage, convective diffusion in the middle position of high-shear regions is the predominant mechanism, denoted by the middle concentrative distribution in figure 10(a). Elongated rods are subject to shear-induced alignment with the ambient flow, leading to a slight decay in the peak values. Pre-asymptotic evolution of ellipsoidal gyrotactic swimmers uncovers richer interactions of gyrotaxis and shear, in contrast with the wall accumulation of spherical active particles for 
$t\gg 1$. The elongated shape draws remarkable differences, as illustrated in figure 10(b) on the axis of the pipe where the maximum velocity of Poiseuille flow is sampled, in contrast with the obvious streamwise alignment (rheotaxis) of infinitely thin rods. The diversity of peak values with various shapes of active particles is presented in figures 10(c,d) during the dispersion regime, with a sharp decrease in the concentration at the periphery. This is owing to the competitive physics behind the shear-induced alignment and directed locomotion. Shear alignment and vortical rotation weaken the concentrated behaviour, and hence cause the postponement of the equilibrium state for ellipsoidal self-propelled particles. For asymptotically large times, spherical active particles exhibit biased upswimming, while the swimming directions of the elongated rods are disrupted by the shear-induced alignment. The shear alignment not only influences the asymptotic angular distribution, but exerts a significant impact on the pre-asymptotic active dispersion of three-dimensional swimmers under the interactions of Brownian rotational diffusion, gyrotaxis-driven locomotion, wall accumulation and Taylor dispersion mechanisms.
Figure 10. Effect of particle anisotropy in an upwelling flow. Shear alignment and vortical rotation are the main mechanisms of ellipsoidal active particles, giving rise to the appreciable phenomena of enhanced upstream rheotaxis, effective trapping and dispersion. Times: (a) 
$t=1$, (b) 
$t=2$, (c) 
$t=4$, and (d) 
$t=10$. Pointwise gyrotactic swimmers are initially oriented vertically upwards and displaced uniformly at the origin of the pipe. Common parameters: 
${Pe} _f=1$, 
$D_t=0$, 
${Pe} _s=1$ and 
$\lambda =1$.
We now investigate further the distributions of number density, polar order and nematic order parameters to give a more mechanistic understanding of the dispersion phenomena. Effective trapping is observed for elongated ellipsoidal particles, while spheroids form wall accumulation at 
$t=2$, as shown in figure 11(a). Streamwise components of the polarisation vector 
$P^z=-\int _0^{2{\rm \pi} } \mathrm {d}\phi \int _0^{{\rm \pi} } \mathrm {d} \theta \, P_0 \sin ^2\theta \cos \phi / \langle P_0 \rangle _O$ could well depict the axial swimming flux of active particles. Additionally, the streamwise nematic order parameter 
$Q^{zz}=\int _0^{2{\rm \pi} } \mathrm {d}\phi \int _0^{{\rm \pi} } \mathrm {d} \theta \, P_0 \sin \theta \,(\sin ^2\theta \cos ^2\phi -\frac {1}{3}) / \langle P_0 \rangle _O$ helps to analyse the streamwise alignment of active particles irrespective of swimming orientation. For spherical swimmers within the asymptotic dispersion regimes, hydrodynamic focusing in an upwelling Poiseuille flow forms the concentration peak at the periphery, while shear alignment makes the particles turn parallel to the flow direction. As a result, absolute values of 
$P^z$ and 
$Q^{rz}$ located near the wall where 
$P_0$ reaches the maximum indicate that 
$\theta$ approaches 
${\rm \pi} /2$, and 
$\phi$ tends to be 
$0$ or 
${\rm \pi}$ on average. On the other hand, unique effective trapping is observed for gyrotactic ellipsoids. For extremely elongated particles, this transient phenomenon is more appreciable, and the maximum in the middle high-shear regions results in the non-monotonic orientation distribution shown in figures 11(b,c). Particle anisotropy, as a critical factor, influences the evolution of hydrodynamic focusing through the orientational configuration of elongated active particles. As 
$\alpha _0$ increases, the elongated swimmers preferentially turn parallel to the flow direction, where strong polarisation in the middle regions is weakened, as seen in figure 11(b). Upstream swimming of active rods is present near the wall, whilst spheroids are convected downstream and moving outwards as a result of gyrotaxis, as shown in figure 11(c). Strong shear alignment is accompanied by effective trapping, compared to isotropy of spheroids in figure 11(d). In summary, the extremely elongated gyrotactic swimmers are revealed to sample the high-shear regions followed by radial accumulation with strong shear alignment and upstream rheotaxis during the pre-asymptotic regime.
Figure 11. Distributions in an upwelling Poiseuille flow as functions of particle anisotropy: (a) number density; (b) wall-normal polarisation vector; (c) wall-normal nematic order parameter; and (d) shear nematic alignment parameter. The time specified corresponds to the effective trapping of elongated active rods observed in the pre-asymptotic regime at 
$t=4$. Pointwise active particles are initially oriented randomly and displaced uniformly at the origin of the pipe. Common parameters: 
${Pe} _f=1$, 
$D_t=0$, 
${Pe} _s=1$ and 
$\lambda =1$.
5.2. Macroscopic transport coefficients
Although the steady-state gyrotactic focusing and active dispersion could be well predicted through a normal distribution (Kessler Reference Kessler1985b; Pedley & Kessler Reference Pedley and Kessler1992), the observed distributions of both spherical and ellipsoidal active particles during the pre-asymptotic are normally not Gaussian. This indicates a need to understand the various perceptions of pre-asymptotic non-Gaussian cumulants of the probability density distribution that persist long from the initial stage to the pre-asymptotic stage. Pre-asymptotic variation of macroscopic transport characteristics of ellipsoidal active particles is henceforth derived in upwelling shear flows through a vertical pipe, as shown in figure 12.
Figure 12. Temporal evolution of macroscopic transport characteristics as a function of particle anisotropy in upwelling shear flows: (a) drift, (b) dispersivity, and (c) skewness. Different aspect ratios (
$\alpha _0=0$, 0.5 and 
$1$) of spherical and elongated gyrotactic swimmers are investigated under the effect of shear alignment and vortical rotation. The dominance of shear-induced alignment over vortical rotation leads to upstream swimming during the pre-asymptotic regime. The persistent effective trapping and enhanced dispersion of elongated particles could contribute to cell separation and cultivation. All pointwise gyrotactic swimmers are initially oriented randomly and discharged with uniform distribution at the origin of the pipe. Common parameters: 
${Pe} _f=1$, 
$D_t=0$, 
${Pe} _s=1$ and 
$\lambda =1$.
According to the generalised Taylor dispersion theory (Frankel & Brenner Reference Frankel and Brenner1993; Hill & Bees Reference Hill and Bees2002; Jiang & Chen Reference Jiang and Chen2020), explicit expression of the drift could be deduced rigorously as 
$U_d=\overline {[{Pe} _f\,U(r)- {Pe} _s \sin \theta \cos \phi ] P_0^\infty }$. The drift is the total product of the whole longitudinal speed and the number density distribution. Particle anisotropy is embodied implicitly in the number density since the shape factor would significantly influence the local distribution. The value of the drift shows the relation between the average speed of the fluid particles and active particles over the cross-section. The asymptotic constant value of 
$U_d$ is non-trivial because a stable phase deviation exists for the synchronous motion of a fluid particle in an Eulerian sense and active particles from the Lagrangian perspective. In other words, the active particles behave as if there were no biased force in arbitrary radial positions, as opposed to the actual situation. This could be interpreted as the active particles ‘forgetting’ some of their initial positional information when sampling all the regions for various times and dispersing for asymptotically large times. The temporal evolution of the drift 
$U_d$ is plotted for different shapes of gyrotactic swimmers in figure 12(a). In the case of upwelling flows, active particles exhibit biased upswimming for most of the time and reach the equilibrium state with non-zero drift values for different aspect ratios. With the particles elongated, the pre-asymptotic absolute value of maximum drift (
$|U_d|$) decreases as a result of the streamwise alignment against the transient gyrotaxis-relevant locomotion in the preferential direction. On the other hand, the asymptotic values of 
$|U_d|$ during the dispersion regime increase with the aspect ratios, showing that elongated rods eventually move faster aligned with the streamwise flow. Although effects of the shape have a minor impact on the qualitative tendency of the drift curves, we observe transiently positive 
$U_d$ for the extremely elongated ellipsoidal active particles (e.g. 
$\alpha _0=0.75$ and 
$1$) even under the combined action of upwelling convection and gyrotaxis. The positive value corresponds to the effective trapping phenomena in which particles are likely to sample lower flow speeds near the wall. The non-monotonic variation of drift with time is caused by the transition from the gyrotaxis-driven upswimming towards the dispersion mechanisms. The non-zero negative drift physically means that the moving speed of the cloud of gyrotactic swimmers is greater than the mean speed of the upwelling ambient flow. It is interesting to note that when particles accumulate at the wall in the dispersion regime, where it experiences almost no flow convection, they move even faster than the average flow due to the comparable motility on their own.
In the light of the generalised Taylor dispersion theory, explicit expression of the asymptotic dispersivity could be deduced rigorously as 
$D_T=D_t+\overline { [{Pe} _f\,U(r)-}$ 
$\overline{{Pe} _s \sin \theta \cos \phi ] b }$, wherein 
$b(r,\theta,\phi )$ is the solution of 
$\mathcal {L}b=P_0^\infty [{Pe} _f\, U(r)- {Pe} _s \sin \theta$ 
$\cos \phi -U_d]$. It represents the eventual dispersivity as a sum of molecular diffusivity and convective dispersivity. The influence of particle anisotropy is embedded into the initial information of probability (the integration of 
$b$), which will not be lost with time like other initial conditions, however. Figure 12(b) presents the pre-asymptotic variation of dispersivity of gyrotactic active particles with different shapes. During the initial stage, the dispersivity increases as the gyrotactic particles swim along the upwelling ambient flow. Unlike the monotonic increase in dispersivity with time for the spherical particles without a background flow in figure 5(b), or with a downwelling flow in figure 9(b), the dispersivity of gyrotactic active particles in an upwelling flow falls rapidly when the swimmers are trapped in the high-shear regions, and increases to a steady-state enhanced value, as denoted by the stagnation point corresponding to the concave-up 
$D_T$ during the pre-asymptotic dispersion regime. The transition from pre-asymptotic near-wall peaks to wall accumulation is sensitive to aspect ratios of active particles. Over the whole temporal course, elongated swimmers exhibit a monotonically enhanced dispersivity in contrast with the spherical ones owing to the shear-induced alignment. The shear alignment not only influences the asymptotic angular distribution, but exerts a significant impact on the pre-asymptotic active dispersion of three-dimensional swimmers under the interactions of biased swimming, hydrodynamic focusing and Taylor dispersion mechanisms. The persistent effective trapping is appreciable for the ellipsoidal micro-swimmers and enhanced with the elongation, shedding new light on applications of cell separation for various micro-organisms with different aspect ratios. Note that the time scale to reach the equilibrium state is relevant to the shape of active particles. The elongation of ellipsoidal particles is revealed to be beneficial for dispersion due to the shear alignment and vortical rotation mechanisms. The skewness of spherical and ellipsoidal active particles is plotted in figure 12(c). The negative skewness during the initial regime represents a left-skewed distribution of mean probability density distribution under the convection of ambient flow. The variation of skewness is much more complicated for elongated rods when the shear alignment takes turns to dominate during the convective regime. Particle anisotropy exerts remarkable impacts on 
$\gamma _1$ during the pre-asymptotic stage, whereas all the curves overlap approximately for 
$t>0.5$. The curves first decrease and then increase to zero, resulting in an inflexion point owing to the gyrotaxis-relevant directed locomotion in the high-shear regions. The mean concentration distribution of elongated rods approaches normality more slowly than the spherical swimmers as shown in figure 10, since the shear-induced alignment rotates them to swim parallel to the flow on average against the biased locomotion during the dispersion regime. The revealed characteristics could contribute to cell separation and purification of elongated gyrotactic micro-organisms, for example, by adjusting the external field.
6. Conclusion
A vertical pipe flow inherently introduces complicated multi-scale effects corresponding to smooth transitions among gyrotaxis-relevant mechanisms. Based on an intensive investigation of pre-asymptotic dispersion of active particles through a vertical pipe, the present study is devised to explore the evolutional origin of hydrodynamic focusing by effects of active swimming, external flow and particle anisotropy.
With the assumption of axisymmetric rigid spheroids and rods, model gyrotactic active particles are exemplified in this paper. The main results are as follows.
(i) In the absence of background flow, we deduce theoretically an explicit benchmark solution of transient dispersion for self-propelled Brownian particles with three-dimensional motions without arbitrary approximations made to the moment equations. The ballistic and diffusive scaling laws are deduced with the predominance of self-propelled swimming, and computation is verified with the benchmark solution and Lagrangian particle tracking.
(ii) With weak shear, when swimming dispersion dominates, hydrodynamic focusing of gyrotactic active particles and the approach to the normality of non-gyrotactic particles are reproduced. Specifically, hydrodynamic focusing phenomena could indicate the wall accumulation in an upwelling flow or gyrotactic focusing in a downwelling flow, which will eventually replace the transient near-wall peaks.
(iii) With moderate shear, when double dispersion mechanisms exist, particle anisotropy matters profoundly concerning the swimming behaviour, collective phenomena and time scale to reach an equilibrium status. The rotational dynamics of active particles is affected by the flow and swimming Péclet numbers: the orientational moments vary non-monotonically with
${Pe} _f$ and ${Pe} _s$. Additionally, the active swimming, external flow and particle anisotropy bring persistent asymmetries into the particle dispersion process, giving rise to the appreciable phenomena of enhanced upstream rheotaxis, effective trapping and dispersion.(iv) With strong shear, when Taylor dispersion dominates, flow strength and direction are revealed to shape directly the concentration distribution. Note that a unique gyrotactic concentration layer named transient effective trapping is observed as a new dispersion regime. This indicates that different strategies of biomass transport and flow control should be adopted for various pressure-driven active flows.
(v) Interestingly, inherent transitions from effective trapping to hydrodynamic focusing are reported for the first time. After the convection-dominated regime, the active particles have a tendency to peak in the middle regions away from the pipe wall induced by the intrinsic feature of linearly varying shear in Poiseuille pipe flow during the pre-asymptotic regime, with a smooth transition towards hydrodynamic focusing during the dispersion regime. Transient effective trapping of gyrotactic active particles through a vertical pipe is uniquely different from the commonly observed concentrative phenomena in horizontal flows.
The findings in this paper are subject to at least three limitations. First, we assume a dilute suspension of model active particles that has no impact on the ambient flow. Potential discrepancies with experimental results could be due to the local counter-effect of micro-swimmers to shear, as tested by Croze et al. (Reference Croze, Bearon and Bees2017). Bio-convection could be present for the density differences between the self-propelled particles and the background flow (Hill & Pedley Reference Hill and Pedley2005; Bees Reference Bees2020). Dynamics of swimmers exposed to an external flow with hydrodynamic interactions exerts remarkable differences on surface accumulation (Chilukuri, Collins & Underhill Reference Chilukuri, Collins and Underhill2015). Extensions to include the hydrodynamic interactions (Takatori & Brady Reference Takatori and Brady2017), bed-load particle motions (Wu et al. Reference Wu, Jiang, Zeng and Fu2023) and high-concentration suspensions (Zade et al. Reference Zade, Costa, Fornari, Lundell and Brandt2018) with the current theoretical framework are desired for future work. Second, a basic parallel flow through a pipe is applied for mathematical convenience. A natural progression of the present work is to explore the dispersion in a turbulent channel flow (Guo & Chen Reference Guo and Chen2022; Li et al. Reference Li, Gong, Jiang, Zhan, Wang, Fu, Xu and Wu2023) or three-dimensional pipe turbulence. Richer rheological properties and macroscopic patterns (Saintillan Reference Saintillan2018; Morris Reference Morris2020) are expected under the three-dimensional flow structure. Third, idealised reflective boundary conditions are adopted, which guarantees only a zero probability flux normal to walls. For realistic cell–wall interaction, more complicated behaviour is observed for micro-organisms near the wall (Zeng, Jiang & Pedley Reference Zeng, Jiang and Pedley2022). Appropriate boundary conditions and dynamical description for the cell–cell, cell–wall and cell–fluid interactions are still open questions.
In conclusion, the contributions of the present work are: (a) intensively elaborating the multi-scale effects of rheological factors of pressure-driven flows through a vertical pipe on the dynamics of active Brownian particles; and (b) a first analytical attempt to show an overall picture of active transport regimes induced by swimming and shear dispersion within a broader range of parameter space. In addition to recovering the hydrodynamic focusing, upstream rheotaxis and non-monotonic dispersion behaviour explored by previous work, we investigate systematically the flow strength, activity, gyrotaxis, particle anisotropy and Brownian diffusion on the dynamics of active particles in Poiseuille pipe flows. In particular, we summarise the transient effective trapping mechanism as an origin of hydrodynamic focusing through a vertical pipe, which has not been analysed in detail to the best of our knowledge. The results of this work might enrich our understanding of the peculiar collective phenomena and dispersive dynamics of non-spherical active particles in pipe flows. From the applied perspective, the present work could potentially shed light on cell separation and cultivation related to micro-algae reactor technology as well as flow control and biomass transfer in confined environments.
Funding
This work is supported by the National Natural Science Foundation of China (grant nos 52079001 and 52109093).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Method of bi-orthogonal expansion
For completeness, we give analytical solution procedures of moments of the p.d.f. through a vertical pipe with the bi-orthogonal expansion technique. Following Jiang & Chen (Reference Jiang and Chen2020), the unbounded axial coordinate 
$z$ is chosen as the global space variable 
${\boldsymbol {Q}}= (z)$, while the local space is 
${\boldsymbol {q}}= (r, \psi, \theta, \phi )$. According to (2.13), the governing moment equation of the p.d.f. (2.4) reads (Aris Reference Aris1956)
\begin{equation} \frac{\partial P_n}{\partial t} +\mathcal{L}P_n = n (n - 1) D_t P_{n - 2} + n \left[ {Pe}_{{f}}\,U (r) -{Pe}_{{s}} \sin \theta \cos \phi \right] P_{n - 1},\quad n = 0, 1, \ldots, \end{equation}
where 
$\mathcal {L}$ is an operator defined in the local space 
${\boldsymbol {q}}$, explicitly as
\begin{align} \mathcal{L} P_n &\triangleq - \frac{{Pe}_{{s}} \sin \theta \sin \phi}{r}\, \frac{\partial (r P_n)}{\partial r} + \frac{{Pe}_{{s}} \cos \theta}{r}\, \frac{\partial P_n}{\partial \psi} - D_{{t}}\left[\frac{1}{r}\,\frac{\partial}{\partial r} \left(r\,\frac{\partial P_n}{\partial r} \right) + \frac{1}{r^2}\, \frac{\partial^2 P_n}{\partial \psi^2} \right] \nonumber\\ &\quad + \frac{1}{\sin \theta}\,\frac{\partial (\dot{\theta} \sin \theta P_n)}{\partial \theta} + \frac{\partial (\dot{\phi} P_n)}{\partial \phi} - \frac{1}{\sin \theta}\,\frac{\partial}{\partial \theta} \left( \sin \theta\, \frac{\partial P_n}{\partial \theta} \right) - \frac{1}{\sin^2 \theta}\, \frac{\partial^2 P_n}{\partial \phi^2}. \end{align}
Here, 
$P_{- 1} = P_{- 2} = 0$ is introduced to write the governing equation for all the moments in a consistent form.
The boundary conditions of 
$P_n$ (
$n = 0, 1, \ldots$) are in the same form as those of 
$P$, i.e. (2.5) and (2.6). For reflective conditions,
$$\begin{gather} P_n \left( r,\psi,\theta ,\phi ,t \right) = P_n \left( r,\psi,\theta ,-\phi ,t \right),\quad\mathrm{at}\ r = 1, \end{gather}$$
$$\begin{gather}\frac{\partial P_n}{\partial r} \left( r,\psi,\theta ,\phi ,t \right) =- \frac{\partial P_n}{\partial r} \left( r,\psi,\theta ,-\phi ,t \right),\quad \mathrm{at}\ r = 1. \end{gather}$$
For 
$\psi$, periodic conditions are imposed at the boundary:
\begin{equation} \left.\begin{gathered} P_n |_{\psi = 0} = P_n |_{\psi = 2{\rm \pi}},\\ \left. \frac{\partial P_n}{\partial \psi}\right|_{\psi = 0} = \left. \frac{\partial P_n}{\partial \psi} \right|_{\psi =2 {\rm \pi}}. \end{gathered}\right\} \end{equation}
In the orientation space, periodic conditions for 
$\phi$ read
\begin{equation} \left.\begin{gathered} P_n |_{\phi =-{\rm \pi}} = P_n |_{\phi = {\rm \pi}},\\ \left. \frac{\partial P_n}{\partial \phi} \right|_{\phi =-{\rm \pi}} = \left. \frac{\partial P_n}{\partial \phi} \right|_{\phi = {\rm \pi}}, \end{gathered}\right\} \end{equation}
and finite conditions for 
$\theta$ are given as
\begin{equation} P_n |_{\theta = {\rm \pi}/2} \neq \infty, \quad P_n |_{\theta =-{\rm \pi}/2} \neq \infty. \end{equation}The initial conditions according to (2.12) are
\begin{equation} \left.\begin{gathered} P_0 |_{t = 0} = \frac{1}{2 {\rm \pi}}\,P^{(0)}(r, \phi), \\ P_n |_{t = 0} = 0, \quad n = 1, 2, \ldots. \end{gathered}\right\} \end{equation}The governing formulation of the global moments can also be obtained through similar procedures. With the integration by parts formula, we deduce
\begin{equation} \overline{\mathcal{L}P_n} = 0, \quad n = 0, 1, \ldots, \end{equation}under the reflective condition (A4). With simple algebra, we have
\begin{equation} \frac{\mathrm{d} M_n}{\mathrm{d} t} = n (n - 1) D_t M_{n - 2} + n\, \overline{\left( {Pe}_f\,U(r) - {Pe}_s \sin \theta \cos \phi \right) P_{n - 1}}, \quad n = 1,2 \ldots. \end{equation}In particular, the drift is exactly the local-weighted average of the axial velocity component
\begin{equation} U_d = \overline{\left( {Pe}_f\,U(r) - {Pe}_s \sin \theta \cos \phi \right) P_0}. \end{equation}Although extensive research has been carried out on solutions of moments of a solvent with the method of separation of variables and/or integral transform method (Barton Reference Barton1983), the existing accounts fail to solve the auxiliary Sturm–Liouville problem in the case of active particles due to the complexity of function basis for the expansion. The auxiliary eigenvalue formulation of (A1) in terms of moments of the p.d.f. could be expressed as
\begin{equation} \mathcal{L} f_i = \lambda_i f_i, \end{equation}
where 
$\lambda _i$ is the eigenvalue (
$i=1,2,\ldots$) and 
$f_i$ the associated eigenfunction with the reflective boundary conditions of 
$P_n$. In the case of 
$\lambda _1=0$, 
$f_1$ asymptotes to the long-time limit of 
$P_0$, as has been discussed fully in Jiang & Chen (Reference Jiang and Chen2019, Reference Jiang and Chen2020).
It is almost technically intractable to express the solution of the eigenfunction 
$f_i$ in an explicit form, as a result of the complicated operator 
$\mathcal {L}$. Employing the Galerkin method, we approximate the eigenvalues and eigenfunction (
$\lambda _i$ and 
$f_i$) with proper precision. Concrete expressions of the bases with respect to the reflective condition are shown hereinafter. Once we have obtained the basis functions, written as 
$\{e_i\}_{i=1}^{\infty }$, corresponding to the specific boundary conditions, we can express the eigenfunction 
$f_i$ as
\begin{equation} f_i = \sum_{j = 1}^{\infty} \phi_{i j} e_j. \end{equation}
Here, 
$\phi _{ij}$ denotes the expansion coefficient. Corresponding to the local operator 
$\mathcal {L}$, we can derive the bilinear form 
$A({\cdot }, {\cdot })$ (or matrix form denoted by 
$\boldsymbol{\mathsf{A}}$) in the light of the bases. The elements read
\begin{equation} A_{i j} = A(e_i, e_j) = \langle e_i, \mathcal{L} e_j \rangle, \quad i=1,2,\ldots,\enspace j=1,2,\ldots. \end{equation}
Here, 
$\langle {\cdot }, {\cdot } \rangle$ presents the associated inner product, which is defined for specified functions 
$f$ and 
$g$ as
\begin{equation} \langle \,f, g\rangle \triangleq \int_{0}^{1} r \,\mathrm{d} r \int_{0}^{\rm \pi} \sin \theta \,\mathrm{d} \theta \int_{0}^{2 {\rm \pi}} \mathrm{d} \phi\, f(r, \theta, \phi)\,g(r, \theta, \phi). \end{equation}A weak expression of (A12) can be given as
\begin{equation} \boldsymbol{\mathsf{A}} \boldsymbol{\phi}_i= \lambda_i \boldsymbol{\phi}_i. \end{equation}
The vector of the eigenfunction coefficients 
$f_i$ is denoted as 
${\boldsymbol {\phi }}_i =( \phi _{i 1}, \phi _{i 2}, \ldots ) ^{\mathrm {T}}$. With proper truncations on the expansion (A13), a Galerkin solution could be derived in terms of the eigenfunction 
$f_i$.
Now that the eigenvalue 
$\lambda _i$ and eigenfunction 
$f_i$ have been found, one may naturally follow the seminal framework of Barton (Reference Barton1983) to expand the local moments of the p.d.f. for active particles as
\begin{equation} P_n (r, \theta, \phi, t) = \sum_{i = 1}^{\infty} p_{n i} (t) \exp({-\lambda_i t})\,f_i (r, \theta, \phi), \quad n=0,1,\ldots, \end{equation}
where 
$p_{n i}$ gives the corresponding expansion coefficient. On the method of moments for solute dispersion, Barton (Reference Barton1983, § 3) resolved the technical difficulties of deriving the expansion coefficients 
$p_{n i}$ with elements in the bilinear form.
In contrast with the case of passive particles, the local operator 
$\mathcal {L}$ related with the boundary conditions may be non-self-adjoint under the combined action of swimming and rotational motion of active particles. Moreover, 
$\boldsymbol{\mathsf{A}}$ can be indeed non-symmetric so that eigenvalues and eigenvectors are complicated. Indeed, one will fail to apply directly the integral transform method of Barton for the non-orthogonal set of functions 
$\{\,f_i\}_{i=1}^{\infty }$, that is, the inner product 
$\langle \,f_i, f_j\rangle \neq 0$ for arbitrary 
$i \neq j$.
From the above discussion, we finally resort to a bi-orthogonality relation with 
$\{\,f_i\}_{i=1}^{\infty }$ (Strand et al. Reference Strand, Kim and Karrila1987; Brezinski Reference Brezinski1991; Jiang & Chen Reference Jiang and Chen2021). We denote the dual counterpart with the superscript ‘
$\star$’ and deduce a bi-orthogonality relation
\begin{equation} \langle \,f^{{\star}}_i, f_j\rangle = \delta_{i j},\end{equation}
where 
$\{\,f^{\star }_i\}_{i=1}^{\infty }$ is the corresponding eigenfunction of the adjoint operator 
$\mathcal {L}^{\star }$), and 
$\delta$ is the Kronecker delta.
In the light of the Galerkin method with the base 
$\{e_i\}_{i=1}^{\infty }$, we can solve 
$f^{\star }_i$ as
\begin{equation} \boldsymbol{\mathsf{A}}^{{\star}} {\boldsymbol{\phi}}_i^{{\star}} = \lambda_i {\boldsymbol{\phi}}_i^{{\star}}, \end{equation}
where 
$\boldsymbol{\mathsf{A}}^{\star }$ represents the corresponding matrix of 
$\mathcal {L}^{\star }$ (which equals the transpose of 
$\boldsymbol{\mathsf{A}}$), and 
${\boldsymbol {\phi }}_i^{\star }$ is the coefficient vector of 
$f_i^{\star }$. We apply the series expansion with the same basis by analogy with (A13):
\begin{equation} f^{{\star}}_i = \sum_{i = 1}^{\infty} \phi^{{\star}}_{i j} e_j, \end{equation}
and 
${\boldsymbol {\phi }}_i^{\star } =(\phi ^{\star }_{i 1}, \phi ^{\star }_{i 2}, \ldots )^{\mathrm {T}}$. The eigenvalues of 
$\boldsymbol{\mathsf{A}}^{\star }$ are consistent with those of 
$\boldsymbol{\mathsf{A}}$ (Strand et al. Reference Strand, Kim and Karrila1987; Jiang & Chen Reference Jiang and Chen2021). With the bi-orthogonal group 
$\{\,f_i\}_{i=1}^{\infty }$ and 
$\{\,f^{\star }_i\}_{i=1}^{\infty }$, we can now utilise the general expressions of Barton (Reference Barton1983) for expansion coefficients of moments in (A17) using the bi-orthogonal relation (A18).
To close the problem, we hereinafter derive the specific basis functions 
$\{e_i\}_{i=1}^{\infty }$ satisfying the reflective boundary conditions of moments. The eigenfunctions of the Laplace operator in a spherical coordinate system read
\begin{equation} \varDelta_q = \frac{1}{r}\,\frac{\partial}{\partial r} \left( r\, \frac{\partial}{\partial r} \right) + \frac{1}{\sin \theta}\, \frac{\partial}{\partial \theta} \left( \sin \theta\, \frac{\partial}{\partial \theta} \right) + \frac{1}{\sin^2 \theta}\,\frac{\partial^2}{\partial \phi^2}. \end{equation}
The eigenfunctions 
$\varDelta _q$ can be devised with the polyadic surface spherical harmonics and Bessel functions (Brenner Reference Brenner1964)
$$\begin{gather} R^{c}_{0}\,\sqrt{\frac{2l + 1}{4 {\rm \pi}} } \mathrm{P}_l (\cos \theta), \end{gather}$$
$$\begin{gather} R^{c}_{m}\,\sqrt{2}\,\sqrt{\frac{2l + 1}{4 {\rm \pi}}\,\frac{(l - m) !}{(l + m) !}} \cos (m \phi)\,\mathrm{P}_l^m (\cos \theta), \quad m = 1, 2,\ldots, l, \end{gather}$$
$$\begin{gather} R^{c}_{m}\,\sqrt{2}\,\sqrt{\frac{2l + 1}{4 {\rm \pi}}\,\frac{(l - m)!}{(l + m) !}} \sin (m \phi)\,\mathrm{P}_l^m (\cos \theta)\, \quad m = 1, 2, \ldots, l. \end{gather}$$
Here, 
$l = 0, 1, 2, \ldots$, 
$\mathrm {P}_l$ is the Legendre polynomial, and 
$\mathrm {P}^m_l$ are the associated Legendre polynomials (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010)
\begin{equation} \mathrm{P}^m_l (x) = (- 1)^m (1 - x^2)^{m / 2}\, \frac{\mathrm{d}^m}{\mathrm{d}\kern 0.06em x^m} \mathrm{P}_l (x), \end{equation}
with the Condon–Shortley phase 
$(- 1)^m$. The solutions of
\begin{equation} R_{n}^{{c}}(r)=\sqrt{2}\,\frac{\mathrm{J}_{0}\left(\sqrt{-\beta_{n}}\,r\right)}{\mathrm{J}_{0} \left(\sqrt{-\beta_{n}}\right)}, \quad n=0,1,2, \ldots, \end{equation}and
\begin{equation} R^{s}_{n}=\sqrt{2}\,\frac{\mathrm{J}_{0}\left(\sqrt{-\gamma_{n}}\,r\right)}{\mathrm{J}_{1} \left(\sqrt{-\gamma_{n}}\right)}, \quad n=1,2, \ldots, \end{equation}
are consistent with Jiang & Chen (Reference Jiang and Chen2020). Here, 
$\mathrm {J}_{0}$ is the first-kind Bessel function, with 
$\beta _n$ denoting the 
$n$th-order non-negative zeros of
\begin{equation} \frac{\mathrm{d}\mathrm{J}_{0}}{\mathrm{~d} r}\left(\sqrt{-\beta_{n}}\right)= \mathrm{J}_{1}\left(\sqrt{-\beta_{n}}\right)=0, \quad n=0,1,2, \ldots. \end{equation}
Note that 
$\beta _0=0$ and 
$R_{0}^{{c}}(r)=\sqrt {2}$. Here, 
$\gamma _n$ are the 
$n$th-order zeros of 
$\mathrm {J}_{0}(\sqrt {-\gamma _{n}})=0$ for 
$n=1,2, \ldots$. The reflection bases are constructed with the sequence of spherical harmonics and Bessel functions (A22)–(A24) in the light of (A26) and (A27), i.e. orthonormal with the inner product (A15).
Appendix B. Lagrangian particle simulation of three-dimensional micro-swimmers through a pipe
We conduct Lagrangian particle simulation of three-dimensional micro-swimmers in vertical pipe flows. The equivalent stochastic differential equations corresponding to the fundamental Smoluchowski equation (Thorn & Bearon Reference Thorn and Bearon2010) in a Cartesian coordinate system read
$$\begin{gather} \mathrm{d}\kern0.06em \boldsymbol{x} = \boldsymbol{U}(\boldsymbol{x})\,\mathrm{d} t+V_{s} \boldsymbol{p}\,\mathrm{d} t, \end{gather}$$
$$\begin{gather}\mathrm{d}\boldsymbol{p} = \left\{\frac{1}{2 B}\,[\boldsymbol{k}-(\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{p}) \boldsymbol{p}]+\frac{1}{2}\,\boldsymbol{\omega} \wedge \boldsymbol{p}\right\} \mathrm{d} t+\sqrt{2 D_{r}}\,\mathrm{d} \boldsymbol{W}, \end{gather}$$
with the rotational diffusion denoted by 
$\sqrt {2 D_{r}}\,\mathrm {d} \boldsymbol {W}$. (B1) and (B2) can be discretised numerically in the light of the two-step Heun method.
For a two-dimensional case, one may take it for granted that the distribution of 
$\theta$ near 
$\theta =0$ follows a Gaussian distribution (
$\phi$ is uniformly distributed) that satisfies
\begin{equation} \frac{\partial P}{\partial t} = D_r\,\frac{\partial^2 P}{\partial \theta^2}, \end{equation}
and thus 
$\langle \theta ^2\rangle =2 D_r t$. However, this is not exactly true for three-dimensional micro-swimmers. In fact, for Brownian motion on a spherical surface, the correct equation should be
\begin{equation} \frac{\partial P}{\partial t} = D_r\,\frac{1}{\sin \theta}\, \frac{\partial}{\partial \theta} \left( \sin \theta\,\frac{\partial P}{\partial \theta} \right). \end{equation}According to Brillinger (Reference Brillinger1997, (3.5)) and Chirikjian (Reference Chirikjian2009), the corresponding stochastic differential equation is
\begin{equation} \mathrm{d} \theta = \frac{\sigma^2}{2 \tan \theta}\,\mathrm{d} t + \sigma\, \mathrm{d} W, \end{equation}
where 
$\sigma ^2=2D_r$.
The stochastic noise of orientation diffusion vectors should be rotated corresponding to the real-time direction of 
$\boldsymbol {p}$, as pointed out by Thorn & Bearon (Reference Thorn and Bearon2010). Note that the rotational Brownian motion can be approximated by a two-dimensional diffusion process in a plane (Berg Reference Berg1993, p. 84):
\begin{equation} \theta^2=\theta_x^2+\theta_y^2. \end{equation}
Note that both 
$\theta _x$ and 
$\theta _y$ follow 
$N(0, 2 D_r\,{\rm \Delta} t)$. Here, 
$\theta \sim \mathrm {Rayleigh} (\sigma )$ is Rayleigh distributed if 
$\theta ={\sqrt {\theta _x^{2}+\theta _y^{2}}}$, where 
$\theta _x \sim N(0,\sigma ^{2})$ and 
$\theta _y \sim N(0,\sigma ^{2})$ are independent normal random variables. This gives motivation to denote the symbol ‘
$\sigma$’ in the above parametrisation of the Rayleigh distribution. Rigorously, 
$\theta$ in fact follows the Rayleigh distribution (
$\theta \sim \mathrm {Rayleigh} (\sigma )$).
Reflecting boundary conditions are implemented as (Volpe et al. Reference Volpe, Gigan and Volpe2014; Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016)
\begin{equation} r_n \rightarrow 2-r_n,\quad \theta_n \rightarrow {\rm \pi}-\theta_n,\quad \phi_n \rightarrow 2{\rm \pi}-\phi_n,\quad \text{if } r_n >1. \end{equation}The previous and prediction information of swimming direction and positional distribution is stored for the purpose of reflections. When the swimmers escape from the pipe boundary, the particle position is reflected on the tangent line at the contact point of the circle, and the new orientation distribution is denoted by the normalised vector joining the intact point and final position. In particular, we observe the exceptional cases of multiple reflections within a single random walk. The positions of particles for each step should be checked to ensure that they are finally inside the periphery of the pipe.
As shown in figure 13, this Lagrangian particle simulation collaborates perfectly with our solutions in the light of the bi-orthogonal expansion method. That is, by series expansions of spherical harmonics and Bessel functions, we have conducted appropriate computation on the macro-scale transport coefficients and local distributions.
Figure 13. Comparison of Lagrangian particle simulation (symbols) and theoretical (lines) first moment 
$M_1$, and MSD 
$\sigma ^2$, with red and blue for non-gyrotactic (
$\lambda =0$) and gyrotactic (
$\lambda =1$) active particles, respectively. In this downwelling Poiseuille flow, pointwise swimmers are initially oriented vertically upwards and displaced uniformly at the origin of the pipe. A total of 
$10^5$ particles are released for each run, with time step 
$10^{-3}$. Common parameters: 
${Pe} _f=1$, 
${Pe}_s=1$ and 
$D_t=0$.
