# Emergent rheotaxis of shape-changing swimmers in Poiseuille flow

B.J. Walker, K. Ishimoto, Clément Moreau, E.A. Gaffney, M.P. Dalwadi

► To cite this version:

B.J. Walker, K. Ishimoto, Clément Moreau, E.A. Gaffney, M.P. Dalwadi. Emergent rheotaxis of shape-changing swimmers in Poiseuille flow. Journal of Fluid Mechanics, 2022, 944, pp.R2. 10.1017/jfm.2022.474. hal-05050468

HAL Id: hal-05050468

https://hal.science/hal-05050468v1

Submitted on 29 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
J M

RAPIDS

<figure><img src="image_2.png" /></figure>

Emergent rheotaxis of shape-changing swimmers
in Poiseuille flow

B.J. Walker<fnref n="1" />, †, K. Ishimoto<fnref n="2" />, C. Moreau<fnref n="2" />, E.A. Gaffney<fnref n="3" /> and M.P. Dalwadi<fnref n="1" />

(Received 16 March 2022; revised 26 April 2022; accepted 26 May 2022)

A simple model for the motion of shape-changing swimmers in Poiseuille flow was recently proposed and numerically explored by Omori *et al.* (*J. Fluid Mech.*, vol. 930, 2022, A30). These explorations hinted that a small number of interacting mechanics can drive long-time behaviours in this model, cast in the context of the well-studied alga *Chlamydomonas* and its rheotactic behaviours in such flows. Here, we explore this model analytically via a multiple-scale asymptotic analysis, seeking to formally identify the causal factors that shape the behaviour of these swimmers in Poiseuille flow. By capturing the evolution of a Hamiltonian-like quantity, we reveal the origins of the long-term drift in a single swimmer-dependent constant, whose sign determines the eventual behaviour of the swimmer. This constant captures the nonlinear interaction between the oscillatory speed and effective hydrodynamic shape of deforming swimmers, driving drift either towards or away from rheotaxis.

**Key words:** micro-organism dynamics

**1. Introduction**

The behaviours of microswimmers in flows have long been a topic of broad theoretical and experimental study. Recently, Omori *et al.* (2022) numerically explored a model of a shape-changing swimmer in Poiseuille flow, posed in the context of the alga *Chlamydomonas* for comparison with their experimental findings. Their investigations suggested an interesting and subtle connection between the long-time behaviours of the microswimmer and the details of its changing speed and shape, with certain conditions

† Email address for correspondence: benjamin.walker@ucl.ac.uk

<footnote n="1">¹Department of Mathematics, University College London, London WC1H 0AY, UK</footnote>

<footnote n="2">²Research Institute for Mathematical Sciences, Kyoto University, Kyoto 606-8502, Japan</footnote>

<footnote n="3">³Wolfson Centre for Mathematical Biology, Mathematical Institute, University of Oxford,
Oxford OX2 6GG, UK</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. Notation and set-up. We illustrate a model swimmer in Poiseuille flow, located at a transverse displacement y from the midline of the parabolic flow profile. The swimming direction θ is measured from the midline, with θ = 0 corresponding to downstream swimming.</figcaption></figure>

apparently necessary for long-time upstream-facing swimming in the flow, referred to as
rheotaxis by Omori *et al.* (2022).

Their ordinary differential equation (ODE) model may be simply stated in terms of a
transverse coordinate y and the swimmer orientation θ as

$$
\frac{dy}{dt} = \omega u(\omega t) \sin \theta, \quad (1.1a)
$$

$$
\frac{d\theta}{dt} = \gamma y(1 - B(\omega t) \cos 2\theta), \quad (1.1b)
$$

with given initial conditions and with reference to the set-up of figure 1, where all
quantities are considered dimensionless. We refer the interested reader to the work of
Omori *et al.* (2022) for a full derivation of the model. The functions *u* and *B* capture
the time-dependent active swimming speed and shape-capturing Bretherton constant
(Bretherton 1962), respectively. These prescribed functions are assumed to be oscillatory
with a shared high frequency *ω* ≫ 1. In particular, the swimming speed naturally scales
with *ω* in (1.1a), as the Stokes equations are both linear and independent of time except
via boundary conditions (Kim & Karrila 2005). However, such a velocity scaling is
absent from the explicitly stated equations of Omori *et al.* (2022), although it is present
in their explored parameter regimes. The equations of Stokes flow are appropriate in
the low-Reynolds-number and low-frequency-Reynolds-number regimes associated with
many microswimmers, including *Chlamydomonas* (Guasto, Rusconi & Stocker 2012),
and we will restrict our analysis to such regimes of practical interest. Here, *γ* is a fixed
characteristic shear rate of the flow, non-negative without loss of generality. This model
neglects any interactions of the swimmer with solid boundaries typically associated with
Poiseuille flow, and we will proceed without additional consideration of boundary effects.

