3.2.1 Reflection basis functions

One can use the eigenfunctions of the Laplace operator,

(3.6)$$\begin{eqnarray}\unicode[STIX]{x1D6E5}_{q}=\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\right)+\frac{1}{\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\left(\sin \unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)+\frac{1}{\sin ^{2}\unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{2}},\end{eqnarray}$$

on the local space to find the required basis functions. The inner product for given functions $f$ and $g$ is

(3.7)$$\begin{eqnarray}\langle f,g\rangle \triangleq \int _{0}^{1}r\text{d}r\int _{0}^{\unicode[STIX]{x03C0}}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D719}\,f(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})g(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}).\end{eqnarray}$$

Note that $\unicode[STIX]{x1D6E5}_{q}$ is self-adjoint with the reflective, finite and periodic boundary conditions (2.18)–(2.24). One needs to check the radial term. Integration by parts gives

$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{0}^{1}r\,\text{d}r\int _{0}^{\unicode[STIX]{x03C0}}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D719}\,g\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}r}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =\left.\int _{0}^{\unicode[STIX]{x03C0}}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D719}\left(gr\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}r}-rf\frac{\unicode[STIX]{x2202}g}{\unicode[STIX]{x2202}r}\right)\right|_{r=0}^{r=1}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\int _{0}^{1}r\,\text{d}r\int _{0}^{\unicode[STIX]{x03C0}}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D719}\,f\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}g}{\unicode[STIX]{x2202}r}\right).\nonumber\end{eqnarray}$$

The reflection conditions (2.18) and (2.19) ensure that all the boundary terms are equal to zero. Namely,

(3.8)$$\begin{eqnarray}\int _{0}^{\unicode[STIX]{x03C0}}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D719}\left.\left(rg\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}r}\right)\right|_{r=1}=\int _{0}^{\unicode[STIX]{x03C0}}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D719}\left.\left(rf\frac{\unicode[STIX]{x2202}g}{\unicode[STIX]{x2202}r}\right)\right|_{r=1}=0\end{eqnarray}$$

because both $(\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}r(1,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}))$ and $(\unicode[STIX]{x2202}g/\unicode[STIX]{x2202}r(1,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}))$ are odd functions with respect to $\unicode[STIX]{x1D719}$.

To construct the eigenfunctions for $\unicode[STIX]{x1D6E5}_{q}$, one can use the spherical harmonics,

(3.9)$$\begin{eqnarray}\displaystyle & \displaystyle R^{\text{c}}(r)\sqrt{\frac{2l+1}{4\unicode[STIX]{x03C0}}}\text{P}_{l}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(3.10)$$\begin{eqnarray}\displaystyle & \displaystyle R^{\text{c}}(r)\sqrt{2}\sqrt{\frac{2l+1}{4\unicode[STIX]{x03C0}}\frac{(l-m)!}{(l+m)!}}\cos (m\unicode[STIX]{x1D719})\text{P}_{l}^{m}(\cos \unicode[STIX]{x1D703}),\quad m=1,2,\ldots ,l, & \displaystyle\end{eqnarray}$$

(3.11)$$\begin{eqnarray}\displaystyle & \displaystyle R^{\text{s}}(r)\sqrt{2}\sqrt{\frac{2l+1}{4\unicode[STIX]{x03C0}}\frac{(l-m)!}{(l+m)!}}\sin (m\unicode[STIX]{x1D719})\text{P}_{l}^{m}(\cos \unicode[STIX]{x1D703}),\quad m=1,2,\ldots ,l, & \displaystyle\end{eqnarray}$$

where $l=0,1,2,\ldots ,$ $\text{P}_{l}$ are the Legendre polynomials and $\text{P}_{l}^{m}$ are the associated Legendre polynomials (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010),

(3.12)$$\begin{eqnarray}\text{P}_{l}^{m}(x)=(-1)^{m}(1-x^{2})^{m/2}\frac{\text{d}^{m}}{\text{d}x^{m}}\text{P}_{l}(x),\end{eqnarray}$$

with the Condon–Shortley phase $(-1)^{m}$.

The undetermined functions $R^{\text{c}}(r)$ and $R^{\text{s}}(r)$ satisfy

(3.13)$$\begin{eqnarray}\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}R(r)}{\unicode[STIX]{x2202}r}\right)-l(l+1)R(r)=\unicode[STIX]{x1D712}R(r),\end{eqnarray}$$

where $\unicode[STIX]{x1D712}$ is the eigenvalue for $\unicode[STIX]{x1D6E5}_{q}$. For $R^{\text{c}}(r)$, the reflection conditions (2.18) and (2.19) require

