1. Introduction
In a seminal effort, Taylor (Reference Taylor1953) examined the axial dispersion of a solute in pipe Poiseuille flow. Using a physically motivated version of the multiple scales analysis employed here (see § 2), Taylor showed that, for sufficiently long times, the cross-sectionally averaged solute concentration, while being convected with the mean speed (
$U_m/2$
), also spreads diffusively along the flow direction with an effective diffusivity given by
$D({{\it{Pe}}^2}/{192})$
; here,
$D$
is the diffusivity of the solute molecules, and
${\it{Pe}}=U_mR/D$
is the Péclet number based on the pipe radius
$R$
and centreline speed
$U_m$
. In a later effort, Aris (Reference Aris1956) derived the equations governing the temporal evolution of the longitudinal moments of the cross-sectionally averaged concentration profile, from the convection–diffusion equation governing the solute concentration field. Based on the relation between the second- and higher-order moments, he showed that the solute concentration profile, as a function of the flow coordinate, approached a Gaussian form for long times. Further, the profile variance increased linearly in time, with this increase driven by an effective diffusivity
$D_{\it{eff}}=D [1+({{\it{Pe}}^2}/{192}) ]$
. Thus, the solute diffusivity
$D$
was shown by Aris (Reference Aris1956) to contribute additively to Taylor’s original result. For laminar pressure-driven channel flow, corresponding to the plane Poiseuille velocity profile, a similar analysis gives
$D_{\it{eff}}=D [1+({8Pe^2}/{945}) ]$
, where
${\it{Pe}}=U_mH/D$
, with
$H$
being the channel half-width (Wooding Reference Wooding1960; Dorfman & Brenner Reference Dorfman and Brenner2001).
The
$O(Pe^2)$
shear-enhanced contribution to the effective longitudinal diffusivity, in the above examples, arises from an interplay of the background shearing flow and Brownian diffusion. It may be understood via scaling arguments for the flow-aligned component of the variance,
$\langle x^2 \rangle - \langle x\rangle ^2$
, of a Brownian particle (the solute molecule) subject to the pertinent shearing flow. Considering the channel flow case, for instance, a particle starting at the centreline (gradient coordinate
$z=0$
), is convected by the spatially varying component of the flow,
$u\!\sim \!U_m({z^2}/{H^2})$
, for a time
$t$
, leading to
$x\!\sim \!(U_mt)( {z^2}/{H^2})$
, so
$\langle x^2 \rangle \!\sim \!({(U_mt)^2}/{H^4})\langle z^4 \rangle$
. Using
$\langle z^4 \rangle \!\sim \!D^2t^2$
for the diffusively spreading Gaussian profile in the gradient direction, gives
$\langle x^2 \rangle \!\sim \!({(U_m D)^2}/{H^4}) t^4$
(Latini & Bernoff Reference Latini and Bernoff2001). This superdiffusive growth changes to a diffusive one when the particle has sampled the largest velocity variations of
$O(U_m)$
for
$t \sim H^2/D$
. Writing
$t^4$
as
$t^3\times t$
in
$\langle x^2 \rangle$
above, replacing
$t^3$
by
$H^6/D^3$
, and comparing with
$\langle x^2 \rangle - \langle x\rangle ^2 \sim D_{\it{eff}}t$
, gives the
${\it{Pe}}$
-dependent term in
$D_{\it{eff}}$
in the previous paragraph. As mentioned above, the solute diffusivity
$D$
enters as an independent additive contribution (Aris Reference Aris1956). Note that, unlike a linear flow (Foister & Van De Ven Reference Foister and Van De Ven1980; Subramanian & Brady Reference Subramanian and Brady2004), the mean position of a particle starting at the centreline also evolves with time, in a reference frame moving with
$U_m$
. However, arguments analogous to the one above give
$\langle x \rangle \sim -({U_m D}/{H^2}) t^2$
for short times, implying that the variance,
$\langle x^2 \rangle -\langle x \rangle ^2$
, scales in the same manner as
$\langle x^2 \rangle$
.
Moving beyond passive solutes, in this paper, we examine the Taylor dispersion of anisotropic microswimmers in pressure-driven channel flow. Understanding the spread of microbial populations subject to an ambient flow is of fundamental interest in both the natural habitats of microswimmers (Stamm Reference Stamm1991; Kurz et al. Reference Kurz, Secchi, Stocker and Jimenez-Martinez2023; Markale et al. Reference Markale, Carrel, Kurz, Morales, Holzner and Jiménez-Martínez2023), as well as for biotechnological engineering and applications (Fernández et al. Reference Fernández, Camacho and Chisti1999; Reid Reference Reid1999; Donlan Reference Donlan2001; Anderson et al. Reference Anderson, Clarke, Arkin and Voigt2006; Scott et al. Reference Scott, Davey, Dennis, Horst, Howe, Lea-Smith and Smith2010; Croze et al. Reference Croze, Sardina, Ahmed, Bees and Brandt2013; Felfoul et al. Reference Felfoul2016). Consequently, the interplay between swimming and a background shearing flow has been extensively studied in recent years; see, for instance, the review by Wheeler et al. (Reference Wheeler, Secchi, Rusconi and Stocker2019).
The solution of the active dispersion problem is rendered non-trivial by the range of behaviour exhibited by microswimmers even in simple unidirectional shearing flows (Durham, Kessler & Stocker Reference Durham, Kessler and Stocker2009; Zöttl & Stark Reference Zöttl and Stark2013; Rusconi, Guasto & Stocker Reference Rusconi, Guasto and Stocker2014; Berman et al. Reference Berman, Ferguson, Bizzak, Solomon and Mitchell2022; Vennamneni et al. Reference Vennamneni, Nambiar and Subramanian2020b
). While the passive solute in the channel flow example above spreads uniformly across the cross-section for
$t\!\gt \!H^2/D$
, non-spherical swimmers, depending on
${\it{{\it{Pe}}}}_r$
and
$\kappa$
, exhibit either high-shear trapping (HST) due to a wallward drift (Durham et al. Reference Durham, Kessler and Stocker2009; Rusconi et al. Reference Rusconi, Guasto and Stocker2014; Bearon & Hazel Reference Bearon and Hazel2015; Ezhilan & Saintillan Reference Ezhilan and Saintillan2015; Vennamneni et al. Reference Vennamneni, Nambiar and Subramanian2020b
), or low-shear trapping (LST) due to a drift towards the channel centreline (Barry et al. Reference Barry, Rusconi, Guasto and Stocker2015; Vennamneni et al. Reference Vennamneni, Nambiar and Subramanian2020b
; Berman et al. Reference Berman, Ferguson, Bizzak, Solomon and Mitchell2022; Borgnino et al. Reference Borgnino, Boffetta, Cencini, De Lillo and Gustavsson2022; Ganesh et al. Reference Ganesh, Douarche, Dentz and Auradou2023). Here,
$\kappa$
is the swimmer aspect ratio and
${\it{{\it{Pe}}}}_r= U_m/(HD_r)$
is the rotary Péclet number measuring the relative importance of shear-induced rotation at the rate
$U_m/H$
, and relaxation due to rotary diffusion with diffusivity
$D_r$
, for active Brownian swimmers (ABPs); ABPs constitute one of two canonical swimmer types, the other being run-and-tumble particles (RTP). For large
${\it{{\it{Pe}}}}_r$
and slender swimmers (
$\kappa \gg 1$
), the aforementioned trapping regimes are demarcated by
${\it{{\it{Pe}}}}_r \approx 0.07\kappa ^3$
, with LST occurring for
${\it{{\it{Pe}}}}_r \gt 0.07\kappa ^3$
, see figure 1(a) (Vennamneni et al. Reference Vennamneni, Nambiar and Subramanian2020b
). It is worth adding that an anisotropic passive solute has also been predicted to exhibit mild HST in a plane Poiseuille flow (Nitsche & Hinch Reference Nitsche and Hinch1997), although this has not been observed in experiments (Strednak et al. Reference Strednak, Shaikh, Butler and Guazzelli2018).
(a) Migration diagram for ABPs in channel flow (Vennamneni et al. Reference Vennamneni, Nambiar and Subramanian2020b
). (b) Schematic of the active Taylor dispersion problem;
$U_s D_r^{-1}$
is the persistence length associated with an ABP, with
$U_s$
being the swimming speed and
$D_r$
the rotary diffusivity.
Of particular significance to the active dispersion problem is the onset of a ‘collapse’ for swimmers with
$\kappa \gtrsim 2$
, within the LST regime, first identified in Vennamneni et al. (Reference Vennamneni, Nambiar and Subramanian2020b
). As shown in figure 1(a), for such swimmers, the drift towards the centreline is strong enough that the swimmer concentration approaches a delta function, localized at the centreline, for
${\it{{\it{Pe}}}}_r \rightarrow \infty$
. The delta-function limiting form is an artefact of the assumption underlying the multiple time scales framework in Vennamneni et al. (Reference Vennamneni, Nambiar and Subramanian2020b
), which regards the swimmer mean free path as being negligibly small – real swimmers are expected to localize in a finite thickness ‘Knudsen layer’ with a thickness of
$O(U_s/D_r)$
or smaller,
$U_s$
being the swimming speed. Nevertheless, the spatial localization accompanying the onset of the centreline-collapse regime restricts the velocity variations sampled by the swimmers with increasing
${\it{{\it{Pe}}}}_r$
, potentially leading to reduced dispersion. Herein, we map out classical and anomalous scaling regimes for the dispersion of prolate spheroidal swimmers (
$\kappa \gt 1$
) in pressure-driven channel flow. ‘Classical’ is used here in a broader sense to refer to both the passive solute, and active swimmers where the
$D_{\it{eff}}$
-scaling may be obtained from the passive one, merely by replacing the solute diffusivity with the
${\it{{\it{Pe}}}}_r$
-dependent swimmer translational diffusivity in the scaling arguments above. ‘Anomalous’ refers to cases where centreline collapse leads instead to relatively rare tail events controlling the large-
${\it{{\it{Pe}}}}_r$
dispersion.
The outline of the paper is as follows. In § 2, we present the mathematical framework for deriving the longitudinal diffusivity that determines the linear-in-time growth of the variance of a microswimmer population of a finite spatial extent, in plane Poiseuille flow, for sufficiently long times. Section 2.1 introduces the kinetic equation that governs the time evolution of the swimmer probability density in position–orientation space; the swimmers (i.e. ABPs) being assumed to execute active Brownian motion. Following this, the method of multiple scales is used to obtain a formal expression for the longitudinal diffusivity (
$D_{\it{eff}}$
), as a function of
${\it{{\it{Pe}}}}_r$
and
$\kappa$
, in sufficiently wide channels. The multiple scales analysis is spread across §§ 2.2.1 and 2.2.2. In the former subsection, the small ratio of the swimmer mean-free-path to the channel width (the swimmer Knudsen number
$\epsilon$
) allows one to separate the time scale characterizing orientation relaxation, from that pertaining to drift-cum-diffusion in physical space. This leads to a two-dimensional drift–diffusion equation governing the swimmer concentration that is a function of the flow and gradient coordinates. In § 2.2.2, we further separate the slow evolution along the gradient direction, from the slower one along the flow direction, leading to a one-dimensional diffusion equation for the swimmer concentration, as a function of the flow coordinate, along with the closed form expression for the diffusivity
${D}_{\it{eff}}$
valid for
$\epsilon \ll 1$
,
${\it{Pe}}\gg 1$
. In § 2.2.3, we derive the generalized longitudinal diffusion expression
$D_{\it{eff}}$
valid for arbitrary
${\it{Pe}}$
or
${\it{{\it{Pe}}}}_r$
. The results are presented in § 3. Section 3.1 focuses on the nature of the swimmer concentration profiles resulting from migration along the gradient direction, and characterizes the transition from high-shear to LST with increasing
${\it{{\it{Pe}}}}_r$
. Section 3.2 presents the results for
$D_{\it{eff}}$
, obtained from a numerical evaluation of the above expression, as a function of
$\kappa$
and
${\it{{\it{Pe}}}}_r$
. In all cases, the computed results are supported by analysis and/or scaling arguments. The emphasis is on the anomalous scaling regime of
$D_{\it{eff}}$
, which arises due to the centreline collapse of swimmers with
$\kappa \gtrsim 2$
. We present our main conclusions in § 4, while pointing to the presence of anomalous scaling in a recently published simulation effort (Ganesh et al. Reference Ganesh, Douarche, Dentz and Auradou2023). The scope for future extensions of the present effort is also discussed. Appendix A contains numerical details regarding the evaluation of
$D_{\it{eff}}$
, pertaining to computing both the underlying swimmer orientation distributions and concentration profiles; herein, we also present results for
$D_{\it{eff}}$
, as a function of
${\it{{\it{Pe}}}}_r$
, for swimmers below the centreline-collapse threshold. Appendix B contains the results of Langevin simulations that have been primarily used to validate the
$D_{\it{eff}}$
values for spherical swimmers, obtained from the multiple scales analysis, over the entire range of
${\it{{\it{Pe}}}}_r$
.
2. Mathematical formulation
In § 2.1, we begin by writing down the equation governing the time evolution of the swimmer probability density in position–orientation space, and defining all the terms that appear therein. In § 2.2, we employ a multiple scales analysis to derive the effective longitudinal diffusivity that characterizes the rate of spread of a microswimmer population in a wide channel; ‘wide’ here corresponds to the limit of the channel width being much larger than the swimmer mean free path (
$U_s/D_r \ll H$
in figure 1
b).
2.1. Governing equations
The governing equation for the time dependent probability density of the position (
$\boldsymbol{x}$
) and orientation (
$\boldsymbol{p}$
) of an ABP,
$\varOmega (\boldsymbol{x},\boldsymbol{p},t)$
, is given by (Subramanian & Koch Reference Subramanian and Koch2009; Saintillan & Shelley Reference Saintillan and Shelley2015; Vennamneni et al. Reference Vennamneni, Garg and Subramanian2020a
,
Reference Vennamneni, Nambiar and Subramanianb
)
\begin{equation} \frac {\partial \varOmega }{\partial t} + \boldsymbol{\nabla} _{\!\boldsymbol{x}} \!\boldsymbol{\cdot }\! [( U_{s}\boldsymbol{p}+\boldsymbol{u} ) \varOmega ]+\boldsymbol{\nabla}_{\!p} \!\boldsymbol{\cdot }\! (\dot {\boldsymbol{p}}\varOmega ) -D_r {\nabla} ^{2}_{p}\varOmega =0. \end{equation}
The probability density
$\varOmega$
is convected in physical space both due to intrinsic swimming with speed
$U_s$
along the swimmer orientation, and by the ambient plane Poiseuille flow defined by
$\boldsymbol{u} = U_m[1-(z/H)^2]\boldsymbol{1}_1$
. Here,
$\boldsymbol{1}_1$
is the unit vector along the flow direction (the
$x$
-axis), with
$z\in [-H,H]$
being the transverse coordinate; as mentioned in § 1,
$H$
is the channel half-width and
$U_m$
the maximum speed at the channel centreline. The third term in (2.1) denotes the convection of
$\varOmega$
in orientation space due to rotation by the ambient shear, with the rotation rate given by the Jeffery equation,
$\dot {\boldsymbol{p}}=B [ \boldsymbol{E}\!\boldsymbol{\cdot }\! \boldsymbol{p}\!-\!\boldsymbol{p}(\boldsymbol{E}\!:\!\boldsymbol{p}\boldsymbol{p}) ]+\boldsymbol{W}\!\boldsymbol{\cdot }\!\boldsymbol{p}$
,
$B = (\kappa ^2-1)/(\kappa ^2+1)$
(Jeffery Reference Jeffery1922; Leal & Hinch Reference Leal and Hinch1971) being the Bretherton constant for the spheroidal swimmers of aspect ratio
$\kappa$
examined here. Here
$\boldsymbol{E} = -z(U_m/H)(\boldsymbol{1}_1\boldsymbol{1}_3 + \boldsymbol{1}_3\boldsymbol{1}_1)$
,
$\boldsymbol{W} = -z(U_m/H)(\boldsymbol{1}_1\boldsymbol{1}_3 - \boldsymbol{1}_3\boldsymbol{1}_1)$
represent the rate-of-strain and vorticity tensors of the plane Poiseuille flow, with
$\boldsymbol{1}_3$
being the unit vector in the gradient direction. The final term in (2.1) denotes the relaxation of
$\varOmega$
to isotropy (
$\varOmega = 1/4\pi$
) due to rotary diffusion,
$D_r$
being the rotary diffusivity.
To implement the multiple scales analysis in § 2.2 below, we find it convenient to choose different length scales along the gradient and flow directions. While
$H$
is evidently the relevant scale in the gradient direction, we use
$L$
to denote the characteristic size of the swimmer population along the flow direction. There is an arbitrariness in the choice of
$L$
– for a channel of an infinite length, for instance,
$L$
may be made as large as possible, in order for the separation of the gradient and flow-directed time evolutions to occur. A physical significance may nevertheless be attributed to the minimum estimate of this longitudinal scale, denoted
$L_{\it{min}}$
, by enforcing that, on scales larger than
$L_{\it{min}}$
, the evolution along the flow direction becomes slower than that along the gradient direction. An estimate for
$L_{\it{min}}$
may then be obtained from equating the corresponding diffusion times,
$H^2/D \sim L_{\it{min}}^2/D_{\it{eff}}$
. While
$L_{\it{min}}$
for active swimmers will have to await the expression for
$D_{\it{eff}}$
obtained below, that for the passive solute may be obtained using
$D_{\it{eff}} \sim Pe^2 D$
for
${\it{Pe}} \gg 1$
, as given in the introductory paragraph, which gives
$L_{\it{min}} \sim Pe.H$
. Thus, while a stronger ambient shear does lead to a larger effective diffusivity along the flow direction, the onset of this enhanced diffusion regime takes a commensurately longer time. On scales larger than
$L_{\it{min}}$
, the evolution of the solute or swimmer concentration can be described using a one-dimensional diffusion equation with diffusivity
$D_{\it{eff}}$
, in a reference frame moving with the mean speed. In cases where the transiently evolving concentration field, as a function of the transverse (gradient or radial) coordinate, may be obtained in closed form, a more precise estimate of
$L_{\it{min}}$
is obtained from the time required to approach a steady state. This is how Taylor (Reference Taylor1953) obtained an estimate for
$L_{\it{min}}$
, for pipe Poiseuille flow. The estimate is the same as the one above, but contains an additional numerical prefactor related to the zero of the Bessel function that characterizes the radial variation of the solute concentration (see (16) in Taylor (Reference Taylor1953)).
Given the above preamble on the choice of length scales, we use
$D_r^{-1}$
as the characteristic time scale, to obtain the kinetic equation in the following non-dimensional form:
\begin{equation} \frac {\partial \varOmega }{\partial t} +\varLambda ( \epsilon \,p_1 +{\it{{\it{Pe}}}}_r u_1)\frac {\partial \varOmega }{\partial x} + \epsilon \,p_3 \frac {\partial \varOmega }{\partial z}+{\it{{\it{Pe}}}}_r \dot {\gamma } \boldsymbol{\nabla}_{\!p} \!\boldsymbol{\cdot }\! (\tilde {\dot {\boldsymbol{p}}}\varOmega ) - {\nabla} ^{2}_{p}\varOmega =0, \end{equation}
where
$\tilde {\dot {\boldsymbol{p}}}=\dot {\boldsymbol{p}}/({\rm d}u/{\rm d}z)$
represents the non-dimensional swimmer rotation rate and
$\dot {\gamma }(z)=-2z$
is the linearly varying shear rate associated with the plane Poiseuille profile. In (2.2),
$\epsilon =U_{s}/HD_r$
and
$\varLambda \epsilon =U_{s}/LD_r$
denote the swimmer Knudsen numbers with
$H$
and
$L$
as the choices for the macroscopic scale. Here
${\it{{\it{Pe}}}}_r = U_m/HD_r$
is the rotary Péclet number introduced in § 1, and governs the relative magnitudes of shear-induced rotation and stochastic relaxation due to rotary diffusion, with
${\it{{\it{Pe}}}}_r \gg 1$
denoting the shear-dominant limit. In (2.2),
$\varOmega$
satisfies the following normalization constraint:
\begin{equation} \int {\rm d}\boldsymbol{p} {\rm d}x {\rm d}z\, \varOmega =1, \end{equation}
along with no-flux conditions at the channel walls (
$z = \pm 1$
) which may be written as
\begin{equation} \int {\rm d}\boldsymbol{p}\, (U_{s}\boldsymbol{p}\varOmega ) \boldsymbol{\cdot }\boldsymbol{\hat {n}}=0, \end{equation}
where
$\boldsymbol{\hat {n}} = \pm \boldsymbol{1}_z$
is the unit normal to the channel walls.
2.2. Multiple scales analysis for wide channels
The multiple scales analysis is performed in two steps. In § 2.2.1, we separate the fast time scale associated with orientation relaxation due to rotary diffusion, from the slow time scale that characterizes spatial drift-cum-diffusion in the flow and gradient directions. Next, in § 2.2.2, adopting a reference frame moving with the mean speed of the swimmer population, we separate the spatial time scales into a slow time scale characterizing drift-cum-diffusion in the gradient direction, and a slower one characterizing diffusion along the flow direction. In involving positional degrees of freedom alone, the analysis in § 2.2.2 bears a resemblance to the original analysis for a passive solute devoid of internal degrees of freedom.
2.2.1. Separating the orientation and physical space time scales
For
$\epsilon \ll 1$
, the orientation distribution relaxes to a quasisteady state commensurate with the local shear rate, in a time of
$O(D_r^{-1})$
, while the spatial redistribution of swimmers requires a time of
$O(H^2/D_t)$
or longer. The latter time may equivalently be written as
$O(D_r^{-1}\epsilon ^{-2})$
, on using the (nominal) estimate
$D_t \sim U_s^2/D_r$
for the swimmer translational diffusivity. We therefore define
$t_1 = t$
and
$t_2 = \epsilon ^2 t$
as the fast and slow time variables, and rewrite the time derivative in the form
\begin{equation} \frac {\partial }{\partial t} = \frac {\partial }{\partial t_1} + \epsilon ^2 \frac {\partial }{\partial t_2}, \end{equation}
as per the multiple scales formalism (Subramanian & Brady Reference Subramanian and Brady2004; Kevorkian & Cole Reference Kevorkian and Cole2012), while expanding the swimmer probability density as
\begin{align} \varOmega (\boldsymbol{x},\boldsymbol{p},t)\!&=\varOmega ^{(0)}(x,z,\boldsymbol{p},t_1,t_2;{\it{{\it{Pe}}}}_r) + \epsilon \, \varOmega ^{(1)}(x,z,\boldsymbol{p},t_1,t_2;{\it{{\it{Pe}}}}_r) \nonumber \\ &\quad +\epsilon ^2\varOmega ^{(2)}(x,z,\boldsymbol{p},t_1,t_2;{\it{{\it{Pe}}}}_r) + O(\epsilon ^3). \end{align}
Substituting (2.5) and (2.6) into the governing kinetic equation (2.2), one obtains the following equations at successive orders in
$\epsilon$
:
\begin{align} {O}(1)\frac {\partial \varOmega ^{(0)}}{\partial t_1} + {\it{{\it{Pe}}}}_r \dot {\gamma } \boldsymbol{\nabla}_{\!p} \!\boldsymbol{\cdot }\! (\tilde {\dot {\boldsymbol{p}}}\varOmega _{0}) - {\nabla} ^{2}_{p}\varOmega ^{(0)}&=0, \end{align}
\begin{align} {O}(\epsilon )\frac {\partial \varOmega ^{(1)}}{\partial t_1} + {\it{{\it{Pe}}}}_r \dot {\gamma } \boldsymbol{\nabla}_{\!p} \!\boldsymbol{\cdot }\! \big(\tilde {\dot {\boldsymbol{p}}}\varOmega ^{(1)}\big)- {\nabla} ^{2}_{p}\varOmega ^{(1)}&=- p_3 \frac {\partial \varOmega ^{(0)}}{\partial z}, \end{align}
\begin{align} {O}\big(\epsilon ^2\big)\frac {\partial \varOmega ^{(2)}}{\partial t_1} + {\it{{\it{Pe}}}}_r \dot {\gamma } \boldsymbol{\nabla}_{\!p} \!\boldsymbol{\cdot } \big(\tilde {\dot {\boldsymbol{p}}}\varOmega ^{(2)}\big)- {\nabla} ^{2}_{p}\varOmega ^{(2)}&=- p_3 \frac {\partial \varOmega ^{(1)}}{\partial z} - (\varLambda /\epsilon ^2) {\it{{\it{Pe}}}}_{r}\,u_1 \frac {\partial \varOmega ^{(0)}}{\partial x}\nonumber\\&\quad - \frac {\partial \varOmega ^{(0)}}{\partial t_2}. \end{align}
In including convection by the plane Poiseuille flow only in (2.9), at
$O(\epsilon ^2)$
, we have assumed
$\varLambda /\epsilon ^2 \lesssim O(1)$
. That this is true for
${\it{{\it{Pe}}}}_r$
fixed will be shown a posteriori. The above choice for
$\varLambda$
implies that swimming-induced convection along
$z$
appears at
$O(\epsilon )$
, while that along
$x$
appears only at
$O(\epsilon ^3)$
, and therefore, is not included in (2.7)–(2.9).
The absence of spatial derivatives in (2.7) points to the
$\boldsymbol{x}$
and
$\boldsymbol{p}$
dependencies in
$\varOmega ^{(0)}$
being separable, so that
\begin{equation} \varOmega ^{(0)}= I(x,z,t_2)G(\boldsymbol{p},t_1;{\it{{\it{Pe}}}}_r), \end{equation}
with
$G(\boldsymbol{p},t_1;{\it{{\it{Pe}}}}_r)$
for
$t_1 \gg 1$
being the quasisteady orientation distribution commensurate with the local shear rate, achieved for dimensional times much greater than
$O(D_r^{-1})$
. The function
$I(x,z,t_2)$
denotes the swimmer concentration field that varies along the flow and gradient directions, with spatial homogeneity being assumed along the spanwise (vorticity) direction. The orientation distribution
$G$
satisfies
\begin{equation} {\it{{\it{Pe}}}}_r \dot {\gamma } \boldsymbol{\nabla}_{\!p} \!\boldsymbol{\cdot }\! (\tilde {\dot {\boldsymbol{p}}}G) - {\nabla} ^{2}_{p}G=0, \end{equation}
for large
$t_1$
, corresponding to the quasisteady state above. Without loss of generality, one may impose the constraint
\begin{equation} \int {\rm d}\boldsymbol{p}\,G(\boldsymbol{p},t_1;{\it{{\it{Pe}}}}_r)= 1. \end{equation}
Use of (2.6), (2.10) and (2.12) in (2.3) gives
\begin{equation} \int {\rm d}x {\rm d}z \,I = 1, \end{equation}
at leading order, and
\begin{equation} \int {\rm d}\boldsymbol{p}\,\varOmega ^{(1)}= \int {\rm d}\boldsymbol{p}\,\varOmega ^{(2)} = 0, \end{equation}
at
$O(\epsilon )$
and
$O(\epsilon ^2)$
, respectively. Substituting the separable form of
$\varOmega ^{(0)}$
from (2.10) into (2.8), one obtains
\begin{equation} \frac {\partial \varOmega ^{(1)}}{\partial t_1} + {\it{{\it{Pe}}}}_r \dot {\gamma } \boldsymbol{\nabla}_{\!p} \!\boldsymbol{\cdot }\! (\tilde {\dot {\boldsymbol{p}}}\varOmega ^{(1)})- {\nabla} ^{2}_{p}\varOmega ^{(1)}= - p_3 G \frac {\partial I}{\partial z}- p_3\frac {\partial G}{\partial z}I, \end{equation}
where
$G$
depends on
$z$
only via the
$z$
-dependent shear rate of the plane Poiseuille flow. On account of linearity, one may write
\begin{equation} \varOmega ^{(1)} = \varOmega _{11} \frac {\partial I}{\partial z} + \varOmega _{12} I, \end{equation}
where
$\varOmega _{11}$
and
$\varOmega _{12}$
satisfy equations similar to (2.15), but with the individual forcing functions appearing on the right-hand side instead (
$-p_3 G$
for
$\varOmega _{11}$
and
$-p_3({\partial G}/{\partial z})$
for
$\varOmega _{12}$
). The
$O(\epsilon )$
normalization constraint in (2.14) implies
\begin{equation} \int {\rm d}\boldsymbol{p} \varOmega _{11} = \int {\rm d}\boldsymbol{p} \varOmega _{12} = 0. \end{equation}
Using the above form for
$\varOmega ^{(1)}$
in (2.9), and integrating the resulting equation over orientation space leads to the following drift–diffusion equation governing
$I$
:
\begin{align} \frac {\partial I}{\partial t_2} + \varLambda {\it{{\it{Pe}}}}\,u_{1} \frac {\partial I}{\partial x} = \frac {\partial }{\partial z}\left (D_{zz} \frac {\partial I}{\partial z} -V_{z}I\right )\!, \end{align}
where we have replaced
${\it{{\it{Pe}}}}_{r}/\epsilon ^2$
, in the convective term on the left-hand side by
${\it{Pe}}$
, this being the active analogue of the solute Péclet number based on the swimmer translational diffusivity (
$D_t \sim U_s^2/D_r$
). The drift and diffusion coefficients on the right-hand side of (2.18) are given by
\begin{align} V_z(z;{\it{{\it{Pe}}}}_r,\kappa ) &= \displaystyle \int p_3\varOmega _{12}\,{\rm d}\boldsymbol{p}, \end{align}
\begin{align} D_{zz}(z;{\it{{\it{Pe}}}}_r,\kappa ) &= - \displaystyle \int p_3\varOmega _{11}\,{\rm d}\boldsymbol{p}, \end{align}
in non-dimensional form, and as indicated by the arguments, are functions of both
${\it{{\it{Pe}}}}_r$
and
$\kappa$
.
2.2.2. Separating the time scales corresponding to spatial diffusion along the gradient and flow directions
Prior to separating the physical space time scales, we transform to a reference frame moving with the long-time average speed of the swimmer population. Denoting this mean speed as
$\bar {u}_1$
, defining a new flow coordinate
$\bar {x} = x - \bar {u}_1t_2$
and writing
$u_1 = \bar {u}_1 + u_1^{\prime}$
,
$u_1^{\prime}$
being the fluctuating component of the streamwise velocity, (2.18) takes the form
\begin{equation} \frac {\partial I}{\partial t_2} =- \varLambda Pe \,u^{\prime }_{1} \frac {\partial I}{\partial \bar {x}}+ \frac {\partial }{\partial z}\left (D_{zz} \frac {\partial I}{\partial z} - V_z I\right )\!, \end{equation}
where the mean convection has now been subsumed by the chosen reference frame. It will be shown below, in § 2.2.3, that
$\varLambda Pe \ll 1$
provided
${\it{{\it{Pe}}}}_r$
is large. In this limit, one may again follow the multiple scales framework, and write the time derivative in (2.21) as
\begin{equation} \frac {\partial I}{\partial t_2}= \frac {\partial I}{\partial t_2} + \varLambda ^2 \frac {\partial I}{\partial t_3}, \end{equation}
where
$\varLambda$
rather than
$\varLambda Pe$
is used as the small parameter for convenience. In (2.22),
$t_2$
is the slow time variable defined in § 2.2.1, while
$t_3 = \varLambda ^2t_2$
is the slower time variable pertaining to evolution along the flow direction. Expanding
$I$
in the form
\begin{equation} I(\bar {x},z,t_2;\varLambda ) = I^{(0)}(\bar {x},z,t_2,t_3) + \varLambda I^{(1)}(\bar {x},z,t_2,t_3) + \varLambda ^2 I^{(2)}(\bar {x},z,t_2,t_3), \end{equation}
and substituting (2.22) and (2.23) in (2.21), one has the following equation governing
$I^{(0)}$
:
\begin{align} {O}(1)\qquad&\frac {\partial I^{(0)}}{\partial t_2} =\frac {\partial }{\partial z}\left (D_{zz} \frac {\partial I^{(0)}}{\partial z} - V_z I^{(0)} \right )\!, \end{align}
at leading order. The absence of
$\bar {x}$
-derivatives in (2.24) implies a separable dependence on
$\bar {x}$
and
$z$
, so one may write
\begin{equation} I^{(0)} = K(\bar {x},t_3)F(z,t_2;{\it{{\it{Pe}}}}_r,\kappa ). \end{equation}
Here,
$F(z,t_2;{\it{{\it{Pe}}}}_r,\kappa )$
determines the dependence of the swimmer concentration on the transverse coordinate
$z$
, and evolves on a time scale of
$O(H^2/D_t)$
. The amplitude function
$K(\bar {x},t_3)$
governs the longitudinal spread of the swimmer population for times much greater than
$O(H^2/D_t)$
, in a reference frame moving with the mean speed
$\bar {u}_1$
. Again, without loss of generality, one may use the normalization condition
\begin{equation} \int ^{1}_{-1} {\rm d}z \, I(\bar {x},z,t_2)=K(\bar {x},t_{3}), \end{equation}
leading to
\begin{equation} \int _{-1}^1 F {\rm d}z = 1 \end{equation}
at leading order, and to
\begin{equation} \int _{-1}^1 I^{(1)} {\rm d}z = \int _{-1}^1 I^{(2)} {\rm d}z =0 \end{equation}
at
$O(\varLambda )$
and
$O(\varLambda ^2)$
, respectively.
For times much longer than
$O(H^2/D_t)$
,
$F$
in (2.25) attains a steady state,
$F_s(z;{\it{{\it{Pe}}}}_r,\kappa )$
, which satisfies
\begin{equation} \frac {\rm d}{{\rm d}z}\!\!\left (D_{zz} \frac {{\rm d}F_s}{{\rm d}z} - V_z F_s \right ) = 0. \end{equation}
Integrating once, and setting the integration constant to be zero to satisfy the no-flux boundary conditions at the channel walls, given by
\begin{eqnarray} \left (D_{zz} \frac {{\rm d}F_s}{{\rm d}z} - V_z F_s \right )\biggr \rvert _{z=\pm 1}= 0, \end{eqnarray}
we obtain
\begin{equation} F_s(z;{\it{{\it{Pe}}}}_r,\kappa ) = {\mathcal N} \exp \left [ \displaystyle \int _{-1}^z \frac {V_z(z')}{D_{zz}(z')}{\rm d}z' \right ]\!. \end{equation}
Use of (2.27) and (2.31) gives the normalization constant as
\begin{equation} \mathcal{N}=\frac {1}{\int _{-1}^{1} {\rm d}z \exp \left (\int _{-1}^{z}{\rm d}t\,\frac {V_{z}(t)}{D_{zz}(t)} \right )}. \end{equation}
Note that the one-dimensional transport scenario in (2.29) implies that the flux associated with (2.31) is, in fact, zero at every
$z$
-station, and not merely at the channel walls. Also note that the no-flux condition in (2.30) is a coarse-grained version of the integral boundary condition in position–orientation space, given in (2.4), and is expected to only be valid outside of wall layers. We comment further on this aspect in § 4.
The expression (2.31) serves to define the mean speed of the swimmer population, which is now given by
$\bar {u}_1 = \textstyle \int _{-1}^1 u_1F_s(z) {\rm d}z$
; note that
$\textstyle \int _{-1}^1 u_1^{\prime}(z)F_s(z) {\rm d}z = 0$
. For the classical case of a uniformly distributed passive solute,
$F_s(z) = {1}/{2}$
and
$\bar {u}_1= {2}/{3}$
. In contrast, the mean speed of the swimmer population is in general a function of
$\it{Pe}_r$
and
$\kappa$
. As will be seen in figure 2(b), for
$\kappa \lt 2$
, the swimmer concentration approaches a smooth limiting profile for
$Pe_r\,\rm {to}\,\infty$
, and
$\bar{u}_1$
tends to a limiting value less than unity, that is a function of
$\kappa$
alone. On the other hand, in figure 3(b) below,
$\lim _{\it{{\it{Pe}}}_r\rightarrow \infty }\bar {u}_1 = 1$
for active swimmers with
$\kappa \gtrsim 2$
, which exhibit a centreline collapse (Vennamneni et al. Reference Vennamneni, Nambiar and Subramanian2020b
), so
$F_s(z) \rightarrow \delta (z)$
for
${\it{{\it{Pe}}}}_r \rightarrow \infty$
.
(a) Half-channel swimmer concentration profiles (
$F_{s}(z;{\it{{\it{Pe}}}}_r,\kappa )$
) as a function of
${\it{{\it{Pe}}}}_r$
, for
$\kappa =1.6$
; horizontal dashed line (black) corresponds to the uniform profile at
${\it{{\it{Pe}}}}_r = 0$
. The inset shows the initial depletion of swimmer concentration near the centreline for small
${\it{{\it{Pe}}}}_r$
. (b) An enlarged view of the large-
${\it{{\it{Pe}}}}_r$
profiles near the channel centreline, illustrating convergence to the infinite-
${\it{{\it{Pe}}}}_r$
limiting form,
$F_s^\infty = 0.2|z|^{-0.6}$
, which appears as a dashed red curve in both (a) and (b). (c) The profile variance,
$\langle z^{2} \rangle =\int _{-1}^{1}{\rm d}z\ z^{2}F_{s}(z)$
, plotted as a function of
${\it{{\it{Pe}}}}_{r}$
, for
$\kappa =1.6$
; the accompanying insets show examples of swimmer concentration profiles in the two trapping regimes.
Half-channel swimmer concentration profiles (
$F_{s}(z;{\it{{\it{Pe}}}}_r,\kappa )$
) for different
${\it{{\it{Pe}}}}_r$
: (a)
$\kappa = 1.2$
and (b)
$\kappa = 3.6$
. Horizontal dashed lines (black) correspond to the uniform profile at
${\it{{\it{Pe}}}}_r = 0$
, with the infinite-
$P\!e_r$
profile,
$F_s^\infty (z)=0.4|z|^{-0.19}$
, shown as a dashed red curve in (a); black arrows in (b) denote direction of profile movement with increasing
${\it{{\it{Pe}}}}_r$
. (c) Here
$\langle z^2 \rangle$
versus
${\it{{\it{Pe}}}}_r$
for
$\kappa = 1.2$
and
$3.6$
. (d) Unscaled and (e) scaled concentration profiles, over the interval
${\it{{\it{Pe}}}}_r\in (100,150\,000)$
, for
$\kappa =3.6$
.
We now derive the equation governing the amplitude function
$K(\bar {x},t_3)$
by considering the equations that arise, at
$O(\varLambda )$
and
$O(\varLambda ^2)$
, from using (2.22) and (2.23) in (2.21):
\begin{align} {O}(\varLambda )\qquad\frac {\partial I^{(1)}}{\partial t_2} - \frac {\partial }{\partial z}\!\!\left (D_{zz} \frac {\partial I^{(1)}}{\partial z} - V_z I^{(1)} \right ) =&- {\it{{\it{Pe}}}}\, u^{\prime}_1 \frac {\partial I^{(0)}}{\partial \bar {x}}, \end{align}
\begin{align} {O}(\varLambda ^2)\qquad\frac {\partial I^{(2)}}{\partial t_2} - \frac {\partial }{\partial z}\!\!\left (D_{zz} \frac {\partial I^{(2)}}{\partial z} - V_z I^{(2)} \right ) =& - {\it{{\it{Pe}}}}\, u^{\prime}_1 \frac {\partial I^{(1)}}{\partial \bar {x}} - \frac {\partial I^{(0)}}{\partial t_3}. \end{align}
Substituting the separable form of
$I^{(0)}$
from (2.25) above, in (2.33), one obtains, at
$O(\varLambda )$
,
\begin{equation} \frac {\partial I^{(1)}}{\partial t_2} - \frac {\partial }{\partial z}\!\!\left (D_{zz} \frac {\partial I^{(1)}}{\partial z} - V_z I^{(1)} \right ) = - {\it{{\it{Pe}}}}\,u^{\prime}_1F_s\frac {\partial K}{\partial \bar {x}}. \end{equation}
For large
$t_2$
, the time derivative in (2.35) may be neglected, and (2.35) is then readily integrated once to yield the following inhomogeneous first-order ordinary differential equation for
$I^{(1)}$
:
\begin{equation} \frac {{\rm d}I^{(1)}}{{\rm d}z} - \frac {V_z}{D_{zz}} I^{(1)} = Pe \frac {\partial K}{\partial \bar {x}} \frac {1}{D_{zz}(z)}\displaystyle \int _{-1}^z u^{\prime}_1(z')F_s(z')\,{\rm d}z', \end{equation}
on using the no-flux condition at the lower channel boundary
\begin{equation} \frac {{\rm d}I^{(1)}}{{\rm d}z} - \frac {V_z}{D_{zz}} I^{(1)} = 0\,\,\, \text{at}\,\,\, z = -1. \end{equation}
One may solve (2.36) using an integrating factor (
$\mathcal{N}/F_s(z)$
), whence one obtains
\begin{equation} I^{(1)}(\bar {x},z) = Pe \frac {\partial K}{\partial \bar {x}} F_s(z) \displaystyle \int _{-1}^z {\rm d}z' \frac {1}{F_s(z')D_{zz}(z')}\left (\displaystyle \int _{-1}^{z'} {\rm d}z^{\prime\prime} u^{\prime}_1(z^{\prime\prime})F_s(z^{\prime\prime}) \right ) + J(\bar {x})\frac {F_s(z)}{\mathcal{N}}. \end{equation}
The above expression already satisfies the no-flux condition at the upper boundary (
$z=1$
) regardless of
$J(\bar {x})$
. The latter function may be determined from the normalization constraint at
$O(\varLambda ^2)$
in (2.28), although the particular expression for
$J(\bar {x})$
does not affect the expression for the effective longitudinal diffusivity,
$\bar {D}_{\it{eff}}$
, obtained below.
Substituting (2.38) in (2.34), integrating with respect to
$z$
over
$[-1,1]$
, and using
$\int _{-1}^{1}{\rm d}z F_{s}(z)u^{\prime }_{1}(z)=0$
, one obtains the following equation governing
$K$
:
\begin{equation} \frac {\partial K}{\partial t_3} = \bar {D}_{\it{eff}}\frac {\partial ^2 K}{\partial \bar {x}^2}, \end{equation}
in scaled variables, with
\begin{equation} \bar {D}_{\it{eff}} = - Pe^2 \displaystyle \int _{-1}^1 {\rm d}z\, u^{\prime}_1(z) F_s(z) \displaystyle \int _{-1}^z {\rm d}z' \frac {1}{F_s(z')D_{zz}(z')}\left (\displaystyle \int _{-1}^{z'} {\rm d}z^{\prime\prime} u^{\prime}_1(z^{\prime\prime})F_s(z^{\prime\prime}) \right )\!. \end{equation}
One may simplify (2.40) using an integration by parts, while also using the symmetry of the integrand about the centreline (
$z=0$
), leading to the following simpler expression for the diffusivity:
\begin{equation} \bar {D}_{\it{eff}} = 2Pe^2 \displaystyle \int _{0}^1 {\rm d}z \frac {1}{F_s(z)D_{zz}(z)}\left (\displaystyle \int _{0}^{z} {\rm d}z' u^{\prime}_1(z')F_s(z') \right )^2 , \end{equation}
which also makes explicit its positive definite character. One may rewrite (2.41) further by dividing and multiplying by
$F_s$
, leading to
\begin{equation} \bar {D}_{\it{eff}} = 2Pe^2 \displaystyle \int _{0}^1 {\rm d}z \frac {F_s(z)}{D_{zz}(z)}\left (\frac {1}{F_s(z)}\displaystyle \int _{0}^{z} {\rm d}z' u^{\prime}_1(z')F_s(z') \right )^2 , \end{equation}
which is broadly consistent with an interpretation of the longitudinal diffusivity of the swimmer population arising as the product of the variance of the streamwise velocity fluctuation weighted by the swimmer concentration, and a correlation time for these fluctuations that scales inversely with the diffusivity along the gradient direction.
On reverting to dimensional variables via
$t_3 = \varLambda ^2 \epsilon ^2 D_r t$
and
$x = L\bar {x}$
, and going back to the original laboratory reference frame, the longitudinal length scale
$L$
(introduced solely for scaling purposes) cancels out, and (2.39) reduces to the convection–diffusion equation
\begin{equation} \frac {\partial K}{\partial t} + \bar {u}_1\frac {\partial K}{\partial x} = \bar {D}_{\it{eff}}\frac {\partial ^2 K}{\partial {x^2}}, \end{equation}
with the longitudinal diffusivity, in dimensional form, being given by
\begin{equation} {D}_{\it{eff}} = \frac {2U_s^2}{D_r} Pe^2 \displaystyle \int _{0}^1 {\rm d}z \frac {F_s(z)}{D_{zz}(z)}\left (\frac {1}{F_s(z)}\displaystyle \int _{0}^{z} {\rm d}z' u^{\prime}_1(z')F_s(z') \right )^2. \end{equation}
The assumptions made during the course of the above derivation imply that (2.44) is valid for
$\epsilon \ll 1$
,
$\varLambda Pe \ll 1$
or, equivalently,
$\epsilon \ll 1, {\it{{\it{Pe}}}}_r \gg 1$
. The equivalence of the small-
$\varLambda Pe$
and large-
${\it{{\it{Pe}}}}_r$
assumptions is shown after having clarified the different scaling regimes for
$D_{\it{eff}}$
in § 3.2.
2.2.3. Longitudinal diffusivity (
$D_{\it{eff}}$
) for arbitrary Péclet numbers
For small
${\it{{\it{Pe}}}}_r$
, the dominant rotary diffusion renders the swimmer orientation distribution isotropic. The isotropy of the orientation distribution implies a vanishing drift (
$V_z = 0$
), and therefore, a uniform swimmer concentration profile (
$F_s = {1}/{2}$
). Thus, in this limit, one expects an exact equivalence between an active swimmer and a solute molecule with diffusivity
$D_t = U_s^2/6D_r$
; in dimensionless terms, the latter corresponds to
$D_{zz} = {1}/{6}$
in (2.44). Using
$F_s = {1}/{2}$
,
$D_{zz} = {1}/{6}$
in (2.44) gives
$D_{\it{eff}} = D_t({8Pe^2}/{945})$
, which is the active analogue of the
$O(Pe^2)$
shear-enhanced diffusion in ambient plane Poiseuille flow, quoted in § 1 (Wooding Reference Wooding1960), the only difference being that the swimmer diffusivity now replaces the solute diffusivity. Thus, although (2.44) is supposed to be valid only for large
${\it{{\it{Pe}}}}_r$
, it gives the correct result even for
${\it{{\it{Pe}}}}_r \ll 1$
provided
${\it{Pe}} \gg 1$
. The restriction of (2.44) to
${\it{Pe}}\gg 1$
is largely an apparent one, since for
${\it{Pe}} \sim O(1)$
,
${\it{{\it{Pe}}}}_r \sim O(\epsilon ^2) \ll 1$
. Rotary diffusion is rapid enough for the swimmers to diffuse in a nearly isotropic fashion and the aforementioned analogy between an active swimmer and a solute with translational diffusivity
$D_t$
continues to hold. Hence, in a manner similar to the passive scenario (Aris Reference Aris1956), one expects to circumvent the restriction to large
${\it{Pe}}$
in (2.44) by including the intrinsic translational diffusivity of the swimmer as an independent additive contribution, leading to
\begin{align} \frac {D_{\it{eff}}}{D_{t}} &= 1 + 12 Pe^2 \displaystyle \int _{0}^1 {\rm d}z \frac {F_s(z)}{D_{zz}(z)}\left (\frac {1}{F_s(z)}\displaystyle \int _{0}^{z} {\rm d}z' u^{\prime}_1(z')F_s(z') \right )^2 \nonumber \\ &= 1 + \frac {12 Pe^{2}_{r}}{\epsilon ^{4}} \displaystyle \int _{0}^1 {\rm d}z \frac {F_s(z)}{D_{zz}(z)}\left (\frac {1}{F_s(z)}\displaystyle \int _{0}^{z} {\rm d}z' u^{\prime}_1(z')F_s(z') \right )^2, \\[9pt] \nonumber \end{align}
as the final expression for
$D_{\it{eff}}$
. This gives the correct results for
$D_{\it{eff}}$
in both the limits
${\it{{\it{Pe}}}}_r \ll 1, Pe \gg 1$
(
$({D_{\it{eff}}}/{D_t}) = ({8Pe^2}/{945})$
, as mentioned above) and
${\it{Pe}}, {\it{{\it{Pe}}}}_r \ll 1$
(
$({D_{\it{eff}}}/{D_t}) = 1$
), and it is therefore reasonable to expect (2.45) to provide a good approximation for arbitrary
${\it{Pe}}$
and
${\it{{\it{Pe}}}}_r$
, provided only that
$\epsilon \ll 1$
. Using Langevin simulations, we show this to be the case for spherical swimmers in Appendix B for
$\epsilon$
up until at least
$0.2$
; see table 1 in Appendix B. A limited set of results in table 2, in Appendix B, show that (2.45) remains accurate even for non-spherical swimmers for
${\it{{\it{Pe}}}}_r \sim O(1)$
. The expression (2.45) also allows one to justify the assumption, of
$\varLambda /\epsilon ^2$
being order unity, made in § 2.2.1 when separating the orientation and spatial degrees of freedom. Since the integral in (2.45) is only a function of
${\it{{\it{Pe}}}}_r$
and
$\kappa$
, for these parameters being fixed and with
$\epsilon \ll 1$
, one has
$D_{zz} \sim D_t$
,
$D_{\it{eff}} \sim D_t\epsilon ^{-4}$
. Using these in the relation
$H^2/D_{zz} \sim L_{\it{min}}^2/D_{\it{eff}}$
gives
$\varLambda ^2 \sim H^2/L_{\it{min}}^2 \sim \epsilon ^4$
, and therefore
$\varLambda /\epsilon ^2 \sim O(1)$
.
The results for the effective longitudinal diffusivity, presented in the following section, are based on a numerical evaluation of the integral in (2.45) for swimmers of various aspect ratios (
$\kappa$
), and as a function of
${\it{{\it{Pe}}}}_r$
. This involves the numerical evaluation of the drift and diffusivity expressions that enter the expression for
$F_s$
in (2.31), using truncated expansions in spherical harmonics for the orientation probability density
$G$
. Analogous expansions are used for the other probability densities involved (
$\varOmega _{11}$
and
$\varOmega _{12}$
). Details of these calculations are given in Appendix A.
3. Results and discussion
3.1. Migration along the gradient direction
We first discuss swimmer migration in the
$z$
-direction, based on the expression for the transverse swimmer concentration profile
$F_s$
, given by (2.31), with the drift (
$V_z$
) and diffusivity (
$D_{zz}$
) entering this expression being given by (2.19) and (2.20), respectively. This is because transverse shear-induced migration directly affects dispersion along the flow direction. Note that the above expressions were originally obtained in Vennamneni et al. (Reference Vennamneni, Nambiar and Subramanian2020b
) using a version of the multiple scales analysis detailed in § 2, but without the amplitude function
$K$
and the associated slower time variable, this being equivalent to assuming the swimmer population to be infinite in extent along the flow direction. The expression (2.31) was then used to understand the high-shear and LST regimes, and specifically, the transition from the former to the latter for any finite
$\kappa$
with increasing
${\it{{\it{Pe}}}}_r$
. It was found that
$D_{zz}$
and
$V_z$
transition from a small-
${\it{{\it{Pe}}}}_r$
plateau to a
${\it{{\it{Pe}}}}_r^{-2}$
-scaling regime, with increasing
${\it{{\it{Pe}}}}_r$
. This change in scaling is accompanied by a reversal in drift (a change in sign of
$V_z$
), from initially pointing wallward to eventually pointing towards the centreline, which leads to the aforementioned change in the trapping regime across a threshold
${\it{{\it{Pe}}}}_r$
. For large
$\kappa$
,
$D_{zz}$
and
$V_z$
exhibit an intermediate asymptotic scaling of
$O({\it{{\it{Pe}}}}_r^{-({4}/{3})})$
in the interval
$1 \ll {\it{{\it{Pe}}}}_r \ll \kappa ^3$
, prior to the eventual
${\it{{\it{Pe}}}}_r^{-2}$
-scaling regime, and the transition
${\it{{\it{Pe}}}}_r$
is therefore of
$O(\kappa ^3)$
in this limit; as shown in figure 1(a), the large-
$\kappa$
asymptote (depicted by a dashed line) of the boundary separating the high-shear and LST regimes is given by
${\it{{\it{Pe}}}}_r \approx 0.07 \kappa ^3$
.
Since the drift and diffusivity are both
$O({\it{{\it{Pe}}}}_r^{-2})$
for sufficiently large
${\it{{\it{Pe}}}}_r$
, the drift-to-diffusivity ratio asymptotes to a
${\it{{\it{Pe}}}}_r$
-independent form away from the centreline, given by
$V_{z}/D_{zz}=-C(\kappa )/\vert z \vert$
. For
$\kappa$
below the centreline-collapse threshold (
$\kappa \approx 2$
), this leads to the limiting concentration profile,
$F_s^\infty (z;\kappa ) = ({1}/{2})(1-C)|z|^{-C}$
with
$C(\kappa ) \lt 1$
being an aspect-ratio-dependent constant. Figure 2(a) illustrates the approach of the transverse swimmer concentration profiles, at a fixed
$x$
and for
$t_2 \gg 1$
, towards the aforesaid limiting form (shown as a dashed red curve), with increasing
${\it{{\it{Pe}}}}_r$
. The profiles correspond to
$\kappa = 1.6$
(
$C \approx 0.6$
for this case), and the accompanying inset highlights the dip near the centreline (
$F_s(0) \lt 0.5$
) for the small-
${\it{{\it{Pe}}}}_r$
profiles, a signature of HST, that morphs into a peak (
$F_s(0) \gt 0.5$
) with increasing
${\it{{\it{Pe}}}}_r$
. Figure 2(b) shows a zoomed-in version of the large-
${\it{{\it{Pe}}}}_r$
profiles near the channel centreline to highlight the continued accumulation of swimmers in this region with increasing
${\it{{\it{Pe}}}}_r$
, which leads to an continued increase in the amplitude of the said peak;
$F_s(0)$
diverges in the limit
${\it{{\it{Pe}}}}_r \rightarrow \infty$
. Figure 2(c) shows the corresponding plot for the profile variance,
$\langle z^2\rangle =\textstyle \int _{-1}^1 z^2 F_s(z;{\it{{\it{Pe}}}}_r,\kappa ) {\rm d}z$
, a measure of the spread of the profile in the gradient direction. The variance is seen to exhibit a non-monotonic transition from a zero-
${\it{{\it{Pe}}}}_r$
plateau (corresponding to a uniform profile, and equalling
$1/3$
), to a lower infinite-
${\it{{\it{Pe}}}}_r$
one (corresponding to the limiting profile above, and equalling
$(1-C)/(3-C)=0.166$
for
$C = 0.6$
). The non-monotonicity arises from the initial increase in the variance in the HST regime, due to swimmers migrating towards the walls, being followed by the eventual decrease driven by swimmers drifting back towards the centreline. It must be mentioned that, within the framework used here, spherical swimmers remain uniformly distributed (
$C = 0$
for
$\kappa = 1$
) regardless of
${\it{{\it{Pe}}}}_r$
, with the latter parameter only determining the rate of approach to the uniform state. The uniform state results due to the boundaries only serving to enforce the no-flux constraint. Thus, specific swimmer–boundary interactions, hydrodynamic or otherwise, are neglected. Such interactions have been shown to lead to accumulation within wall layers in narrow channels with
$\epsilon \sim O(1)$
, see Elgeti & Gompper (Reference Elgeti and Gompper2009, Reference Elgeti and Gompper2013) and Chilukuri, Collins & Underhill (Reference Chilukuri, Collins and Underhill2015). We further comment on the nature of such interactions in § 4.
Figures 3(a) and 3(b) plot the swimmer concentration profiles corresponding to a fixed
$x$
and for
$t_2 \gg 1$
, as a function of
${\it{{\it{Pe}}}}_r$
, for
$\kappa = 1.2$
and
$3.6$
. The aim here is to contrast the nature of the transition from high-shear to LST, for
$\kappa$
below (
$\kappa = 1.2$
) and above (
$\kappa = 3.6$
) the centreline-collapse threshold (
$\kappa \approx 2$
). Figure 3(c) contrasts the behaviour of the corresponding profile variances. The behaviour for the smaller aspect ratio,
$\kappa = 1.2$
, is similar to that already seen in figure 2. Thus, in figure 3(a), the concentration profiles approach a
${\it{{\it{Pe}}}}_r$
-independent form,
$F_s^\infty (z) = 0.4|z|^{-0.19}$
, for
${\it{{\it{Pe}}}}_r \rightarrow \infty$
; this limiting profile corresponds to
$C = 0.19$
, and is again shown as a dashed red curve. Similarly, the profile variance,
$\langle z^2\rangle \!$
in figure 3(c) (the blue curve), after an initial increase, asymptotes to a lower infinite-
${\it{{\it{Pe}}}}_r$
plateau. Despite this similarity, a comparison of figures 2(c) and 3(c) shows that the profile variance for
$\kappa = 1.6$
exhibits a greater maximum and a smaller infinite-
${\it{{\it{Pe}}}}_r$
plateau, than for
$\kappa = 1.2$
; further, the variance is seen to dip below
$1/3$
at a larger
${\it{{\it{Pe}}}}_r$
, indicative of a delayed transition from HST to LST. These are, in fact, already indicative of the singular nature of the
$\kappa = \infty$
limit examined below, when the transition in trapping regime is postponed to
${\it{{\it{Pe}}}}_r = \infty$
, so the swimmer concentration profiles only exhibit HST.
For
$\kappa \gtrsim 2$
,
$C(\kappa )\gt 1$
and
$F_s^\infty \propto |z|^{-C}$
becomes non-normalizable due to the non-integrable divergence at
$z = 0$
. The implication is that swimmers with
$\kappa \gtrsim 2$
collapse onto the centreline, with
$F_s(z;{\it{{\it{Pe}}}}_r,\kappa ) \rightarrow \delta (z)$
for
${\it{{\it{Pe}}}}_r \rightarrow \infty$
. This is seen to be the case for
$\kappa =3.6$
in figure 3(b) where, for
${\it{{\it{Pe}}}}_r\rightarrow \infty$
, the swimmer concentration approaches zero away from the centreline, while still diverging at
$z = 0$
(figure 3
b, inset);
$C= 2.41$
for this aspect ratio. Mirroring this trend,
$\langle z^2 \rangle$
in figure 3(c) (the red curve), following an initial increase, approaches zero as
${\it{{\it{Pe}}}}_r^{1-C} \equiv {\it{{\it{Pe}}}}_r^{-1.41}$
for
${\it{{\it{Pe}}}}_r \rightarrow \infty$
. Figures 3(d) and 3(e) further illustrate the profile localization in the neighbourhood of the centreline, and the approach towards the singular limiting form (
$\delta (z)$
). The scaled concentration profiles in figure 3(e) show that
$F_s\!\propto \!{\it{{\it{Pe}}}}_r^{1-C}|z|^{-C}$
outside a shrinking
$O({\it{{\it{Pe}}}}_r^{-1})$
central core, where
$F_s$
saturates at an
$O({\it{{\it{Pe}}}}_r)$
maximum.
It is worth remarking on the non-monotonicity of the concentration profiles in figures 3(d) and 3(e). As is seen from figures 1 and 2, this non-monotonicity is, in fact, a generic feature of profiles in the low-shear-trapping regime. For the plane Poiseuille profile being examined, the local rotary Péclet number, defined based on the local shear rate (rather than using
$U_m$
and
$H$
), is
$2z{\it{{\it{Pe}}}}_r$
. This Péclet number decreases from a (absolute) maximum of
$2{\it{{\it{Pe}}}}_r$
at either wall, to zero at the centreline. Thus, for any
${\it{{\it{Pe}}}}_r$
, however large, there is always a neighbourhood of the centreline where the local Péclet number is small enough that swimmers in this region exhibit a wallward migration. Therefore, for an order unity aspect ratio, swimmers in the region
$O({\it{{\it{Pe}}}}_r^{-1}) \lt |z| \lt 1$
migrate towards the centreline, while those in the region
$|z| \lt O({\it{{\it{Pe}}}}_r^{-1})$
migrate away from it (that is, wallward). This change in migration direction causes the swimmers to converge onto, and the swimmer concentration to attain a maximum at, a location an
$O({\it{{\it{Pe}}}}_r^{-1})$
distance away from the centreline. We end this subsection by noting that, right at the centreline-collapse threshold,
$F_z(z) \propto 1/z$
for
$z \gt O({\it{{\it{Pe}}}}_r^{-1})$
, which leads to
$\langle z^2 \rangle \sim O(\ln {\it{{\it{Pe}}}}_r)^{-1}$
for
${\it{{\it{Pe}}}}_r \rightarrow \infty$
.
3.2. Dispersion along the flow direction
We now examine how the migration behaviour described above, and that leads to centreline localization at large
${\it{{\it{Pe}}}}_r$
, affects dispersion. Using the functional form for
$F^\infty _{s}$
given in § 3.1, and the large-
${\it{{\it{Pe}}}}_r$
asymptote,
$D_{zz} \approx C_1(\kappa )/(-2z{\it{{\it{Pe}}}}_r)^2$
,
$-2z{\it{{\it{Pe}}}}_r$
being the local
$z$
-dependent Péclet number defined above (Vennamneni et al. Reference Vennamneni, Nambiar and Subramanian2020b
), in the second integral term in (2.45), one obtains the limiting form of the longitudinal diffusivity as
\begin{equation} \lim _{\it{{\it{Pe}}}_r \rightarrow \infty } \frac {D_{\it{eff}}}{D_t}= \epsilon ^{-4}{\it{{\it{Pe}}}}_r^4\frac {192(1-C)}{C_1(C-3)^2\left(315-143C+21C^2-C^3\right)}, \end{equation}
for
$\kappa \lesssim 2$
. Figure 4 plots
$D_{\it{eff}}/D_t$
obtained from numerically evaluating (2.45), as a function of
${\it{{\it{Pe}}}}_r$
, for
$\kappa = 1.2$
and
$3.6$
; the large-
${\it{{\it{Pe}}}}_r$
asymptote, given by (3.1), appears as a dashed red line for
$\kappa = 1.2$
. For this aspect ratio,
$D_{\it{eff}}/D_t$
transitions from the initial zero-
${\it{{\it{Pe}}}}_r$
plateau to a passive
${\it{{\it{Pe}}}}_r^2$
-regime that extends over the interval
$\epsilon ^2 \ll {\it{{\it{Pe}}}}_r \ll 1$
. There is then a second transition to a
${\it{{\it{Pe}}}}_r^4$
-scaling regime across
${\it{{\it{Pe}}}}_{r} \approx 3$
. The
${\it{{\it{Pe}}}}_r^4$
-scaling for large
${\it{{\it{Pe}}}}_r$
is readily explained using the classical expression for
$D_{\it{eff}}$
(see § 1), with the only change being that one uses
$D = D_{zz} \sim {\it{{\it{Pe}}}}_r^{-2}$
for
${\it{{\it{Pe}}}}_r \gg 1$
, instead of the originally constant solute diffusivity. Thus, the increased dispersion arises from the reduction in
$D$
at large
${\it{{\it{Pe}}}}_r$
. This reduction, which owes its origin to shear-induced rotation, leads to an increase in the time required for the swimmer to sample a new transverse location, in turn increasing the decorrelation time between convective displacements along the flow direction. While the
${\it{{\it{Pe}}}}_r^4$
-scaling, arising from vorticity-induced rotation, has been identified earlier as providing an upper bound for microswimmer dispersion (Dehkharghani et al. Reference Dehkharghani, Waisbord, Dunkel and Guasto2019; Peng & Brady Reference Peng and Brady2020), we show that this scaling applies exactly to spheroidal swimmers with
$1 \leqslant \kappa \lesssim 2$
. The prefactor in the
${\it{{\it{Pe}}}}_r^4$
-asymptote, obtained using (3.1) while also accounting for the aspect-ratio dependence of both
$C$
and
$C_1$
, turns out to be a non-monotonic function of
$\kappa$
, with a maximum at
$\kappa \approx 1.4$
(see Appendix A and figure 8
a). Swimmers with
$\kappa \approx 1.4$
are thus the fastest dispersing prolate swimmers.
Here
$D_{\it{eff}}/D_t$
versus
${\it{{\it{Pe}}}}_r$
for
$\kappa =1.2$
and
$3.6$
;
$\epsilon =0.1$
. The dashed black curve is the passive asymptote obtained from (2.45) using
$D_{zz}=1/6, F_s=1/2$
; the dashed red line is the large-
${\it{{\it{Pe}}}}_r$
-asymptote for
$\kappa = 1.2$
, obtained from (3.1) using
$C = 0.19$
and
$C_1 = 2.12$
; the black dash–dotted line is the
${\it{{\it{Pe}}}}_r^{5-C}$
-anomalous asymptote (
$C=2.41$
) for
$\kappa = 3.6$
. The vertical blue band marks the onset of transition from classical to anomalous scaling, and includes Péclet numbers realized in recent experiments (see § 4).
Centreline collapse leads to a novel
$D_{\it{eff}}$
-scaling regime for
$\kappa \gtrsim 2$
. This regime may, in fact, be anticipated from (3.1) – the coefficient of
${\it{{\it{Pe}}}}_r^4$
in this expression approaches zero for
$\kappa \rightarrow 2^-$
since
$C\rightarrow 1^-$
, suggesting that, for
$\kappa \gtrsim 2$
,
$D_{\it{eff}}$
must eventually increase as a smaller power of
${\it{{\it{Pe}}}}_r$
. The suggestion is confirmed by the plot for
$\kappa = 3.6$
in figure 4. For this aspect ratio, following the passive
${\it{{\it{Pe}}}}_r^2$
-regime, the subsequent approach to the classical
${\it{{\it{Pe}}}}_r^4$
-regime described above, is interrupted by
$D_{\it{eff}}/D_t$
transitioning to an anomalous scaling regime characterized by a lower exponent (
$2.59$
); furthermore, this exponent decreases with increasing
$\kappa$
(see figure 5, and figure 8
b in Appendix A).
(a) Here
$D_{\it{eff}}/D_t$
versus
${\it{{\it{Pe}}}}_r$
for
$\kappa \in [6,12]$
, and for
$\kappa = \infty$
;
$\epsilon =0.1$
. The passive
${\it{Pe}}^{2}_{r}$
-asymptote and the classical
${\it{{\it{Pe}}}}_r^{10/3}$
-asymptote appear as dashed black and blue lines, respectively; as in figure 4, vertical lines at the
${\it{{\it{Pe}}}}_r$
values indicated mark the transitions between different scaling regimes. Panel (b) shows a magnified view of the transition from the
${\it{{\it{Pe}}}}_r^{10/3}$
-scaling regime, to the anomalous
${\it{{\it{Pe}}}}_r^{5-C}$
-scaling regime; here,
$C = 3.75, 4.49, 5.04, 5.48$
for
$\kappa = 6, 8, 10$
and 12, respectively. Anomalous asymptotes appear as dash–dotted black lines, and black stars corresponding to the intersection of these asymptotes with the
${\it{{\it{Pe}}}}_r^{10/3}$
-asymptote, mark the transition from classical to anomalous scaling; the inset view highlights the rapid increase in
$C$
as a function of
$\kappa$
. (c) The dispersion exponent (
$\alpha$
), in the (truncated) intermediate asymptotic regime, plotted as a function of
$\kappa$
; the (highlighted) exponent for
$\kappa = 4$
pertains to the simulations of Ganesh et al. (Reference Ganesh, Douarche, Dentz and Auradou2023), and is obtained from the fit shown in the right-hand side plot of figure 6(b) below.
The anomalous scaling regime above does not conform to the classical scenario. Thus, scaling arguments along the lines of those in § 1, but based on a channel width of
$O(H{\it{{\it{Pe}}}}_r^{-1})$
, representative of the central core of collapsed swimmers (see figure 3
e), do not work. Instead, dispersion is controlled by the small fraction of swimmers outside this core, that sample the largest velocity variations. This can be seen by estimating the contributions to the integral in (2.45) from the core (
$|z| \sim O({\it{{\it{Pe}}}}_r^{-1})$
) and tail (
$z \sim O(1)$
) regions. The majority of the swimmers are in the core, and
$\bar {u}_{1}=1$
at leading order. Accounting for the
$O({\it{{\it{Pe}}}}_r^{1-C})$
fraction of swimmers outside, with
$O(1)$
velocity deficits, leads to the more accurate estimate
$\bar {u}_1\!\sim \!1-O({\it{{\it{Pe}}}}_r^{1-C})$
for
$1\lt C \lt 3$
. In the core region, one therefore has
$z' \sim O({\it{{\it{Pe}}}}_r^{-1}), F_s \sim O({\it{{\it{Pe}}}}_r), u^{\prime}_1 \sim O({\it{{\it{Pe}}}}_r^{1-C})$
, leading to the inner integral
$\textstyle \int _0^{O({\it{{\it{Pe}}}}_r^{-1})}\!u^{\prime}_1F_s {\rm d}z'$
in (2.45) being
$O({\it{{\it{Pe}}}}_r^{1-C})$
; using
$D_{zz} \sim O(1)$
owing to the vanishingly small shear, the core contribution,
$\textstyle \int _0^{O({\it{{\it{Pe}}}}_r^{-1})}({{\rm d}z}/{F_sD_{zz}}).\,O\,[{\it{{\it{Pe}}}}_r^{2(1-C)}]$
, comes out to be
$O({\it{{\it{Pe}}}}_r^{-2C})$
for
$1 \lt C \lt 3$
. For
$C\gt 3$
, the fraction of swimmers outside the core is small enough for
$\bar {u}_1$
to differ from the centreline maximum only due to the
$O({\it{{\it{Pe}}}}_r^{-2})$
parabolic variation of the Poiseuille profile within the core. This leads to
$u_1^{\prime} \sim O({\it{{\it{Pe}}}}_r^{-2})$
, and the resulting core-region contribution to the integral is
$ O({\it{{\it{Pe}}}}_r^{-6})$
. In the tail region,
$z' \sim O(1), F_s \sim O({\it{{\it{Pe}}}}_r^{1-C}), u^{\prime}_1 \sim O(1)$
. The inner integral is again
$O({\it{{\it{Pe}}}}_r^{1-C})$
, since the net swimmer flux is zero in the chosen reference frame, implying equal and opposite core and tail region contributions. Using
$D_{zz} \sim O({\it{{\it{Pe}}}}_r^{-2}), z \sim O(1)$
, the tail-region contribution comes out to be
$O({\it{{\it{Pe}}}}_r^{3-C})$
, with this estimate unaffected by the transition in the core-region scaling across
$C = 3$
. For
$C \gt 1$
, the tail region contribution is dominant. Accounting for the prefactor
${\it{{\it{Pe}}}}_r^2/\epsilon ^4$
in front of the integral in (2.45),
$D_{\it{eff}}/D_t \propto {\it{{\it{Pe}}}}_r^{5-C(\kappa )}\!\epsilon ^{-4}$
for
$\kappa \!\gtrsim \!2$
. With
$C = 2.41$
for
$\kappa = 3.6$
, this estimate (shown by the black dash–dotted line) agrees well with the numerics in figure 4. For the marginal case of
$\kappa \approx 2$
, as already mentioned in § 3.1, the profile variance exhibits a logarithmic decrease with
${\it{{\it{Pe}}}}_r$
, and this leads to
$D_{\it{eff}}/D_t \propto {\it{{\it{Pe}}}}_r^4\epsilon ^{-4}(\ln {\it{{\it{Pe}}}}_r)^{-1}$
for
${\it{{\it{Pe}}}}_r \gg 1$
; we have verified the onset of this logarithmic scaling, starting at
${\it{{\it{Pe}}}}_r \approx 100$
, for
$\kappa = 1.99$
.
Figure 5(a) plots
$D_{\it{eff}}/D_t$
versus
${\it{{\it{Pe}}}}_r$
for larger
$\kappa$
upto 12, and helps highlight the emergence of a new intermediate asymptotic regime with
$D_{\it{eff}}/D_t \propto {\it{{\it{Pe}}}}_r^{{10}/{3}}$
(the blue dashed line), instead of the
${\it{{\it{Pe}}}}_r^4$
-scaling regime for
$\kappa = 3.6$
seen in figure 4. This regime pertains to sufficiently slender swimmers that exhibit HST for
$1\ll {\it{{\it{Pe}}}}_r \ll \kappa ^3$
with
$D_{zz}, V_z \sim O({\it{{\it{Pe}}}}_r^{-{4}/{3}})$
in this
${\it{{\it{Pe}}}}_r$
-interval (Vennamneni et al. Reference Vennamneni, Nambiar and Subramanian2020b
) – this
$D_{zz}$
-scaling, when used in the scaling argument of § 1, yields the exponent (
$10/3$
) above, attesting to the classical nature of this regime. Figure 5(a) also shows the transition to the anomalous scaling regime being postponed to progressively larger
${\it{{\it{Pe}}}}_r$
, with increasing
$\kappa$
. The magnified view in figure 5(b) better highlights this aspect. The black stars in this figure mark the intersection of the classical (
${\it{{\it{Pe}}}}_r^{{10}/{3}}$
) and anomalous (
${\it{{\it{Pe}}}}_r^{5-C}$
) asymptotes, and show the transition
${\it{{\it{Pe}}}}_r$
increasing from
$748$
for
$\kappa =6$
, to
$3142$
for
$\kappa = 12$
. Since the onset of anomalous scaling correlates to LST and the associated centreline collapse, the transition
${\it{{\it{Pe}}}}_r$
must eventually increase as
$\kappa ^3$
, as already shown in figure 1(a); although, one needs to go to
$\kappa$
larger than
$12$
to observe this scaling. Figure 5(c) plots the dispersion exponent in the intermediate asymptotic regime, that precedes the anomalous one, as a function of
$\kappa$
. As already pointed out in the context of figure 4, the theoretically predicted exponent is not realized owing to onset of anomalous scaling. This is again seen in figure 5(c), where the dispersion exponent is seen to equal
$4$
only for
$\kappa \lt 2$
, in the absence of the anomalous regime. For larger
$\kappa$
, the exponent decreases from
$4$
, slowly approaching
$10/3$
for
$\kappa \gg 1$
. The slow approach is because the appearance of the
${\it{{\it{Pe}}}}_r^{-{4}/{3}}$
-scaling regime for
$D_{zz}$
, over a substantial
${\it{{\it{Pe}}}}_r$
-interval, requires going to
$\kappa$
substantially greater than
$10$
. It is worth mentioning that the exponents between
$4$
and
$10/3$
in figure 5(c), representative of a truncated intermediate aysmptotic regime, must be treated as empirical fits, obtained using an adjustable
${\it{{\it{Pe}}}}_r$
-interval. They have no fundamental basis, and the hope is that they nevertheless serve as a guide for future experiments or computations that might attempt to fit dispersion results over a finite
${\it{{\it{Pe}}}}_r$
-interval. On the whole, figure 5(a–c) show that
$\kappa = \infty$
is a singular limit when the anomalous regime recedes to infinity, and
$D_{\it{eff}}/D_t$
transitions from an initial plateau to the passive
${\it{Pe}}^{2}_{r}$
-scaling regime, and finally to the
${\it{{\it{Pe}}}}_r^{{10}/{3}}$
-scaling regime. Thus, the infinitely slender approximation, often used in the literature (Bearon & Hazel Reference Bearon and Hazel2015; Ezhilan & Saintillan Reference Ezhilan and Saintillan2015; Borgnino et al. Reference Borgnino, Boffetta, Cencini, De Lillo and Gustavsson2022), fails to capture anomalous dispersion for large but finite
$\kappa$
.
There are a couple of additional noteworthy points in figure 5. First, the increase in
$C(\kappa )$
, shown in the inset of figure 5(c), causes the anomalous exponent (
$5-C$
) to dip below
$2$
for
$\kappa \gtrsim 4.5$
, implying that the large-
${\it{{\it{Pe}}}}_r$
dispersion for such swimmers is slower than passive particles with the same
$\kappa$
. A second more remarkable feature is that
$C \approx 5$
for
$\kappa \approx 10$
, increasing further and appearing to eventually scale as
$\ln \kappa$
(not shown). This implies a near-zero anomalous exponent for
$\kappa \approx 10$
, and a negative one for larger
$\kappa$
, implying an eventual decrease in
$D_{\it{eff}}/D_t$
with increasing
${\it{{\it{Pe}}}}_r$
for
$\kappa \gt 10$
. Thus, for
$\kappa \gt 10$
,
$D_{\it{eff}}/D_t$
starts off on the low-
${\it{{\it{Pe}}}}_r$
plateau, increases as
${\it{Pe}}^{2}_{r}$
for
$\epsilon ^{2} \ll {\it{{\it{Pe}}}}_r \ll 1$
, transitions to a
${\it{{\it{Pe}}}}_r^{ {10}/{3}}$
-scaling across
${\it{{\it{Pe}}}}_r \sim O(1)$
, before eventually decreasing as
${\it{Pe}}^{-(C-5)}$
for
${\it{{\it{Pe}}}}_r \gg 1$
.
$D_{\it{eff}}$
attains a maximum at the transition between the latter two regimes. For slender swimmers, the transition
${\it{{\it{Pe}}}}_r$
is
$O(\kappa ^3)$
, implying a maximum
$D_{\it{eff}}$
of
$O(\kappa ^{10}D_t)$
, that serves as a flow-independent bound on slender-swimmer dispersion. As mentioned earlier, achieving this asymptotic scaling requires the transition
${\it{{\it{Pe}}}}_r$
to scale as
$\kappa ^3$
, and this requires going to
$\kappa$
values greater than
$12$
.
It is also worth commenting on the effect of an independent source of translational diffusion on the dispersion scalings above. Such a source would correspond to a term,
$\hat {D}_t {\nabla} ^2_{\boldsymbol{x}} \varOmega$
, in the governing kinetic equation (2.2), with
$\hat {D}_t$
assumed to be independent of
$D_r$
and the ambient shear flow. The effect of such a term on migration dynamics was analysed in Vennamneni et al. (Reference Vennamneni, Nambiar and Subramanian2020b
), where its inclusion was shown to lead to a uniform swimmer concentration profile for sufficiently large
${\it{{\it{Pe}}}}_r$
(see figures 15 and 16 in Vennamneni et al. (Reference Vennamneni, Nambiar and Subramanian2020b
)). This may be understood from the fact that, for large
${\it{{\it{Pe}}}}_r$
, the translational diffusivity (
$D_{zz}$
), arising from swimming combined with rotary diffusion, decreases with
${\it{{\it{Pe}}}}_r$
regardless of the trapping regime – as
$D_t {\it{{\it{Pe}}}}_r^{-{4}/{3}}$
in the HST regime, and as
$D_t {\it{{\it{Pe}}}}_r^{-2}$
in the LST one. Therefore, the ratio
$\hat {D}_t/D_{zz}$
is either
$O(\hat {D}_tD_r/U_s^2\times {{\it{Pe}}}_r^{ {4}/{3}})$
or
$O(\hat {D}_t.D_r/U_s^2\times{{\it{Pe}}}_r^{2})$
, depending on the trapping regime, and becomes asymptotically large for any non-zero
$\hat {D}_t$
. This implies that it is passive (rather than active) diffusion that will resist the drift (
$V_z$
) driving swimmer migration for large
${\it{{\it{Pe}}}}_r$
, and since
$V_z$
also decreases in a manner analogous to
$D_{zz}$
, with increasing
${\it{{\it{Pe}}}}_r$
, the passive diffusive flux will eventually render the swimmer concentration uniform in the transverse direction. This uniformity is analogous to the passive solute, and implies that, for any non-zero
$\hat {D}_t$
,
$D_{\it{eff}}/D_t$
must revert to a passive
${\it{{\it{Pe}}}}_r^2$
-scaling regime for sufficiently large
${\it{{\it{Pe}}}}_r$
(Peng & Brady Reference Peng and Brady2020). Observation of the departure from classical scaling for swimmer dispersion – either the
${\it{{\it{Pe}}}}_r^4$
-regime for
$\kappa \sim O(1)$
, or the
${\it{{\it{Pe}}}}_r^{10/3}$
-regime for large
$\kappa$
-requires that
$\hat {D}_tD_r/U_s^2\times{{\it{Pe}}}_r^2 \lesssim 1$
for
$\kappa \sim O(1)$
, or
$\hat {D}_tD_r/U_s^2\times{{\it{Pe}}}_r^{ {4}/{3}} \lesssim 1$
for
$\kappa \gg 1$
at the
${\it{{\it{Pe}}}}_r$
corresponding to the onset of anomalous scaling regime.
Having summarized the scaling regimes of
$D_{\it{eff}}$
above, we are in a position to examine the assumption of small
$\varLambda Pe$
, made in the second stage of the multiple scales analysis, in § 2.2.2. For swimmers with
$\kappa \lesssim 2$
, one has
$D_{zz} \sim D_t {\it{{\it{Pe}}}}_r^{-2}$
,
$D_{\it{eff}} \sim D_t({\it{{\it{Pe}}}_r^4}/{\epsilon ^4})$
that, when used in
$H^2/D_{zz} \sim L_{\it{min}}^2/D_{\it{eff}}$
, gives
$\varLambda Pe \sim O({\it{{\it{Pe}}}}_r^{-2})$
. For swimmers with
$1 \ll {\it{{\it{Pe}}}}_r \ll \kappa ^3$
(infinitely slender swimmers being included in this category), one has
$D_{zz} \sim D_t {\it{{\it{Pe}}}}_r^{- {4}/{3}}$
,
$D_{\it{eff}} \sim D_t({\it{{\it{Pe}}}_r^{ {10}/{3}}}/{\epsilon ^4})$
, leading to
$\varLambda Pe \sim O({\it{{\it{Pe}}}}_r^{- {4}/{3}})$
. Finally, for swimmers with
$1 \ll \kappa ^3 \ll {\it{{\it{Pe}}}}_r$
that conform to the anomalous scaling regime, one has
$D_{zz} \sim D_t\kappa ^2 {\it{{\it{Pe}}}}_r^{-2}$
. In this case,
$L_{\it{min}}$
may be directly estimated as the distance travelled along the flow direction with speed
$U_m$
in a time of
$O(H^2/D_{zz})$
, which gives
$\varLambda \sim \epsilon ^2\kappa ^2/{\it{{\it{Pe}}}}_r^3$
and
$\varLambda Pe \sim \kappa ^2/{\it{{\it{Pe}}}}_r^2$
; note that the
$L_{\min }$
estimates for the above cases are
$\it{HPe Pe}^{2}_{r}$
,
$\it{HPePe}^{4/3}_{r}$
and
$HPe {\it{{\it{Pe}}}}_r^2/\kappa ^2$
, respectively, and may be compared with the estimate for the passive solute obtained in § 2.1. Importantly, for all regimes examined,
$\varLambda Pe \ll 1$
for
${\it{{\it{Pe}}}}_r \gg 1$
, thereby validating the separation of the flow and gradient-direction evolutions in § 2.2.2.
4. Conclusions and discussion
In this work, we have characterized the influence of activity and anisotropy on Taylor dispersion of prolate spheroidal swimmers in plane Poiseuille flow. Activity enhances dispersion but importantly, can also reduce it above a threshold aspect ratio, with the reduction being driven by the singular collapse of swimmers towards the channel centreline for large
${\it{{\it{Pe}}}}_r$
. This collapse manifests as an anomalous dispersion exponent, determined by the swimmers experiencing large velocity fluctuations during rare excursions away from the main population that is confined in an
$O({\it{{\it{Pe}}}}_r^{-1})$
central core. The collapse becomes sufficiently intense for slender swimmers, leading to
$D_{\it{eff}}$
attaining a maximum at a finite
${\it{{\it{Pe}}}}_r$
, and whose magnitude is determined by
$\kappa$
and the intrinsic swimmer diffusivity
$D_t$
; dispersion in such cases may be regarded as being limited by the swimmer geometry, rather than the maximum attainable flow rate. While the importance of tail events being responsible for the anomalous scaling of higher-order moments is known for passive Taylor dispersion (Haynes & Vanneste Reference Haynes and Vanneste2014), and for other problems in non-equilibrium statistical physics (Sreenivasan & Antonia Reference Sreenivasan and Antonia1997), their role in determining the second moment (via
$D_{\it{eff}}$
) itself, for active swimmers, is new.
For
$\kappa \sim O(1)$
, departure from the classical
${\it{{\it{Pe}}}}_r^4$
-scaling regime sets in at
${\it{{\it{Pe}}}}_r \lesssim 100$
, as indicated by the blue vertical band in figure 4, and is achievable in experiments. The experiments of Rusconi et al. (Reference Rusconi, Guasto and Stocker2014) and Dehkharghani et al. (Reference Dehkharghani, Waisbord, Dunkel and Guasto2019) have accessed maximum
${\it{{\it{Pe}}}}_r$
values of
$50$
and
$100$
, respectively, while those of Barry et al. (Reference Barry, Rusconi, Guasto and Stocker2015) extend to a
${\it{{\it{Pe}}}}_r$
of
$167$
. Although attaining the anomalous scaling regime, characterized by a fixed but reduced dispersion exponent, requires one to go to higher
${\it{{\it{Pe}}}}_r$
, the reduction in
$D_{\it{eff}}$
is already significant at the smaller
${\it{{\it{Pe}}}}_r$
quoted above, that correspond to the transition region between the classical and anomalous regimes (the aforesaid blue band). Referring again to figure 4, one finds the onset of centreline collapse to lead to a
$55 \, \%$
reduction in
$D_{\it{eff}}$
at
${\it{{\it{Pe}}}}_r=80$
, and a
$70\, \%$
reduction at
${\it{{\it{Pe}}}}_r=167$
, relative to the (hypothetical) values obtained by assuming the classical
${\it{{\it{Pe}}}}_r^4$
-regime to continue to indefinitely large
${\it{{\it{Pe}}}}_r$
.
Langevin simulation results from Ganesh et al. (Reference Ganesh, Douarche, Dentz and Auradou2023) depicting: (a) the normalized flow-velocity variance (
$\sigma ^{2}_{U_{p}}=\textstyle \int _{-1}^1 {u_1^{\prime}}^2(z)F_s(z) {\rm d}z$
),
$U$
being the average flow velocity, and (b) the longitudinal diffusivity (denoted as
$D_{\parallel }$
in (b i), with the intrinsic swimmer translational diffusivity denoted as
$D_0$
), as functions of
$P\!e$
for ABPs of various aspect ratios (denoted
$q$
here);
${\it{Pe}}=25 {\it{{\it{Pe}}}}_r$
. RBC in the legend refers to reflecting boundary conditions. Panel (b ii) shows the longitudinal effective diffusivity
$D_{\it{eff}}/D_{t}$
plotted as a function of
${\it{{\it{Pe}}}}_r$
, for
$\kappa =4$
, with the anomalous asymptote (
${\it{Pe}}^{5-C}_r$
) as a dash–dotted black line. The dashed blue line represents the exponent (
$3.75$
), observed by Ganesh et al. above, and obtained from fitting the numerical results over the interval
${\it{{\it{Pe}}}}_r \sim 10{-}50$
(marked by vertical dashed lines) accessed in the simulations.
Interestingly, the recent simulations of Ganesh et al. (Reference Ganesh, Douarche, Dentz and Auradou2023), of spheroidal ABPs in the pressure-driven channel flow configuration, appear to confirm the centreline collapse of swimmers with
$\kappa \gtrsim 2$
, and the resulting anomalous dispersion. We reproduce two of the figures of Ganesh et al. (Reference Ganesh, Douarche, Dentz and Auradou2023) as figures 6(a) and 6(b). The former plots the flow-velocity variance,
$\textstyle \int _{-1}^1 {u_1^{\prime}}^2(z)F_s(z) {\rm d}z$
, as a function of
${\it{{\it{Pe}}}}_r$
for swimmers of different aspect ratios. The velocity variance is closely related to the profile variance,
$\langle z^2 \rangle$
, plotted in figures 2(c) and 3(c). In particular, the centreline-collapse of swimmers should manifest as a continuous decrease of both the profile and flow-velocity variances with increasing
${\it{{\it{Pe}}}}_r$
. As highlighted in figure 6(a), a slow decrease of the velocity variance (not recognized by the authors) is indeed seen for
$\kappa =4$
, starting from
${\it{Pe}} \lesssim 500$
, which translates to
${\it{{\it{Pe}}}}_r \lesssim 20$
(using the relation
${\it{Pe}} = 25{\it{{\it{Pe}}}}_r$
stated on pages 5 and 9 of Ganesh et al. (Reference Ganesh, Douarche, Dentz and Auradou2023)). Also evident from figure 6(a) is the exceedingly slow approach, of the velocity variance for
$\kappa =2$
, to the infinite-
${\it{{\it{Pe}}}}_r$
plateau, a behaviour that is expected for aspect ratios in the vicinity of the collapse threshold; recall that the profile variance is predicted to only decay as
$(\ln {\it{{\it{Pe}}}}_r)^{-1}$
at the threshold aspect ratio. The aforementioned decrease in the velocity variance, for
$\kappa = 4$
, leads to
$D_{\it{eff}}$
increasing with an exponent smaller than
$4$
in figure 6(bi). The reduced
$D_{\it{eff}}$
values in this figure, for
$\kappa = 10$
, are also consistent with the lower exponent of
$10/3$
predicted here for sufficiently slender swimmers, see figure 5(a). It is worth noting that, despite the
$D_{\it{eff}}$
in figure 6(b) growing at a rate slower than
${\it{{\it{Pe}}}}_r^4$
for
$\kappa = 10$
, the corresponding velocity variance in figure 6(a) exhibits a clear large-
${\it{{\it{Pe}}}}_r$
plateau, confirming the absence of a centreline-collapse. It is this difference between
$\kappa = 4$
and
$10$
that constitutes the clearest evidence for anomalous scaling in the former case. In fact, the reduced exponent (
$3.75$
) for
$\kappa =4$
in figure 6(b) corresponds to the intermediate asymptotic regime highlighted via figure 5(b). Although less than
$4$
, it does not represent the actual anomalous exponent, owing to
${\it{{\it{Pe}}}}_r$
in the above simulations not being large enough to access the true asymptotic regime. The true anomalous exponent for
$\kappa = 4$
is much lower, being given by
$5-C \approx 2.31$
; see figure 6(bii).
Several earlier efforts have focused on near-wall accumulation of swimmers, and the associated upstream swimming (Hill et al. Reference Hill, Kalkanci, McMurry and Koser2007; Elgeti & Gompper Reference Elgeti and Gompper2009; Costanzo et al. Reference Costanzo, Di Leonardo, Ruocco and Angelani2012; Kaya & Koser Reference Kaya and Koser2012; Elgeti & Gompper Reference Elgeti and Gompper2013; Chilukuri et al. Reference Chilukuri, Collins and Underhill2015; Ezhilan & Saintillan Reference Ezhilan and Saintillan2015; Figueroa-Morales et al. Reference Figueroa-Morales, Mino, Rivera, Caballero, Clément, Altshuler and Lindner2015). Recent work has served to further reinforce the subtle orientation dependence of swimmer–boundary interactions, with both hydrodynamic and non-hydrodynamic factors playing important roles (Molaei et al. Reference Molaei, Barry, Stocker and Sheng2014; Contino et al. Reference Contino, Lushi, Tuval, Kantsler and Polin2015; Buchner et al. Reference Buchner, Muller, Mehmood and Tam2021; Zeng, Jiang & Pedley Reference Zeng, Jiang and Pedley2022). While it is desirable to incorporate the results of such investigations, in the form of more detailed position–orientation space boundary conditions for the governing kinetic equation (Jiang & Chen Reference Jiang and Chen2019; Wang et al. Reference Wang, Jiang, Chen and Tao2022), our work, entailing simple no-flux boundary conditions for a drift–diffusion equation in physical space, nevertheless highlights independent bulk mechanisms that drive migration and influence dispersion. Migration away from boundaries, driven by bulk mechanisms, can potentially affect incipient biofilm formation, and the resulting dispersion can affect the spread of microbial populations in wide channels. In general, one expects both wall-induced and bulk mechanisms to be relevant, with the former confining a fraction of the swimmer population to thin wall layers. Importantly, for a fixed swimmer concentration, the fraction of near-wall swimmers is expected to decrease with increasing channel width and also with increasing flow rate; swimmers outside these layers must disperse in the manner described here. The above expectations are broadly consistent with the findings of Maretvadakethope et al. (Reference Maretvadakethope, Hazel, Vasiev and Bearon2023) who have compared swimmer concentration profiles obtained from a continuum formulation, using a range of boundary conditions, with those obtained from simulations of individual swimmers. Among other things, the authors find the swimmer concentration in the bulk to remain insensitive to a change from a specular-reflection (equal incident and reflected angles) to a diffuse-reflection (uniform distribution of reflected angles) boundary condition; this being especially true for large Péclet numbers.
In addition to the aspect-ratio dependence, the dispersion behaviour also depends sensitively on the overall swimmer geometry – that is, on whether swimmer has a rod-like or a disk-like morphology. While discoid swimmers (modelled as oblate spheroids, with
$\kappa \lt 1$
, that swim along the symmetry axis) exhibit a reverse transition from low to HST with increasing
${\it{{\it{Pe}}}}_r$
(Vennamneni & Subramanian Reference Vennamneni and Subramanian2025), there is not an analogous wall-collapse for
${\it{{\it{Pe}}}}_r \rightarrow \infty$
. As a result, in contrast to the dispersion of slender-rod swimmers shown in figure 5, the
$D_{\it{eff}}$
for flat-disk swimmers with
$\kappa \ll 1$
transitions from a zero-
${\it{{\it{Pe}}}}_r$
plateau, via an intermediate asymptotic regime with
$D_{\it{eff}}/D_t\!\propto \!{\it{{\it{Pe}}}}_r^{ {8}/{3}}$
, to an eventual
${\it{{\it{Pe}}}}_r^4$
-scaling regime – importantly, all of these exponents conform to classical scaling arguments (Vennamneni & Subramanian Reference Vennamneni and Subramanian2025). The above is also consistent with the results for oblate swimmers with
$\kappa = 0.5$
in figure 6(a), which show the emergence of a velocity-variance plateau at
${\it{{\it{Pe}}}}_r$
values much lower than the complementary prolate swimmer with
$\kappa = 2$
.
While we have only examined rotary diffusing microswimmers (i.e. ABPs) here, naturally occurring swimmers exhibit a combination of run-and-tumble dynamics and rotary diffusion (Berg Reference Berg1993). In contrast to the small-amplitude orientation changes characterizing rotary diffusion, tumbles lead to instantaneous and impulsive changes in swimmer orientation, with the distribution of tumble events being well approximated by a Poisson process (Berg Reference Berg1993; Subramanian & Koch Reference Subramanian and Koch2009; Subramanian & Nott Reference Subramanian and Nott2011). The importance of tumbling vis-à-vis rotary diffusion may be characterized by the parameter
$\tau D_r$
, with
$\tau$
being the average time interval between successive tumbles (the average duration of an individual run). Here
$\tau D_r \rightarrow 0$
corresponds to the ABPs examined here, while
$\tau D_r = 0$
corresponds to RTP where tumbles alone control the orientation decorrelation. Our earlier work has shown that shear-induced migration exhibits a transition from HST to LST, with increasing
${\it{{\it{Pe}}}}_r$
, regardless of
$\tau D_r$
. The difference in the intrinsic swimming mechanism only manifests as a change in the
$\kappa$
-scaling of the boundary between the trapping regimes; in contrast to that shown in figure 1(a), the transition
${\it{{\it{Pe}}}}_r$
scales as
$O(\kappa )$
, for large
$\kappa$
, for RTP. Importantly, our preliminary results also show a centreline collapse behaviour for RTP, although the threshold aspect ratio for the collapse to occur appears much higher than that for ABPs. Nevertheless, it appears reasonable to conclude that the dispersion regimes examined here are robust to the incorporation of run-and-tumble dynamics (
$\tau D_r \sim O(1)$
) in the swimmer motion.
The dispersion behaviour found here has implications beyond plane Poiseuille flow. Onset of anomalous scaling, and geometrically limited dispersion for large
$\kappa$
, should occur in other unidirectional shearing flows, although the exponents may differ. To the extent that the plane Poiseuille profile is a local approximation of a velocity extremum, dispersion of rod-like swimmers in a general unidirectional shearing flow might be limited by collapse onto velocity extrema, as illustrated in figure 7. A consequence of this collapse, especially for large
${\it{{\it{Pe}}}}_r$
, is the emergence of hysteresis in migration and dispersion. For instance, in figure 7, one expects the number of swimmers localized in the neighbourhood of a given extremum to depend on the number of swimmers initially present within the basin of attraction of the extremum, with this basin being defined as the region between the locations corresponding to the maximum slope (magnitude) on either side of the extremum. The number of swimmers in this basin clearly depends on the initial distribution, implying that the population distribution among the different swimmer clusters that form after the collapse, and thence, the resulting dispersion of the swimmer population, is a function of the initial swimmer distribution. However, this is in contrast to the expression for the steady swimmer concentration profile, given in (2.31), that only depends on the drift and diffusivity coefficients, each of which is determined only by the quasisteady local equilibrium that develops at a given position (
$z$
) in response to the local shear, and is therefore independent of initial conditions. The implication of the above is that the initial-condition-dependent localization, driven by collapse onto the nearest extremum, only leads to a quasisteady state, and that this must be followed for long times (scaling exponentially with
${\it{{\it{Pe}}}}_r$
) by barrier-hopping dynamics, on a potential given by
$-\textstyle \int ^z ({V_z(z')}/{D_{zz}(z')})\,{\rm d}z'$
, that acts to equilibrate the populations in the different swimmer clusters, which finally leads to initial-condition-independent dynamics for sufficiently long times.
A schematic depicting slender-swimmer migration, at large
${\it{{\it{Pe}}}}_r$
, in a general unidirectional shearing flow. The original population is apportioned into multiple clusters localized in the neighbourhood of the velocity extrema. The subsequent dispersion behaviour will be influenced by barrier-hopping dynamics between clusters.
The non-monotonic dependence of slender swimmer dispersion on
${\it{{\it{Pe}}}}_r$
, illustrated in figure 5, also occurs for the dispersion of gyrotactic swimmers in pressure-driven channel flow (Bees & Croze Reference Bees and Croze2010; Croze et al. Reference Croze, Sardina, Ahmed, Bees and Brandt2013), although there has been no detailed analysis of scaling exponents. Unlike our findings, gyrotactic dispersion depends sensitively on channel orientation, with non-monotonicity being most pronounced for vertically down-welling flow, owing to the centreline-focusing that occurs in this configuration (Kessler Reference Kessler1985). The
$\kappa$
-dependence of the dispersion differs, with elongated swimmers exhibiting a weaker non-monotonicity. It would be of interest to examine dispersion under the combined effects of gyrotaxis and centreline collapse.
Finally, Taylor, in a later paper, also examined dispersion of a solute in turbulent pipe flow (Taylor Reference Taylor1954). For swimmer suspensions, increasing swimmer concentration and/or activity leads to ‘active turbulence’ or ‘collective motion’ – a novel spatiotemporally chaotic state with negligible fluid inertia (Saintillan & Shelley Reference Saintillan and Shelley2007; Subramanian & Koch Reference Subramanian and Koch2009; Koch & Subramanian Reference Koch and Subramanian2011; Krishnamurthy & Subramanian Reference Krishnamurthy and Subramanian2015; Martinez et al. Reference Martinez2020; Alert, Casademunt & Joanny Reference Alert, Casademunt and Joanny2022). It would evidently be of interest to examine how the anomalous dispersion found here is affected by the onset of active turbulence. Also of interest is to extend the present formulation to account for the role of swimmer-handedness on orientation dynamics identified in recent efforts (Jing et al. Reference Jing, Zöttl, Clément and Lindner2020; Zöttl et al. Reference Zöttl, Tesser, Matsunaga, Laurent, Du Roure and Lindner2023).
Funding
L.V. would like to thank the financial support of GITAM University (grant no. 2023/0237). G.S. would like to acknowledge financial support from Department of Biotechnology (DBT) grant in Life Science Research, Education and Training at JNCASR (grant no. BT/INF/22/SP27679/2018).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical evaluation of the longitudinal dispersion coefficient
All of the plots in §§ 2 and 3, including the concentration profiles and longitudinal diffusivities, have been determined from a numerical evaluation of the underlying orientation distributions (
$G(\boldsymbol{p})$
,
$\varOmega _{11}$
and
$\varOmega _{12}$
). In (2.11), for instance, the unknown probability density
$G(\boldsymbol{p})$
is expressed as a truncated expansion in spherical harmonics in the form
\begin{equation} G(p)=\sum _{l=0}^{N}\sum _{m=-l}^{l}a_{l,m}Y_{l}^{m}(\boldsymbol{p}), \end{equation}
where
$Y_{l}^{m}(\boldsymbol{p})= {\sqrt {(({(2l+1)}/{4\pi })({(l-m)!}/{{(l+m)!}}))}P^m_{l}(\cos \theta ) \exp (\mathrm{i}m\phi )}$
represents the spherical harmonic of order
$l$
and degree
$m$
;
$\theta$
and
$\phi$
being the polar and azimuthal angles of the spherical coordinate system (Abramowitz & Stegun Reference Abramowitz and Stegun1965). Substitution of this expansion, and use of the orthogonality property of the spherical harmonics (
$\textstyle \int {\rm d}\boldsymbol{p} Y_{l}^{m}(\boldsymbol{p}) Y_{l^{\prime }}^{{{m}^{\prime }}*}(\boldsymbol{p}) = \delta _{ll'}\delta _{mm'}$
, where
$Y_{l^{\prime }}^{{{m}^{\prime }}*}(\boldsymbol{p})$
represents the conjugate spherical harmonic) leads to a system of linear equations for the unknown coefficients
$a_{l,m}$
(Vennamneni et al. Reference Vennamneni, Nambiar and Subramanian2020b
). Along similar lines, using the expansions
\begin{equation} \varOmega _{11}(p)=\displaystyle \sum _{l=0}^{N}\sum _{m=-l}^{l}b_{l,m}Y_{l}^{m}(\boldsymbol{p}) \end{equation}
and
\begin{equation} \varOmega _{12}(p)=\displaystyle \sum _{l=0}^{N}\sum _{m=-l}^{l}c_{l,m}Y_{l}^{m}(\boldsymbol{p}) \end{equation}
in (2.8), again leads to linear systems that may be solved to obtain
$b_{l,m}$
and
$c_{l,m}$
. The resulting expansions for
$\varOmega _{11}$
and
$\varOmega _{12}$
are then used in (2.20) and (2.19) to calculate
$D_{zz}$
and
$V_z$
. We choose
$N=80$
that ensures that the above spherical harmonics expansions converge, in turn leading to converged results for the integrals defining
$D_{zz}$
and
$V_z$
. We then evaluate the integral in the expression for the steady state concentration profile
$F_s(z;{\it{{\it{Pe}}}}_r,\kappa )$
, given by (2.31), using Gauss–Legendre quadrature, taking care to increase the number of quadrature points while also concentrating them in the vicinity of the centreline, so as to obtain converged results with increasing
${\it{{\it{Pe}}}}_r$
. For
$\epsilon =0.1$
, the integrals involved in
$D_{\it{eff}}$
, as given by (2.45), are again evaluated using Gauss–Legendre quadrature, taking particular care to ensure adequate resolution of the shrinking (collapsed) region near the centreline with increasing
${\it{{\it{Pe}}}}_r$
.
Figure 8(a) depicts the longitudinal diffusivities for multiple aspect ratios (
$\kappa = 1$
,
$1.4$
and
$1.9$
) below
$\kappa \approx 2$
. The main plot shows the expected transition to a
${\it{{\it{Pe}}}}_r^4$
-scaling regime for sufficiently large
${\it{{\it{Pe}}}}_r$
, with the curves for different aspect ratios seemingly coincident. The separation between the different curves, owing to
$\kappa$
-dependence of the prefactor in the
${\it{Pe}}^{4}_r$
-asymptote (which determines the vertical position of the asymptote) is seen in the top inset, however, which provides a magnified view of a subset of the large-
${\it{{\it{Pe}}}}_r$
interval; of interest here is the continued discrepancy between the numerics and the large-
${\it{{\it{Pe}}}}_r$
asymptote, for
$\kappa = 1.9$
, which is evidence for the
${\it{{\it{Pe}}}}_r^4$
-regime receding to infinity with approach towards the critical aspect ratio; the intermediate asymptotic regime is expected to scale as
${\it{{\it{Pe}}}}_r^4/\ln {\it{{\it{Pe}}}}_r$
. The bottom inset shows the prefactor above to be a non-monotonic function of
$\kappa$
. It starts at 253.9 for
$\kappa =1$
, increases to
$321.9$
for
$\kappa = 1.2$
, attains a maximum of
$354$
for
$\kappa \approx 1.4$
, and then decreases as
$\kappa$
continues to increase, eventually reaching zero at
$\kappa \approx 2$
. Figure 8(b) highlights the onset of anomalous scaling for aspect ratios immediately above
$2$
, as evident from the deviation of
$D_{\it{eff}}/D_t$
from the classical
${\it{Pe}}^{4}_{r}$
asymptote. As indicated in the main manuscript, this deviation begins at
${\it{{\it{Pe}}}}_r \approx 100$
.
(a) Here
$D_{\it{eff}}/D_t$
versus
${\it{{\it{Pe}}}}_r$
, with
$\epsilon =0.1$
, for
$\kappa =1$
,
$1.4$
and
$1.9$
. The black dashed curve denotes the passive asymptote, obtained from (2.45) with
$D_{zz}=1/6$
,
$F_{s}=1/2$
; the top inset compares the large-
${\it{{\it{Pe}}}}_r$
variation, for the said aspect ratios, with the corresponding asymptotes (red, blue and green dashed lines) which exhibit the classical
${\it{{\it{Pe}}}}_r^4$
-scaling – these are obtained from (3.1) using
$C=0,C_1\!\approx \!8/3$
for
$\kappa =1$
;
$C\! \approx \!0.39,C_1\!\approx \!1.85$
for
$\kappa =1.4$
;
$C\!\approx \!0.9,C_1\!\approx \!1.69$
for
$\kappa =1.9$
; the bottom inset shows variation of the asymptotic prefactor,
$D^{*}=(D_{\it{eff}}\vert _{\it{{\it{Pe}}}_{r}\rightarrow \infty }/D_{t})Pe^{-4}_{r}$
, with
$\kappa$
. (b)
$D_{\it{eff}}/D_t$
versus
${\it{{\it{Pe}}}}_r$
with
$\epsilon =0.1$
, for
$\kappa =2.4$
, 3 and 3.6. The ‘classical’
${\it{{\it{Pe}}}}_r^4$
and the passive
${\it{Pe}}^{2}_r$
-asymptotes appear as a dashed red line and a dashed black curve, respectively; the anomalous
$\kappa$
-dependent asymptotes (
$\propto {\it{{\it{Pe}}}}_r^{5-C}$
) appear as dash–dotted black lines, with the
$C$
values indicated in the legend.
Appendix B. Langevin simulations of ABP swimmers
In order to verify the theoretical results presented in § 3, we have carried out stochastic simulations to independently determine the longitudinal diffusivity of a population of swimmers, and its variation with
${\it{{\it{Pe}}}}_r$
. Owing to the computational expense involved, particularly at large
${\it{{\it{Pe}}}}_r$
, these simulations have only been performed for a limited number of cases, the focus being largely on spherical swimmers (
$\kappa =1$
) with
$\epsilon = U_s/(H D_r) = 0.2$
. The simulations involve time-integration of the Langevin equations governing the positions (
$\boldsymbol{x}_i$
) and orientations (
$\boldsymbol{p}_i$
) of
$N$
microswimmers. Interactions between swimmers are neglected, in which case the Langevin equations for the
$i^{\mathrm{th}}$
microswimmer are given by
\begin{align} \frac {{\rm d}\boldsymbol{x}_i}{{\rm d}t} &= U_s \boldsymbol{p} + \boldsymbol{u}, \end{align}
\begin{align} \frac {{\rm d}\boldsymbol{p}_i}{{\rm d}t} &= B\left [ \boldsymbol{E}\boldsymbol{\cdot } \boldsymbol{p}-\boldsymbol{p}(\boldsymbol{E}:\boldsymbol{p}\boldsymbol{p})\right ]+\boldsymbol{W}\boldsymbol{\cdot }\boldsymbol{p} + D_r \boldsymbol{T} \times \boldsymbol{p} \\[9pt] \nonumber \end{align}
where
$\boldsymbol{T}$
is the delta-correlated fluctuating torque acting to change the swimmer orientation (Dhont Reference Dhont1996); these equations are consistent with the convection and rotary diffusion terms in (2.1). We chose
$N = 1000$
for the results presented here, but have verified that the results are independent of
$N$
. The time evolution is done using a second-order Runge–Kutta method, with the time step being changed depending on the specific parameter values, in particular,
${\it{{\it{Pe}}}}_r$
. More details with regard to these stochastic simulations are given in Vennamneni et al. (Reference Vennamneni, Nambiar and Subramanian2020b
).
In addition to confirming the theoretical predictions pertaining to long-time evolution, the simulations also allow access to the transient behaviour prior to the long-time diffusive regime. The transient dynamics offer an interesting contrast to the passive case, although this difference has not been recognized in the literature; for instance, the transient dynamics is fitted with an exponential ansatz in Ganesh et al. (Reference Ganesh, Douarche, Dentz and Auradou2023). The nature of both the short-time (
$\langle z^2 \rangle \sim Dt, \langle x^2 \rangle \sim ({(U_mD)^2}/{H^4}t^4)$
for
$t \ll H^2/D$
) and long-time (
$\langle z^2 \rangle \sim O(H^2), \langle x^2 \rangle \sim (Pe^2 D)t$
for
$t \gg H^2/D$
) evolutions of the gradient and flow-aligned variance components, in the passive case, have already been described in § 1. In order to examine the transient evolution of these variance components in the active case, the initial swimmer population is chosen to be localized in the neighbourhood of the centreline: thus, at
$t= 0$
, the
$1000$
swimmers are uniformly distributed in a cuboidal volume defined by
$(0,H)$
in the axial direction,
$(0.98 H, 1.02 H)$
in the gradient direction and
$(0, H)$
in the vorticity direction;
$H$
being the channel half-width. Periodic boundary conditions are used in both gradient and vorticity directions, while there are no boundaries in the flow direction since we aim to determine the long-term dispersion. We have also verified that the aforementioned transient scalings in the passive limit are achieved by setting
$D_r \gg 1$
for
$H$
and
$U_s$
fixed.
Figures 9(a) and 9(b) show the swimmer variances along the gradient (
$\langle z^2 \rangle$
) and flow (
$\langle x^2 \rangle$
) directions as a function of
$t$
, for the highest
${\it{{\it{Pe}}}}_r\,(= 100)$
, and for different choices of the time step. The latter serves to illustrate numerical convergence; since the simulations are for
$\kappa = 1$
, the approach of
$\langle z^2 \rangle$
towards the value (
$1/3$
), corresponding to a uniform distribution in the gradient direction, is a good indication of convergence. Next, we compare the numerical results for
$D_{\it{eff}}$
obtained from the long-time diffusive regime in figure 9(b), with the theoretical value given by (2.45), for varying
${\it{{\it{Pe}}}}_r$
. The numerical and analytical results for the axial diffusivity (
$D_{\it{eff}}/D_{t}$
), in table 1, show good agreement for all
${\it{{\it{Pe}}}}_r$
examined. Table 2 shows similarly good agreement between the theoretical and numerical axial diffusivities for non-spherical swimmers (
$\kappa$
varying from
$2-20$
) at
${\it{{\it{Pe}}}}_r=1$
. It is worth emphasizing here that the axial diffusivity results quoted are for
$\epsilon = 0.2$
, as opposed to those in the main manuscript which are for
$\epsilon = 0.1$
, and show that the theory for wide channels remains valid at least until
$\epsilon = 0.2$
.
The variances characterizing the spherical swimmer distribution in the (a) flow (
$x$
) and (b) gradient (
$z$
) directions, as a function of time, for
${\it{{\it{Pe}}}}_r = 100$
; the different curves correspond to different choices for the time step used in the numerical integration.
Non-dimensional axial diffusivity (
$D_{\it{eff}}/D_t$
) results of numerical (Langevin simulations) and theory ((2.45) in the manuscript) for spherical swimmers (
$\kappa =1$
) obtained at different values of
${\it{{\it{Pe}}}}_r$
with
$\epsilon =0.2$
.
Non-dimensional axial diffusivity (
$D_{\it{eff}}/D_t$
) results of numerical (Langevin simulations) and theory ((2.45) in the manuscript) for different values of
$\kappa$
at
${\it{{\it{Pe}}}}_r=1$
with
$\epsilon =0.2$
.
Finally, it is of interest to note that, at intermediate times (that is, times much longer than the initial ballistic regime but shorter than those corresponding to the diffusive regime),
$\langle x^2 \rangle \sim t^3$
and
$\langle z^2 \rangle \sim t^{ {1}/{2}}$
for
${\it{{\it{Pe}}}}_r =100$
(see figures 9
a and 9
b). These intermediate asymptotic scalings arise for all large
${\it{{\it{Pe}}}}_r$
, provided one starts with an initial condition localized around the centreline, and owe their origin to the shear-limited diffusion (
$D_{zz} \propto D_t {\it{{\it{Pe}}}}_r^{-2}$
) along the gradient direction at large
${\it{{\it{Pe}}}}_r$
. To see this, note that
$ z^2 \propto D_{zz}t$
with
$D_{zz} \sim D_t ({\it{{\it{Pe}}}}_r z/H)^{-2}$
,
${\it{{\it{Pe}}}}_r z/H$
being the local Péclet number, which finally leads to
$\langle z^2 \rangle \propto H {\it{{\it{Pe}}}}_r^{-1} (D_tt)^{ {1}/{2}}$
. Using this scaling (
$z \sim H^{ {1}/{2}}(D_t t)^{ {1}/{4}}{\it{{\it{Pe}}}}_r^{- {1}/{2}}$
) in
$x \sim U_m(z/H)^2 t$
leads to
$\langle x^2 \rangle \sim {\it{{\it{Pe}}}}_r^{-2}(U_m^2D_t/H^2)t^3$
. These intermediate-time asymptotes ought to be incorporated in any empirical forms used to fit simulation results for the mean-squared displacements; this, as opposed to the exponential ansatz that is normally used in such fits, an example of which was quoted above.