Via the numerical explorations of Omori *et al.* (2022), this model is noted to give rise
to a range of complex, long-time behaviours, perhaps the most remarkable of which is
conditional convergence towards a central upstream-facing configuration. In this study,
we will aim to analytically uncover the driving factors behind these long-time dynamics.
Via a multiple-scale asymptotic analysis (Bender & Orszag 1999), as recently applied to
similar models of swimming (Gaffney *et al.* 2022; Ma, Pujara & Thiffeault 2022; Walker
*et al.* 2022), we will show how the effective swimmer behaviour can be captured by a
Hamiltonian-like quantity, whose slow evolution accurately encodes the long-time trends
of behaviour noted by Omori *et al.* (2022). Further, we will identify a markedly simple
relation between the eventual behaviour of the swimmer and its oscillating speed and
shape, enabling the deduction of long-time dynamics through the calculation of a single
swimmer-dependent constant.
## 2. Direct asymptotic analysis

The time scales present in the model of Omori et al. (2022) are best identified through a change of variable. Defining $z(t) := y(t)/\omega^{1/2}$, the system reads

$$ \frac{dz}{dt} = \omega^{1/2}u(\omega t) \sin \theta, \qquad (2.1a) $$

$$ \frac{d\theta}{dt} = \gamma\omega^{1/2}z(1 - B(\omega t) \cos 2\theta). \qquad (2.1b) $$

This suggests a natural fast time scale $T := \omega t$, that of the oscillating swimmer speed and shape. Additionally, $O(1)$ oscillations of $u$ and $B$ in (2.1) drive $O(1)$ changes in $z$ and $\theta$ over an intermediate time scale of $t = O(\omega^{-1/2})$. We will later see that these changes correspond to quasi-periodic orbits of the swimmer in the flow, quantifying motion on this time scale via $\tau := \omega^{1/2}t$. In addition to the two time scales $T$ and $\tau$ evident from this system of equations, Omori et al. (2022) observed behavioural changes over longer time scales, with $t = O(1)$. Hence, we expect the system to evolve on three separated time scales, corresponding to $T$, $\tau$ and $t$ each being $O(1)$. Our overarching aim is to characterise and understand the behaviours of the system over the time scale $t = O(1)$.

To this end, we conduct a multiple-scale analysis and formally write $z(t) = z(T, \tau, t)$ and $\theta(t) = \theta(T, \tau, t)$, treating each time variable as independent (Bender & Orszag 1999). This transforms the proper time derivative via

$$ \frac{d}{dt} \mapsto \omega \frac{\partial}{\partial T} + \omega^{1/2} \frac{\partial}{\partial \tau} + \frac{\partial}{\partial t}, \qquad (2.2) $$

which transforms (2.1) into the system of partial differential equations

$$ \omega z_T + \omega^{1/2} z_{\tau} + z_t = \omega^{1/2} u(T) \sin \theta, \qquad (2.3a) $$

$$ \omega\theta_T + \omega^{1/2}\theta_{\tau} + \theta_t = \omega^{1/2}\gamma z(1 - B(T) \cos 2\theta). \qquad (2.3b) $$

Here and hereafter, subscripts for $t$, $\tau$ and $T$ denote partial derivatives. We will later remove the extra degrees of freedom that we have introduced by imposing periodicity of the dynamics in the intermediate and fast variables $\tau$ and $T$. Expanding $z$ and $\theta$ in powers of $\omega^{-1/2}$ as $z \sim z_0 + \omega^{-1/2}z_1 + \omega^{-1}z_2 + \dots$ and $\theta \sim \theta_0 + \omega^{-1/2}\theta_1 + \omega^{-1}\theta_2 + \dots$, we obtain the $O(\omega)$ balance

$$ z_{0T} = 0, \quad \theta_{0T} = 0, \qquad (2.4a,b) $$

so that $z_0 = z_0(\tau, t)$ and $\theta_0 = \theta_0(\tau, t)$ are independent of $T$. To determine how $z_0$ and $\theta_0$ depend on $\tau$ and $t$, we must proceed to higher asymptotic orders.

We next consider the balance of $O(\omega^{1/2})$ terms in (2.3), which reads

$$ z_{1T} + z_{0\tau} = u(T) \sin \theta_0, \qquad (2.5a) $$

$$ \theta_{1T} + \theta_{0\tau} = \gamma z_0 (1 - B(T) \cos 2\theta_0). \qquad (2.5b) $$

The Fredholm solvability condition for (2.5) is equivalent to averaging over a period in $T$ and enforcing $T$-periodicity of $z_1$ and $\theta_1$. Introducing the averaging operator $\langle \cdot \rangle$, defined
via its action on functions $v(T, \tau, t)$ via

$$ \langle v \rangle(\tau, t) := \int_0^1 v(T, \tau, t) \, dT, \qquad (2.6) $$

we obtain the averaged equations

$$ z_{0\tau} = \langle u \rangle \sin \theta_0, \qquad (2.7a) $$

$$ \theta_{0\tau} = \gamma z_0 (1 - \langle B \rangle \cos 2\theta_0), \qquad (2.7b) $$