(3.14)$$\begin{eqnarray}\displaystyle & \displaystyle \cos (m\unicode[STIX]{x1D719})R^{\text{c}}(1)=\cos (-m\unicode[STIX]{x1D719})R^{\text{c}}(1),\quad \forall ~\unicode[STIX]{x1D719}\in [0,2\unicode[STIX]{x03C0}], & \displaystyle\end{eqnarray}$$

(3.15)$$\begin{eqnarray}\displaystyle & \displaystyle \cos (m\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}R^{\text{c}}(1)}{\unicode[STIX]{x2202}r}=-\cos (-m\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}R^{\text{c}}(1)}{\unicode[STIX]{x2202}r},\quad \forall ~\unicode[STIX]{x1D719}\in [0,2\unicode[STIX]{x03C0}], & \displaystyle\end{eqnarray}$$

namely,

(3.16)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}R^{\text{c}}(1)}{\unicode[STIX]{x2202}r}=0.\end{eqnarray}$$

Therefore, the solution for $R^{\text{c}}(r)$ is

(3.17)$$\begin{eqnarray}R_{n}^{\text{c}}(r)=\sqrt{2}\,\frac{\text{J}_{0}(\sqrt{-\unicode[STIX]{x1D6FD}_{n}}r)}{\text{J}_{0}(\sqrt{-\unicode[STIX]{x1D6FD}_{n}})},\quad n=0,1,2,\ldots ,\end{eqnarray}$$

where $\text{J}_{0}$ is the Bessel function of the first kind and $\unicode[STIX]{x1D6FD}_{n}$ is the $n$th non-negative zero of

(3.18)$$\begin{eqnarray}\frac{\text{d}\text{J}_{0}}{\text{d}r}(\sqrt{-\unicode[STIX]{x1D6FD}_{n}})=\text{J}_{1}(\sqrt{-\unicode[STIX]{x1D6FD}_{n}})=0,\quad n=0,1,2,\ldots .\end{eqnarray}$$

Notice that $\unicode[STIX]{x1D6FD}_{0}=0$ and $R_{0}^{\text{c}}(r)=\sqrt{2}$. The corresponding eigenvalue is

(3.19)$$\begin{eqnarray}\unicode[STIX]{x1D712}_{nl}^{\text{c}}=-l(l+1)-\unicode[STIX]{x1D6FD}_{n}^{2},\quad n=0,1,2,\ldots .\end{eqnarray}$$

Similarly, for $R^{\text{s}}(r)$, the reflection conditions (2.18) and (2.19) require

(3.20)$$\begin{eqnarray}\displaystyle & \displaystyle \sin (m\unicode[STIX]{x1D719})R^{\text{s}}(1)=\sin (-m\unicode[STIX]{x1D719})R^{\text{s}}(1),\quad \forall ~\unicode[STIX]{x1D719}\in [0,2\unicode[STIX]{x03C0}], & \displaystyle\end{eqnarray}$$

(3.21)$$\begin{eqnarray}\displaystyle & \displaystyle \sin (m\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}R^{\text{s}}(1)}{\unicode[STIX]{x2202}r}=-\sin (-m\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}R^{\text{s}}(1)}{\unicode[STIX]{x2202}r},\quad \forall ~\unicode[STIX]{x1D719}\in [0,2\unicode[STIX]{x03C0}], & \displaystyle\end{eqnarray}$$

namely,

(3.22)$$\begin{eqnarray}R^{\text{s}}(1)=0.\end{eqnarray}$$

Therefore, the solution for $R^{\text{s}}(r)$ is

(3.23)$$\begin{eqnarray}R^{\text{s}}(r)=\sqrt{2}\,\frac{\text{J}_{0}(\sqrt{-\unicode[STIX]{x1D6FE}_{n}}r)}{\text{J}_{1}(\sqrt{-\unicode[STIX]{x1D6FE}_{n}})},\quad n=1,2,\ldots ,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}_{n}$ is the $n$th zero of

(3.24)$$\begin{eqnarray}\text{J}_{0}(\sqrt{-\unicode[STIX]{x1D6FE}_{n}})=0,\quad n=1,2,\ldots .\end{eqnarray}$$

Note that $R^{\text{s}}(r)$ is not orthogonal to $R^{\text{c}}(r)$ with respect to $r$, namely

(3.25)$$\begin{eqnarray}\int _{0}^{1}R^{\text{s}}(r)R^{\text{c}}(r)\,r\,\text{d}r\neq 0.\end{eqnarray}$$

