1. Introduction
Oscillatory flows occur in various aquatic and biological systems, including tidal and wind-driven currents (Bars, Cébron & Gal Reference Bars, Cébron and Gal2015; Wang, Li & Dong Reference Wang, Li and Dong2021), tributary inflows subject to diurnal discharge regulation (Long et al. Reference Long, Ji, Yang, Cheng, Yang, Liu, Liu and Lorke2020), and biological flows (Secomb Reference Secomb2017; Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019; Alaminos-Quesada et al. Reference Alaminos-Quesada, Gutiérrez-Montes, Coenen and Sánchez2024), as well as in engineered fluidic systems (Reis et al. Reference Reis, Gonçalves, Aguedo, Gomes, Teixeira and Vicente2006; Vedel, Olesen & Bruus Reference Vedel, Olesen and Bruus2010; Hacioglu & Narayanan Reference Hacioglu and Narayanan2016). Compared with steady flows, mass transfer in oscillatory flows is more complex. A substance patch can alternately expand and contract as the advective velocity profile develops, decays and reverses. The extent of streamwise dispersion in this process depends on cross-sectional passive diffusion or active migration, which determines the ability to exploit differential advection between streamlines. It is therefore of interest to investigate how active particles, whose cross-sectional migration is jointly influenced by shear and motility (Bees & Croze Reference Bees and Croze2010; Bearon, Bees & Croze Reference Bearon, Bees and Croze2012; Croze et al. Reference Croze, Sardina, Ahmed, Bees and Brandt2013, Reference Croze, Bearon and Bees2017), disperse differently from solute particles, whose cross-sectional migration is purely diffusive and independent of the flow profile under laminar conditions.
Research on solute dispersion in oscillatory flows began shortly after the seminal work of Taylor (Reference Taylor1953) on dispersion in steady pipe flows. A significant milestone came when Aris (Reference Aris1960) applied his method of concentration moments (Aris Reference Aris1956) to analyse the period-averaged asymptotic dispersion of solute in oscillatory pipe flows. A key finding from previous studies (Bowden Reference Bowden1965; Holley, Harleman & Fischer Reference Holley, Harleman and Fischer1970; Fischer et al. Reference Fischer, List, Koh, Imberger and Brooks1979) is that in a purely oscillatory flow, the effective period-averaged dispersivity is significantly lower than in a steady flow of the same amplitude, due to concentration contraction when the flow reverses. Building on Aris’s (Reference Aris1956) method of moments, Brenner & Edwards (Reference Brenner and Edwards1993) systematically developed generalised Taylor dispersion (GTD) theory – a versatile framework capable of addressing both time- and space-dependent transport processes. The core formulation of GTD involves the so-called a posteriori trial forms of the first- and second-order concentration moments, which, however, rely on intuitive judgement.
In contrast to the aforementioned studies on period-averaged asymptotic dispersion, Chatwin (Reference Chatwin1975) focused on the variation within the oscillation period. Through statistical analysis, he inferred that as the solute tracers ultimately sample the entire cross-section, their time-dependent mean speed equals the instantaneous cross-sectional average of the streamwise advection velocity. Based on this inference, he further conjectured that the time-dependent dispersivity oscillates at twice the flow frequency, due to the multiplicative interaction of time-dependent terms in the definition of mean square displacement. Although Chatwin’s (Reference Chatwin1975) formal proof of these conclusions relied on the simplifying assumption that the concentration varies linearly in the streamwise direction, his inferences were subsequently validated by other researchers (Yasuda Reference Yasuda1984; Mukherjee & Mazumder Reference Mukherjee and Mazumder1988).
Researchers have identified two key parameters governing solute dispersion in oscillatory flows: the oscillation period and the cross-sectional diffusion time scale. Consider solute dispersion in an oscillatory flow driven by a time-periodic pressure gradient, with an initial uniform line source in the cross-section, released when the mean flow is zero. During the first half of the oscillation cycle, the solute patch spreads in one streamwise direction. If cross-sectional diffusion is relatively slow, then the patch largely returns to its original position during the second half of the cycle, producing only limited net dispersion in the streamwise direction. Therefore, transient negative dispersivity may arise in the latter half of this expansion–contraction process, even though the period-averaged dispersivity remains positive. This phenomenon motivated Smith (Reference Smith1982) to develop a more robust delayed-diffusion model to address the singularities caused by transient negative dispersivity. Yasuda (Reference Yasuda1984) later clarified the seemingly paradoxical transient negative dispersivity by proposing an alternative procedure for its calculation. He compared the commonly used approach – averaging across the cross-section before calculating the mean square displacement – with the newly proposed method, where the mean square displacement is calculated prior to averaging. While the first method reproduces the previously observed transient negative dispersivity, the second method, which he regarded as more reasonable, consistently yields a positive transient dispersivity.
In light of the aforementioned research on dispersion in oscillatory flows, most studies assume simplified cross-sectional transport processes characterised by a constant and steady diffusivity. Under this setting, the cross-sectional distribution eventually becomes uniform, i.e. the zeroth-order longitudinal concentration moment (i.e. the marginal distribution of solute across the cross-section) approaches a uniform state, significantly simplifying the transient and asymptotic calculations. For instance, based on an initial condition of uniform line release in the cross-section, Ding et al. (Reference Ding, Hunt, McLaughlin and Woodie2021) derived explicit expressions for the concentration moments up to third order. While this setting holds for solute dispersion in laminar oscillatory flows, it is not suitable for either active particles with non-diffusive migration (Fung, Bearon & Hwang Reference Fung, Bearon and Hwang2022), or turbulent conditions with non-uniform cross-sectional diffusivity (Bowden Reference Bowden1965).
Beyond the special case of the ultimate uniform zeroth-order streamwise concentration moment in the cross-section, the GTD theory for time-periodic processes, proposed by Brenner & Edwards (Reference Brenner and Edwards1993, ch. 6), addresses the dispersion of Brownian particles subjected to a time-periodic external force in the cross-section (Haber, Brenner & Shapira Reference Haber, Brenner and Shapira1990; Shapiro & Brenner Reference Shapiro and Brenner1990a , Reference Shapiro and Brennerb ). In these circumstances, the zeroth-order streamwise concentration moment becomes time-periodic as well. However, although the GTD theory for time-periodic processes of Brenner & Edwards (Reference Brenner and Edwards1993) has solved the time-dependent zeroth-order streamwise concentration moment, which adequately reflects the asymptotic transport mechanisms in the cross-section, it does not extend to how the drift and dispersivity behave asymptotically within the oscillation period, as only period-averaged solutions are obtained.
In this study, we consider the dispersion of active particles in oscillatory channel flows between parallel plates driven by a time-periodic pressure gradient. The active particles propel at a constant speed, referred to as swimmers hereafter, mimicking micro-organisms such as motile micro-algae and bacteria. The swimming directions are governed by the Jeffery equation (Jeffery Reference Jeffery1922), with a possible modification from gyrotaxis, typically induced by bottom-heaviness (Pedley & Kessler Reference Pedley and Kessler1992). The original Jeffery equation describes how the angular velocity of a spheroid depends on flow shear and rate-of-strain. For elongated swimmers suspended in a steady laminar pressure-driven flow, the steady concentration distribution shows a non-trivial response to the mean shear rate (Vennamneni, Nambiar & Subramanian Reference Vennamneni, Nambiar and Subramanian2020), as reported for both strongly elongated bacteria cells (Rusconi, Guasto & Stocker Reference Rusconi, Guasto and Stocker2014; Bearon & Hazel Reference Bearon and Hazel2015) and less elongated algal cells (Barry et al. Reference Barry, Rusconi, Guasto and Stocker2015). On the other hand, gyrotactic swimmers also possess a shear-sensitive feature. They are efficiently guided by the flow to accumulate near regions with the fastest downwelling speed (Kessler Reference Kessler1985a , Reference Kesslerb ), which can lead to a modification of the flow profile and even hydrodynamic instabilities in relatively concentrated suspensions (Kessler Reference Kessler1986; Hwang & Pedley Reference Hwang and Pedley2014b ; Bees Reference Bees2020; Fung, Bearon & Hwang Reference Fung, Bearon and Hwang2020; Fung & Hwang Reference Fung and Hwang2020; Ishikawa, Dang & Lauga Reference Ishikawa, Dang and Lauga2022; Fung Reference Fung2023; Wang, Jiang & Chen Reference Wang, Jiang and Chen2023; Ishikawa, Brumley & Pedley Reference Ishikawa, Brumley and Pedley2025). Thus both the elongated shape and gyrotaxis introduce a dependence of the swimmers’ cross-sectional migration on the local time-periodic shear. As a result, the existing transient solutions for the concentration moment equations (Mukherjee & Mazumder Reference Mukherjee and Mazumder1988) cannot be directly applied due to the inherent shear-dependent transport in the orientation and position space. Instead, a reformulation of the solutions for the moment equations, derived from the Smoluchowski equation for the transport of orientable swimmers in oscillatory flows, should be pursued to better understand the associated physics. It is important to note that several orientational-averaged models have been developed to characterise swimmer transport in position space, including the Pedley–Kessler model (Pedley & Kessler Reference Pedley and Kessler1990), the two-step GTD model (Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003), and a more recent new model by Fung et al. (Reference Fung, Bearon and Hwang2022). However, these models are generally considered inapplicable when the flow shear varies rapidly, such as in high-frequency oscillatory flows (Caldag & Bees Reference Caldag and Bees2025). Therefore, these models are not employed here.
For the long-time asymptotic dispersion regime, we revisit the classical GTD theory by Brenner & Edwards (Reference Brenner and Edwards1993), with an extension to the asymptotic oscillatory behaviours of the drift and dispersivity during an oscillation period. Such an extension is particularly useful for investigating the dispersion mechanism at a specific phase of the oscillation. For the transient dispersion processes, we adopt the method of moments proposed by Aris (Reference Aris1956). It is important to note that the introduction of self-propulsion results in the non-self-adjointness of the eigenvalue problem, which poses challenges for the solution technique of Barton (Reference Barton1983). This challenge is partially addressed by employing the bi-orthogonal expansion method (Strand, Kim & Karrila Reference Strand, Kim and Karrila1987; Nambiar et al. Reference Nambiar, Phanikanth, Nott and Subramanian2019; Jiang & Chen Reference Jiang and Chen2021). However, the time-dependent eigenvalue problem becomes significantly more complex to solve, especially when an additional orientation space is involved. Even for solutes, related works addressing the time-dependent eigenvalue problem, or equally complex Green’s function problem, have only derived the concentration moments up to second order (Yasuda Reference Yasuda1984; Mukherjee & Mazumder Reference Mukherjee and Mazumder1988; Wu et al. Reference Wu, Zeng, Chen, Li, Shao, Wang and Jiang2012). Therefore, while it is theoretically possible to solve for the transient dispersion of orientable swimmers in oscillatory flows, the dispersivity still presents significant technical challenges, let alone the calculation of higher-order statistics.
Recently, Jiang & Chen (Reference Jiang and Chen2025) applied a two-time-variable expansion to the concentration moment equations before solving them with the eigenfunction expansion method (Barton Reference Barton1983), which is particularly efficient for higher-order concentration moments. They introduced an auxiliary oscillation time variable that characterises the inherent oscillation in the dispersion due to the continuous expansion–contraction processes of the concentration. The concept underlying this new method actually shares a fundamental principle with the GTD theory for asymptotic dispersion – both feature an inherent periodic process governing the dispersion, while other transient modes decay gradually. The present work extends the two-time-variable expansion method by Jiang & Chen (Reference Jiang and Chen2025) to analyse the transient dispersion of swimmers in oscillatory flows, while also advancing the GTD theory to predict long-time asymptotic periodic behaviour of drift and dispersivity, rather than only their period-averaged values.
The remainder of this paper is organised as follows. In § 2, we describe the problem set-up for swimmer dispersion in oscillatory flows between parallel plates. In § 3, we develop the general solution procedure for the transient moment equations. In § 4, we revisit the GTD theory for long-time asymptotic dispersion, and extend it to account for phase-resolved drift and dispersivity. In § 5, we present and discuss the results. Finally, in § 6, we provide concluding remarks.
2. Problem set-up
2.1. Governing equation
As shown in figure 1, we consider the dispersion of a patch of swimmers released into an oscillatory flow between two vertical, parallel plates separated at distance
$W^{\ast }$
. Assuming a dilute suspension of swimmers, and thereby neglecting the forces exerted by the swimmer, the flow is driven by a pressure gradient composed of a steady part (which includes the contribution from gravitational acceleration) and a zero-mean oscillatory part:
\begin{equation} G_p^{\ast }(z^{\ast }, t^{\ast }) = -P_0^{\ast }\, \boldsymbol{e}_x - Q_0^{\ast } \cos \omega ^{\ast } t^{\ast }\, \boldsymbol{e}_x. \end{equation}
The velocity profile is given by (Von Kerczek Reference Von Kerczek1982)
\begin{gather} U^{\ast }(z^{\ast },t^{\ast }) = \frac {P_0^{\ast } (W^{\ast }-z^{\ast }) z^{\ast }}{2\nu ^{\ast }} + \frac {Q_0^{\ast }}{\omega ^{\ast }}\, \textrm{Im}\!\left \{ \left [\! 1- \frac { \cosh \left [ {(1+\textrm{i}) (2z^{\ast }-W^{\ast })}/{(2 \delta ^{\ast })} \right ]}{\cosh \left [ {(1+\textrm{i}) W^{\ast }}/({2\delta ^{\ast }})\right ] } \right ]\! \textrm{e}^{\textrm{i} \omega ^{\ast } t^{\ast }} \right \}\!, \end{gather}
where
$\nu ^{\ast }$
is the kinematic viscosity of water,
$\omega ^{\ast }$
is the angular frequency of oscillation, and
$\delta ^{\ast } = \sqrt {2\nu ^{\ast }/\omega ^{\ast }}$
denotes the thickness of the Stokes layer.
Figure 1. Schematic illustration of swimmer dispersion in a vertical oscillatory channel flow. (
$a$
) A gyrotactic swimmer in the oscillatory channel flow experiences a viscous torque and a gravitational torque, in addition to rotational diffusion. (
$b$
) Time evolution of oscillatory velocity profile for the case with
$P_0^{\ast } = 0$
,
$Q_0^{\ast } = 8 \nu ^{\ast } D_r^{\ast }/W^{\ast }$
,
$\delta ^{\ast } = 0.86 W^{\ast }$
and
$\omega ^{\ast } = D_r^{\ast }$
.
The swimmer propels itself at a constant speed
$V_s^{\ast }$
along an unsteady direction described by a unit vector
$\boldsymbol{p}$
. We note that the swimmers may exhibit more complex behaviours if they are allowed to rotate out of the
$x$
–
$z$
plane, even when the flow velocity gradient exists solely in the
$z$
-direction. An example is the resonate alignment of helical swimmers in oscillatory channel flows, as observed in Hope et al. (Reference Hope, Croze, Poon, Bees and Haw2016). As a first step in investigating oscillatory active dispersion, we therefore restrict our attention to particles whose rotation is constrained to the
$x$
–
$z$
plane, for mathematical convenience. In this case, the orientation of the swimmer is characterised by the angle
$\theta$
, defined as the anticlockwise angle between
$\boldsymbol{e}_x$
and
$\boldsymbol{p}$
, such that
$\boldsymbol{p} = \cos \theta\, \boldsymbol{e}_x + \sin \theta\, \boldsymbol{e}_z$
.
Apart from the random rotational diffusion, the deterministic rate of change of the swimmer’s direction is influenced by the vorticity, rate-of-strain and gravitational reorientation, as described by the modified Jeffery’s orbit (Jeffery Reference Jeffery1922; Pedley & Kessler Reference Pedley and Kessler1992):
\begin{equation} \dot {\boldsymbol p}^{\ast } = \frac {1}{2} \boldsymbol{\omega }^{\ast } \times \boldsymbol{p} + \alpha _0 \boldsymbol{p} \boldsymbol{\cdot } \unicode{x1D640}^{\ast } \boldsymbol{\cdot } (\unicode{x1D644}- \boldsymbol{pp}) + \frac {1}{2B^{\ast }} \boldsymbol{p}\times \boldsymbol{k} \times \boldsymbol{p}. \end{equation}
Here,
\begin{equation} \boldsymbol{\omega} ^{\ast } = \frac {\partial U^{\ast }}{\partial z^{\ast }} \boldsymbol{e}_z \times \boldsymbol{e}_x \end{equation}
denotes the vorticity,
$\alpha _0$
is the Bretherton parameter (with
$\alpha _0=0$
for a sphere, and
$\alpha _0=1$
for a slender rod),
\begin{equation} \unicode{x1D640}^{\ast } = \frac {1}{2} \frac {\partial U^{\ast }}{\partial z^{\ast }}\left ( \boldsymbol{e}_x \boldsymbol{e}_z + \boldsymbol{e}_z \boldsymbol{e}_z \right ) \end{equation}
is the rate-of-strain tensor,
$\unicode{x1D644}$
is the identity tensor,
$B^{\ast }$
is the gravitational reorientation time scale, and
$\boldsymbol{k} = -\boldsymbol{e}_x$
is the unit vector directed opposite to gravity.
The governing equation of the probability density function (p.d.f.)
$P$
of the swimmers takes the form of the Smoluchowski equation, with the dimensional form written as
\begin{gather} \frac {\partial P}{\partial t^{\ast }}+ \left (\frac {\partial }{\partial x^{\ast }} \boldsymbol{e}_x + \frac {\partial }{\partial z^{\ast }} \boldsymbol{e}_z \right ) \boldsymbol{\cdot } \boldsymbol{J}_p^{\ast }+ \left (\frac {\partial }{\partial \theta } \boldsymbol{e}_{\theta } \right ) \boldsymbol{\cdot } \boldsymbol{j}_r^{\ast }=0, \end{gather}
where
\begin{equation} \boldsymbol{J}_p^{\ast } = [V_s^{\ast } \boldsymbol{p} + U^{\ast }(z^{\ast },t^{\ast }) \boldsymbol{e}_x] P- D_t^{\ast } \left (\frac {\partial P}{\partial x^{\ast }} \boldsymbol{e}_x + \frac {\partial P}{\partial z^{\ast }} \boldsymbol{e}_z \right ) \end{equation}
and
\begin{equation} \boldsymbol{j}_r^{\ast } = \dot {\boldsymbol p}^{\ast } P-D_r^{\ast } \frac {\partial P}{\partial \theta } \boldsymbol{e}_{\theta } \end{equation}
are the fluxes in the position and orientation spaces, respectively. Here,
$D_r^{\ast }$
and
$D_t^{\ast }$
represent the rotational and translational diffusivity, respectively.
We non-dimensionalise the problem by
\begin{align} \left . \begin{array}{c} \displaystyle t=t^{\ast }D_r^{\ast },\quad \omega = \frac {\omega ^{\ast }}{D_r^{\ast }}, \quad z = \frac {z^{\ast }}{W^{\ast }}, \quad x=\frac {x^{\ast }}{W^{\ast }},\quad \textit {Wo} = \frac {W^{\ast }}{\delta ^{\ast }}, \quad U=\frac {U^{\ast }}{D_r^{\ast } W^{\ast }},\\[12pt] \displaystyle {\textit{Pe}}_s=\frac {V_s^{\ast }}{D_r^{\ast }W^{\ast }},\quad {\textit{Pe}}_{\!f}^{s} = \frac {U_s^{\ast }}{D_r^{\ast }W^{\ast }},\quad {\textit{Pe}}_{\!f}^{o} = \frac {U_o^{\ast }}{D_r^{\ast }W^{\ast }}, \quad \lambda = \frac {1}{2B^{\ast } D_r^{\ast }} ,\quad D_t=\frac { D_t^{\ast }}{D_r^{\ast }{W^{\ast }}^2}. \end{array}\!\! \right \} \end{align}
Here, the characteristic velocities
$U_s^{\ast } = {W^{\ast }}^2 P_0^{\ast }/(8\nu ^{\ast })$
and
$U_o^{\ast } = {W^{\ast }}^2 Q_0^{\ast }/(8\nu ^{\ast })$
are defined based on the steady and oscillatory pressure gradients, respectively. In (2.8),
$t$
denotes the dimensionless time, rescaled with the rotational diffusion time scale;
$\omega$
is the angular frequency of oscillation non-dimensionalised by the rotational diffusivity;
$z$
and
$x$
are the cross-sectional and streamwise coordinates, respectively, rescaled with the channel width;
$\textit {Wo}$
is the Womersley number;
$U$
is the dimensionless flow velocity, normalised by the characteristic velocity required to cross the channel width over one rotational diffusion time scale;
${\textit{Pe}}_s$
is the dimensionless swimming Péclet number;
${\textit{Pe}}_{\!f}^s$
and
${\textit{Pe}}_{\!f}^o$
correspond to the steady and oscillatory flow Péclet numbers, respectively;
$D_t$
is the dimensionless translational diffusivity rescaled with the characteristic diffusivity associated with diffusing across the channel width within one rotational diffusion time scale; and
$\lambda$
is the dimensionless gravitactic bias parameter.
The dimensionless Smoluchowski equation is given by
\begin{gather} \frac {\partial P}{\partial t} + (U + {\textit{Pe}}_s \cos \theta ) \frac {\partial P}{\partial x} - D_t \frac {\partial ^2 P}{\partial x^2} + {\textit{Pe}}_s \sin \theta\, \frac {\partial P}{\partial z} - D_t \frac {\partial ^2 P}{\partial z^2} + \frac {\partial (\dot {\theta } P)}{\partial \theta } - \frac {\partial ^2 P}{\partial \theta ^2} = 0, \end{gather}
where the dimensionless rate of change of
$\theta$
is
\begin{gather} \dot {\theta }(z,\theta ,t) = \frac {1}{2} \frac {\partial U}{\partial z} (-1 + \alpha _0 \cos 2 \theta ) + \lambda \sin \theta , \end{gather}
and the dimensionless flow velocity is
\begin{gather} U(z,t) = 4{\textit{Pe}}_{\!f}^{s}\,(1-z)z + \frac {4\, {\textit{Pe}}_{\!f}^{o}}{\textit {Wo}^2}\, \textrm{Im}\left \{ \left [ 1- \frac {\cosh \left [\textit {Wo}\,(1+\textrm{i}) (2z-1)/2\right ]}{\cosh [\textit {Wo}\, (1+\textrm{i})/2]} \right ] \textrm{e}^{\textrm{i} \omega t} \right \}. \end{gather}
We further define
$T = 2\unicode{x03C0} /\omega$
as the dimensionless oscillation period.
2.2. Periodic, boundary and initial conditions
The periodic conditions in the orientation variable
$\theta$
are naturally satisfied:
\begin{align} P|_{\theta =0} &= P|_{\theta =2\unicode{x03C0} }, \end{align}
\begin{align} \frac {\partial P}{\partial \theta }|_{\theta =0}&= \frac {\partial P}{\partial \theta }|_{\theta =2\unicode{x03C0} }. \end{align}
The interaction between swimmer and wall is inherently complex (Maretvadakethope et al. Reference Maretvadakethope, Hazel, Vasiev and Bearon2023; Zeng et al. Reference Zeng, Jiang, Guan, Lee and Chen2025), involving potential wall-induced modifications to swimming speed, angular velocity and rotational diffusivity (Kantsler et al. Reference Kantsler, Dunkel, Polin and Goldstein2013; Zeng, Jiang & Pedley Reference Zeng, Jiang and Pedley2022). However, since the primary aim of this work is to elucidate the dispersion mechanism governed predominantly by swimming and oscillatory shear, we adopt one of the most idealised and widely used boundary conditions in the continuum models – the reflective boundary conditions:
\begin{align} P|_{\theta =\theta _0} &= P|_{\theta =2\unicode{x03C0} -\theta _0}\quad \mbox{on} \ z=0, 1, \end{align}
\begin{align} \frac {\partial P}{\partial z}\bigg|_{\theta =\theta _0}&= - \frac {\partial P}{\partial z}\bigg|_{\theta =2\unicode{x03C0} -\theta _0}\quad \mbox{on} \ z=0, 1. \end{align}
The reflective boundary conditions ensure the no-flux condition at the channel walls,
\begin{equation} \int _0^{2\unicode{x03C0} } \left ( {\textit{Pe}}_s \sin \theta\, \frac {\partial P}{\partial z} - D_t \frac {\partial ^2 P}{\partial z^2} \right ) \textrm{d} \theta = 0\quad \mbox{on} \ z=0, 1, \end{equation}
thereby guaranteeing the conservation of the total number of swimmers.
The initial condition for
$P(x,z,\theta ,t)$
is formally prescribed as a source located in
$x=0$
:
\begin{equation} P|_{t=0} = I_{ini}(z, \theta )\, \delta (x), \end{equation}
where
$I_{ini}(z,\theta )$
represents the initial distribution in the cross-section-orientation space and satisfies the normalisation condition
\begin{equation} \int _0^1 \int _0^{2 \unicode{x03C0} } I_{ini}(z, \theta ) \, \textrm{d}\theta\, \textrm{d} z = 1. \end{equation}
3. Two-time-variable expansion for the transient moments
3.1. Smoluchowski equation in the two-time-variable system
Following Jiang & Chen (Reference Jiang and Chen2025), we introduce a two-time-variable expansion into the problem – the basic time variable and an auxiliary oscillatory time variable:
\begin{equation} t_0 \triangleq t, \quad t_1 \triangleq \omega t. \end{equation}
The auxiliary oscillatory time variable
$t_1$
is introduced to characterise the intrinsic periodic oscillation of the flow, and is applied exclusively to the time-periodic oscillatory advection term, thus the advection velocity becomes
\begin{equation} U(z,t_1) =4\,{\textit{Pe}}_{\!f}^{s}\,(1-z)z + \frac {4\, {\textit{Pe}}_{\!f}^{o}}{\textit {Wo}^2}\, \textrm{Im}\left \{ \left [ 1- \frac {\cosh \left [\textit {Wo}\,(1+\textrm{i}) (2z-1)/2 \right ]}{\cosh [\textit {Wo}\, (1+\textrm{i})/2]} \right ] \textrm{e}^{\textrm{i} t_1} \right \}. \end{equation}
As the original time variable
$t$
splits into two new time variables,
$t_0$
and
$t_1$
, the p.d.f. must be redefined accordingly to reflect its dependence on both time scales:
\begin{equation} \hat {P}(x,z,\theta ,t_0,t_1) \triangleq P(x,z,\theta ,t=t_0), \end{equation}
where the hat notation is used to distinguish the new two-time-variable p.d.f. from the original p.d.f.
Since
$t_1$
is introduced to characterise the time-periodic behaviour of the system, we impose periodicity of
$\hat {P}$
in
$t_1$
as
\begin{equation} \hat {P}(x,z,\theta ,t_0,t_1) = \hat {P}(x,z,\theta ,t_0,t_1+2\unicode{x03C0} ). \end{equation}
Accordingly, we map
$t_1$
to the interval
$[0,2\unicode{x03C0} )$
via
\begin{equation} t_1 = t_1 \bmod 2\unicode{x03C0} . \end{equation}
The use of two time variables offers two main advantages. First, it circumvents the complexity associated with solving the time-dependent eigenvalue problem, such as in Mukherjee & Mazumder (Reference Mukherjee and Mazumder1988). Second, it clarifies the underlying physical interpretation of the dispersion problem:
$t_0$
represents the time scale associated with the evolving dispersion from the initial release, while
$t_1$
captures the time-periodic dispersion process due to the flow oscillation. As will be further discussed in § 4, this approach aligns naturally with the GTD theory, wherein
$t_1$
emerges as the sole relevant time variable governing the asymptotic dispersion regime.
Under the two-time-variable expansion, the time derivative transforms as
\begin{equation} \frac {\partial }{\partial t} \longrightarrow \frac {\partial }{\partial t_0} + \omega \frac {\partial }{\partial t_1}. \end{equation}
It is worth noting that splitting of the time derivative has also been employed in several studies investigating the emergent dynamics of single swimmers with periodically varying shape and/or speed (Gaffney et al. Reference Gaffney, Dalwadi, Moreau, Ishimoto and Walker2022; Walker et al. Reference Walker, Ishimoto, Gaffney, Moreau and Dalwadi2022a
,
Reference Walker, Ishimoto, Moreau, Gaffney and Dalwadib
, Reference Walker, Ishimoto and Gaffney2023; Dalwadi et al. Reference Dalwadi, Moreau, Gaffney, Ishimoto and Walker2024a
,
Reference Dalwadi, Moreau, Gaffney, Walker and Ishimotob
). However, the methodologies diverge from this point on. In the multiple-time-scale approach, the periodic variation in shape or speed is assumed to be fast (
$\omega \gg 1$
), which naturally defines
$t_1$
as a fast time scale to facilitate perturbation analysis. In contrast, the two-time-variable method employed here does not assume any separation of time scale between
$t_0$
and
$t_1$
.
Another less directly related category of works focuses on deriving effective evolution equations for the concentration of swimmer populations subject to rotational and/or translational noises. Representative works include Bearon & Hazel (Reference Bearon and Hazel2015), Vennamneni et al. (Reference Vennamneni, Nambiar and Subramanian2020) and Fung et al. (Reference Fung, Bearon and Hwang2022), all of which also employ multiple-time-scale perturbation analysis. These studies typically introduce a slow time scale to derive effective transport equations tailored to different physical scenarios. However, none explicitly considers time-periodic transport processes. Instead, their scale hierarchy is based on the small swimming Péclet number, which compares the mean straight swimming length to the confinement scale.
Substituting (3.6) into the original transport equation (2.9) yields the two-time-variable formulation of the Smoluchowski equation:
\begin{gather} \frac {\partial \hat {P}}{\partial t_0} + \omega \frac {\partial \hat {P}}{\partial t_1}+ U_x \frac {\partial \hat {P}}{\partial x} - D_t \frac {\partial ^2 \hat {P}}{\partial x^2} + {\textit{Pe}}_s \sin \theta\, \frac {\partial \hat {P}}{\partial z} - D_t \frac {\partial ^2 \hat {P}}{\partial z^2} + \frac {\partial (\dot {\theta } \hat {P})}{\partial \theta } - \frac {\partial ^2 \hat {P}}{\partial \theta ^2} = 0, \end{gather}
where
\begin{equation} U_x(z,\theta ,t_1) = U(z,t_1) + {\textit{Pe}}_s \cos \theta \end{equation}
is the streamwise translational velocity, comprising both oscillatory advection and autonomous swimming.
It is important to note that by implementing the two-time-variable expansion and imposing periodic conditions on the oscillatory time variable
$t_1$
, this variable can effectively be regarded and treated as a pseudo-spatial variable defined on the interval
$[0,2\unicode{x03C0} )$
, associated with a convection term
$\omega (\partial ({\cdot })/\partial t_1)$
. The absence of a diffusion term in (3.7) for
$t_1$
actually is a direct consequence of intrinsic relationship between
$t_1$
and
$t_0$
: specifically,
$t_1$
evolves proportionally with
$t_0$
at a rate governed by
$\omega$
, which manifests solely through the convective term in the transport equation. We also note that the definition of
$t_0$
relates bidirectionally with the original time variable
$t$
, whereas
$t_1$
relates unidirectionally from
$t$
. Consequently, specifying
$t_0$
uniquely determines the original single time variable
$t$
, while specifying
$t_1$
alone does not.
In summary, the application of the two-time-variable expansion appears to increase the dimensionality of the Smoluchowski equation by 1. However, this transformation facilitates the subsequent solution of moment equations. Equation (3.7) also underscores a key distinction between passive and active dispersion: the motility of the swimmers causes the streamwise translation velocity and cross-sectional migration velocity to depend not on only the position but also on the orientation, which is simultaneously influenced by the background shear.
3.2. Moment equations in the two-time-variable system
The
$n$
th-order local p.d.f. moment is defined as
\begin{equation} \hat {P}_n(z,\theta ,t_0,t_1) \triangleq \int _{-\infty }^{\infty } x^n \hat {P} \,\textrm{d}x. \end{equation}
The governing equations for the first three local p.d.f. moments,
$\hat {P}_0$
,
$\hat {P}_1$
and
$\hat {P}_2$
, can be derived from (3.7), under the assumptions that
$\hat {P} \to 0$
and
$\partial \hat {P}/\partial x \to 0$
as
$|x| \to \infty$
:
\begin{align}&\qquad\,\,\, \frac {\partial \hat {P}_0}{\partial t_0} + \omega \frac {\partial \hat {P}_0}{\partial t_1} + \mathcal{L}_{co} \hat {P}_0 = 0, \end{align}
\begin{align}&\qquad \frac {\partial \hat {P}_1}{\partial t_0} + \omega \frac {\partial \hat {P}_1}{\partial t_1}+\mathcal{L}_{co} \hat {P}_1 = U_x \hat {P}_0, \end{align}
\begin{align}& \frac {\partial \hat {P}_2}{\partial t_0} + \omega \frac {\partial \hat {P}_2}{\partial t_1}+\mathcal{L}_{co} \hat {P}_2 = 2D_t\hat {P}_0 + 2U_x \hat {P}_1, \end{align}
where
$\mathcal{L}_{co}$
denotes the flux operator in the cross-section-orientation
$(z,\theta )$
space:
\begin{equation} \mathcal{L}_{co}({\cdot }) \triangleq {\textit{Pe}}_s \sin \theta\, \frac {\partial ({\cdot })}{\partial z} - D_t \frac {\partial ^2 ({\cdot })}{\partial z^2} + \frac {\partial ({\dot {\theta } ({\cdot })})}{\partial \theta } - \frac {\partial ^2 ({\cdot })}{\partial \theta ^2}. \end{equation}
The periodic conditions satisfied by
$P$
in
$\theta$
-space are inherited by
$\hat {P}_n$
:
\begin{align} \hat {P}_n|_{\theta =0} &= \hat {P}_n|_{\theta =2\unicode{x03C0} }, \end{align}
\begin{align} \frac {\partial \hat {P}_n}{\partial \theta }|_{\theta =0}&= \frac {\partial \hat {P}_n}{\partial \theta }|_{\theta =2\unicode{x03C0} }. \end{align}
Similarly,
$\hat {P}_n$
satisfies the periodicity in
$t_1$
:
\begin{equation} \hat {P}_n|_{t_1=0} = \hat {P}_n|_{t_1=2\unicode{x03C0} }. \end{equation}
The reflective boundary conditions imposed on
$P$
in the
$(z,\theta )$
space also apply to
$\hat {P}_n$
:
\begin{align} \hat {P}_n|_{\theta =\theta _0} &= \hat {P}_n|_{\theta =2\unicode{x03C0} -\theta _0}\quad \mbox{on} \ z=0, 1, \end{align}
\begin{align} \frac {\partial \hat {P}_n}{\partial z}|_{\theta =\theta _0} &= - \frac {\partial \hat {P}_n}{\partial z}|_{\theta =2\unicode{x03C0} -\theta _0}\quad \mbox{on} \ z=0, 1. \end{align}
For the initial condition of
$\hat {P}$
, the following relation can be deduced based on (3.1) and (3.3):
\begin{equation} \hat {P}|_{t_0=0,t_1=0} = P|_{t=0}. \end{equation}
Thus we use
\begin{equation} \hat {P}|_{t_0\,=\,0} = I_{ini}(z,\theta )\, \delta (x) \end{equation}
to satisfy (3.15), with the normalisation condition
\begin{equation} \frac {1}{2\unicode{x03C0} } \int _{-\infty }^{\infty } \int _{0}^{1} \int _{0}^{2\unicode{x03C0} } \int _{0}^{2\unicode{x03C0} } \hat {P}|_{t_0\,=\,0} \,\textrm{d}t_1\, \textrm{d}\theta \, \textrm{d}z\,\textrm{d}x = 1. \end{equation}
The initial conditions for the moments are
\begin{align}& \hat {P}_0|_{t_0\,=\,0} = I_{ini}(z,\theta ), \end{align}
\begin{align}&\quad\,\, \hat {P}_1|_{t_0\,=\,0} = 0, \end{align}
\begin{align}&\quad\,\, \hat {P}_2|_{t_0\,=\,0} = 0. \end{align}
The total moments integrated over
$(z,\theta )$
are defined as
\begin{gather} \big\langle \hat {P}_n\big \rangle _{z,\theta }(t_0,t_1)= \int _0^1 \int _0^{2\unicode{x03C0} } \hat {P}_n \,\textrm{d}\theta \,\textrm{d}z, \end{gather}
with the corresponding governing equations
\begin{equation} \frac {\partial \big \langle \hat {P}_n\big \rangle _{z,\theta }}{\partial t_0} + \omega \frac {\partial \big\langle \hat {P}_n\big\rangle _{z,\theta }}{\partial t_1} = n(n-1) D_t \big\langle \hat {P}_{n-2} \big\rangle _{z,\theta } + n \big\langle U_x \hat {P}_{n-1} \big\rangle _{z,\theta }. \end{equation}
Note that we have implicitly used the relations
$\langle \mathcal{L}_{co} \hat {P}_n \rangle _{z,\theta } = 0$
, based on the no-flux condition at the boundary and the periodicity in
$\theta$
. It is noted that
$\langle \hat {P}_n\rangle _{z,\theta }$
can be interpreted as the total moments in the single-time-variable dispersion problem, which are also used to calculate the transient evolution of drift and dispersivity.
A further integration of (3.20) over
$t_1 \in [0,2\unicode{x03C0} ]$
, followed by division by
$2\unicode{x03C0}$
, yields
\begin{equation} \frac {\partial \hat {M}_n}{\partial t_0} = n(n-1) D_t \hat {M}_{n-2} + n\, \overline {\left \langle U_x \hat {P}_{n-1} \right \rangle _{z,\theta }}, \end{equation}
where
\begin{equation} \overline {({\cdot })} \triangleq \frac {1}{2\unicode{x03C0} } \int _0^{2\unicode{x03C0} } ({\cdot }) \,\textrm{d}t_1 \end{equation}
represents the mean operation over
$t_1$
, and
\begin{equation} \hat {M}_n = \overline {\left \langle \hat {P}_n \right \rangle _{z,\theta }} \end{equation}
can be interpreted as the total moments in the two-time-variable system viewing
$t_1$
as a periodic variable. Equation (3.21) describes how the moments evolve from a period-averaged perspective.
3.3. Eigenfunction expansions and bi-orthogonality relation
Up to this point, we can solve (3.10a
)–(3.10c
) successively, subject to the appropriate boundary and initial conditions. Next, we define the new effective flux operator
$\mathcal{L}_{coo}$
in the cross-section-orientation-oscillation space
$(z,\theta ,t_1)$
, incorporating the advection term in
$t_1$
:
\begin{equation} \mathcal{L}_{coo}({\cdot }) \triangleq \mathcal{L}_{co}({\cdot }) + \omega \frac {\partial ({\cdot })}{\partial t_1}. \end{equation}
To solve for the moment equations using the Barton (Reference Barton1983) eigenfunction expansion method, we seek to solve the eigenvalue problem
\begin{equation} \mathcal{L}_{coo} f_i = \lambda _i f_i, \end{equation}
where
$f_i$
are the eigenfunctions, and
$\lambda _i$
are the corresponding eigenvalues. Due to the orientable motility of swimmers,
$\mathcal{L}_{coo}$
is non-self-adjoint; therefore, we employ a Galerkin method to numerically solve (3.25). The main steps involved in this approach are as follows.
-
(i) Find the appropriate basis functions in the
$(z,\theta ,t_1)$
space that satisfy the reflective boundary conditions (2.13) in
$(z,\theta )$
space, and the periodicity in
$t_1$
space. These basis functions are denoted as
$\{e_i\}_{i=1}^{\infty }$
, with the normalisation condition(3.26)where
\begin{equation} \int _{0}^{1} \int _{0}^{2\unicode{x03C0} } \int _{0}^{2\unicode{x03C0} } e_i e_{\!j} \,\textrm{d}t_1\, \textrm{d}\theta \, \textrm{d}z = \delta _{ij}, \end{equation}
$\delta _{ij}$
is the Kronecker delta function. The basis functions satisfying the reflective boundary conditions, as given in Wang, Jiang & Chen (Reference Wang, Jiang and Chen2022a
) for steady channel flow, are modified here by multiplying them by a normalised Fourier series:(3.27)
\begin{equation} \frac {1}{\sqrt {2\unicode{x03C0} }},\quad \frac {\sin (n t_1)}{\sqrt {\unicode{x03C0} }},\quad \frac {\cos (n t_1)}{\sqrt {\unicode{x03C0} }}, \quad n = 1, 2, 3, \ldots . \end{equation}
-
(ii) Construct the inner product matrix
$\unicode{x1D63C}$
, whose elements are expressed as(3.28)
\begin{equation} A_{ij} = A(e_i,e_{\!j}) = \int _{0}^{1} \int _{0}^{2\unicode{x03C0} } \int _{0}^{2\unicode{x03C0} } e_i (\mathcal{L}_{coo} e_{\!j}) \,\textrm{d}t_1\, \textrm{d}\theta \, \textrm{d}z. \end{equation}
-
(iii) Solve the weak formulation of the eigenvalue problem (3.25) at a truncation degree:
(3.29)where
\begin{equation} \unicode{x1D63C} \boldsymbol{\phi }_i = \lambda _i \boldsymbol{\phi }_i, \end{equation}
$\boldsymbol{\phi }_i$
is the vector of the coefficients of
$f_i$
.
Due to non-self-adjointness of
$\mathcal{L}_{coo}$
, the dual eigenfunctions
$\{f_i^{\star }\}_{i=1}^{\infty }$
are introduced for eigenfunction expansion, which satisfy the bi-orthogonality relation with the eigenfunctions
$\{f_i\}_{i=1}^{\infty }$
(Jiang & Chen Reference Jiang and Chen2021). The bi-orthogonality relation is given by
\begin{equation} \int _{0}^{1} \int _{0}^{2\unicode{x03C0} } \int _{0}^{2\unicode{x03C0} } f_i f_{\!j}^{\star } \,\textrm{d}t_1\, \textrm{d}\theta \, \textrm{d}z = \delta _{ij}. \end{equation}
This condition ensures mutual orthogonality between the eigenfunctions and their dual counterparts. The dual eigenfunctions satisfy the following eigenvalue problem for the adjoint operator
$\mathcal{L}_{coo}^{\star }$
:
\begin{equation} \mathcal{L}_{coo}^{\star } f_i^{\star } = \lambda _i f_i^{\star }. \end{equation}
With the obtained eigenfunctions, dual eigenfunctions and eigenvalues, the local p.d.f. moments can be expanded into the series (Barton Reference Barton1983; Jiang & Chen Reference Jiang and Chen2021)
\begin{equation} \hat {P}_n(z,\theta ,t_1,t_0) = \sum _{i=1}^{\infty } p_{ni}\, \textrm{e}^{-\lambda _i t_0}\, f_i(z,\theta ,t_1),\quad n=0,1,\ldots , \end{equation}
where
$p_{ni}$
are the expansion coefficients for local p.d.f. moments of each order.
For the first step in solving for the zeroth local p.d.f. moment
$\hat {P}_0$
, we find the expansion coefficients
$p_{0i}$
using the initial condition (3.18a
):
\begin{align} p_{0i} = \int _{0}^{1} \int _{0}^{2\unicode{x03C0} } \int _{0}^{2\unicode{x03C0} } I_{ini}(z,\theta )\, f_i^{\star } \,\textrm{d}t_1\, \textrm{d}\theta \, \textrm{d}z, \end{align}
where
$N$
denotes the truncation degree. The higher-order local p.d.f. moments
$\hat {P}_1$
and
$\hat {P}_2$
can be successively solved with the series representations given in Jiang & Chen (Reference Jiang and Chen2021, App. A).
With the solved transient moments in the two-time-variable system, we return to the moments in the single-time-variable system using
\begin{equation} P_n(z,\theta ,t) = \hat {P}_n(z,\theta ,t_0=t,t_1=\omega t) \end{equation}
for straightforward definitions of the drift and dispersivity:
\begin{align}& U_d \triangleq \frac {\textrm{d} \left \langle P_1 \right \rangle _{z,\theta }}{\textrm{d} t}, \end{align}
\begin{align}&\,\, D_T \triangleq \frac {1}{2}\frac {\textrm{d} \sigma ^2}{\textrm{d} t}, \end{align}
where
\begin{equation} \sigma ^2 = \left \langle P_2 \right \rangle _{z,\theta } - \left \langle P_1 \right \rangle _{z,\theta }^2 \end{equation}
is the mean square displacement of the cross-section-averaged concentration.
We validate our transient moments solutions by comparing with Brownian dynamics (BD) simulations, as described in Appendix A. Figure 2 presents a typical case of gyrotactic swimmers released as a uniform line source. As shown, the results of method of moments and BD simulations are in good agreement.
Figure 2. Comparison of the results obtained from moments equations and BD simulations over the first two oscillation periods. (
$a$
) First-order total moment
$\langle P_1 \rangle _{z,\theta }$
. (
$b$
) Mean square displacement of the cross-section-averaged concentration
$\sigma ^2$
. Note that quantities obtained with moments equations are expressed with the single original time variable
$t$
using the substitutions
$t_0 \longrightarrow t$
and
$t_1 \longrightarrow \omega t$
, and the hat symbols are simultaneously removed. Parameters are
${\textit{Pe}}_s=0.1$
,
${\textit{Pe}}_{\!f}^s=0$
,
${\textit{Pe}}_{\!f}^o=1$
,
$\alpha _0=0$
,
$\lambda =2.19$
,
$\omega =1$
,
$\textit {Wo}=1.72$
,
$I_{ini} =1/(2\unicode{x03C0} )$
.
In the method of moments, the components of the basis functions in each variable are truncated at the following degrees:
$50$
for
$z$
,
$12$
for
$\theta$
, and
$4$
for
$t_1$
. Additionally,
$500$
pairs of eigenvalues and eigenfunctions are retained. These truncation degrees are determined through independence tests, comparing the moments and cross-sectional concentration.
4. The GTD theory for the long-time asymptotic periodic dispersion
According to GTD theory for time-periodic dispersion problems (Brenner & Edwards Reference Brenner and Edwards1993), the long-time asymptotic solution of zeroth-order local moment
$P_0$
, denoted by
$P_0^{\infty }$
, becomes periodic over the oscillation cycle. Consequently, it is natural to seek the governing equation for the redefined zeroth-order local moment
$\hat {P}_0^{\infty }$
within the two-time-variable formulation by neglecting the time derivative with respect to
$t_0$
, while retaining the derivative with respect to
$t_1$
:
\begin{equation} \mathcal{L}_{coo} \hat {P}_0^{\infty } = 0. \end{equation}
Equation (4.1) can be solved using a Galerkin method with the same basis functions as described in § 3:
\begin{equation} \hat {P}_0^{\infty }(z,\theta ,t_1) = \sum _{i=1}^{\infty } a_i\, e_i(z,\theta ,t_1), \end{equation}
which is equivalent to the expression in (3.32) for
$n=0$
in the limit
$t_0 \to \infty$
. We expect
$\hat {P}_0^{\infty }$
to exhibit oscillatory behaviour for gyrotactic and elongated swimmers under the reflection boundary conditions, since their cross-sectional migration is predominantly influenced by the oscillatory shear.
4.1. Long-time periodic drift
The long-time asymptotic drift velocity can be directly evaluated as
\begin{equation} U_d^{\infty }(t_1) = \big\langle \hat {P}_0^{\infty } U_x \big\rangle _{z,\theta }, \end{equation}
where the averaging is performed over the
$(z,\theta )$
phase space only. As a result,
$U_d^{\infty }$
remains time-periodic. To characterise the net transport over one oscillation cycle, we define the period-averaged asymptotic drift velocity as
\begin{equation} \overline {U_d^{\infty }} = \frac {1}{2\unicode{x03C0} } \int _0^{2\unicode{x03C0} } U_d^{\infty } \,\textrm{d}t_1. \end{equation}
4.2. Long-time period-averaged dispersivity
We assume that the solution for the first-order total p.d.f. moment, further averaged over the oscillatory time variable
$t_1$
, takes the form
\begin{equation} \hat {M}_1(t_0) \sim \overline {U_d^{\infty }} t_0 + \overline {\left \langle b \right \rangle _{z,\theta }}+ \text{e.d.t.} \end{equation}
Equation (4.5) can be interpreted as follows: from a period-averaged perspective, the centroid of the swimmer distribution moves with the period-averaged asymptotic drift velocity
$\overline {U_d^{\infty }}$
, offset by a constant term
$\overline {\left \langle b \right \rangle _{z,\theta }}$
, and accompanied by an exponential decay term (e.d.t.) in
$t_0$
.
Based on (4.5) and the relationship between
$\hat {P}_1$
and
$\hat {M}_1$
given in (3.23), we assume that the solution for the first-order local p.d.f. moment takes the form
\begin{equation} \hat {P}_1(z,\theta ,t_1,t_0) \sim \hat {P}_0^{\infty }\, \overline {U_d^{\infty }} t_0 + b(z,\theta ,t_1) + \text{e.d.t.}, \end{equation}
where we define the scalar field
\begin{equation} B(z,\theta ,t_1) = \frac {b}{\hat {P}_0^{\infty }} = \lim _{t_0 \to \infty } \left ( \frac {\hat {P}_1}{\hat {P}_0^{\infty }} - \hat {M}_1 \right ) + \overline {\left \langle b \right \rangle _{z,\theta }}. \end{equation}
Inspection of (4.7) reveals that
$B(z,\theta ,t_1)$
can be interpreted as the deviation between the mean streamwise position of swimmers located at
$(z,\theta ,t_1)$
and the overall mean streamwise position of all swimmers, with the addition of a constant
$\overline {\left \langle b \right \rangle _{z,\theta }}$
. The scalar field
$B(z,\theta ,t_1)$
, which is often called the ‘Brenner field’, thus encapsulates all relevant information regarding the dispersion mechanisms (Brenner & Edwards Reference Brenner and Edwards1993; Haugerud, Linga & Flekkøy Reference Haugerud, Linga and Flekkøy2022; Wang et al. Reference Wang, Jiang, Zeng and Chen2025).
Substituting (4.6) into the governing equation for
$\hat {P}_1$
, (3.10b
), we obtain the equation governing
$b(z,\theta ,t_1)$
:
\begin{equation} \mathcal{L}_{coo} b = \hat {P}_0^{\infty } \left (U_x - \overline {U_d^{\infty }}\right ). \end{equation}
Rather than solving directly for
$b$
, we introduce a normalised Brenner field
$b_N \triangleq b - \overline {\left \langle b \right \rangle _{z,\theta }} \hat {P}_0^{\infty }$
, which satisfies
\begin{equation} \overline {\left \langle b_N \right \rangle _{z,\theta }} = 0. \end{equation}
The governing equation for
$b_N$
takes the same form as that for
$b$
:
\begin{equation} \mathcal{L}_{coo} b_N = \hat {P}_0^{\infty } \left (U_x - \overline {U_d^{\infty }}\right ). \end{equation}
The period-averaged dispersivity is defined as
\begin{equation} \overline {D_T^{\infty }} \triangleq \lim _{t_0 \to \infty } \frac {1}{2} \frac {\partial }{\partial t_0} \left (\hat {M}_2 - \hat {M}_1^2\right ) . \end{equation}
This expression can be simplified using (3.21), yielding
\begin{equation} \overline {D_T^{\infty }} = D_t + \lim _{t_0 \to \infty } \left (\overline {\left \langle U_x \hat {P}_1 \right \rangle _{z,\theta }} - \hat {M}_1 \frac {\partial \hat {M}_1}{\partial t_0} \right ). \end{equation}
Substituting (4.5) and (4.6) into (4.12), we arrive at the final expression for the period-averaged dispersivity:
\begin{equation} \overline {D_T^{\infty }} = D_t + \overline { \left \langle U_x b_N \right \rangle _{z,\theta } } . \end{equation}
The preceding derivation can be viewed as an extension of the framework presented by Brenner & Edwards (Reference Brenner and Edwards1993, ch. 6), which addresses scenarios without coupling between fluxes in the global and local spaces. In contrast, our formulation incorporates the coupling between the global flux in the streamwise direction (
$x$
) and the local flux in the cross-section-orientation space
$(z,\theta )$
, arising from the swimmers’ self-propulsion. In the absence of such coupling, the analysis simplifies considerably, as
$\hat {P}_0^{\infty }$
becomes independent of
$t_1$
.
The long-time asymptotic period-averaged drift and dispersivity predicted by the GTD theory can be alternatively derived using a two-time-scale homogenisation method, as detailed in Appendix B. However, neither approach has, to the best of our knowledge, been extended to characterise the long-time asymptotic periodic (phase-dependent) dispersivity.
4.3. Long-time periodic dispersivity
To proceed further, we integrate (4.6) over the
$(z,\theta )$
phase space to obtain
\begin{equation} \big\langle \hat {P}_1 \big\rangle _{z,\theta } \sim \big\langle \hat {P}_0^{\infty } \big\rangle _{z,\theta } \overline {U_d^{\infty }} t_0 + \left \langle b \right \rangle _{z,\theta } + \text{e.d.t.} \end{equation}
Thus the long-time asymptotic periodic dispersivity can be calculated using (3.20) and (4.14):
\begin{align} D_T^{\infty }(t_1) &\triangleq \frac {1}{2} \lim _{t_0 \to \infty }\frac {\textrm{d}}{\textrm{d} t}\left (\left \langle P_2 \right \rangle _{z,\theta } - \left \langle P_1 \right \rangle _{z,\theta }^2\right ) \notag \\ & = \frac {1}{2} \left ( \frac {\partial \left \langle \hat {P}_2 \right \rangle _{z,\theta }}{\partial t_0} + \omega \frac {\partial \left \langle \hat {P}_2 \right \rangle _{z,\theta }}{\partial t_1}\right ) -\left \langle \hat {P}_1 \right \rangle _{z,\theta } \left \langle U_x \hat {P}_0^{\infty } \right \rangle _{z,\theta } \notag \\ & = D_t + \left \langle \left (U_x - U_d^{\infty }\right ) b_N \right \rangle _{z,\theta }, \end{align}
where we have used the conservation condition in the
$(z,\theta)$
space:
\begin{equation} \left \langle \hat {P}_0^{\infty }\right \rangle _{z,\theta } = 1. \end{equation}
The above equation is deduced as follows. First, we have
\begin{equation} \frac {\partial \hat {P}_0^{\infty }}{\partial t_0} = 0. \end{equation}
Integrating (4.17) over
$(z,\theta )$
phase space yields
\begin{equation} \frac {\partial \left \langle P_0^{\infty } \right \rangle _{z,\theta }}{\partial t_0} = 0. \end{equation}
Substituting (4.18) into (3.20) with
$n=0$
results in
\begin{equation} \frac {\partial \left \langle P_0^{\infty } \right \rangle _{z,\theta }}{\partial t_1} = 0. \end{equation}
Thus
$\left \langle P_0^{\infty } \right \rangle _{z,\theta }$
must be a constant. Applying the normalisation condition (3.17), we arrive at (4.16).
We validate our GTD solutions for the periodic drift and dispersivity by comparing them with BD simulations and method of moments. Figure 3 presents the same case as in figure 2, but with a time far from the initial release for the method of moments and BD simulations, as the GTD solutions are intended for long times. Once again, we find good agreement between the results from GTD theory, the method of moments, and BD simulations.
Figure 3. Comparison of the results obtained from the moments equations, GTD and BD simulations over an oscillation period long after the initial release (
$t \in [14T,15T]$
): (
$a$
) drift
$U_d$
, (
$b$
) dispersivity
$D_T$
. Note that quantities obtained with moment equations and GTD are expressed with the single original time variable
$t$
using the substitutions
$t_0 \longrightarrow t$
and
$t_1 \longrightarrow \omega t$
. The parameters are consistent with those used in figure 2.
5. Results and discussion
We use a practical parameter space, as listed in table 1, referencing the model organisms Chlamydomonas augustae with relaxation on the shape (Bretherton parameter
$\alpha _0$
) to encompass a broader range of micro-organisms, such as chain-forming micro-algae (Lovecchio et al. Reference Lovecchio, Climent, Stocker and Durham2019) and bacteria (Ran & Arratia Reference Ran and Arratia2024). To clarify how oscillatory shear modulates the dispersion of swimmers by influencing their local rotation, we also present the results for the solute. The solute diffusivity is set equal to the steady effective diffusivity of two-dimensional non-gyrotactic swimmers in an unbounded quiescent fluid, given by
${\textit{Pe}}_s^2/(n^2-n)+ D_t$
(Cates & Tailleur Reference Cates and Tailleur2013), where
$n=2$
denotes the number of spatial dimensions. For simplicity, we adopt a uniform line release:
Table 1. Parameters used for swimmers in this work. The values of
$V_s^{\ast }$
,
$B^{\ast }$
and
$D_r^{\ast }$
are based on the model organisms Chlamydomonas augustae (data primarily sourced from Pedley & Kessler (Reference Pedley and Kessler1990) and Hwang & Pedley (Reference Hwang and Pedley2014a
,Reference Hwang and Pedley
b
)). Note that since constant values of
$W^{\ast }$
,
$\nu ^{\ast }$
and
$D_r^{\ast }$
are used, the Womersley number
$\textit {Wo}$
is uniquely determined by the relation
$\textit {Wo} = 1.72 \sqrt {\omega }$
.
\begin{equation} I_{ini} = \frac {1}{2\unicode{x03C0} }. \end{equation}
In the following, we present all results using the original single-time-variable system. Quantities derived using the method of moments, which formally depend on both
$t_0$
and
$t_1$
, are now expressed with the original single time variable
$t$
. Quantities obtained from the GTD theory, which periodically depend on
$t_1$
, are presented as functions of the original time variable modulo the oscillation period, i.e.
$t \bmod T$
.
5.1. Effects of swimming ability and oscillatory flow strengths
In this subsection, we investigate the effects of swimming ability,
${\textit{Pe}}_s$
, and oscillatory flow strengths,
${\textit{Pe}}_{\!f}^o$
, on both the transient and long-time asymptotic dispersion. We focus on spherical non-gyrotactic swimmers (SNS), where the distribution in
$(z,\theta )$
space remains uniform for the uniform line release (5.1) under reflective boundary conditions, i.e.
$\hat {P}_0(z,\theta ,t_1) = 1/(2\unicode{x03C0} )$
(Jiang & Chen Reference Jiang and Chen2021). As a result, the transient drift
$U_d$
and the long-time asymptotic periodic drift
$U_d^{\infty }$
are identical for solute and SNS, perfectly tracking the instantaneous mean flow speed, as seen in figures 4(
$a$
) and 4(
$c$
).
Figure 4. (
$a$
,
$b$
) Transient drift
$U_d$
and dispersivity
$D_T$
of solute and SNS over the first three periods following a uniform line release for several oscillatory flow Péclet numbers
${\textit{Pe}}_{\!f}^o$
. (
$c$
,
$d$
) Long-time asymptotic periodic drift
$U_d^{\infty }$
and dispersivity
$D_T^{\infty }$
of solute and SNS over one period for several oscillatory flow Péclet numbers
${\textit{Pe}}_{\!f}^o$
. Parameters for flow:
${\textit{Pe}}_{\!f}^s=0$
,
$\omega =1$
,
$\textit {Wo}=1.72$
. Parameters for solute:
${\textit{Pe}}_s=0$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0.005$
. Parameters for SNS:
${\textit{Pe}}_s=0.1$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0$
.
As shear disrupts the cross-sectional migration of swimmers in an oscillatory manner, their dispersivity behaves differently from the solute, which undergoes simple molecular diffusion across the cross-section. Figures 4(
$b$
) and 4(
$d$
) show the subtle differences in the time evolution of dispersivity of SNS compared to the solute. Initially, following release,
$D_T$
of the solute starts at its molecular diffusivity
$D_t=0.005$
, whereas
$D_T$
of SNS begins at zero and gradually approaches that of the solute. This can be attributed to the fact that the mean square displacement of SNS exhibits Brownian diffusive behaviour with diffusivity
$D_t$
(set to zero in our work) at very short time scales in the absence of flow (Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016). Beyond this initial stage, the difference in
$D_T$
between solute and SNS becomes most pronounced at
${\textit{Pe}}_{\!f}^o=1$
, where the flow dominates swimming (
${\textit{Pe}}_{\!f}^o \gg {\textit{Pe}}_s$
); in other cases, the oscillatory behaviour of
$D_T$
remains relatively weak. The most notable differences appear at the peaks and troughs of the oscillatory
$D_T$
– i.e.
$D_T$
of SNS exceeds that of the solute at the peaks, but falls below it at the troughs. This phenomenon can be interpreted as follows. When the flow develops strongly in one direction, SNS are dispersed more effectively due to suppressed cross-streamline migration under Jeffery orbits, which is reminiscent of the inverse dependence of Taylor dispersivity on molecular diffusivity in steady flows. However, when the flow reverses direction, SNS tend to return more closely to their initial position compared to the solute, as their ability to diffuse across the cross-section is limited. For both solute and SNS, the emergence of negative
$D_T$
arises from our use of the widely adopted definition of
$\sigma ^2$
, which involves averaging across the cross-section before computing the mean square displacement, rather than the alternative definition proposed by Yasuda (Reference Yasuda1984), which calculates mean square displacement at each streamline first.
At long times,
$D_T^{\infty }$
becomes periodic for both solute and SNS. Furthermore, figure 4(
$d$
) reveals that the actual period of
$D_T^{\infty }$
is
$T/2$
, due to the two symmetrical dispersion cycles that occur when the mean flow velocity points in opposite directions, which is in line with the preliminary conjecture of Chatwin (Reference Chatwin1975) and the more formal results given by Yasuda (Reference Yasuda1984). The two dispersion cycles are identical in the variations of
$D_T^{\infty }$
, which is independent of phase of the initial release. Each cycle can be further separated into two intervals: one with positive
$D_T^{\infty }$
, corresponding to the attenuation of mean flow velocity in one direction, and one with negative
$D_T^{\infty }$
, corresponding to the strengthening of mean flow velocity in the opposite direction. Moreover, the integration of
$D_T^{\infty }$
over time during the positive interval outweighs the integration during the negative interval, leading to a positive period-averaged dispersivity.
The long-time asymptotic period-averaged dispersivity
$\overline {D_T^{\infty }}$
is plotted in figure 5 as functions of
${\textit{Pe}}_s$
and
${\textit{Pe}}_{\!f}^o$
. Note that the period-averaged drift
$\overline {U_d^{\infty }}$
is always zero for both solute and SNS due to the uniform distribution in the
$(z,\theta )$
space, so it is not shown here for conciseness. Increasing trends of
$\overline {D_T^{\infty }}$
with
${\textit{Pe}}_s$
and
${\textit{Pe}}_{\!f}^o$
are observed, indicating the positive influences of swimming ability and oscillatory flow strength on the dispersion at the period-averaged level over asymptotically long times. Furthermore, in figure 5(
$b$
), we observe that as
${\textit{Pe}}_{\!f}^o$
increases,
$\overline {D_T^{\infty }}$
of the solute exceeds that of SNS, indicating a net weakening effect of motility on dispersion in oscillatory flow. However, we speculate that this trend may not hold across a broader range of oscillation frequency
$\omega$
.
Figure 5. Long-time asymptotic period-averaged dispersivity
$\overline {D_T^{\infty }}$
as functions of (
$a$
) the swimming Péclet number
${\textit{Pe}}_s$
and (
$b$
) the oscillatory flow Péclet number
${\textit{Pe}}_{\!f}^o$
. Parameters for flow:
${\textit{Pe}}_{\!f}^s=0$
,
$\omega =1$
,
$\textit {Wo}=1.72$
. Parameters in (
$a$
):
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0$
. Parameters in (
$b$
):
${\textit{Pe}}_s=0$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0.005$
for solute,
${\textit{Pe}}_s=0.1$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0$
for SNS.
5.2. Effects of oscillation frequency
In this subsection, we examine the previous conjecture regarding the influence of the oscillation frequency
$\omega$
on active dispersion in oscillatory flows. Figure 6 presents the transient and long-time asymptotic periodic drift and dispersivity for several values of
$\omega$
, with comparisons between solute and SNS. While both
$U_d$
and
$U_d^{\infty }$
continue to follow the instantaneous mean flow velocity, their magnitudes decrease and their phases shift rightwards as
$\omega$
increases. This behaviour reflects the increasing influence of fluid inertia, characterised by the Womersley number
$\textit {Wo}$
.
Figure 6. (
$a$
,
$b$
) Transient drift
$U_d$
and dispersivity
$D_T$
of solute and SNS over the first three periods following a uniform line release for several oscillation frequencies
$\omega$
. (
$c$
,
$d$
) Long-time asymptotic periodic drift
$U_d^{\infty }$
and dispersivity
$D_T^{\infty }$
of solute and SNS over one period for several oscillation frequencies
$\omega$
. Parameters for solute:
${\textit{Pe}}_s=0$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0.005$
. Parameters for SNS:
${\textit{Pe}}_s=0.1$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0$
.
The amplitude of dispersivity oscillations diminishes more noticeably with increasing
$\omega$
, due to the quadratic dependence of dispersivity on the flow velocity. At the lowest oscillation frequency considered (
$\omega =0.1$
), the transient dispersivity of SNS is significantly larger than that of the solute, except at the troughs, where the two values nearly converge. At
$\omega =1$
, the dispersivity curves of SNS and solute nearly overlap, with only minor differences at the peaks and troughs. At the highest oscillation frequency (
$\omega =10$
), the dispersivity of SNS is consistently smaller than that of the solute. The complex influence of oscillation frequency is illustrated more clearly from a long-time asymptotic period-averaged perspective. As shown in figure 7, although the dispersivity
$\overline {D_T^{\infty }}$
for both SNS and solute declines rapidly with increasing
$\omega$
,
$\overline {D_T^{\infty }}$
for SNS exceeds that of the solute at relatively low oscillation frequencies
$\omega \lt 0.3$
, but becomes smaller at higher oscillation frequencies
$\omega \gt 0.3$
. Since
$\omega$
is defined as the ratio of dimensional oscillation frequency to the rotational diffusivity, this observation underscores the intricate coupling between swimming dynamics and oscillatory shear in determining the dispersion of SNS.
Figure 7. Long-time asymptotic period-averaged dispersivity
$\overline {D_T^{\infty }}$
as a function of oscillation frequency
$\omega$
. The parameters for solute and SNS are consistent with those used in figure 6.
5.3. Effects of superimposed steady component on dispersion in oscillatory flows
This subsection investigates the effects of a superimposed steady component on dispersion in oscillatory flows. Figure 8 presents the transient and long-time asymptotic variations of drift and dispersivity for several steady flow Péclet numbers
${\textit{Pe}}_{\!f}^s$
, with a fixed oscillatory flow Péclet number
${\textit{Pe}}_{\!f}^o=1$
. The combination of oscillatory flow and steady flow does not alter the uniform distribution in the
$(z,\theta )$
space; therefore, the drift remains equivalent to the instantaneous mean flow velocity, contributed by both the steady and oscillatory components. In terms of dispersivity, the presence of a steady flow leads to an increased oscillation amplitude – specifically, the peaks rise more than the troughs fall. Additionally, only a single dispersion cycle occurs over an asymptotic oscillation period, unlike the two cycles observed previously when
${\textit{Pe}}_{\!f}^s = 0$
.
Figure 8. (
$a$
,
$b$
) Transient drift
$U_d$
and dispersivity
$D_T$
of solute and SNS over the first three periods following a uniform line release for several steady flow Péclet numbers
${\textit{Pe}}_{\!f}^s$
. (
$c$
,
$d$
) Long-time asymptotic periodic drift
$U_d^{\infty }$
and dispersivity
$D_T^{\infty }$
of solute and SNS over one period for several steady flow Péclet numbers
${\textit{Pe}}_{\!f}^s$
. Parameters for flow:
${\textit{Pe}}_{\!f}^o=1$
,
$\omega =1$
,
$\textit {Wo}=1.72$
. Parameters for solute:
${\textit{Pe}}_s=0$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0.005$
. Parameters for SNS:
${\textit{Pe}}_s=0.1$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0$
.
The introduction of a steady component also results in a more complex variation in the dispersivity of SNS. As illustrated in figure 9, there is a dual effect of an oscillatory component on the dispersion of SNS when a steady component is present: at small
${\textit{Pe}}_{\!f}^s$
, oscillation enhances dispersion, whereas at large
${\textit{Pe}}_{\!f}^s$
, oscillation inhibits it. In contrast, for the solute, the addition of oscillation consistently enhances dispersion compared with purely steady flow.
Figure 9. Long-time asymptotic period-averaged dispersivity
$\overline {D_T^{\infty }}$
plotted against the steady flow Péclet number
${\textit{Pe}}_{\!f}^s$
. Parameters for flow:
$\omega =1$
,
$\textit {Wo}=1.72$
. The parameters for solute and SNS are consistent with those used in figure 8.
5.4. Effects of gyrotaxis and elongation
This subsection discusses the effects of gyrotaxis and elongation on oscillatory dispersion. Figure 10 compares the transient and long-time asymptotic periodic drift and dispersivity for solute, SNS, spherical gyrotactic swimmers (SGS), and elongated non-gyrotactic swimmers (ENS). While ENS slightly always drift slower than the instantaneous mean flow, SGS exhibit more pronounced variations – drifting faster during upward flow, and slower during downward flow. This behaviour is not attributed to the well-known response of SGS to shear flow: in steady conditions, SGS typically undergo gyrotactic focusing near the centre or walls of the channel for downwelling or upwelling flows, resulting in concentration accumulations significantly exceeding the mean, and causing enhanced drift during downward flow, and reduced drift during upward flow. However, in oscillatory flows, the cross-sectional concentration distribution of SGS cannot immediately follow the time-varying shear profile due to the constraints of weak swimming strength (
${\textit{Pe}}_s=0.1$
) and limited gyrotactic reorientation. As shown in figure 11(
$a$
), the cross-sectional concentration profile of SGS remains relatively flat in the central region over a period, with notable gradients only occurring near the walls – starkly contrasting with the steady-flow case. Figure 11(
$c$
) further presents the long-time asymptotic periodic local mean swimming direction component along the cross-section
$\left \langle p_z^{\infty } \right \rangle _{\theta }$
of SGS over a period. Significant gradients in
$\left \langle p_z^{\infty } \right \rangle _{\theta }$
are again confined to the near-wall region, suggesting corresponding localised concentration variations.
Figure 10. (
$a$
,
$b$
) Transient drift
$U_d$
and dispersivity
$D_T$
of solute, SNS, SGS and ENS over the first three periods following a uniform line release. (
$c$
,
$d$
) Long-time asymptotic periodic drift
$U_d^{\infty }$
and dispersivity
$D_T^{\infty }$
of solute, SNS, SGS and ENS over one period. Parameters for solute:
${\textit{Pe}}_s=0$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0.005$
. Parameters for SNS:
${\textit{Pe}}_s=0.1$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0$
. Parameters for SGS:
${\textit{Pe}}_s=0.1$
,
$\alpha _0=0$
,
$\lambda =2.19$
,
$D_t=0$
. Parameters for ENS:
${\textit{Pe}}_s=0.1$
,
$\alpha _0=1$
,
$\lambda =0$
,
$D_t=0$
.
Figure 11. Long-time asymptotic periodic cross-sectional distribution
$\langle P_0^{\infty } \rangle _{\theta }\triangleq \int _0^{2\unicode{x03C0} } P_0^{\infty } \, \textrm{d} \theta$
, and local mean swimming direction components along the cross-section
$\langle p_z^{\infty } \rangle _{\theta } \triangleq (\int _0^{2\unicode{x03C0} } P_0^{\infty } \sin \theta \, \textrm{d} \theta )/ (\int _0^{2\unicode{x03C0} } P_0^{\infty } \, \textrm{d} \theta )$
and along the streamwise direction
$\langle p_x^{\infty } \rangle _{\theta } \triangleq (\int _0^{2\unicode{x03C0} } P_0^{\infty } \cos \theta \, \textrm{d} \theta )/(\int _0^{2\unicode{x03C0} } P_0^{\infty } \, \textrm{d} \theta )$
over one period: (
$a$
,
$c$
,
$e$
) SGS, (
$b$
,
$d$
,
$f$
) ENS. Flow parameters:
${\textit{Pe}}_{\!f}^s=0$
,
${\textit{Pe}}_{\!f}^o=1$
,
$\omega =1$
,
$\textit {Wo}=1.72$
. Particle parameters are consistent with those used in figure 10.
A full sampling of the non-uniform shear field is crucial to the shear-induced trapping of ENS, which typically requires a long time and streamwise distances (Rusconi et al. Reference Rusconi, Guasto and Stocker2014; Vennamneni et al. Reference Vennamneni, Nambiar and Subramanian2020). Moreover, the concentration enhancement in steady flow remains of the same order as the mean concentration. As shown in figure 11(
$b$
), the cross-sectional concentration profiles of ENS are even more uniform, with only weak gradients near the walls. Figures 11(
$d$
) and 11(
$f$
) present the long-time asymptotic periodic local mean cross-sectional and streamwise swimming direction components over a period,
$\left \langle p_z^{\infty } \right \rangle _{\theta }$
and
$\left \langle p_x^{\infty } \right \rangle _{\theta }$
, respectively, further confirming the weak accumulation and alignment of ENS in oscillatory flows.
The stronger response of SGS in terms of concentration distribution and mean swimming direction to oscillatory flow leads to its distinct dispersivity compared with other types of particle. In contrast, the relatively uniform and weak responses of ENS result in dispersivity characteristics more similar to those of solute and SNS, as shown in figures 10(
$b$
) and 10(
$d$
).
Figure 12 presents the long-time asymptotic period-averaged drift
$\overline {U_d^{\infty }}$
and dispersivity
$\overline {D_T^{\infty }}$
as functions of the gravitactic bias parameter
$\lambda$
and Bretherton parameter
$\alpha _0$
. Gyrotaxis is found to induce a negative mean drift velocity, which can be explained by the steady upward streamwise alignment observed in figure 11(
$e$
). It also weakens the overall dispersivity. Elongation, on the other hand, does not produce a net drift for non-gyrotactic swimmers, but slightly enhances
$\overline {U_d^{\infty }}$
for gyrotactic swimmers. For both swimmer types, elongation leads to a modest increase in
$\overline {D_T^{\infty }}$
.
Figure 12. Long-time asymptotic period-averaged drift
$\overline {U_d^{\infty }}$
and dispersivity
$\overline {D_T^{\infty }}$
as functions of (
$a$
,
$c$
) gravitactic bias parameter
$\lambda$
, and (
$b$
,
$d$
) Bretherton parameter
$\alpha _0$
. Flow parameters are consistent with those used in figure 11.
6. Concluding remarks
This work combines a two-time-variable expansion for the transient dispersion and the GTD theory for the long-time asymptotic periodic dispersion to investigate the Taylor–Aris dispersion of active particles in oscillatory channel flows. The two-time-variable expansion, which is introduced to capture the periodicity in transient evolution to the Taylor regime due to oscillation of the flow, gains deeper interpretation through the lens of GTD theory: in the long-time asymptotic limit, the dispersion problem simplifies to a periodic problem governed solely by the oscillatory time variable
$t_1=\omega t$
.
Traditional approximating models based on orientation–position separation, such as the two-step GTD model (Bearon, Hazel & Thorn Reference Bearon, Hazel and Thorn2011), typically assume a quasi-steady and quasi-uniform shear. These models solve the equilibrium orientation pointwise, then compute the drift velocity and diffusivity tensor in position space. However, such assumptions break down in flows with strong spatial inhomogeneity (Bearon et al. Reference Bearon, Hazel and Thorn2011; Jiang & Chen Reference Jiang and Chen2020; Wang et al. Reference Wang, Jiang and Chen2022a
). Recently, Caldag & Bees (Reference Caldag and Bees2025) also demonstrated the failure of the two-step GTD model at high oscillation frequencies, where flow conditions deviate significantly from quasi-steady assumption. At the other end of the oscillation spectrum, while these models may provide accurate predictions for very low-frequency oscillations, they often require prohibitively long simulations and large computational domains to reach the asymptotic Taylor regime (Caldag & Bees Reference Caldag and Bees2025). In contrast, our transient dispersion framework efficiently captures key statistical features in terms of moments, without resolving the full concentration field explicitly. This makes it particularly suitable for low-frequency oscillatory flows as well. Therefore, compared to traditional approximating models, our method accommodates a broader range of oscillation frequencies
$\omega$
, offering a promising tool for understanding and controlling dispersion of active particles in oscillatory flows.
Employing both approaches, we conduct detailed analyses of swimmer dispersion following release from a uniform line source, considering potential additional effects from gyrotaxis and elongation. For spherical non-gyrotactic swimmers, although their zeroth-order moment remains uniform – similar to that of the solute – the presence of oscillatory flow can either enhance or reduce their dispersion relative to the solute with comparable molecular diffusivity. The enhancement, which occurs at low frequencies, is attributed to the swimmers’ reduced effective ability to migrate across streamlines, as described by the Jeffery equation, which captures the linear dependence of rotation on shear rate for spherical particles. In contrast, the reduction at high frequencies is due to the near resetting of swimmers’ streamwise positions after each short oscillation cycle, a consequence of limited cross-sectional homogenising. The case of superimposing a steady component onto an oscillatory flow is also investigated, revealing a dual effect of oscillation on swimmer dispersion, in contrast to the monotonic enhancement observed for solute. Two potential terms in the Jeffery orbit with bias to gravity and rate-of-strain – gyrotaxis and elongation – are subsequently considered. While both induce non-uniform cross-sectional distributions, only gyrotaxis significantly alters the dispersion characteristics, owing to its higher sensitivity to shear. By contrast, although elongation allows swimmers to respond to shear gradients, their trapping in high-shear regions requires extended time and distance. Consequently, in an oscillatory environment, where such sustained exposure is absent, shear-induced trapping fails to develop appreciably as in steady flows, thus has a limited effect on dispersion.
This work assumes an extremely dilute suspension, neglecting swimmer–swimmer interaction. In practice, for naturally buoyant swimmers (typically 5 % denser than water), buoyancy-driven effects dominate swimmer–swimmer interactions in steady gyrotactic plumes (Fung et al. Reference Fung, Bearon and Hwang2020). In oscillatory flows, the development of gyrotactic plumes is hindered, weakening both buoyancy effects and far-field interactions. If buoyancy–flow coupling becomes non-negligible, then local swimmer accumulation may alter density distributions, resulting in a more pronounced phase lag between flow velocity and pressure gradient. Additionally, considering far-field interactions for puller-type swimmers may increase the effective viscosity (Saintillan Reference Saintillan2018), further amplifying this phase lag.
Swimmer–boundary interactions represent another important factor that can influence active particle dispersion in oscillatory flows. The reflective boundary conditions employed in the current study do not capture wall accumulation effects, which may lead to an underestimation of near-wall concentrations. Previous studies have shown that wall-accumulating swimmers exhibit reduced dispersivity in moderate steady flows (Jiang & Chen Reference Jiang and Chen2019). In contrast, in strong steady flows, the dispersivity approaches that of non-accumulating swimmers due to the suppression of wall polarisation by shear. In oscillatory flows, we speculate that wall accumulation could be more pronounced than in steady flows, owing to the reduced shear rate near the wall. Such enhanced accumulation may further suppress dispersion.
The current work can be extended to other periodic processes beyond oscillatory flows. Recent studies (Omori et al. Reference Omori, Kikuchi, Schmitz, Pavlovic, Chuang and Ishikawa2022; Walker et al. Reference Wang, Jiang, Chen and Tao2022b ) suggest that the periodic variations in swimming speed and body shape may explain the experimentally observed rheotaxis and centreline migration of swimmers. Additionally, based on this work, a multiple-time-variable expansion can be employed to address the situations when multiple time-periodic processes coexist, each with distinct base frequencies. Further development could involve constructing a period-resolved transport model capturing the full regimes of concentration evolution. Although applying the long-time asymptotic period-averaged one-dimensional dispersion equation (B16) avoids the singularity associated with transient negative diffusivity, it resolves none of the cross-sectional concentration, the characteristics within the period, and the transient dynamics. A potential improvement would be to approximate the concentration distribution by fully utilising the transient concentration moments at each streamline, such as applying the Edgeworth expansion (Chatwin Reference Chatwin1970; Wang & Chen Reference Wang and Chen2017; Guo, Jiang & Chen Reference Guo, Jiang and Chen2020; Li et al. Reference Li, Gong, Jiang, Zhan, Wang, Fu, Xu and Wu2023; Guan & Chen Reference Guan and Chen2024) or Gill’s generalised dispersion model (Gill Reference Gill1967a , Reference Gillb ), both of which have been shown to be efficient in solute dispersion.
Acknowledgements
We are grateful to Professor G.Q. Chen for valuable discussions during the initial conceptualisation of this work.
Funding
This work is supported by the National Key R&D Program of China (grant no. 2024YFC3210902), the IWHR Research and Development Support Program (grant no. HY0199A112021), and the Independent Research Project of State Key Laboratory of Water Cycle and Water Security in River Basin (grant nos SKL2024TS12 and SKL2025KYQD02). B.W. acknowledges support from the Youth Talent Lifting Project of the Department of Hydraulics, IWHR (grant no. HY121003A0012025).
Declaration of interests
The authors report no conflict of interest.
Appendix A. The BD simulations
We use BD simulations to validate our solutions for the transient moment equations and the long-time asymptotic GTD theory. The simulations are performed using the Langevin equation for the single-particle motion:
\begin{align}& \textrm{d} x = U\, \textrm{d} t + {\textit{Pe}}_s \cos \theta \textrm{d} t + \sigma _x\, \textrm{d} W_x, \end{align}
\begin{align}&\quad\,\,\, \textrm{d} z = {\textit{Pe}}_s \sin \theta \textrm{d} t + \sigma _z\, \textrm{d} W_z, \end{align}
\begin{align}&\qquad\quad \textrm{d} \theta = \dot {\theta } \textrm{d} t + \sigma _r\, \textrm{d} W_{\theta }. \end{align}
Here
$W_x$
,
$W_z$
and
$W_{\theta }$
are independent standard Brownian motions,
$\sigma _x = \sigma _z = \sqrt {2 D_t}$
and
$\sigma _r = \sqrt {2}$
.
After discretising (A1a
)–(A1c
) using the Euler–Maruyama forward scheme, we simulate
$10^5$
trajectories of swimmers with time step
$\Delta t = 10^{-4} T$
. For swimmers outside the boundary
$z \in [0,1]$
, a reflection operation on both position and orientation is performed with respect to the adjacent wall, ensuring consistency with the reflective boundary conditions given in (2.13). Concentration moments and distributions are then extracted to validate the theoretical models.
Here, we perform more comprehensive validations of the moment equation solutions using BD simulations, covering a wide range of oscillation frequencies (
$\omega = 0.1, 1, 10$
) and various particle types. As shown in figure 13, the results from the BD simulations and the moment equations show good agreement, confirming that the validity of the two-time-variable expansion method is independent of the oscillation frequency.
Figure 13. Comparison of the results obtained from moments equations (ME) and BD simulations over the first two oscillation periods: (a,c,e) first-order total moment
$\langle P_1 \rangle _{z,\theta }$
; (b,d,f) mean square displacement of the cross-section-averaged concentration
$\sigma ^2$
; for (a,b)
$\omega =0.1$
, (c,d)
$\omega =1$
, (e,f)
$\omega =10$
. Note that quantities obtained with moments equations are expressed with the original single time variable
$t$
using the substitutions
$t_0 \longrightarrow t$
and
$t_1 \longrightarrow \omega t$
, and the hat symbols are simultaneously removed. Parameters for flow:
${\textit{Pe}}_{\!f}^s=0$
,
${\textit{Pe}}_{\!f}^o=1$
. Parameters for solute:
${\textit{Pe}}_s=0$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0.005$
. Parameters for SNS:
${\textit{Pe}}_s=0.1$
,
$\alpha _0=0$
,
$\lambda =0$
,
$D_t=0$
. Parameters for SGS:
${\textit{Pe}}_s=0.1$
,
$\alpha _0=0$
,
$\lambda =2.19$
,
$D_t=0$
. Parameters for ENS:
${\textit{Pe}}_s=0.1$
,
$\alpha _0=1$
,
$\lambda =0$
,
$D_t=0$
. The initial condition is
$I_{ini} =1/(2\unicode{x03C0} )$
.
Appendix B. Two-time-scale homogenisation for the long-time asymptotic period-averaged drift and dispersivity
In this appendix, we present a two-time-scale homogenisation method for determining the long-time asymptotic period-averaged drift and dispersivity. The single-time-scale homogenisation method (Pavliotis & Stuart Reference Pavliotis and Stuart2008, ch. 12) has been widely applied to dispersion in steady flows (Pavliotis Reference Pavliotis2005; Chen & Thiffeault Reference Chen and Thiffeault2021; Wang et al. Reference Wang, Jiang and Chen2022b ; Guan et al. Reference Guan, Jiang, Tao, Chen and Lee2024), and is equivalent to the small wavenumber expansion in the complex Fourier space (Peng & Brady Reference Peng and Brady2020; Peng Reference Peng2024). Here, due to the inherent unsteadiness of the forcing flow, a two-time-scale homogenisation is needed.
We first rewrite (2.9) in a co-moving streamwise coordinate,
$x_0 \triangleq x - \overline {U_d^{\infty }}\, t$
, where
$\overline {U_d^{\infty }}$
is the period-averaged drift:
\begin{equation} \frac {\partial P}{\partial t}+ \left (U_x - \overline {U_d^{\infty }}\right ) \frac {\partial P}{\partial x_0} - D_t \frac {\partial ^2 P}{\partial x_0^2} + \mathcal{L}_{co} P = 0. \end{equation}
Note that while we use the same symbol
$\overline {U_d^{\infty }}$
here, its value is yet to be determined from the perturbation equations, distinct from its previous use in § 4.
We then consider the two time scales
\begin{equation} \tau _0 \triangleq t, \quad \tau _1 \triangleq \epsilon ^2 t, \end{equation}
where
$\tau _0$
represents the fast time scale compared with the slow diffusive time scale
$\tau _1$
, since
$\epsilon$
is a small parameter. Equation (B1) is rewritten as
\begin{equation} \frac {\partial \tilde {P}}{\partial \tau _0} + \epsilon ^2 \frac {\partial \tilde {P}}{\partial \tau _1} + \left (U_x - \overline {U_d^{\infty }}\right ) \frac {\partial \tilde {P}}{\partial x} - D_t \frac {\partial ^2 \tilde {P}}{\partial x^2} + \mathcal{L}_{co} \tilde {P} = 0, \end{equation}
where
$\tilde {P}(x,z,\theta ,\tau _0,\tau _1) = P(x,z,\theta ,t)$
.
We next introduce the typical diffusive length scaling
\begin{equation} \xi = \epsilon x_0, \end{equation}
where
$\epsilon$
is a small parameter, which can be interpreted as
$W^{\ast }/L^{\ast }$
, with
$L^{\ast }$
denoting the characteristic streamwise length scale of the dispersing swimmer patch. Under this rescaling, the transport equation becomes
\begin{equation} \frac {\partial \tilde {P}}{\partial \tau _0} + \epsilon ^2 \frac {\partial \tilde {P}}{\partial \tau _1} + \epsilon \left (U_x - \overline {U_d^{\infty }}\right ) \frac {\partial \tilde {P}}{\partial \xi } - \epsilon ^2 D_t \frac {\partial ^2 \tilde {P}}{\partial \xi ^2} + \mathcal{L}_{co} \tilde {P} = 0. \end{equation}
We expand
$\tilde {P}$
as a regular perturbation series in powers of
$\epsilon$
:
\begin{equation} \tilde {P} = \tilde {P}_0 + \epsilon \tilde {P}_1 + \epsilon ^2 \tilde {P}_2 + O(\epsilon ^3). \end{equation}
To ensure solvability at each order, we impose the solvability conditions
\begin{align}& \int _0^1 \int _0^{2\unicode{x03C0} } \tilde {P}_1 \, \textrm{d}\theta \, \textrm{d}z=0, \end{align}
\begin{align}& \int _0^1 \int _0^{2\unicode{x03C0} } \tilde {P}_2 \, \textrm{d}\theta \, \textrm{d}z=0, \end{align}
together with periodicity conditions on
$P_1$
and
$P_2$
in
$\tau _0$
, each with period
$T$
.
The perturbation equations at successive orders of
$\epsilon$
are obtained by substituting (B6) into (B5) and collecting terms of the same order.
\begin{align}&\qquad\qquad\qquad\qquad O(1): \frac {\partial \tilde {P}_0}{\partial \tau _0} + \mathcal{L}_{co} \tilde {P}_0 = 0. \end{align}
\begin{align}&\qquad\qquad O(\epsilon ): \frac {\partial \tilde {P}_1}{\partial \tau _0} + \Big(U_x - \overline {U_d^{\infty }}\Big) \frac {\partial \tilde {P}_0}{\partial \xi } + \mathcal{L}_{co} \tilde {P}_1 = 0. \end{align}
\begin{align}& O\big(\epsilon ^2\big): \frac {\partial \tilde {P}_2}{\partial \tau _0} + \frac {\partial \tilde {P}_0}{\partial \tau _1} + \Big(U_x - \overline {U_d^{\infty }}\Big) \frac {\partial \tilde {P}_1}{\partial \xi } - D_t \frac {\partial ^2 \tilde {P}_0}{\partial \xi ^2} + \mathcal{L}_{co} \tilde {P}_2 = 0. \end{align}
The leading-order solution
$\tilde {P}_0$
is assumed to take the separable form
\begin{equation} \tilde {P}_0(\xi , z, \theta , \tau _0, \tau _1) = g_0(z, \theta , \tau _0)\, C(\xi , \tau _1), \end{equation}
where
$g_0$
represents the conditional probability density in the
$(z,\theta ,\tau _0)$
space, and
$C$
denotes the streamwise, period- and cross-section-averaged concentration. Substituting this ansatz into the
$O(1)$
perturbation problem yields the governing equation for
$g_0$
:
\begin{equation} \frac {\partial g_0}{\partial \tau _0} + \mathcal{L}_{co} g_0 = 0, \end{equation}
with periodic conditions
$g_0|_{\tau _0} = g_0|_{\tau _0+T}$
.
The first-order correction
$\tilde {P}_1$
is assumed to take the form
\begin{equation} \tilde {P}_1(\xi ,z,\theta ,\tau _0,\tau _1) = \chi (z,\theta ,\tau _0)\, \frac {\partial C}{\partial \xi } + g_0(z,\theta ,\tau _0)\, h(\xi ,\tau _1), \end{equation}
where
$\chi$
satisfies the same boundary and periodic conditions as
$g_0$
. Substituting this expression into the
$O(\epsilon )$
perturbation equation yields the cell problem for
$\chi$
:
\begin{equation} \frac {\partial \chi }{\partial \tau _0} + \mathcal{L}_{co} \chi = \left (\overline {U_d^{\infty }} - U_x\right ) g_0 . \end{equation}
Integrating the above equation over
$\tau _0 \in [0,T]$
,
$z\in [0,1]$
and
$\theta \in [0,2\unicode{x03C0} ]$
eliminates the left-hand side by boundary and periodic conditions, leading to the expression for the period-averaged drift:
\begin{equation} \overline {U_d^{\infty }} = \frac {1}{T}\int _{0}^{T} \int _0^1 \int _0^{2\unicode{x03C0} } U_x g_0\, \textrm{d} \theta \, \textrm{d} z \, \textrm{d} \tau _0. \end{equation}
Applying the solvability condition, i.e. integrating the
$O(\epsilon ^2)$
perturbation equation (B8c
) over
$\tau _0 \in [0, T]$
,
$z\in [0,1]$
and
$\theta \in [0,2\unicode{x03C0} ]$
, yields the effective one-dimensional dispersion equation
\begin{equation} \frac {\partial C}{\partial \tau _1} = \overline {D_T^{\infty }} \frac {\partial ^2 C}{\partial \xi ^2}, \end{equation}
where
$\overline {D_T^{\infty }}$
denotes the period-averaged effective dispersivity with the same symbol as in § 4. The expression for
$\overline {D_T^{\infty }}$
is given by
\begin{equation} \overline {D_T^{\infty }} = D_t + \frac {1}{T}\int _{0}^{T} \int _0^1 \int _0^{2\unicode{x03C0} } \left [ (\overline {U_d^{\infty }} - U_x) \chi \right ] \, \textrm{d} \theta \, \textrm{d} z \, \textrm{d} \tau _0. \end{equation}
Comparing the governing equations and resulting expressions (B10) and (B12)–(B14) with their counterparts in the GTD framework, specifically (4.1), (4.4), (4.10) and (4.13), it becomes evident that the two-time-scale homogenisation method is formally equivalent to the GTD theory in capturing the long-time asymptotic period-averaged drift and dispersivity. However, we note that the current two-time-scale homogenisation method cannot capture the temporal variations of drift and dispersivity within a period. Such variations, which are retained in our extended GTD formulation, are averaged out in the homogenisation framework.
It is also important to distinguish the current two-time-scale homogenisation from the two-time-variable expansion introduced in § 3. In the two-time-variable expansion,
$t_0$
and
$t_1$
do not necessarily have different orders of magnitude, no perturbation analysis is required, and the solutions for the transient moments equations are exact. In contrast, the two-time-scale homogenisation method aims to capture the long-time asymptotic period-averaged dispersion behaviour by introducing two time scales of different orders: the normal time scale
$\tau _0$
, and the asymptotic effective diffusive time scale
$\tau _1$
, given that
$\epsilon$
is a small parameter.
Converting back to the original coordinate, (B14) becomes the period-averaged one-dimensional dispersion equation
\begin{equation} \frac {\partial C}{\partial t} = \overline {U_d^{\infty }} \frac {\partial C}{\partial x} + \overline {D_T^{\infty }} \frac {\partial ^2 C}{\partial x^2}, \end{equation}
which is valid in the long-time asymptotic limit.
Open access