where $\langle u \rangle$ and $\langle B \rangle$ are the averages of $u(T)$ and $B(T)$, respectively, representing the average speed and shape of the model swimmer. In particular, $\langle u \rangle$ and $\langle B \rangle$ are constant, with the dynamics being rendered trivial if $\langle u \rangle = 0$; we exclude this case from our analysis and henceforth take $\langle u \rangle > 0$ without further loss of generality (the mapping $\theta \mapsto \theta + \pi$ transforms $\langle u \rangle < 0$ into the positive case). We will also assume that $|\langle B \rangle| < 1$, which imposes only a minimal restriction on the admissible swimmer shapes, since $|B| \ge 1$ is typically associated with objects of exceedingly large aspect ratio (Bretherton 1962).

Of particular note, if viewed as a system of ODEs in $\tau$, the system of (2.7) corresponds precisely to the original dynamical system of (2.1), suitably scaled, but with the time-varying speed and shape parameters replaced by their averages. We will shortly return to these equations and explore the ramifications of this observation in detail, in particular noting the existence of a Hamiltonian-like quantity, but first complete our analysis of the $O(\omega^{1/2})$ problem to determine the form of $z_1$ and $\theta_1$ for later convenience.

Without solving (2.7), we can deduce the form of $z_1$ and $\theta_1$ by substituting (2.7) into (2.5), yielding the simplified system

$$ z_{1T} = [u(T) - \langle u \rangle] \sin \theta_0, \qquad (2.8a) $$

$$ \theta_{1T} = -\gamma z_0 [B(T) - \langle B \rangle] \cos 2\theta_0. \qquad (2.8b) $$

Integrating (2.8) in $T$, recalling that $z_0$ and $\theta_0$ are independent of $T$, yields the solution

$$ z_1 = I_u(T) \sin \theta_0 + \bar{z}_1(\tau, t), \qquad (2.9a) $$

$$ \theta_1 = -\gamma z_0 I_B(T) \cos 2\theta_0 + \bar{\theta}_1(\tau, t), \qquad (2.9b) $$

where $\bar{z}_1$ and $\bar{\theta}_1$ are functions of $\tau$ and $t$, undetermined at this order, and we define

$$ I_u(T) := \int_0^T [u(\tilde{T}) - \langle u \rangle] d\tilde{T}, \quad I_B(T) := \int_0^T [B(\tilde{T}) - \langle B \rangle] d\tilde{T}. \qquad (2.10a,b) $$

Noting that $I_u(T)$ and $I_B(T)$ are $T$-periodic with period one, it follows that $z_1$ and $\theta_1$ are $T$-periodic with the same period.

In principle, one could proceed to the next asymptotic order to determine how $z_0$ and $\theta_0$ evolve in $t$ through the derivation of an additional solvability condition. However, here, this procedure would be complicated by the absence of an explicit solution to the nonlinear system of equations (2.7), compounded by the potentially $t$-dependent period of the solution, which would require using the generalised method of Kuzmak (1959). To circumvent this difficulty, we instead turn our attention back to (2.7), seeking further understanding of the leading-order dynamics over the intermediate time scale $\tau$.
Rheotaxis of shape-changing swimmers in Poiseuille flow

<figure><img src="image_4.png" /><figcaption>Figure 2. Phase portrait of motion on the intermediate time scale $\tau$. Solutions of (3.1) are closed orbits in the $z_0-\theta_0$ plane for constant $H_0$, symmetric in both $z_0 = 0$ and $\theta_0 = \pi$. Solutions in the shaded region, where $H_0 > g(0)$, do not cross $z_0 = 0$, corresponding to tumbling motion and monotonic evolution of $\theta_0$. Trajectories with $H_0 < g(0)$ instead exhibit swinging motion, with $\theta_0$ oscillating between two values. The black contour $H_0 = g(0)$ separates these regimes, with the direction of motion in the phase plane indicated by arrowheads, recalling that $\gamma \ge 0$. The point $(z_0, \theta_0) = (0, \pi)$ corresponds to rheotaxis, with $H_0 = g(\pi)$.</figcaption></figure>

### 3. A Hamiltonian-like quantity

If treated as a system of ODEs, we noted that (2.7) closely resembles the original swimming problem, with averages taking the place of oscillatory swimming speeds and swimmer shapes. In fact, the equivalent ODE problem has been extensively studied, with Zöttl & Stark (2013) thoroughly exploring this dynamical system and identifying a Hamiltonian-like constant of motion. Motivated by their study, we identify an analogous first integral of (2.7):

$$H_0(t) := \frac{\gamma}{2\langle u \rangle} z_0^2 + g(\theta_0), \quad (3.1)$$

for $H_0(t) \in [g(\pi), \infty)$, where, for $\langle B \rangle \in (-1, 1)$,

$$g(\theta_0) := \frac{\operatorname{arctanh}\left(\sqrt{\frac{2\langle B \rangle}{1+\langle B \rangle}} \cos \theta_0\right)}{\sqrt{2\langle B \rangle(1+\langle B \rangle)}}, \quad g'(\theta_0) = -\frac{\sin \theta_0}{1-\langle B \rangle \cos 2\theta_0}, \quad (3.2a,b)$$