The corresponding eigenvalue is

(3.26)$$\begin{eqnarray}\unicode[STIX]{x1D712}_{nl}^{\text{s}}=-l(l+1)-\unicode[STIX]{x1D6FE}_{n}^{2},\quad n=1,2,\ldots .\end{eqnarray}$$

Note that the above sequence of spherical harmonics (3.9), (3.10) and (3.11) with (3.17) and (3.23) is orthonormal with respect to the inner product (3.7). It forms a basis for functions on the local space satisfying the reflective, finite and periodic boundary conditions (2.18)–(2.24). Thus, we call it a ‘reflection basis’.

3.2.2 Galerkin method

Now we apply the Galerkin method with the above reflection basis to obtain the solution of $P_{0}^{\infty }$. We use $\{{e_{i}^{\text{r}}\}}_{i=1}^{\infty }$ to denote the basis; then the expression for $P_{0}^{\infty }$ is

(3.27)$$\begin{eqnarray}P_{0}^{\infty }(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\mathop{\sum }_{i=1}^{\infty }q_{i}e_{i}^{\text{r}}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}),\end{eqnarray}$$

where $q_{i}$ is the expansion coefficient to be determined. Note that the BVP for $P_{0}^{\infty }$ is symmetric with respect to the line $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$,

(3.28)$$\begin{eqnarray}P_{0}^{\infty }(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=P_{0}^{\infty }(r,\unicode[STIX]{x03C0}-\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}).\end{eqnarray}$$

Therefore, only even-symmetric associated Legendre polynomials are used in the expansion, i.e. with $l+m$ even.

To obtain the Galerkin equation, i.e. the weak formulation of (3.5) by inner product with the reflection basis $\{{e_{i}^{\text{r}}\}}_{i=1}^{\infty }$, one needs to construct the bilinear form $a(\cdot ,\cdot )$ for the local operator ${\mathcal{L}}$. In matrix form, the element of the corresponding matrix is

(3.29)$$\begin{eqnarray}\unicode[STIX]{x1D608}_{ij}=a(e_{i}^{\text{r}},e_{j}^{\text{r}})=\langle e_{i}^{\text{r}},{\mathcal{L}}e_{j}^{\text{r}}\rangle ,\quad i=1,2,\ldots ,~j=1,2,\ldots .\end{eqnarray}$$

Because $P_{0}^{\infty }$ is independent of $\unicode[STIX]{x1D713}$, equation (3.5) reduces to

(3.30)$$\begin{eqnarray}-\mathit{Pe}_{s}\sin \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719}\frac{1}{r}\frac{\unicode[STIX]{x2202}(rP_{0}^{\infty })}{\unicode[STIX]{x2202}r}+\unicode[STIX]{x1D735}_{p}\boldsymbol{\cdot }(\dot{\boldsymbol{p}}P_{0}^{\infty })=D_{t}\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}P_{0}^{\infty }}{\unicode[STIX]{x2202}r}\right)+(\unicode[STIX]{x1D735}_{p}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{p})P_{0}^{\infty }.\end{eqnarray}$$

The diffusive term (the right-hand side of (3.30)) is Laplacian, and thus it is easy to treat with the reflection basis, resulting in a diagonal matrix for the corresponding bilinear form.

However, the ‘convective’ term (the left-hand side of (3.30)) is laborious to handle. To simplify the bilinear form with spherical harmonics, recursion relations and identities of associated Legendre polynomials can be imposed (Strand & Kim Reference Strand and Kim1992; Bees, Hill & Pedley Reference Bees, Hill and Pedley1998; Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003), but it still leads to lengthy expressions and tedious calculations. We follow Doi & Edwards (Reference Doi and Edwards1978) and treat the local operator ${\mathcal{L}}$ in a more systematic way. The angular momentum operators in quantum mechanics are introduced and related to ${\mathcal{L}}$. Integrals of the spherical harmonics are given by the Wigner $3j$-symbols. More details are shown in appendix A. Finally, truncation of the series (3.27) to some degree $N$ gives a Galerkin solution for $P_{0}^{\infty }$.

4.1.1 Distribution in phase space