taking the appropriate limits and branches where required. As $H_0(t)$ is effectively a constant of motion over the intermediate time scale $\tau$, (3.1) demonstrates that solutions to (2.7) are closed orbits in $z_0-\theta_0$ phase space over this time scale, with $\theta_0$ appropriately understood to be taken modulo $2\pi$.

We illustrate this phase space in figure 2, equivalent to the plot of Zöttl & Stark (2013, figure 2b). In particular, it is helpful to emphasise the two distinct behavioural regimes on this time scale: (1) endless tumbling, where the swimmer does not cross $z_0 = 0$ and $\theta_0$ varies monotonically; and (2) periodic swinging, where the swimmer repeatedly crosses the midline of the flow and $\theta_0$ oscillates between two values, $\theta_0 \in [\theta_{min}, \theta_{max}]$, readily computed from (3.1). The separating trajectory passes through $(z_0, \theta_0) = (0, 0)$ and corresponds to $H_0 = g(0)$. Here, $H_0 > g(0)$ corresponds to tumbling, shaded grey in figure 2, and $H_0 < g(0)$ corresponds to swinging. Of note, the period of these dynamics over the intermediate time scale, which we denote by $P_\tau$, depends non-trivially on $H_0$.

We can identify these dynamics with those reported both numerically and experimentally by Omori et al. (2022). In particular, as $H_0(t)$ approaches its minimum
of $g(\pi)$, the trajectory in the phase space approaches a single point, $(z_0, \theta_0) = (0, \pi)$, corresponding to a swimmer that is directed upstream on the midline of the flow. This is precisely the so-called rheotactic behaviour observed by Omori et al. (2022), which they noted as the long-time behaviour of (1.1) for particular definitions of $u(T)$ and $B(T)$, corresponding also to the experimentally determined behaviours of *Chlamydomonas* in channel flow. This agreement suggests that the long-time dynamics of the full system may be captured by the evolution of the Hamiltonian-like quantity $H_0(t)$.

In order to examine this evolution, we consider the dynamics of the following Hamiltonian-like expression:

$$H(t) := \frac{\gamma}{2\langle u \rangle} z^2 + g(\theta). \qquad (3.3)$$

Taking the proper time derivative of (3.3) and inserting (2.1) yields

$$\frac{dH}{dt} = \gamma\omega^{1/2}z\sin\theta \left[ \frac{u(T)}{\langle u \rangle} - \frac{1 - B(T)\cos 2\theta}{1 - \langle B \rangle \cos 2\theta} \right]. \qquad (3.4)$$

Transforming the time derivative following (2.2), we then insert our expansions of $z$ and $\theta$ into (3.4), noting that $H$ is equal to $H_0$ at leading order. Expanding $H \sim H_0 + \omega^{-1/2}H_1 + \omega^{-1}H_2 + \dots$, the balance at $O(\omega)$ is simply $H_{0T} = 0$, so we deduce that $H_0 = H_0(\tau, t)$, as expected. At $O(\omega^{1/2})$, we have

$$H_{1T} + H_{0\tau} = \gamma z_0 \sin \theta_0 \left[ \frac{u(T)}{\langle u \rangle} - \frac{1 - B(T) \cos 2\theta_0}{1 - \langle B \rangle \cos 2\theta_0} \right]. \qquad (3.5)$$

Averaging over the fast time scale $T$, equivalent to applying the Fredholm solvability condition, we immediately see that the term in square brackets vanishes, so that $H_{0\tau} = 0$ and $H_0 = H_0(t)$ is also independent of $\tau$, as expected.

Finally, we consider the $O(1)$ terms in (3.4), which may be stated as

$$H_{2T} + H_{1\tau} + H_{0t} = h(T, \tau, t), \qquad (3.6)$$

where

$$h(T, \tau, t) := \gamma(z_0\theta_1 + z_1 \sin\theta_0) \left[ \frac{u(T)}{\langle u \rangle} - \frac{1 - B(T) \cos 2\theta_0}{1 - \langle B \rangle \cos 2\theta_0} \right] \\ - \frac{2\gamma z_0\theta_1 \sin\theta_0 \sin 2\theta_0}{(1 - \langle B \rangle \cos 2\theta_0)^2} [B(T) - \langle B \rangle], \qquad (3.7)$$

and we note that the expressions in the square brackets each average to zero over a period in $T$. Inserting our expressions for $z_1$ and $\theta_1$ from (2.9), we have

$$h(T, \tau, t) = \gamma(\sin^2 \theta_0 I_u(T) - \gamma z_0^2 \cos \theta_0 \cos 2\theta_0 I_B(T)) \left[ \frac{u(T)}{\langle u \rangle} - \frac{1 - B(T) \cos 2\theta_0}{1 - \langle B \rangle \cos 2\theta_0} \right] \\ + \frac{\gamma^2 z_0^2 \sin \theta_0 \sin 2\theta_0 \cos 2\theta_0}{(1 - \langle B \rangle \cos 2\theta_0)^2} I_B(T)(B(T) - \langle B \rangle) + [\cdot], \qquad (3.8)$$

with $[\cdot]$ encompassing terms that have zero fast-time-scale average, which here are those involving $\bar{z}_1$ and $\bar{\theta}_1$. It is also helpful to note that $2I_B(T)(B(T) - \langle B \rangle) = (I_B^2)_T$, so that $\langle I_B(B - \langle B \rangle) \rangle = 0$ by the periodicity of $I_B$. Hence, the entire second line of (3.8) will
vanish when averaged over a period in T. Averaging over T and noting the further relations
$\langle I_B B \rangle = \langle I_B \langle B \rangle \rangle$, $\langle I_u u \rangle = \langle I_u \langle u \rangle \rangle$, and

$$
(I_u I_B)_T = I_u(T)(B(T) - \langle B \rangle) + I_B(T)(u(T) - \langle u \rangle), \quad (3.9)
$$

the fast-time-scale average of $h$ is simply

$$
\langle h \rangle = \gamma \cos 2\theta_0 \left( \frac{\gamma}{\langle u \rangle} z_0^2 \cos \theta_0 + \frac{\sin^2 \theta_0}{1 - \langle B \rangle \cos 2\theta_0} \right) \underbrace{\langle I_u (B - \langle B \rangle) \rangle}_{W}, \quad (3.10)
$$

where $W := \langle I_u(B - \langle B \rangle) \rangle$ is constant.

Thus, averaging (3.6) over a period in T and then over a period in τ (here, noting that
H₀ is independent of τ, we can naively average over a single oscillation in τ, despite the t
dependence of Pτ, which would otherwise require the method of Kuzmak (1959)) yields
the long-time evolution equation

$$
\frac{dH_0}{dt} = \gamma f(H_0)W, \tag{3.11}
$$

where

$$
f(H_0) := \frac{1}{P_{\tau}} \int_{0}^{P_{\tau}} \cos 2\theta_0 \left( \frac{\gamma}{\langle u \rangle} z_0^2 \cos \theta_0 + \frac{\sin^2 \theta_0}{1 - \langle B \rangle \cos 2\theta_0} \right) d\tau, \quad (3.12)
$$

recalling that $P_\tau$ denotes the period of the oscillatory dynamics in $\tau$ and $H_0$ is independent of $\tau$ and $T$. Notably, the integrand of (3.12) depends on the swimmer's speed and shape only via their fast-time averages $\langle u \rangle$ and $\langle B \rangle$. Therefore, all of the information encoding the dynamic variation of $u$ and $B$ about their mean is solely contained within the swimmer-dependent constant $W$. Of particular note, if either $u$ or $B$ is constant, then $W=0$ and, hence, $dH_0/dt = 0$, so that there is no long-time drift of $H_0$. This analytically verifies the numerical observations of Omori et al. (2022), who concluded that oscillations in both swimmer speed and shape were required to modify long-time behaviour.

Having reduced the dynamics to the one-dimensional autonomous dynamical system of (3.11), notably independent of $\omega$, it remains to understand $f(H_0)$, the average of a particular function of $z_0$ and $\theta_0$ over a period in $\tau$, which we illustrate in figure 3. In the Appendix, we analytically demonstrate that $f(H_0) \le 0$ for all $H_0$, in agreement with figure 3. Hence, the sign of (3.11) is determined by $W$, which depends only on the dynamics of $u$ and $B$ over a single oscillation. Strictly, there is a higher-order problem to be solved close to $H_0 = g(0)$, evidenced by the cusp-like profile in figure 3, with $P_\tau \to \infty$ and $f(H_0) \to 0$ as $H_0 \to g(0)$. However, noting that $f(H_0) < 0$ either side of $H_0 = g(0)$, this point is half-stable in the context of (3.11) (Strogatz 2018), so that it is unstable in practice and does not materially impact on the evolving dynamics.

Thus, the fixed point at $H_0 = g(\pi)$, corresponding to the rheotactic configuration ($z_0, \theta_0 = (0, \pi)$), is globally stable if $W > 0$ and unstable if $W < 0$. Hence, rheotaxis is the globally emergent behaviour at leading order if $W > 0$, whilst endless tumbling arises if $W < 0$.

4. Summary and conclusions

Our analysis allows us to characterise the long-time behaviour of a swimmer in Poiseuille
flow via the computation of a simple quantity, W, defined in (3.10) and dependent only
on the dynamics of the speed *u* and the shape parameter *B* of the swimmer over a
<figure><img src="image_5.png" /><figcaption>Figure 3. Exemplifying $f(H_0)$. We plot an example $f(H_0)$, as defined in (3.12) and computed numerically, for a range of $H_0$. The non-positivity of $f(H_0)$ is immediately evident, with $f \to 0$ from below as $H_0 \to g(\pi)$ or $H_0 \to g(0)$. As noted in the main text, $f$ is undefined at $H_0 = g(0)$, which we indicate with a hollow circle, but this point is readily seen to be half-stable in the context of the dynamical system of (3.11), so has negligible impact on the dynamics in practice. Here, we have fixed $\gamma = 1$, $\langle u \rangle = 1$ and $\langle B \rangle = 0.5$. The shaded region corresponds to tumbling dynamics.</figcaption></figure>