Examples for the gyrotactic focusing of cells in vertical pipe flows are shown in figure 2. For the case of a downwelling flow, cells accumulate near the centre of the tube (Kessler Reference Kessler and Velarde1984, Reference Kessler1985), with the highest concentration at the centre $(r=0)$. The distributions of cells in the $\unicode[STIX]{x1D719}{-}\unicode[STIX]{x1D703}$ plane at all the radial positions are similar: cells concentrate near the point $\unicode[STIX]{x1D719}=0$, $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$. Note that there is no obvious swimming bias towards the centre ($0<\unicode[STIX]{x1D719}<\unicode[STIX]{x03C0}$), which is significantly different from the biased distribution predicted by the FP equation in both the two-step methods (Bearon et al. Reference Bearon, Hazel and Thorn2011). For the upwelling flow, the results are just reversed: the concentration at the wall ($r=1$) is the highest (Kessler Reference Kessler and Velarde1984, Reference Kessler1985).

Figure 2.

Density plot of local distributions $P_{0}^{\infty }(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ of spherical and ellipsoidal gyrotactic cells at different radial positions: (a,d,g,j) $r=0$; (b,e,h,k) $r=0.5$; (c,f,i,l) $r=1$. (a–c) Spheres in downwelling flow ($\unicode[STIX]{x1D6FC}_{0}=0$); (d–f) ellipsoids in downwelling flow ($\unicode[STIX]{x1D6FC}_{0}=0.31$); (g–i) spheres in upwelling flow ($\unicode[STIX]{x1D6FC}_{0}=0$); (j–l) ellipsoids in upwelling flow ($\unicode[STIX]{x1D6FC}_{0}=0.31$). In all cases, $\unicode[STIX]{x1D706}=2.2$, $\mathit{Pe}_{s}=1$ and $w_{V}=1$.

The structure of the local distribution in the phase space can be captured by the dynamical system theory (Zöttl & Stark Reference Zöttl and Stark2012, Reference Zöttl and Stark2013; Santamaria et al. Reference Santamaria, De Lillo, Cencini and Boffetta2014). A local velocity field for the local transport problem in the phase space $(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ can be defined as (Jiang & Chen Reference Jiang and Chen2019a)

(4.4)$$\begin{eqnarray}\boldsymbol{u}=(-\mathit{Pe}_{s}\sin \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719},\dot{\unicode[STIX]{x1D703}},\dot{\unicode[STIX]{x1D719}}).\end{eqnarray}$$

Swimmers will accumulate near stable critical points (with zero local velocity). For vertical downwelling pipe flow, the critical point $(0,\unicode[STIX]{x03C0}/2,0)$ is stable, where cells swim perfectly upwards, resulting in gyrotactic focusing. In contrast, the critical point $(0,\unicode[STIX]{x03C0}/2,\unicode[STIX]{x03C0})$ is unstable, where cells swim downwards, resulting in an area with the lowest concentration. For the case of the upwelling flow, the cells accumulate near the reflective boundaries due to swinging motions towards the tube wall, similar to the trapping by tumbling motions near high-shear regions (Zöttl & Stark Reference Zöttl and Stark2013; Rusconi, Guasto & Stocker Reference Rusconi, Guasto and Stocker2014). Additionally, for the case of spherical cells without gyrotactic bias ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=0$), note that the local velocity field is divergence-free, as shown in § A.2.1. Therefore, the local transport problem (3.30) with reflective boundary conditions (2.18) and (2.19) has a constant solution $P_{0}^{\infty }=1/(4\unicode[STIX]{x03C0}^{2})$. There is no stable critical point and thus no focusing. Note that the incident probability flux at the wall ($r=1,-\unicode[STIX]{x03C0}<\unicode[STIX]{x1D719}<0$) can be non-zero but is equal to the reflection flux ($0<\unicode[STIX]{x1D719}<\unicode[STIX]{x03C0}$) because the reflective boundary conditions (2.18) and (2.19) for the continuum transport equation (2.13) are a simplified description of the reflection process, assuming that the incidence and reflection of a cell happen in an instant.

For the influence of cell shape, as shown in figure 2(a,d) for the downwelling flow, the main difference between the cases of spheres ($\unicode[STIX]{x1D6FC}_{0}=0$) and ellipsoids ($\unicode[STIX]{x1D6FC}_{0}=0.31$) is due to the effect of shear alignments. The concentrated area of ellipsoids in figure 2(d) is much smaller than that of spheres in figure 2(a) because the swimming direction of prolate cells is aligned closely along the streamlines of the flow by the strain effect. However, the focusing of ellipsoids in the radial direction is weaker than that of spheres. Shear alignments reduce the probability that prolate cells can swim towards the centre of the tube, resulting in a more uniform distribution in the radial direction. The discussion for the upwelling case is similar.

4.1.2 Orientation-space-mean probability density function

The gyrotactic focusing of cells near the centre of the tube can be illustrated more clearly in terms of the orientation-space-mean p.d.f.