single oscillation. In particular, we find that the long-time behaviours take one of three possible forms, given in terms of the leading-order Hamiltonian-like quantity $H_0$ of (3.1): (1) endless tumbling at increasing distance from the midline of the flow ($H_0 \to \infty$); (2) preserved initial behaviour of the swimmer ($H_0 = H_0(0)$); or (3) convergence to upstream rheotaxis, where the swimmer is situated at the midline of the flow ($H_0 \to g(\pi)$). We find that the drift towards these long-time regimes is caused by the delicate higher-order interactions in the system. Specifically, the nonlinear interaction between the small $O(\omega^{-1/2})$ variations from the leading-order system over the fast time scale $t = O(\omega^{-1})$ gives rise to the significant $O(1)$ effect over the slow time scale $t = O(1)$.

In the context of Omori et al. (2022), where $u(T) = \alpha + \beta \sin(2\pi T)$ and $B(T) = \delta + \mu \sin(2\pi T + \lambda)$, we note that in-phase oscillations with $\lambda \in \{n\pi | n \in \mathbb{Z}\}$ immediately lead to $W = 0$, corresponding to case (2) above. Any other values of $\lambda$ lead to $W \neq 0$ (with maximal magnitude for $\lambda \in \{\pi/2 + n\pi | n \in \mathbb{Z}\}$) and long-term evolution of the swimmer behaviour. Notably, shifting $\lambda$ by $\pi$ results in a precise reversal of the sign of $W$ and a corresponding reversal of the sign of $dH_0/dt$. Hence, this shift in phase will precisely flip the fate of the swimmer, with rheotaxis being replaced by tumbling, or vice versa. Concretely, if $\beta\mu > 0$, then $\lambda \in (0, \pi)$ results in tumbling, whilst $\lambda \in (\pi, 2\pi)$ gives rise to rheotaxis. Further, our analysis predicts that swimmers having either $u$ or $B$ constant will not undergo a similar drift over $t = O(1)$ at leading order. In particular, this highlights that rigid externally driven swimmers and Janus particles, associated with constant $B$, would differ fundamentally in their long-time behaviour from shape-changing swimmers with dynamically varying $B$.

In support of our asymptotic analysis, we present three numerical examples in figure 4, approximating $f(H_0)$ with quadrature using the easily obtained numerical solution of the simple ODE system of (2.7) in $\tau$ for fixed $H_0(t)$. In figure 4(a), we compare the asymptotic and full numerical solutions through the evolution of $H$, adopting the sinusoidal model of Omori et al. (2022) and demonstrating excellent agreement between the solutions. This numerical validation spans the distinct dynamical regimes of tumbling and swinging, with different values of the phase shift $\lambda$ giving rise to distinct behaviours from otherwise identical initial conditions. Figure 4(b) captures a transition between behaviours, as observed by Omori et al. (2022) for $\lambda = 6\pi/5$, where we have predicted the bounds of $z$ oscillations by combining the solution of (3.11) with (3.1) evaluated at $\theta_0 = \theta_{min}$ and $\theta_0 = \pi$. We anticipate that similar calculations may be used to predict collisions between swimmers and channel boundaries in related experimental set-ups, though theoretical consideration of boundary effects in future work is warranted and may complicate analysis.
<figure><img src="image_6.png" /><figcaption>Rheotaxis of shape-changing swimmers in Poiseuille flow</figcaption></figure>

Figure 4. Numerical validation. (a) The value of *H*, as computed from the full numerical solution of (2.1) and the approximation of (3.11), shown as blue and black curves, respectively, for three phase shifts $\lambda \in \{4\pi/5, \pi, 6\pi/5\}$. Small, rapid oscillations in the full numerical solution are visible in the inset. (b) The asymptotically predicted bounds of *z* oscillations for $\lambda = 6\pi/5$ are shown as black curves, with the rapidly oscillating full solution shown in blue, highlighting excellent agreement even when the full solution transitions from tumbling dynamics towards rheotactic behaviour. Here, we have taken $(\alpha, \beta, \delta, \mu) = (1, 0.5, 0.32, 0.3)$ and $\lambda \in \{4\pi/5, \pi, 6\pi/5\}$ in the sinusoidal model of Omori *et al.* (2022), fixing $\gamma = 1, \omega = 50$ and $(z, \theta) = (1, \pi/4)$ initially. The shaded regions correspond to tumbling dynamics.

A promising direction for such an exploration is the augmentation of the model of Zöttl & Stark (2012) with the effects of an oscillatory swimming speed and swimmer shape, equivalent to extending the studied model of Omori *et al.* (2022) to a confined geometry.

More generally, the study of long-time swimmer dynamics in other scenarios, such as extensional flows, merits further investigation via multiple-scale analysis, with the behaviours of microswimmers in flow having been the subject of extensive research interest since the early twentieth century, as summarised by Bretherton & Rothschild (1961). Recent examples of this include, but are not limited to, the works of Marcos *et al.* (2012), Miki & Clapham (2013) and Uspal *et al.* (2015), which consider the rheotaxis of bacteria, spermatozoa and active particles, respectively. The study of microswimmers in flow has also recently been considered in the context of theoretical control and guidance (Walker *et al.* 2018; Moreau & Ishimoto 2021; Moreau *et al.* 2021), further motivating the development of our understanding of the potentially subtle interactions between swimmer shape, swimming speed and background flows.