(4.5)$$\begin{eqnarray}\langle P_{0}^{\infty }\rangle _{O}\triangleq \int _{0}^{\unicode[STIX]{x03C0}}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D719}\,P_{0}^{\infty }(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}),\end{eqnarray}$$

i.e. the marginal p.d.f. in the position plane $r$–$\unicode[STIX]{x1D713}$ (here independent of $\unicode[STIX]{x1D713}$). With the orientation space eliminated, results are compared to those by the two-step models for the transport problem in the position space (Bearon et al. Reference Bearon, Bees and Croze2012).

First, we discuss the case of the downwelling flow. As shown in figure 3, the gyrotactic focusing near the axis of the tube is enhanced with increased flow rates (also quantified by the ratio between speeds $w_{V}$ with a fixed swimming speed). However, the distributions of the spherical cells in figure 3(a–c) are roughly the same. Namely, the focusing of spherical cells is nearly dependent only on $w_{V}$, the ratio between the Péclet numbers. Note that the steady gyrotactic focusing can be predicted by a Gaussian distribution (Kessler Reference Kessler and Velarde1984, Reference Kessler1985; Pedley & Kessler Reference Pedley and Kessler1992). Near the centre of the tube with sufficiently small shear, Bearon et al. (Reference Bearon, Bees and Croze2012) gave an asymptotic solution based on the two-step GTD method, in current notation,

(4.6)$$\begin{eqnarray}\langle P_{0}^{\infty }\rangle _{O}\sim \text{e}^{-\unicode[STIX]{x1D706}w_{V}r^{2}},\end{eqnarray}$$

which depends only on $w_{V}$ and $\unicode[STIX]{x1D706}$. We remark that this approximate result is hard to obtain directly from the local transport equation (3.30) in the whole phase space by asymptotic analysis, but can be simply derived from the effective transport equation only in the position space by the two-step method; see Bearon et al. (Reference Bearon, Bees and Croze2012, equation (24)). In the absence of gyrotactic bias ($\unicode[STIX]{x1D706}=0$), there is no focusing and the radial distribution is uniform.

For the influence of cell shape, shear alignments of ellipsoids ($\unicode[STIX]{x1D6FC}_{0}=0.31$) can weaken the gyrotactic focusing, which is consistent with previous studies (Bearon et al. Reference Bearon, Hazel and Thorn2011; Croze et al. Reference Croze, Sardina, Ahmed, Bees and Brandt2013) for channel flows. Compared with the spherical case ($\unicode[STIX]{x1D6FC}_{0}=0$), distributions of ellipsoids are much more uniform in the radial direction. More importantly, the decrease of concentration near the centre of the tube from the spherical case to the ellipsoidal case is greater with both larger flow rates (quantified by $w_{V}$) and larger swimming Péclet numbers $\mathit{Pe}_{s}$, as shown in figure 3. The reason is that shear alignments by the strain effect can inhibit prolate cells from swimming towards the centre of the tube, as discussed in § 4.1.1. The increase of both $\mathit{Pe}_{s}$ and $w_{V}$ corresponds to a larger flow Péclet number $\mathit{Pe}_{f}$ due to the definition of $w_{V}$ (4.3), as the ratio between the Péclet numbers. Thus, the shear rate and the strain also increase, enhancing the inhibiting effect.

Figure 3.

Orientation-space-mean p.d.f.s $\langle P_{0}^{\infty }\rangle _{O}$ of cells in the downwelling flow for different swimming Péclet numbers $\mathit{Pe}_{s}$: (a,d,g) $\mathit{Pe}_{s}=0.1$; (b,e,h) $\mathit{Pe}_{s}=1$; (c,f,i) $\mathit{Pe}_{s}=10$; and different ratios between speeds $w_{V}$: (a–c) $w_{V}=0.1$; (d–f) $w_{V}=1$; (g–i) $w_{V}=10$. ‘Spherical’ denotes the case of gyrotactic spherical cells ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=2.2$). ‘Ellipsoidal’ denotes ellipsoidal cells ($\unicode[STIX]{x1D6FC}_{0}=0.31$, $\unicode[STIX]{x1D706}=2.2$). ‘No gyrotaxis’ denotes spherical cells without gyrotactic bias ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=0$). ‘Two-step, GTD’ and ‘Two-step, PK’ denote the results of spherical cells ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=2.2$) by the two-step method with the GTD and PK models, respectively. ‘Asymptotic’ denotes the associated asymptotic Gaussian distribution (4.6).

For the case of spherical gyrotactic cells ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=2.2$), both the two-step method with the GTD model and that with the PK model (Bearon et al. Reference Bearon, Bees and Croze2012) give successful predictions with small swimming Péclet numbers $\mathit{Pe}_{s}$ and relative variations of shear rates of the Poiseuille flow (quantified by $w_{V}$). As confirmed in figure 3(a,b,d), these three curves are on top of each other. However, outside the required parameter range in the $\mathit{Pe}_{s}{-}w_{V}$ plane, as shown in figure 3(e,f,h,i), there are considerable differences of the results by the PK model from the exact results by the one-step method. With stronger shear, more uniform distributions are predicted in the PK model than in GTD models, due to neglecting the local spatial distribution (thus valid only for weak flows) (Bearon et al. Reference Bearon, Hazel and Thorn2011, Reference Bearon, Bees and Croze2012; Croze et al. Reference Croze, Bearon and Bees2017). On the other hand, the predictions given by the two-step GTD method are still in agreement with the exact results, even in figure 3(i) with large $\mathit{Pe}_{s}$ and $w_{V}$. To some extent, the gyrotactic focusing accumulates cells in the regions with low shear and far from the tube wall, which enlarges the appropriate parameter range for the two-step method (Bearon et al. Reference Bearon, Hazel and Thorn2011).

Figure 4.

Orientation-space-mean p.d.f.s $\langle P_{0}^{\infty }\rangle _{O}$ of cells in the upwelling flow for different swimming Péclet numbers $\mathit{Pe}_{s}$: (a,d,g) $\mathit{Pe}_{s}=0.1$; (b,e,h) $\mathit{Pe}_{s}=1$; (c,f,i) $\mathit{Pe}_{s}=10$; and different ratios between speeds $w_{V}$: (a–c) $w_{V}=0.1$; (d–f) $w_{V}=1$; (g–i) $w_{V}=10$. ‘Spherical’ denotes the case of gyrotactic spherical cells ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=2.2$). ‘Ellipsoidal’ denotes ellipsoidal cells ($\unicode[STIX]{x1D6FC}_{0}=0.31$, $\unicode[STIX]{x1D706}=2.2$). ‘No gyrotaxis’ denotes spherical cells without gyrotactic bias ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=0$). ‘Two-step, GTD’ and ‘Two-step, PK’ denote the results for spherical cells ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=2.2$) by the two-step method with the GTD and PK models, respectively. ‘Asymptotic’ denotes the associated asymptotic Gaussian distribution (4.6) with $\unicode[STIX]{x1D706}=-2.2$.

Now we turn to the case of the upwelling flow. The distinct difference is that the gyrotactic focusing is near the walls (Kessler Reference Kessler and Velarde1984, Reference Kessler1985), as shown in figure 4. As in the upwelling case, the orientation-space-mean p.d.f. can still be matched by the asymptotic solution (4.6) with the reversed sign of shear, although the shear near the wall is much larger than that near the centre of the tube. When the shear rate is large, shown in figure 4(g–i), the focusing is so intense that higher-order basis functions in the radial direction are needed for the series expansion (3.27) to reduce the truncation error. Also, distributions of spherical cells are nearly independent of swimming Péclet number $\mathit{Pe}_{s}$ with a fixed ratio between speeds $w_{V}$, but it is not true for the ellipsoids. The shear alignments can also flatten the distributions in the radial direction. The discussion about the predictions by the two-step methods is analogous: the GTD model gives more successful results than the PK model, even inside the ‘breakdown’ parameter region.

4.2.1 Influence of shear flow

First, we discuss the influence of shear on the overall drift of spherical cells. Note that $U_{d}$ (3.34) is the drift above the mean flow because the frame of reference has been transformed in (2.11). With the swimming Péclet number $\mathit{Pe}_{s}$ fixed, changing the ratio between speeds $w_{V}$ in figure 5 corresponds to changing the flow rate and shear rate.

The drift increases monotonically with the flow rate, as shown in figure 5(a,c,e). At small flow rates, where the swimming effect is dominant, the drift is negative due to the upward orientation of the cells by the gravitational torque (Bees & Croze Reference Bees and Croze2010; Croze et al. Reference Croze, Bearon and Bees2017). As the flow rate increases, the drift becomes positive because the gyrotactic focusing accumulates cells near the centre of the tube where the flow speed $\mathit{Pe}_{f}$ is the fastest. Therefore, at high flow rates, it is observed that $U_{d}\sim \mathit{Pe}_{f}=w_{V}\mathit{Pe}_{s}$, which is consistent with the discussion by Croze et al. (Reference Croze, Bearon and Bees2017). In the absence of gyrotactic bias ($\unicode[STIX]{x1D706}=0$), there is no drift above the mean flow because the local distribution is uniform in the whole phase space.