In summary, an asymptotic (three-time-scale) multiple-scale analysis of the swimming model of Omori *et al.* (2022) has revealed a trichotomy of startlingly simple long-time behaviours, determined only by a single swimmer-specific constant of motion that may be readily computed *a priori*. This analysis, which complements the earlier works of Zöttl & Stark (2012, 2014) and Junot *et al.* (2019), confirms the numerical predictions of Omori *et al.* (2022), in agreement with their experimental observations of *Chlamydomonas*, and formally identifies the interacting oscillatory effects needed to elicit the eventual behaviours of endless tumbling and upstream rheotaxis in this model of swimming.

**Acknowledgements.** The computer code used in this study is available at https://gitlab.com/bjwalker/emergent-rheotaxis-in-Poiseuille-flow.

**Funding.** B.J.W. is supported by the Royal Commission for the Exhibition of 1851. K.I. acknowledges JSPS-KAKENHI for Young Researchers (grant no. 18K13456), JSPS-KAKENHI for Transformative Research Areas (grant no. 21H05309) and JST, PRESTO, Japan (grant no. JPMJPR1921). C.M. is supported by the Research Institute for Mathematical Sciences, an International Joint Usage/Research Center located at Kyoto University. M.P.D. is supported by the UK Engineering and Physical Sciences Research Council (grant no. EP/W032317/1).

**Declaration of interests.** The authors report no conflict of interest.
Author ORCIDs.

B.J. Walker https://orcid.org/0000-0003-0853-267X;

K. Ishimoto https://orcid.org/0000-0003-1900-7643;

C. Moreau https://orcid.org/0000-0002-8557-1149;

E.A. Gaffney https://orcid.org/0000-0002-6888-4362;

M.P. Dalwadi https://orcid.org/0000-0001-5017-2116.

Appendix. Deducing the sign of $f(H_0)$

Consider the integrand of (3.12). Recalling the evolution equations of $z_0$ and $\theta_0$ on the intermediate time scale, given in (2.7), this can be written as

$$ \frac{\cos 2\theta_0}{\langle u \rangle (1 - \langle B \rangle \cos 2\theta_0)} (z_0 \sin \theta_0)_\tau. \quad (A1) $$

Integrating by parts then yields

$$ P_{\tau}f(H_0) = \int_{0}^{P_{\tau}} \frac{2z_0 \sin \theta_0 \sin 2\theta_0}{\langle u \rangle (1 - \langle B \rangle \cos 2\theta_0)^2} \theta_0 \tau \, d\tau = \int \frac{2z_0 \sin \theta_0 \sin 2\theta_0}{\langle u \rangle (1 - \langle B \rangle \cos 2\theta_0)^2} \, d\theta_0, \quad (A2) $$

with boundary terms vanishing due to periodicity. With reference to the phase diagram of figure 2, we note that we need only integrate in $\theta_0$ from its minimum attained value up to $\pi$, with both the integrand and the phase-space trajectory being symmetric about both $\theta_0 = \pi$ and $z_0 = 0$. This corresponds to integrating over the branch of the trajectory in the upper-left quadrant of figure 2, with the true value of $P_{\tau}f(H_0)$ then being recovered by appropriate multiplication by 2 or 4, depending on whether the trajectory is one of tumbling or swinging. Hence, we consider the integral only over the range $\theta_0 \in [\theta_{min}, \pi]$, where $\theta_{min} \in [0, \pi]$ is the minimum value attained by $\theta_0$ over an orbit, as specified solely by $H_0(t)$.

In the upper-left quadrant of the phase plane, $z_0$ is a non-negative increasing function of $\theta_0$, evident from (2.7a) and (3.1). In particular, $z_0(\theta_0 + \pi/2) \ge z_0(\pi/2)$ for $\theta_0 \in [0, \pi/2]$. Further, the remaining combination of trigonometric terms in the integrand, denoted by $I(\theta_0)$ for brevity, satisfies $I(\theta_0) \ge 0$ and $I(\theta_0 + \pi/2) = -I(\theta_0)$ for $\theta_0 \in [0, \pi/2]$. Hence, the integral is trivially negative for $\theta_{min} \in [\pi/2, \pi]$, whilst for $\theta_{min} \in [0, \pi/2]$ we have

$$ \begin{aligned} \int_{\theta_{min}}^{\pi} z_0(\theta_0) I(\theta_0) \, d\theta_0 &= \int_{\theta_{min}}^{\pi/2} z_0(\theta_0) I(\theta_0) \, d\theta_0 + \int_{\pi/2}^{\pi} z(\theta_0) I(\theta_0) \, d\theta_0 \\ &\le \int_{0}^{\pi/2} z_0(\pi/2) I(\theta_0) \, d\theta_0 + \int_{0}^{\pi/2} z_0(\theta_0 + \pi/2) I(\theta_0 + \pi/2) \, d\theta_0 \\ &= \int_{0}^{\pi/2} [z_0(\pi/2) - z_0(\theta_0 + \pi/2)] I(\theta_0) \, d\theta_0 \le 0. \end{aligned} \quad (A3) $$