Next, for the dispersivity, its variation with flow rate is much more complex. Note that the dispersion process is driven by both the swimming diffusion effect and the convection (Jiang & Chen Reference Jiang and Chen2019a). When the swimming Péclet number $\mathit{Pe}_{s}$ is not large, the convection effect plays the central role. As shown in figure 5(b,d), the dispersivity rises as the flow rate increases. However, at high flow rates, the dispersivity appears to ‘saturate’ and then falls with $w_{V}$ (Croze et al. Reference Croze, Sardina, Ahmed, Bees and Brandt2013, Reference Croze, Bearon and Bees2017). The reason is that the gyrotactic focusing of cells near the centre of the tube greatly hinders the effect of convection on the dispersion process. Therefore, without gyrotactic bias ($\unicode[STIX]{x1D706}=0$), the dispersivity increases monotonically with the flow rate and is much larger than that of gyrotactic cells.

When the swimming Péclet number $\mathit{Pe}_{s}$ is large, the strong swimming diffusion drives the macrotransport process. As shown in figure 5(f), there are significant differences in the trends of dispersivity with flow rates compared with those in figure 5(b,d): the dispersivity will first decrease, then it increases and saturates like the convection-effect-dominated cases. More importantly, in the absence of gyrotaxis, there is also a sharp decrease in dispersivity before the rising. Namely, the imposed shear can reduce the swimming diffusion effect. To find the reason, note that the swimming diffusion effect can be quantified by $V_{s}^{2}/D_{r}$ (Bearon et al. Reference Bearon, Hazel and Thorn2011). The imposed shear will enhance the rotation of cells, and thus reduce the swimming range, resulting in the decrease of the dispersivity.

Additionally, gyrotaxis can also reduce the swimming diffusion effect. As shown in figure 5(d,f) at small flow rates, i.e. the swimming-effect-dominated region, it results in a much smaller dispersivity compared with the no-bias case. Because gyrotaxis forces cells to swim upwards together, the tendency of the cells to separate from each other by swimming with rotation diffusion is weakened. However, the perfectly upward orientation of cells can be disturbed by the imposed shear due to the enhanced rotation, resulting in an increase in dispersivity, in contrast to the decrease discussed above. In fact, these enhanced and reduced effects of the shear rotation on the dispersion process will be balanced: the dispersivity of the gyrotactic case is comparable to that without the bias, as shown in figure 5(b) when $w_{V}<2$.

Figure 5.

For the downwelling flow, variations of (a,c,e) overall drift $U_{d}$ (3.34) and (b,d,f) dispersivity $D_{T}$ (3.35) with the ratio between speeds $w_{V}$ subject to different fixed swimming Péclet numbers $\mathit{Pe}_{s}$: (a,b) $\mathit{Pe}_{s}=0.1$; (c,d) $\mathit{Pe}_{s}=1$; (e,f) $\mathit{Pe}_{s}=10$. ‘Spherical’ denotes the case of gyrotactic spherical cells ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=2.2$). ‘Ellipsoidal’ denotes ellipsoidal cells ($\unicode[STIX]{x1D6FC}_{0}=0.31$, $\unicode[STIX]{x1D706}=2.2$). ‘No gyrotaxis’ denotes spherical cells without gyrotactic bias ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=0$). ‘Two-step, GTD’ and ‘Two-step, PK’ denote the results for spherical cells ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=2.2$) by the two-step method with the GTD and PK models, respectively. Note that with $\mathit{Pe}_{s}$ fixed, increasing $w_{V}$ corresponds to increasing the mean speed $U_{m}^{\ast }$.

With respect to the influence of cell shape, shear alignments of ellipsoids ($\unicode[STIX]{x1D6FC}_{0}=0.31$) can slightly reduce the drift and greatly enhance the dispersivity, because the gyrotactic focusing is weakened (Croze et al. Reference Croze, Sardina, Ahmed, Bees and Brandt2013). As discussed in § 4.1.1, the local distribution of ellipsoids is more uniform in the radial direction, compared with that of spheres. Therefore, more cells swim outside the fast flow zone, decreasing the drift above the mean flow. Meanwhile, the convection effect on the dispersion process is strengthened.

The predictions by the two-step methods on the dispersion process are compared with the results by the one-step method for the case of spherical gyrotactic cells. When the swimming Péclet number $\mathit{Pe}_{s}$ is small, as required in the two-step methods, both the GTD and PK models provide perfect results for the drift and dispersivity as expected, also shown in figure 5(a,b), although there is still a small deviation in dispersivity at high flow rates (quantified by $w_{V}$). With larger $\mathit{Pe}_{s}$, considerable discrepancies are observed in figure 5(c–f). Instead of growing with the flow rate, the drift by the PK model tends to be a constant much smaller than the exact one at high flow rates (Croze et al. Reference Croze, Sardina, Ahmed, Bees and Brandt2013, Reference Croze, Bearon and Bees2017). As discussed in § 4.1.1, the PK model neglects the local spatial distribution by shear and predicts a more uniform distribution of cells in the radial direction, resulting in a smaller drift similar to the case without gyrotactic bias. In the GTD model, the curve of drift still adheres to the exact one, but curves of dispersivity diverge. The deviation in dispersivity becomes larger as the variation of shear rates in the radial direction (also quantified by the ratio between Péclet numbers $w_{V}$ for the Poiseuille flow) increases and exceeds the requirements of the first step in obtaining the convection–diffusion equation in the position space.

4.2.2 Influence of boundaries

We have tested the predictions produced by the two-step methods on the dispersion process with the variation of shear rates $w_{V}$ at different swimming Péclet numbers $\mathit{Pe}_{s}$. Here, we give a more direct illustration of the influence of $\mathit{Pe}_{s}$, with $w_{V}$ fixed. Note that for C. nivalis with fixed swimming velocity $V_{s}^{\ast }$ and rotational diffusivity $D_{r}^{\ast }$, increasing $\mathit{Pe}_{s}$ corresponds to decreasing the pipe radius $a^{\ast }$, and thus increasing the influence of boundaries.

Changing $\mathit{Pe}_{s}$ hardly affects the gyrotactic focusing, as discussed in § 4.1.1. However, it significantly influences the drift and dispersivity as functions of $w_{V}$, as shown in figure 5. With a fixed ratio between Péclet numbers $w_{V}$, i.e. the relative importance of the shear effect to the swimming effect, increasing $\mathit{Pe}_{s}$ will enhance both these effects. Therefore, the drift and dispersivity vary monotonically with $\mathit{Pe}_{s}$, for both cases of spherical and ellipsoidal cells, as shown in figure 6.

To obtain the convection–diffusion equation in the position space in the two-step methods, a small $\mathit{Pe}_{s}$ is the key requirement, which means that the influence of boundaries can be neglected (Bearon et al. Reference Bearon, Hazel and Thorn2011). As shown in figure 6, increasing $\mathit{Pe}_{s}$ is destructive to the validity of the two-step methods, for both PK and GTD models. The discrepancies in the drift and dispersivity between the two-step methods and the exact ones grow considerably with $\mathit{Pe}_{s}$. Note that, although the reflective boundary conditions (2.18) and (2.19) considered in the current study are an idealized condition for the population models, it is demonstrated that the existence of boundaries will lead to significant deviations in the two-step methods. In fact, when boundaries play a central role, more complicated effects, e.g. hydrodynamic interactions between cells and walls (Berke et al. Reference Berke, Turner, Berg and Lauga2008), should be taken into consideration in both the one-step and two-step methods for the dispersion process.

Figure 6.

For downwelling flow, variations of (a,c,e) overall drift $U_{d}$ (3.34) and (b,d,f) dispersivity $D_{T}$ (3.35) with the swimming Péclet number $\mathit{Pe}_{s}$ subject to different fixed relative variations of shear rates (quantified by the ratio between Péclet numbers $w_{V}$): (a,b) $w_{V}=0.1$; (c,d) $w_{V}=1$; (e,f) $w_{V}=10$. ‘Spherical’ denotes the case of gyrotactic spherical cells ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=2.2$). ‘Ellipsoidal’ denotes ellipsoidal cells ($\unicode[STIX]{x1D6FC}_{0}=0.31$, $\unicode[STIX]{x1D706}=2.2$). ‘No gyrotaxis’ denotes spherical cells without gyrotactic bias ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=0$). ‘Two-step, GTD’ and ‘Two-step, PK’ denote the results of spherical cells ($\unicode[STIX]{x1D6FC}_{0}=0$, $\unicode[STIX]{x1D706}=2.2$) by the two-step method with the GTD and PK models, respectively. Note that for the studied case, increasing $\mathit{Pe}_{s}$ corresponds to decreasing the pipe radius $a^{\ast }$.