Hence, $P_{\tau}f(H_0) \le 0$, so that $f(H_0) \le 0$ for all $H_0$. In particular, this equality is strict unless $H_0 = g(\pi)$ or $H_0 = g(0)$, which correspond to the degenerate rheotactic trajectory ($z_0, \theta_0) = (0, \pi)$ and the separating trajectory that lies between tumbling and swinging behaviours in figure 2. As discussed in the main text, the case with $H_0 = g(0)$ requires the consideration of a higher-order problem, though has negligible impact on the dynamics in practice.
Rheotaxis of shape-changing swimmers in Poiseuille flow

REFERENCES

BENDER, C.M. & ORSZAG, S.A. 1999 *Advanced Mathematical Methods for Scientists and Engineers I.* Springer.

BRETHERTON, F.P. 1962 The motion of rigid particles in a shear flow at low Reynolds number. *J. Fluid Mech.* **14** (2), 284–304.

BRETHERTON, F.P. & ROTHSCHILD, N.M.V. 1961 Rheotaxis of spermatozoa. *Proc. R. Soc. B: Biol. Sci.* **153** (953), 490–502.

GAFFNEY, E.A., DALWADI, M.P., MOREAU, C., ISHIMOTO, K. & WALKER, B.J. 2022 Canonical orbits for rapidly deforming planar microswimmers in shear flow. *Phys. Rev. Fluids* **7** (2), L022101.

GUASTO, J.S., RUSCONI, R. & STOCKER, R. 2012 Fluid mechanics of planktonic microorganisms. *Annu. Rev. Fluid Mech.* **44** (1), 373–400.

JUNOT, G., FIGUEROA-MORALES, N., DARNIGE, T., LINDNER, A., SOTO, R., AURADOU, H. & CLÉMENT, E. 2019 Swimming bacteria in Poiseuille flow: the quest for active Bretherton-Jeffery trajectories. *Europhys. Lett.* **126** (4), 44003.

KIM, S. & KARRILA, S.J. 2005 *Microhydrodynamics: Principles and Selected Applications*. Dover.

KUZMAK, G. 1959 Asymptotic solutions of nonlinear second order differential equations with variable coefficients. *J. Appl. Math. Mech.* **23** (3), 730–744.

MA, K., PUJARA, N. & THIFFEAULT, J.-L. 2022 Reaching for the surface: spheroidal microswimmers in surface gravity waves. *Phys. Rev. Fluids* **7** (1), 014310.

MARCOS, FU, H.C., POWERS, T.R. & STOCKER, R. 2012 Bacterial rheotaxis. *Proc. Natl Acad. Sci.* **109** (13), 4780–4785.

MIKI, K. & CLAPHAM, D.E. 2013 Rheotaxis guides mammalian sperm. *Curr. Biol.* **23** (6), 443–452.

MOREAU, C. & ISHIMOTO, K. 2021 Driving a microswimmer with wall-induced flow. *Micromachines* **12** (9), 1025.

MOREAU, C., ISHIMOTO, K., GAFFNEY, E.A. & WALKER, B.J. 2021 Control and controllability of microswimmers by a shearing flow. *R. Soc. Open Sci.* **8** (8), 211141.

OMORI, T., KIKUCHI, K., SCHMITZ, M., PAVLOVIC, M., CHUANG, C.H. & ISHIKAWA, T. 2022 Rheotaxis and migration of an unsteady microswimmer. *J. Fluid Mech.* **930**, A30.

STROGATZ, S.H. 2018 *Nonlinear Dynamics and Chaos*. CRC.

USPAL, W.E., POPESCU, M.N., DIETRICH, S. & TASINKEVYCH, M. 2015 Rheotaxis of spherical active particles near a planar wall. *Soft Matt.* **11** (33), 6613–6632.

WALKER, B.J., ISHIMOTO, K., GAFFNEY, E.A., MOREAU, C. & DALWADI, M.P. 2022 Effects of rapid yawing on simple swimmer models and planar Jeffery's orbits. *Phys. Rev. Fluids* **7** (2), 023101.

WALKER, B.J., ISHIMOTO, K., WHEELER, R.J. & GAFFNEY, E.A. 2018 Response of monoflagellate pullers to a shearing flow: a simulation study of microswimmer guidance. *Phys. Rev. E* **98** (6), 063111.

ZÖTTL, A. & STARK, H. 2012 Nonlinear dynamics of a microswimmer in Poiseuille flow. *Phys. Rev. Lett.* **108** (21), 218104.

ZÖTTL, A. & STARK, H. 2013 Periodic and quasiperiodic motion of an elongated microswimmer in Poiseuille flow. *Eur. Phys. J. E* **36** (1), 4.

ZÖTTL, A. & STARK, H. 2014 Hydrodynamics determines collective motion and phase behavior of active colloids in quasi-two-dimensional confinement. *Phys. Rev. Lett.* **112** (11), 118101.