# Effects of rapid yawing on simple swimmer models and planar Jeffery's orbits

Benjamin J Walker, Kenta Ishimoto, Eamonn A Gaffney, Clément Moreau,

Mohit P Dalwadi

► To cite this version:

Benjamin J Walker, Kenta Ishimoto, Eamonn A Gaffney, Clément Moreau, Mohit P Dalwadi. Effects of rapid yawing on simple swimmer models and planar Jeffery's orbits. Physical Review Fluids, 2022, 7 (2), pp.023101. 10.1103/physrevfluids.7.023101. hal-05050439

HAL Id: hal-05050439

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

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
# Effects of rapid yawing on simple swimmer models and planar Jeffery's orbits

Benjamin J. Walker<sup>1,2,*</sup> Kenta Ishimoto,<sup>3,†</sup> Eamonn A. Gaffney,<sup>1,‡</sup>

Clément Moreau<sup>3,§</sup> and Mohit P. Dalwadi<sup>1,2,||</sup>

<sup>1</sup>Mathematical Institute, University of Oxford, Oxford, OX2 6GG, United Kingdom

<sup>2</sup>Department of Mathematics, University College London, London, WC1H 0AY, United Kingdom

<sup>3</sup>Research Institute for Mathematical Sciences, Kyoto University, Kyoto, 606-8502, Japan

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

(Received 3 September 2021; accepted 21 December 2021; published 4 February 2022)

Over a sufficiently long period of time, or from an appropriate distance, the motion of many swimmers can appear smooth, with their trajectories appearing almost ballistic in nature and slowly varying in character. These long-time behaviors, however, often mask more complex dynamics, such as the side-to-side snakelike motion exhibited by spermatotozoa as they swim, propelled by the frequent and periodic beating of their flagellum. Many models of motion neglect these effects in favour of smoother long-term behaviors, which are often of greater practical interest than the small-scale oscillatory motion. While it may be tempting to ignore any yawing motion, simply assuming that any effects of rapid oscillations cancel out over a period, a precise quantification of the impacts of high-frequency yawing is lacking. In this study, we systematically evaluate the long-term effects of general high-frequency oscillations on translational and angular motion, cast in the context of microswimmers but applicable more generally. Via a multiple-scales asymptotic analysis, we show that rapid oscillations can cause a long-term bias in the average direction of progression. We identify sufficient conditions for an unbiased long-term effect of yawing, and we quantify how yawing modifies the speed of propulsion and the effective hydrodynamic shape when in shear flow. Furthermore, we investigate and justify the long-time validity of the derived leading-order solutions and, by direct computational simulation, we evidence the relevance of the presented results to a canonical microswimmer.

DOI: 10.1103/PhysRevFluids.7.023101

## I. INTRODUCTION

Swimming in nature is often driven by undulatory deformations, from the snakelike motion of spermatozoa at low Reynolds number to the inertia-dominated movement of fish and aquatic mammals [1]. When viewed over a sufficiently long time frame, or from an adequate distance, the motion of many of these self-propelled objects can appear smooth, with their trajectories appearing almost ballistic in nature and slowly varying in character. These long-term behaviors, however, often

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>*Corresponding author: benjamin.walker@maths.ox.ac.uk</footnote>

<footnote>†ishimoto@kurims.kyoto-u.ac.jp</footnote>

<footnote>‡gaffney@maths.ox.ac.uk</footnote>

<footnote>§cmoreau@kurims.kyoto-u.ac.jp</footnote>

<footnote>‖dalwadi@maths.ox.ac.uk</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 1. The path of a spermatozoon and a model self-propelled swimmer in flow. (a) The digitally tracked path of a human spermatozoon, shown as a black curve, with the smoothed average path shown in blue. (b) We illustrate an idealized swimmer, moving with speed *V* at an orientation *φ* from the laboratory *e*<sub>*X*</sub> axis. Though we show a model flagellated swimmer in shear flow, with *u* = −*γ*Y*e*<sub>*X*</sub>, this setup is also applicable in the absence of background fluid flow, equivalent to taking *γ* = 0, and to general self-propelled particles. Videomicroscopy from Movie 1 of the Supplemental Material of Ishimoto *et al.* [37].</figcaption></figure>

mask the details of complex short-term dynamics that arise from their relatively high-frequency
undulation. Though it is tempting to neglect the effects of yawing over long enough timescales,
rationalizing this by assuming that the effects of rapid oscillations cancel over a period, a precise
investigation of the impacts of high-frequency yawing on the long-term behavior of a swimmer
or particle is lacking, even for the simplest models. Indeed, the nonlinear angular dependence in
some simple models of motion, such as that of an active particle with constant speed and prescribed
orientation [2], suggests that the effects of rapid yawing may not simply average out to zero, as
might be expected in linear systems. Hence, in this study we seek to systematically evaluate the
long-term effects of simple yet general oscillations on translational and angular motion. While we
will frame much of this exploration in the context of both self-propelled and passive bodies in a
Stokesian fluid, much of the analysis, particularly that of Sec. II, applies more generally to active
particle models.

In the well-explored context of microswimming, there are many examples of yawing bodies,
with a multitude of swimmers being propelled by rapid oscillatory drivers. For instance, canonical
among microswimmers is the aforementioned spermatozoon, with a slender flagellum providing
the propulsive force that gives rise to its net head-first motion. Due to the progression of a
traveling wave of displacement along their flagellum, as is typical (though not ubiquitous) for
spermatozoa in common media [3–5], these cells undergo large side-to-side yawing, giving rise to
complex oscillatory trajectories on a fast timescale [see Fig. 1(a)], though their long-term behaviors
nevertheless appear to be smooth. In the study of this swimmer, along with that of many others, a
broad variety of methods have been employed of varying degrees of complexity, ranging from the
simplest phenomenological models, as we will consider in this work and are particularly common in
the large-scale study of active matter [6–8], to the most computationally and geometrically intricate
[9–14]. While the latter class of methods inherently includes the effects of rapid yawing in their
formulation, it being an emergent property of geometrically faithful study, many simpler models
neglect these effects in favour of capturing the smoother, long-term behaviors, which are often
of greater interest than the rapid transient motion in applications [6, 15–24]. However, despite the
success of simple long-term or averaged models, it is not clear that the often seen neglect of the
yawing motion exhibited by these swimmers is justified. Hence, the primary aim of this study will
be to conduct a detailed investigation into the effects of rapid yawing on simple swimmer models,
motivated by the valuable utility and persisting prevalence of such minimal representations.
A common simplifying feature of these models is to assume symmetric swimmer geometry. These geometrically coarse models have been widely employed in many contexts, from large-scale investigations of active matter to single-swimmer studies of guidance and control [2,25–30]. Often the reduction in geometrical complexity affords many advantages, such as numerical efficiency and analytical tractability. For example, a recent *in silico* study of *Leishmania mexicana* evidenced that an ellipsoid can be an appropriate model for this flagellated swimmer in the presence of a background shear flow [31], enabling rapid simulation in comparison to the computationally expensive boundary element method from which it was derived. Ellipsoid models in particular afford the significant advantage of having been the subject of much classical and recent study, with their behavior in shear flow being the topic of Jeffery’s seminal work of 1922 [32]. Since then Jeffery’s analytical results, which we partially recount in Sec. III, have been extended and generalized [33–36]. In particular, the work of Bretherton [33] and the recent contribution of Ishimoto [36] considered broader classes of rigid bodies, including surfaces of revolution and those with particular symmetries. Though our understanding of the behaviors of such geometries in flow benefits significantly from this mature body of study, the impacts of rapid yawing on even the classical results of Jeffery’s formulation remain unknown, despite the noted prevalence of rapid oscillations in microswimmers and the broad utilization of simple models to study them. Thus, as a further aim of this work, we will investigate how, if at all, high-frequency yawing modifies the motion of particles admitted by Jeffery’s theory and its subsequent generalizations for symmetric bodies.

Hence, we proceed by first considering a paradigm model of motion for a self-propelled particle in Sec. II, rendered nontrivial by a prescribed oscillatory orientation that phenomenologically captures the features of rapid yawing. We will derive asymptotic solutions by exploiting the separated timescales of the problem, which relate to translation and rotation. This will allow us to understand the structure of the multiple-scales analysis that we use throughout this paper and to understand several basic implications of rapid yawing. Next, focusing on the case of planar rotations in a background shear flow, we will extend our analysis by introducing an appropriately coupled equation of motion for the rotational dynamics, seeking leading-order solutions in Sec. III. From numerical validation of these leading-order solutions for both translational and angular motion, we will note that the asymptotic solution is surprisingly accurate, even on very long timescales. Motivated by this, we will further seek to establish sufficient conditions for such persistent validity by conducting a higher-order analysis in Sec. IV.

## II. A SELF-PROPELLED, YAWING SWIMMER

We begin by studying a minimal model of a self-propelled yawing swimmer in the plane, in the absence of any background flow. In particular, we will consider a body moving at a constant speed $V$ in an unbounded domain at a time-dependent angle $\phi(t)$ from a fixed laboratory axis [see Fig. 1(b)], denoting time by $t$. Although we will later generalize several of these simplifications, including an extension to three-dimensional dynamics in Appendix A and the effect of a background shear flow coupled to angular motion in Sec. III, we first present the simple two-dimensional problem in the absence of flow. This will allow us to understand the essential structure of the problem and the approach, including several implications of rapid yawing. Here and hereafter, we consider all quantities to be dimensionless. While this entails that we may take $V = 1$ without loss of generality, we will retain the symbol $V$ so as to clearly identify the role of the swimming speed in the analysis that follows.

With the coordinates of the swimmer position written as $(x(t), y(t))$ relative to laboratory unit vectors $e_X$ and $e_Y$, the swimmer motion is governed by the following ordinary differential equations (ODEs):

$$ \frac{dx}{dt} = V \cos \phi, \quad x(0) = 0, \qquad (1a) $$

$$ \frac{dy}{dt} = V \sin \phi, \quad y(0) = 0, \qquad (1b) $$
measuring the orientation of $\phi(t)$ anticlockwise from the X-axis. We note that the initial conditions are chosen for convenience, and may be readily substituted. Here, we will seek to ascertain the effects of a rapidly oscillating $\phi$, modeling the yawing dynamics exhibited by many microswimmers, such as the flagellated spermatozoon or *Leishmania mexicana*.

### A. Sinusoidal yawing

Though we later consider general yawing motions, it is instructive to first prescribe simple, sinusoidal oscillations in the orientation $\phi$, seeking to develop our intuition for the contributions of rapid angular variation. We denote the angular frequency of the oscillations by $\omega \gg 1$, which we take to be large throughout this study in order to consider rapid yawing. Hence, in this section, we prescribe

$$ \phi(t) = \bar{\phi} + A \sin \omega t, \qquad (2) $$

where A is the observed amplitude of the angular yawing dynamics and $\bar{\phi}$ is a constant average orientation. Thus, we seek to analyze the system

$$ \frac{dx}{dt} = V \cos (\bar{\phi} + A \sin \omega t), \qquad (3a) $$

$$ \frac{dy}{dt} = V \sin (\bar{\phi} + A \sin \omega t), \qquad (3b) $$

where $\omega \gg 1$ separates the timescales of oscillation and translation. We will exploit this separation of timescales and pursue a classical multiple-scales analysis [38], introducing the fast variable $T := \omega t$. In the standard manner, we introduce the formal transformation $t \mapsto (t, T)$ and write $x = x(t, T)$ and $y = y(t, T)$, treating $t$ and $T$ as independent. The additional degree of freedom this affords will be removed later by imposing that our dependent variables are $2\pi$-periodic in $T$. The multiple-scales transformation means that $d/dt \mapsto \partial/\partial t + \omega\partial/\partial T$, which transforms Eq. (3) into the following system of partial differential equations (PDEs):

$$ \frac{\partial x}{\partial t} + \omega \frac{\partial x}{\partial T} = V \cos (\bar{\phi} + A \sin T), \qquad (4a) $$

$$ \frac{\partial y}{\partial t} + \omega \frac{\partial y}{\partial T} = V \sin (\bar{\phi} + A \sin T). \qquad (4b) $$

We proceed by posing asymptotic expansions of our dependent variables of the form

$$ x = x_0 + \frac{1}{\omega} x_1 + \dots, \quad y = y_0 + \frac{1}{\omega} y_1 + \dots \quad \text{as } \omega \to \infty, \qquad (5) $$

where the $x_i, y_i$ ($i = 0, 1, \dots$) are ostensibly functions of both $t$ and $T$. Substituting the expansions of Eq. (5) into Eq. (4) and comparing terms of $O(\omega)$ gives the leading-order problems

$$ \frac{dx_0}{dT} = 0, \qquad (6a) $$

$$ \frac{dy_0}{dT} = 0, \qquad (6b) $$

so that $x_0$ and $y_0$ are independent of $T$. Hence, we write $x_0 = \bar{x}_0(t)$ and $y_0 = \bar{y}_0(t)$, adopting the $\bar{\phantom{x}}$ notation to represent functions independent of $T$ throughout. The dependence of these functions on the slow time $t$, and therefore the leading-order behavior of the swimmer, is not yet known. To determine this, we must proceed to the next order, with the corresponding PDEs being

$$ \frac{d\bar{x}_0}{dt} + \frac{\partial x_1}{\partial T} = V \cos (\bar{\phi} + A \sin T), \qquad (7a) $$

$$ \frac{d\bar{y}_0}{dt} + \frac{\partial y_1}{\partial T} = V \sin (\bar{\phi} + A \sin T). \qquad (7b) $$
To obtain appropriate governing equations for $\bar{x}_0$ and $\bar{y}_0$, we integrate these PDEs over $T \in (0, 2\pi)$ and impose periodicity in $T$, which causes the terms involving $x_1$ and $y_1$ to vanish. The subsequent averaged equations are

$$ \frac{d\bar{x}_0}{dt} = V \langle \cos(\bar{\phi} + A \sin T) \rangle, \qquad (8a) $$

$$ \frac{d\bar{y}_0}{dt} = V \langle \sin(\bar{\phi} + A \sin T) \rangle, \qquad (8b) $$

where we have defined the averaging operator over the fast dynamics as

$$ \langle h \rangle = \frac{1}{2\pi} \int_{0}^{2\pi} h(t, T) \, dT, \qquad (9) $$

for functions $h(t, T)$. We note that the PDE multiple-scales formulation has been transformed back into ODEs for the leading-order solutions in Eq. (8), as we would expect. It is instructive to evaluate the averages in Eq. (8) explicitly. In order to do this, we note the following exact integral results:

$$ \frac{1}{2\pi} \int_{0}^{2\pi} \cos(A \sin T) \, dT = J_{0}(A), \quad \frac{1}{2\pi} \int_{0}^{2\pi} \sin(A \sin T) \, dT = 0, \qquad (10) $$

where $J_0$ is the zeroth-order Bessel function of the first kind. The former follows from standard integral representations of the Bessel function, and the latter follows from antisymmetry. Exploiting the linearity of the averaging operator, the averages on the right-hand side of (8) can be expanded as

$$ \langle \cos(\bar{\phi} + A \sin T) \rangle = \cos \bar{\phi} \langle \cos(A \sin T) \rangle - \sin \bar{\phi} \langle \sin(A \sin T) \rangle, \qquad (11a) $$

$$ \langle \sin(\bar{\phi} + A \sin T) \rangle = \sin \bar{\phi} \langle \cos(A \sin T) \rangle + \cos \bar{\phi} \langle \sin(A \sin T) \rangle. \qquad (11b) $$

Making use of the integral results of Eq. (10) in Eq. (11), we may evaluate the right-hand side of Eq. (8) to obtain the following ODEs for the leading-order translation:

$$ \frac{d\bar{x}_0}{dt} = V_e \cos \bar{\phi}, \qquad (12a) $$

$$ \frac{d\bar{y}_0}{dt} = V_e \sin \bar{\phi}, \qquad (12b) $$

where we define the effective translation speed of the yawing swimmer as $V_e = J_0(A)V$. Comparing these equations with the original system of Eq. (1), we see that the leading-order solution to the yawing system of Eq. (3) is given by the solution to ODEs that are analogous to the nonyawing problem. In particular, the effects of yawing manifest as an effective heading of $\bar{\phi}$ and an effective translation speed of $V_e$, at least to leading order. Notably, these results are independent of the yawing rate $\omega$ when $\omega \gg 1$. We will continue to observe this independence of the yawing rate in the leading-order translation in the more general problems that we consider later.

Analyzing the implications of our results in more detail, since $|J_0(A)| \le 1$, with equality only for $A=0$, the effect of sinusoidal yawing is to reduce the effective translation speed to $V_e = J_0(A)V$, with this reduction being in agreement with intuition. Given the oscillatory nature of the Bessel function, as illustrated in Fig. 2(a), there are amplitude ranges that would result in a net reversal of the translational motion, though these amplitudes may be larger than practically observed in typical microswimmers, requiring $|A|$ greater than approximately 2.40. Moreover, as this Bessel function has nontrivial zeros, a notable implication is that there are specific yawing amplitudes that would result in asymptotically small net motion. In these cases, the rapidly oscillating swimmer trajectory approximately resembles a figure of eight, as shown in Fig. 2(b), achieving very little net motion. An additional feature of this simple sinusoidal yawing is that the effective heading is precisely $\bar{\phi}$, the average of the prescribed $\phi$ over a single yawing oscillation, though we show below in Sec. II B that this result need not hold for more general, nonsinusoidal yawing.
<figure><img src="image_5.png" /><figcaption>FIG. 2. A Bessel function and a closed swimming trajectory. (a) The zeroth-order Bessel function of the first kind, $J_0(A)$, plotted as a function of the yawing amplitude $A$, pertinent to the case of sinusoidal yawing. This oscillatory function has magnitude unity only if $A = 0$, so that yawing serves to strictly reduce the effective speed of the swimmer, recalling $V_e = J_0(A)V$. Noting the existence of the zeros of $J_0$, we remark that particular values of $A$ will result in a reversed swimming direction or even zero effective speed, though these correspond to very large amplitude yawing, as might be expected. (b) For $A \approx 2.40$, which approximately satisfies $J_0(A) = 0$, we simulate the full system of ODEs numerically as described in Sec. II B, realizing the closed trajectory predicted by the asymptotic analysis, with displacements on the order of $1/\omega$. Here $V = 1$, $\omega = 20$, $\phi(T) = A \sin T$, and motion is simulated until $t = 10$ from $(x, y) = (0, 0)$.</figcaption></figure>

## B. General yawing

Though our multiple-scales analysis in the previous subsection was simplified by the choice of sinusoidal yawing, it is straightforward to generalize. Here, we briefly outline the analysis and its implications for the general case.

We consider a general observed angular yawing $G(T)$ in place of $A \sin T$ above, and a slowly varying “average” orientation $\bar{\phi}(t)$. These generalizations replace the prescribed $\phi$ in Eq. (2) with

$$ \phi(t, T) = \bar{\phi}(t) + G(T). \qquad (13) $$

Our only requirements for $G$ are that it is periodic (we assume that it has period $2\pi$ without loss of generality), and that $\langle G \rangle$ is constant, which we take to be zero without loss of generality. Intuitively, this latter condition is equivalent to restricting our consideration to swimmers that experience no net change in their orientation over each period of yawing, which we will assume throughout this study. In Appendix A, when considering swimmer motion in three dimensions, we discuss a physically relevant relaxation of this assumption, which is straightforward to accommodate. This condition is commonly satisfied by even complex shape deforming swimmers, such as planar-swimming spermatozoa and the biflagellated alga *Chlamydomonas reinhardtii* [22,39], an example of which we explore in detail in Appendix B.

The analysis of this general case proceeds in the same manner as for Sec. II A, with $A \sin T$ replaced by $G(T)$. Therefore, the $O(1)$ solvability condition is the equivalent of Eq. (8) and the key integral results that we require are the averages of $\cos G$ and $\sin G$. These may be succinctly encoded as the real and imaginary parts of the constant complex number

$$ z_P := \langle \exp iG \rangle = R_P \exp i\theta_P \in \mathbb{C}, \qquad (14) $$

with $R_P$ being the non-negative modulus of $z_P$. This notation allows us to rewrite the drift ODEs for the general case in the compact form

$$ \frac{d\bar{x}_0}{dt} = V_e \cos(\bar{\phi} + \theta_P), \qquad \frac{d\bar{y}_0}{dt} = V_e \sin(\bar{\phi} + \theta_P), \qquad (15a) $$

where

$$ V_e := R_P V, \qquad (15b) $$

mirroring the unperturbed swimmer system. Hence, we may define the general effective speed of propulsion as $V_e = R_P V$, and the complex argument $\theta_P$ as the directional bias introduced by the rapid yawing.
<figure><img src="image_6.png" /></figure>

FIG. 3. Exemplar yawing functions $G(T)$ and the associated translational motion. (a) A sinusoidal yawing function, which is odd about $T_o^* = \pi$, shown dashed. (b) An antiperiodic yawing function, satisfying $G(T + \pi) = -G(T)$, with $G$ on $[0, \pi]$ being the reflection of $G$ on $[\pi, 2\pi]$ along the $T$-axis. (c) A biased yawing function, satisfying neither of the presented antisymmetry conditions and with $\langle \sin G \rangle \neq 0$. (d) The numerical solutions of the full yawing problem, taking $\phi = G(T)$ so that the direction of motion in the absence of yaw is along the $X$-axis, shown for each of the considered yawing functions and offset vertically for visual clarity. Having simulated the motion in each case for the same amount of time, the dependence of the effective speed on the form of the yawing function is apparent, with case (b) giving rise to the largest $V_e$ and, hence, the greatest displacement from the initial configuration. Cases (a) and (b), which satisfy $\langle \sin G \rangle = 0$, give rise to unbiased motion, while the yawing function of panel (c) yields biased motion that is not aligned along the axis, despite the average swimmer orientation $\phi$ being zero. In blue, we plot the leading-order solution of the multiple-scales analysis in each case, demonstrating excellent agreement with the full yawing behavior, and have made use of the asymptotically corrected initial condition of Sec. IV. All considered yawing functions satisfy $\langle G \rangle = 0$, and we have fixed $\omega = 40\pi$ for simulation.

The effective speed of propulsion $V_e$ is directly comparable with the specific case of sinusoidal yawing that we considered in Sec. II A, where $R_p \le 1$ is a more general version of the specific Bessel function result for sinusoidal yawing and inherits the accompanying intuition and remarks. Since $\bar{\phi}(t)$ represents the heading of the swimmer in the absence of yaw, we may interpret the complex argument $\theta_p$ as the directional bias introduced by the rapid yawing. This general bias was not present in the case of sinusoidal yaw (with the exception of the potential reversal of the effective swimmer direction). The reason for this is that the directional bias will vanish if $\langle \sin G \rangle = 0$, which is the case for $G(T) = \sin(T)$. Intuitively, this condition is equivalent to requiring that the average of the component of the velocity perpendicular to the mean heading is zero, an interpretation that is most easily seen if we consider $\phi = 0$, wherein $V \sin G$ is precisely the $Y$-component of the motion. Similarly, $\langle \cos G \rangle$ quantifies the average component of the velocity in the direction of the average heading, along the $X$-axis if $\bar{\phi} = 0$, and, hence, modifies the swimming speed.

It is interesting to consider forms of generalized yawing that lead to unbiased swimming, through our derived condition $\langle \sin G \rangle = 0$. In particular, we may identify two general classes of antisymmetric dynamics that give rise to unbiased swimming. The first category is simply that of all odd functions, satisfying $G(T_o^* + T) = -G(T_o^* - T)$ for some $T_o^*$, which includes $G(T) = A \sin T$ with $T_o^* = 0$, for instance. The second category is all functions that satisfy $G(T + \pi) = -G(T)$, also known as *antiperiodic functions*, a condition that captures the notion of mirror symmetry of the yawing over a period. In Fig. 3 we present exemplar $G(T)$ that satisfy each of these conditions, and note that, for any given $G(T)$, the quantity $z_p$ can be easily computed numerically, yielding the effective propulsive velocity and the directional bias. Here and throughout, systems of ordinary differential equations are solved viaode45, a standard routine in MATLAB [40]. Further, for any
swimmer whose geometrical configuration in a swimmer-fixed reference frame at time $T + \pi$ is simply a reflection of the swimmer configuration at time $T$ in a swimmer-fixed plane, such as the model spermatozoon showcased in Appendix B, we can define the swimmer orientation such that the yawing function $G(T)$ will be antiperiodic. This follows from symmetry and the time independence of Stokes equations, as described in more detail in Appendix C. Hence, we may immediately conclude that the yawing of such swimmers does not contribute any directional bias to their motion at leading order, at least in the context of this simplified model of motion.

As is common in the method of multiple scales, it is pertinent to understand the effects of longer-term drift of the full solution away from the derived leading-order expansion. We investigate this numerically (not shown), finding that the leading-order solution remains accurate on significant timescales, over tens of thousands of yawing oscillations, explored for $\omega = 40\pi$, $V = 1$, and $A = 1$ in the case of sinusoidal yawing. Given that we have only calculated solutions up to leading order, we would expect an accumulating error at $O(1/\omega)$. However, we do not observe this. That is, our asymptotic analysis appears to be more accurate than we might expect it to be. We will later return to this point in Sec. IV and in Appendix A, in which we justify this surprising improved accuracy by consideration of the higher-order problem.

Having set up and illustrated the general framework for considering rapid yawing motions, we next extend our analysis to a class of bodies in the presence of a canonical background flow. This will mean that the angular motion of the swimmer, $\phi$, is coupled to the external flow, such that we will no longer be able to directly prescribe the yawing rotational dynamics, as we did above.

### III. YAWING IN SHEAR FLOW

Described in the seminal work of 1922, Jeffery's orbits capture the motion of a rigid torqu-free ellipsoid in a background shear flow with a no-slip condition on its surface [32]. Through the solution of the governing Stokes equations, Jeffery derived ordinary differential equations in terms of an Euler angle formulation, and in some cases noted their exact solution. The relevant nonlinear equation for the evolution of the azimuthal angle $\phi$ decouples from the other Euler angles and is

$$ \frac{d\phi}{dt} = \frac{\gamma}{2}(1 - B \cos 2\phi) := f(\phi; B). \qquad (16) $$

Here, $\gamma$ denotes the constant shear rate of the background flow, $B$ is a constant derived from the geometrical properties of the ellipsoid, and we will impose an initial condition of $\phi(0) = 0$ for notational simplicity, though this may be readily substituted. The geometrical parameter $B$ differs from the notation of Jeffery [32], instead reflecting the significant generalization of Jeffery's orbits by Bretherton [33], which extended the classical description to include bodies of revolution, among additional contributions. As such, we will refer to $B$ as the Bretherton constant, which, while it may assume any real value, has magnitude less than unity for many well-studied bodies. For $|B| < 1$, which includes both prolate ($B > 0$) and oblate ($B < 0$) spheroids, the dynamics of Eq. (16) are precisely those of the classic Jeffery's orbits, with periodic rotation occurring ad infinitum; for $|B| > 1$, associated most often to objects with very large aspect ratio, the solution approaches one of two stable steady states, as discussed and exemplified by Bretherton [33] and more recently by Ishimoto [35]. This latter work, and its recent generalization [36], further extended the Jeffery's orbit description to bodies with particular symmetries, far removed from the original ellipsoidal considerations of Jeffery. Despite this, the solution to Eq. (16) can be written in the simple form found by Jeffery, which we write as

$$ \phi(t) = \phi_J(t; B) = \arctan \left( \frac{1}{r} \tan \frac{\gamma t}{r + \frac{1}{r}} \right), \quad r := \sqrt{\frac{1+B}{1-B}} \in \mathbb{C}, \qquad (17) $$

where $r$ is precisely the aspect ratio of the rigid body if the object in question is a spheroid. From Eq. (17), the distinction between the cases of $|B| < 1$ and $|B| > 1$ is clear, corresponding to $r$ being real or imaginary, respectively. With reference to its historical origin, we will refer to this solution
as a Jeffery's orbit irrespective of the magnitude of B, with the understanding that this "orbit" may
in fact lead to a steady solution for bodies with Bretherton constants of magnitude larger than unity.

In the case of planar motion, this solution completely describes the angular evolution of the
rigid torque-free object in flow. Owing to the linearity of Stokes equations, we may readily include
consideration of these angular dynamics into the yawing self-propelled swimmer system of the
previous section, leading to the nonlinear dynamical system

$$
\frac{dx}{dt} = V \cos \phi - \gamma y, \qquad (18a)
$$

$$
\frac{dy}{dt} = V \sin \phi, \tag{18b}
$$

$$
\frac{d\phi}{dt} = \frac{\gamma}{2}(1 - B \cos 2\phi) + \omega G'(\omega t). \quad (18c)
$$

Here we have included the effects of a background shear flow, with velocity **u** = −γ**Y**e_X at a point (X, Y) in the plane, driving the rotation and contributing to the translation. We emphasize that the motion of a passive object in a shear flow can be obtained by taking V = 0 in what follows. Similarly, the results of the previous section without a background flow can be recovered by taking γ = 0 in what follows. In full generality, the swimming direction need not coincide with a context-dependent definition of the swimmer orientation, though this is assumed in Eq. (18) to afford notational simplicity. Consistent with this, in what follows, we will assume that the orientation and direction of the swimmer are both captured by φ, noting that a constant difference between these quantities can be trivially included in the dynamics, mapping φ ↦→ φ + C for constant discrepancy C in any equations of rotational motion, for instance, with no significant consequence.

The general yawing dynamics of Sec. II B are included in Eq. (18c) via the derivative of the yawing function, denoted by $G'$. As above, we impose vanishing initial conditions for convenience, but note that these may be readily generalized. Without the yawing, the augmented dynamical system of Eq. (18) represents a minimal model of an active particle in a background shear flow, and has previously been used to model the average motion of complex, flagellated swimmers. One such example is in the context of the monoflagellated parasite *Leishmania mexicana* [31], where a similar model was used to investigate the long-time behaviors of a cell swimming near a boundary. Further, in Appendix B, we evidence the numerical validity of such a minimal model when seeking to study human spermatozoa by comparison with a direct numerical simulation of the angular microswimmer dynamics. We do note, however, that a rigorous theoretical exploration of the surprising general efficacy of such reduced representations remains an interesting avenue for future work. Acknowledging this, we henceforth assume the applicability of Eq. (16) in modeling the angular dynamics of a yaw-free microswimmer in a background shear flow.

In a similar manner to the previous section, we seek a multiple-scales solution to the yawing
system of Eq. (18). Under the same transformation, Eq. (18) becomes

$$
\frac{\partial x}{\partial t} + \omega \frac{\partial x}{\partial T} = V \cos \phi - \gamma y, \quad (19a)
$$

$$
\frac{\partial y}{\partial t} + \omega \frac{\partial y}{\partial T} = V \sin \phi, \qquad (19b)
$$

$$
\frac{\partial \phi}{\partial t} + \omega \frac{\partial \phi}{\partial T} = \frac{\gamma}{2} (1 - B \cos 2\phi) + \omega G'(T). \quad (19c)
$$

Here, noting the decoupling of the angular evolution equation from those of translation, we will
focus on presenting an analysis of the angular dynamics, with the leading-order translational motion
being constructed analogously to Sec. II B with only minimal modification. We pose an asymptotic
expansion for $\phi$ of the form

$$
\phi = \phi_0 + \frac{1}{\omega}\phi_1 + \dots \quad \text{as } \omega \to \infty, \qquad (20)
$$
and equate coefficients of $\omega$ in Eq. (19c). The leading-order angular problem is simply

$$ \frac{\partial \phi_0}{\partial T} = G'(T), \qquad (21) $$

which has general solution

$$ \phi_0 = \bar{\phi}_0(t) + G(T), \qquad (22) $$

precisely mirroring the form of $\phi$ considered in Sec. II B. However, the key difference is that the slow dynamics, $\bar{\phi}_0(t)$, are currently *undetermined* rather than being prescribed, necessitating an analysis of the higher-order problem to determine the angular evolution.

At the next order in $\omega$, the appropriate terms in Eq. (19c) are

$$ \frac{d\bar{\phi}_0}{dt} + \frac{\partial\phi_1}{\partial T} = f(\phi_0; B), \qquad (23) $$

recalling our shorthand $f(\phi; B)$ for the forcing of Jeffery's ODE (16) in the absence of yaw. Averaging Eq. (23) over a period of yawing and imposing periodicity in $T$, we arrive at the averaged ordinary differential equation

$$ \frac{d\bar{\phi}_0}{dt} = \langle f(\phi_0; B) \rangle = \frac{\gamma}{2}(1 - B\langle\cos 2\phi_0\rangle), \qquad (24) $$

which governs $\bar{\phi}_0$. We encountered the average of a slightly different cosine term in Sec. II B; following a similar analysis allows us to encode the key integral results that we require as the real and imaginary parts of the complex number

$$ z_A := \langle \exp 2iG \rangle = R_A \exp i\theta_A \in \mathbb{C}, \qquad (25) $$

similar to the constant $z_P$ from Eq. (14). Such an analysis allows us to rewrite Eq. (24) in terms of $z_A$ as

$$ \frac{d\bar{\phi}_0}{dt} = \frac{\gamma}{2}[1 - B_e \cos(2\bar{\phi}_0 + \theta_A)] = f(\bar{\phi}_0 + \theta_A/2; B_e), \qquad (26a) $$

where

$$ B_e := R_A B \qquad (26b) $$

is the effective Bretherton constant and $\theta_A/2$ is an effective phase shift. Comparing Eq. (26) to that derived by Jeffery, recounted here in Eq. (16), we see that this leading-order solution corresponds to the Jeffery's orbit of a biased angular variable $\psi := \bar{\phi}_0 + \theta_A/2$:

$$ \frac{d\psi}{dt} = f(\psi; B_e). \qquad (27) $$

Therefore, we have shown that any prescribed yawing function $G$ within a Jeffery's orbit will lead, on average, to Jeffery's orbit dynamics at leading order with an effective Bretherton constant $B_e$, where $|B_e| \le |B|$, and with a phase shift of $\theta_A/2$, which is nonzero if $\langle\sin 2G\rangle \ne 0$.

This has particular implications given the markedly different character of Jeffery's orbits associated with $|B| > 1$ and $|B| < 1$. Our results suggest that the introduction of yaw can transition an object from one regime to another. In particular, since yaw strictly reduces the effective magnitude of the Bretherton constant, yawing can elicit nontrivial periodic orbits from bodies that may otherwise have approached a steady state in the background flow. Such a transition is dependent only on the original Bretherton constant and the details of the yawing gait, and is independent of the yawing frequency at leading order. In special cases, for gaits that satisfy $\langle\sin 2G\rangle = 0$, there is no phase shift and we recover standard Jeffery's orbit dynamics, with the Bretherton constant being scaled by $\langle\cos 2G\rangle$. With reference to the discussion of Sec. II B, such a condition is automatically satisfied by yawing functions that satisfy either $G(T_o^* + T) = -G(T_o^* - T)$ or $G(T + \pi) = -G(T)$, for some
<figure><img src="image_7.png" /><figcaption>FIG. 4. Driving qualitative behavioral changes when $|B| > 1$. Considering an object with $|B| > 1$ and adopting simple sinusoidal yawing of amplitude $A$, we plot the numerically computed solutions of Eq. (18c) for $A \in \{0, 0.1\}$, shown as black and gray curves, respectively. In line with the presented asymptotic analysis, the evolution to steady state in the case of $A = 0$ is replaced by a periodic orbit when $A = 0.1$, noting that this choice of $A$ is sufficient to render $|B_e| < 1$ in this case. Here we have fixed $\omega = 10$ and taken $B = -1.002$, giving $B_e \approx -0.99$ for $A = 0.1$.</figcaption></figure>

fixed $T_o^*$. In Fig. 4 we exemplify this behavioral switch by numerically solving the exact angular evolution equation of Eq. (18c) for $G(T) = A \sin T$ and $B = -1.002$, with this Bretherton constant approximately corresponding to the ringlike bodies described by Singh et al. [41]. In line with the predictions of our analysis, the introduction of small-amplitude yawing gives rise to a qualitative change in the behavior of the object in flow, with evolution to a steady state being replaced by periodic dynamics in the presence of sufficient yaw.

In summary, the leading-order system of ODEs for the translational and angular yawing problem in the presence of a shear flow is given by

$$ \frac{d\bar{x}_0}{dt} = V_e \cos(\bar{\phi}_0 + \theta_p) - \gamma \bar{y}_0, \qquad (28a) $$

$$ \frac{d\bar{y}_0}{dt} = V_e \sin(\bar{\phi}_0 + \theta_p), \qquad (28b) $$

$$ \frac{d\bar{\phi}_0}{dt} = \frac{\gamma}{2} [1 - B_e \cos(2\bar{\phi}_0 + \theta_A)], \qquad (28c) $$

identical to the dynamics in the absence of yawing except for the introduction of (1) an effective propulsive speed, $V_e$; (2) a bias to the direction of propulsion, $\theta_p$; (3) an effective geometry via the effective Bretherton constant, $B_e$; and (4) a bias in the Jeffery's orbit dynamics, $\theta_A/2$. The angular biases for the translational and rotational dynamics vanish if the yawing gait satisfies $\langle \sin G \rangle = 0$ and $\langle \sin 2G \rangle = 0$, respectively. We note that these conditions are automatically satisfied for any antiperiodic $G$.

For the sinusoidal case, taking $G(T) = A \sin T$, it is straightforward to compute the following effective coefficients: $V_e = J_0(A)V$, $\theta_p = 0$, $B_e = J_0(2A)B$, and $\theta_A = 0$. Thus, the corresponding leading-order system is given by

$$ \frac{d\bar{x}_0}{dt} = V_e \cos \phi_J(t; B_e) - \gamma \bar{y}_0, \qquad (29a) $$

$$ \frac{d\bar{y}_0}{dt} = V_e \sin \phi_J(t; B_e), \qquad (29b) $$

$$ \phi_0(t) = \phi_J(t; B_e) + A \sin \omega t, \qquad (29c) $$
<figure><img src="image_8.png" /><figcaption>FIG. 5. The angular evolution of a yawing body. (a–c) The leading-order multiple-scales solution $\phi_0$ is shown over time alongside the leading-order solution averaged over a single oscillation $\bar{\phi}_0$, though is not visible at the resolution of (c). Also plotted is the solution to the unperturbed dynamics, $\phi_J$, from which we observe significant differences between the yawing and nonyawing systems. The periods of the slow dynamics are shown as dotted vertical lines in (b), with the yawing giving rise to a shorter period since $|B_e| < |B|$. (d–f) The envelope of the difference between the numerical solution to the full perturbed evolution equation, denoted by $\phi$, and the leading-order multiple-scales solution $\phi_0$. We note that this error is on the order of $10^{-3}$ across the timescales of each of these plots, validating the accuracy of the leading-order multiple-scales solution, with the accumulation of error over time only noticeable in the very longest comparison, which corresponds to approximately 16 full revolutions in the background flow, suggesting a long-term efficacy of the leading-order solution. Here, we have fixed $\gamma = 1$, $B = 0.5$, $A = 1$, and $\omega = 40\pi$.</figcaption></figure>

where $\phi_0$ corresponds to a Jeffery’s orbit with modified Bretherton constant and additive high-frequency yawing. In this context, a lower bound for the magnitude of the yawing needed to switch the sign of the Bretherton constant is given by half the first zero of $J_0$, giving the approximate bound $|A| > 1.2$ as a necessary condition for $B_e/B < 0$. For $|B| < 1$, such a change in sign would correspond to a switch between the behaviors exhibited by prolate and oblate spheroids, so that a yawing oblate spheroid would behave as a prolate spheroid in shear flow, for instance. However, we expect that this lower bound is too large to be of relevance to many biological swimmers.

We compare our asymptotic results to numerical solutions of the full system of Eq. (18), adopting sinusoidal yawing. We fix the parameters to be $\gamma = 1$, $B = 0.5$, $A = 1$, and $\omega = 40\pi$, the latter corresponding to a frequency of 20 oscillations per unit time. In Figs. 5(a)–(c), which display the dynamics on different timescales, respectively, we plot both the leading-order solution $\phi_0$, given by Eq. (29c), and the averaged solution $\bar{\phi}_0 = \phi_J(t; B_e)$. For comparison, we also plot the solution to the nonyawing problem, which is simply $\phi_J(t; B)$, i.e. a Jeffery’s orbit with differing Bretherton constant, with Fig. 5(b) illustrating the significant effects that rapid yawing can have on the dynamics over even a single long-time rotation in the background shear flow, including the period of rotation. At the resolution of these plots, the leading-order solution and the numerical solution of the yawing problem (18c) would be indistinguishable, a property captured by Figs. 5(d)–(f), which display the error $\phi - \phi_0$. We systematically evaluate this error for various $\omega$ in Appendix D, validating the asymptotic analysis and observing the expected $O(1/\omega)$ scaling. For the parameters of the exploration of Fig. 5, the error in the leading-order approximation is on the order of $10^{-3}$, a feature that is retained even over the longest timescale, which constitutes approximately 16 full revolutions in the background flow and 4000 yawing oscillations in this case. With the persistence of this low-magnitude error being somewhat unusual, we now move to investigate the higher-order problem, seeking to understand the source of this prolonged accuracy.
IV. HIGHER-ORDER DRIFT AND BIAS

In the previous sections, we constructed leading-order problems for the translational and rota-
tional motion of a yawing swimmer. When comparing our leading-order asymptotic solutions to
numerical results of the full dynamics, we found that our asymptotic solutions were significantly
more accurate than one might expect. In order to explore this observation, we now seek to calculate
the slow drift in the system at O(1/ω), the next asymptotic order in the system.

We consider the most general system of Eq. (19), with leading-order solutions constructed in
Sec. III. To calculate the higher-order solutions, we begin with the O(1) terms in Eq. (19), which
we recount as

$$
\frac{\partial x_0}{\partial t} + \frac{\partial x_1}{\partial T} = V \cos \phi_0 - \gamma y_0, \quad (30a)
$$

$$
\frac{\partial y_0}{\partial t} + \frac{\partial y_1}{\partial T} = V \sin \phi_0, \qquad (30b)
$$

$$
\frac{\partial \phi_0}{\partial t} + \frac{\partial \phi_1}{\partial T} = \frac{\gamma}{2}(1 - B \cos 2\phi_0). \quad (30c)
$$

Previously, we integrated these equations over a period in T to yield leading-order solvability
conditions and the corresponding solutions for x₀, y₀, and φ₀. In order to proceed to higher order,
we instead solve these equations to obtain

$$
x_1(t, T) = I_x(t, T) - \langle I_x \rangle + \bar{x}_1(t), \quad (31a)
$$

$$
y_1(t, T) = I_y(t, T) - \langle I_y \rangle + \bar{y}_1(t), \quad (31b)
$$

$$
\phi_1(t, T) = I_{\phi}(t, T) - \langle I_{\phi} \rangle + \bar{\phi}_1(t), \quad (31c)
$$

where

$$
I_x(t, T) := V \int_0^T [\cos \phi_0 - \langle \cos \phi_0 \rangle] ds - \gamma \int_0^T [y_0 - \langle y_0 \rangle] ds, \quad (32a)
$$

$$
I_y(t, T) := V \int_0^T [\sin \phi_0 - \langle \sin \phi_0 \rangle] ds, \quad (32b)
$$

$$
I_{\phi}(t, T) := \frac{\gamma B}{2} \int_{0}^{T} [\langle \cos 2\phi_{0} \rangle - \cos 2\phi_{0}] ds. \quad (32c)
$$

Here $\bar{x}_1$, $\bar{y}_1$, and $\bar{\phi}_1$ represent the second-order drifts, and are all as of yet undetermined functions of the slow time.

There are several ways to determine these higher-order drifts. A standard method would be
to repeat the approach employed earlier, collecting the O(1/ω) terms in Eq. (19) and deriving
solvability conditions by integrating over T and imposing periodicity. An alternative approach,
appllicable to this particular problem and affording significant notational simplicity over the previous
method, is to derive ODEs for the combined drift, at both leading and subleading orders, from the
outset. To do this, we integrate the full dynamical system of Eq. (19) over T, before considering any
particular asymptotic expansions, and impose T-periodicity, obtaining

$$
\frac{d\langle x \rangle}{dt} = V \langle \cos \phi \rangle - \gamma \langle y \rangle, \quad (33a)
$$

$$
\frac{d\langle y \rangle}{dt} = V \langle \sin \phi \rangle, \qquad (33b)
$$

$$
\frac{d\langle\phi\rangle}{dt} = \frac{\gamma}{2}(1 - B\langle\cos 2\phi\rangle). \quad (33c)
$$
These differential equations will govern the evolution of the full average dynamics, into which we can substitute the previously identified asymptotic expansions for x, y, and $\phi$. In particular, defining the convenient notation

$$ \bar{x}(t) = \bar{x}_0(t) + \frac{1}{\omega}\bar{x}_1(t), \quad \bar{y}(t) = \bar{y}_0(t) + \frac{1}{\omega}\bar{y}_1(t), \quad \bar{\phi}(t) = \bar{\phi}_0(t) + \frac{1}{\omega}\bar{\phi}_1(t) \quad (34) $$

for the slowly evolving drift of the dynamics, recalling the asymptotic solutions of the previous section, we write

$$ \begin{aligned} x &\sim \bar{x}(t) + \frac{1}{\omega}[I_x(t, T) - \langle I_x \rangle], & y &\sim \bar{y}(t) + \frac{1}{\omega}[I_y(t, T) - \langle I_y \rangle], \\ \phi &\sim \bar{\phi}(t) + G(T) + \frac{1}{\omega}[I_{\phi}(t, T) - \langle I_{\phi} \rangle], & & \end{aligned} \quad (35) $$

with asymptotic errors of $O(1/\omega^2)$. Inserting these asymptotic expansions into the differential equations governing the average dynamics, Eq. (33), and retaining terms of both $O(1)$ and $O(1/\omega)$, we find

$$ \frac{d\bar{x}}{dt} = V\tilde{R}_P \cos(\bar{\phi} + \tilde{\theta}_P) - \gamma\bar{y}, \quad (36a) $$

$$ \frac{d\bar{y}}{dt} = V\tilde{R}_P \sin(\bar{\phi} + \tilde{\theta}_P), \quad (36b) $$

$$ \frac{d\bar{\phi}}{dt} = \frac{\gamma}{2}[1 - B\tilde{R}_A \cos(2\bar{\phi} + \tilde{\theta}_A)], \quad (36c) $$

where

$$ \tilde{z}_P(t) := \left\langle \exp \left( i \left\{ G(T) + \frac{1}{\omega} [I_{\phi}(t, T) - \langle I_{\phi} \rangle] \right\} \right) \right\rangle = \tilde{R}_P \exp i \tilde{\theta}_P \in \mathbb{C}, \quad (37a) $$

$$ \tilde{z}_A(t) := \left\langle \exp \left( 2i \left\{ G(T) + \frac{1}{\omega} [I_{\phi}(t, T) - \langle I_{\phi} \rangle] \right\} \right) \right\rangle = \tilde{R}_A \exp i \tilde{\theta}_A \in \mathbb{C}. \quad (37b) $$

Consistent initial conditions for the ODEs of Eq. (36) can be determined by combining Eqs. (31), (34), and (35), yielding

$$ (\bar{x}(0), \bar{y}(0), \bar{\phi}(0)) = (x(0), y(0), \phi(0)) + \frac{1}{\omega} (\langle I_x \rangle, \langle I_y \rangle, \langle I_{\phi} \rangle)|_{t=0}. \quad (38) $$

Having arrived at these equations for the evolution of the average dynamics, we are able to conclude that the long-time drift, captured by $\bar{x}$, $\bar{y}$, and $\bar{\phi}$ and including terms of order $O(1/\omega)$, satisfies equations of precisely the same form as the leading-order system of Eq. (28). Hence, the conclusions of the previous section, that yawing modifies the average dynamics by modifying the swimmer shape, reducing the swimming speed, and potentially introducing orientational biases, hold even at $O(1/\omega)$. Moreover, any possible deviations away from this must appear at $O(1/\omega^2)$ or higher.

Further, we are able to more precisely quantify the effects of the yawing, with biases and reductions in both swimming and shape parameters now captured at both leading and subleading asymptotic orders. It is straightforward to show that $\tilde{z}_P \sim z_P$ and $\tilde{z}_A \sim z_A$ as $\omega \to \infty$, as we would expect. The more interesting question is the form of the next-order terms of these modified coefficients, that is, at $O(1/\omega)$. A detailed analysis of the next-order terms of Eq. (37) is presented in Appendix E, and we now summarize the key conclusions of this calculation.

First, we find that $\tilde{z}_A - z_A = o(1/\omega)$ if $G$ is antiperiodic. Thus, the average leading-order solution for the angular dynamics $\bar{\phi}_0$ obeys the same ODE as the second-order-accurate solution $\bar{\phi}$, and these solutions are equal subject to imposing the higher-order initial condition of Eq. (38) on $\bar{\phi}_0$. We further note that even if one were to simply employ leading-order initial conditions, the $O(1/\omega)$
error introduced would be periodic and thus bounded in time, so that any long-term growth in error would still occur at $O(1/\omega^2)$ or higher, in accordance with the numerical explorations shown in Fig. 5. Additionally, recalling the conclusions of Sec. III, the imposed antiperiodicity condition is also sufficient to give rise to unbiased angular dynamics, with $\theta_A = 0$ and $\dot{\theta}_A = o(1/\omega)$.

Second, now with reference to the translational dynamics, we find that antiperiodicity alone is not sufficient to yield higher-order agreement between the leading-order and full solutions. However, imposing the additional requirement of $G$ being even about some $T_e^*$ is sufficient to align the average leading-order solution and the higher-order averaged solution, again subject to the prescription of the corrected initial conditions of Eq. (38). Indeed, if the initial condition of the rotational dynamics is specified to next order, but the initial condition of the translational dynamics is only given at leading order, then the resulting error in the translational dynamics at $O(1/\omega)$ will only remain constant in time if $\gamma = 0$, i.e. if there is no background flow. This is a direct result of the advective $\gamma y$ term in Eq. (18a), which gives rise to accumulating errors in $x$ if there is error present in $y$.

Finally, in the absence of flow and with a prescribed orientation, such as in the case of pure translation in Sec. II, we may readily conclude from Eq. (33) that there is no error in the dynamical system at any order, with the averages of $\cos G$, $\sin G$, and $\cos 2G$ being known exactly. These conclusions hold at all asymptotic orders and are not reliant on any assumptions on the yawing function, as described in the extended context of three-dimensional motion in Appendix A. However, due to the potential appearance of additional slower timescales from using the method of multiple scales, we refrain from concluding that the leading-order solution for pure translation, subject to the second-order initial conditions of Eq. (38), is accurate to all orders.

V. DISCUSSION

In this work, we have theoretically and numerically explored the effects of high-frequency yawing oscillations on both objects in Stokes flow and the motion of active particles, from the dynamics of a self-propelled particle to the canonical Jeffery's orbit. In this latter example, employing the asymptotic method of multiple scales as we have done throughout, we have found, perhaps surprisingly, that the slow angular dynamics in the presence of rapid yawing are once again governed by Jeffery's equation in the plane. This supports, though does not justify, the use of Jeffery's orbits for the modeling of yawing bodies, as introducing oscillations to a nonyawing body undergoing a Jeffery's orbit preserves the structure of the dynamical system, a feature that also holds for the translational dynamics. In particular, the general effect of the yawing on the slow dynamics is to (1) modify the effective shape of the particle, (2) introduce a phase shift, and (3) strictly reduce the magnitude of the Bretherton constant and any nonzero swimming speed. Notably, the precise frequency of the rapid oscillation does not affect the slow dynamics at leading order. The potential phase shift in the dynamics can be seen to vanish if the yawing is antiperiodic, with the effects of rapid oscillations otherwise motivating an alternative definition of body orientation that is applicable to yawing objects and explicitly accounts for this phase shift.

In the context of spheroids, which have Bretherton constant of magnitude less than unity, we have seen that high-frequency yawing gives rise to an effective shape that is less elongated, with the extent of this change being governed by the quantity $R_A$, defined in Eq. (25). For bodies with Bretherton constant of magnitude greater than one, which can involve fairly unusual shapes (as illustrated by Bretherton [33]), the effects of yawing can result in significant qualitative differences. In particular, rapid yawing can alter the long-term behavior of the body, since the solutions of Jeffery's ordinary differential equation are of markedly different character for $|B| < 1$ and $|B| > 1$; they are nontrivially periodic in the former and approach a constant in the latter [41,42]. Hence, rapid oscillations can drive qualitative changes in the dynamics of an object in shear flow, and otherwise serve to reduce their period of oscillation. This investigation into a perturbation of a Jeffery's orbit complements the large body of research that utilizes or analyzes the motion of particles in flow, with particular reference to the generalizing studies of Bretherton [33] and Ishimoto [36], among others [43–45].
The condition of *antiperiodic* yawing has been a key property several times during this manuscript. We find that antiperiodicity removes orientational biases from translational and angular motion and bounds the dynamics at next order, resulting in a long-time accuracy of the leading-order dynamics derived in this study, a nonstandard property in general multiple-scales approaches. As noted at the conclusion of Sec. II B, this condition arises from stereotypical beating patterns of canonical swimmers, such as the amplitude-modulated traveling waves of displacement occurring in spermatozoa, and was indeed found to be satisfied in the computational exploration of Appendix B for a model spermatozoon. Hence, while not a necessary condition for removing orientational bias or affording additional accuracy to our approximations, this constraint of antiperiodicity appears to feature in numerous swimming gaits. In particular, as discussed in Appendix C, antiperiodicity of yawing arises from bodies with symmetric swimming gaits, though only if the swimmer orientation is defined such that the average of the yawing function vanishes. In turn, this suggests that there is some notion of a natural definition of swimmer orientation, dependent on the geometry of the swimmer and the details of the swimming gait. With a universal concept of swimmer orientation lacking, leading to definitions that vary from one study to the next, seeking to satisfy this constraint outside the context of the models considered in this work may serve to reduce the freedom in defining swimmer orientation while affording desirable properties, as it has done here, and will be a subject of future exploration.

In our consideration of a model spermatozoon, the computed yawing function appeared to be approximately sinusoidal, similar to that studied in Sec. II A. Serendipitously, the yawing function $G(T) = A \sin T$ is seen to satisfy the conditions of antiperiodicity and evenness. Hence, the analysis of Sec. IV suggests that the effects of rapid yawing dynamics are well captured by the leading-order solution on long timescales, and that the approximations derived in this study may be relied upon for significant durations. This simple yawing form also gave rise to unbiased effective dynamics, preserving the direction of swimming and the swimmer's orientation in the context of a Jeffery's orbit, with the sought parity arguments in the general case stemming from our analysis of this minimal example. Further, in this case, the modifications to the swimming speed in the plane and to the Bretherton constant are given analytically in terms of the Bessel function $J_0$, with conversion between shape and swimming parameters in general also being straightforward, given the details of the yawing function. This example additionally highlights the potential for rapid oscillations to significantly alter motion, with certain amplitudes of sinusoidal yawing giving rise to vanishing displacement or even a reversed swimming direction, though we note that large yawing amplitudes are required in order to achieve such drastic changes.

Throughout this work, we have considered simple, phenomenological models of microswimmers in common use, abstracted away from the microscale details of an individual swimmer. Starting with a study of simple planar motion, rendered nontrivial by the presence of rapid yawing, we proceeded to consider scenarios of increasing complexity, incorporating the effects of geometry in flow in Sec. III and extending our treatment of translation to three dimensions in Appendix A. This latter formulation allows for the exploration of three-dimensional yawing, such as the rapid oscillations generated by a model monotrichous bacterium, which would be axisymmetric and, thus, nonyawing were it not for the finite length of its flagellum. Subject to the applicability of a simple model for swimmer motion, this set-up also allows for the consideration of three-dimensional beating of spermatozoa, for instance, a mode of flagellar motion that is frequently exhibited by these canonical swimmers. Returning to planar swimming, as was the focus of the main text and with particular reference to angular dynamics, we have performed high-resolution computational simulations of a model planar-beating spermatozoon in a shear flow in Appendix B, evidencing the applicability of an effective Jeffery's orbit model for their rotation in this flow. Such a model has already been employed in other contexts, such as to *Leishmania mexicana* by Walker *et al.* [31], though a precise and analytical justification of the suitability of such models remains lacking, and will form the subject of future enquiry. Commensurate with the reports of Walker *et al.* [31], we remark that the planar Jeffery's orbits explored in Sec. III are robust to perturbations out of the plane of the shear flow, with numerical exploration of the three-dimensional dynamics also suggesting a neutral
stability of these perturbations, though this requires analytic assessment in future studies of three-dimensional Jeffery's orbits.

In summary, we have considered the motion of swimmers and active particles that are subject to rapid oscillations in their orientation. Through asymptotic analysis and supporting numerical explorations, we have quantified how swimming speeds, swimming directions, and even the effective geometries of swimmers and microparticles can be influenced by the presence of high-frequency oscillation, both at leading and subleading asymptotic orders. In doing so, we have identified sufficient conditions for oscillations to not introduce directional biases into the emergent dynamics, with the recurring condition of antiperiodic yawing being an observed feature of the motion of a canonical and well-studied swimmer. Finally, in the context of Stokesian fluids, we have seen how the Jeffery's orbits of ellipsoidal and symmetric particles are altered by kinematic oscillations, with the Jeffery's orbit structure of the dynamical system being retained for sinusoidal yawing up to two asymptotic orders, and have highlighted the potential for rapid yawing to enact qualitative long-term changes to the dynamics of bodies in flow.

The data used and generated in this paper are freely available from the University of Oxford Research Archive [46].

## ACKNOWLEDGMENTS

The authors are grateful to Prof. E. Lauga, University of Cambridge, for interesting and motivating discussions on the separated timescales of motion found in many microswimming problems. B.J.W. is supported by the U.K. Engineering and Physical Sciences Research Council (EPSRC), Grant No. EP/R513295/1. K.I. acknowledges JSPS-KAKENHI for Young Researchers (Grant No. 18K13456) and JST, PRESTO, Japan (Grant No. JPMJPR1921). C.M. is a JSPS International Research Fellow (PE20021). K.I. and C.M. were partially supported by the Research Institute for Mathematical Sciences, an International Joint Usage/Research Center located at Kyoto University.

## APPENDIX A: TRANSLATION IN THREE DIMENSIONS

In Sec. II, we considered the motion of a self-propelled particle in the plane under the influence of prescribed yawing. Here, we extend this leading-order analysis to three dimensions and higher orders, with prescribed average orientation and general yawing that need not be confined to a plane. To do so, we will employ a standard spherical coordinate system to describe the orientation of the swimmer, with $\zeta \in [0, \pi]$ denoting the polar angle and $\phi \in [0, 2\pi)$ denoting the azimuthal angle, consistent with the notation of the main text. Explicitly, denoting the third spatial coordinate of the swimmer by $z(t)$, the dynamics are governed by the extended dynamical system

$$ \frac{dx}{dt} = V \sin \zeta \cos \phi, \qquad (A1a) $$

$$ \frac{dy}{dt} = V \sin \zeta \sin \phi, \qquad (A1b) $$

$$ \frac{dz}{dt} = V \cos \zeta, \qquad (A1c) $$

for fixed speed $V$, noting that $\zeta = \pi/2$ corresponds to the planar system of the main text. As in the general analysis of Sec. II B, we will prescribe the evolution of the swimmer orientation, writing $\zeta = \pi/2 + G_{\zeta}(T)$ and $\phi = G_{\phi}(T)$ without loss of generality, for known functions $G_{\zeta}$ and $G_{\phi}$ of the fast yawing timescale $T$. As before, these functions completely describe the evolution of the instantaneous swimming direction and may be deduced from simulation or measured from videomicroscopy, for instance. We will assume that the swimmer angular dynamics is drift free, whereby, without loss of generality, both $G_{\zeta}$ and $G_{\phi}$ have zero mean over a period in $T$, consistent with choosing a coordinate system in which the average swimmer direction is along the X-axis, in
line with the forms of $\zeta$ and $\phi$ as written above. Of note, subject to a simple time-dependent change of coordinates, the analysis of this Appendix may be readily applied if the average heading of the swimmer varies on the slow timescale $t$, though the details of this are notationally cumbersome and are omitted here. Given the reduced setting with a stationary average heading, the dynamical system reads

$$ \frac{dx}{dt} = V \cos G_{\zeta} \cos G_{\phi}, \qquad (A2a) $$

$$ \frac{dy}{dt} = V \cos G_{\zeta} \sin G_{\phi}, \qquad (A2b) $$

$$ \frac{dz}{dt} = -V \sin G_{\zeta}, \qquad (A2c) $$

from which it is clear that $G_{\zeta} = G_{\phi} = 0$, i.e., motion with no yawing, corresponds to swimming along the X-axis. As in the two-dimensional case, we form the corresponding multiple-scales PDEs and seek expansions for the spatial coordinates in inverse powers of $\omega$. In this case, as in Sec. II, the leading-order solutions are independent of the fast timescale $T$, so that we may write, e.g., $x_0 = \bar{x}_0(t)$. These unknown functions correspond to the average motion of the swimmer, with the solvability conditions of the second asymptotic order of the PDEs needed to determine them. These read

$$ \frac{d\bar{x}_0}{dt} + \frac{\partial x_1}{\partial T} = V \cos G_{\zeta} \cos G_{\phi}, \qquad (A3a) $$

$$ \frac{d\bar{y}_0}{dt} + \frac{\partial y_1}{\partial T} = V \cos G_{\zeta} \sin G_{\phi}, \qquad (A3b) $$

$$ \frac{d\bar{z}_0}{dt} + \frac{\partial z_1}{\partial T} = -V \sin G_{\zeta}. \qquad (A3c) $$

Seeking periodicity in $T$ and averaging over the fast timescale, we arrive at the leading-order ODEs

$$ \frac{d\bar{x}_0}{dt} = V \langle \cos G_{\zeta} \cos G_{\phi} \rangle, \qquad (A4a) $$

$$ \frac{d\bar{y}_0}{dt} = V \langle \cos G_{\zeta} \sin G_{\phi} \rangle, \qquad (A4b) $$

$$ \frac{d\bar{z}_0}{dt} = -V \langle \sin G_{\zeta} \rangle. \qquad (A4c) $$

Hence, in order for there to be no directional bias away from the yaw-free heading, which in this coordinate system is along the X-axis, we require both

$$ \langle \cos G_{\zeta} \sin G_{\phi} \rangle = 0 \quad \text{and} \quad \langle \sin G_{\zeta} \rangle = 0. \qquad (A5) $$

These conditions are readily found to be satisfied if $G_{\zeta}$ and $G_{\phi}$ are both antiperiodic, which naturally extends the conditions of the planar case. As noted in the planar scenario and discussed in more detail in Appendix C, these conditions are automatically satisfied by swimmers with antiperiodic gaits given appropriate definitions of the swimmer orientation, now interpreted to be in three dimensions and including the entire geometry of the body. Should these conditions hold, the effective swimming speed of the self-propelled body is given by $V_e = V \langle \cos G_{\zeta} \cos G_{\phi} \rangle$, lower in magnitude than the original speed $V$, as is intuitive. For instance, if both $G_{\zeta}$ and $G_{\phi}$ were sinusoidal, such as $G_{\zeta}(T) = G_{\phi}(T) = A \sin(T)$, then the above conditions are satisfied and the effective velocity is $V_e = V(J_0(2A) + 1)/2$.
With the angular dynamics prescribed, the solution at the next order can be easily obtained, as in the case for planar motion. Given the average dynamics of Eq. (A4), we solve Eq. (A3) to give

$$x_1 = I_x(T) - \langle I_x \rangle + \bar{x}_1(t), \qquad (A6a)$$

$$y_1 = I_y(T) - \langle I_y \rangle + \bar{y}_1(t), \qquad (A6b)$$

$$z_1 = I_z(T) - \langle I_z \rangle + \bar{z}_1(t), \qquad (A6c)$$

here defining

$$I_x = V \int_0^T [\cos G_\zeta(s) \cos G_\phi(s) - \langle \cos G_\zeta \cos G_\phi \rangle] ds, \qquad (A7a)$$

$$I_y = V \int_0^T [\cos G_\zeta(s) \sin G_\phi(s) - \langle \cos G_\zeta \sin G_\phi \rangle] ds, \qquad (A7b)$$

$$I_z = -V \int_0^T [\sin G_\zeta(s) - \langle \sin G_\zeta \rangle] ds. \qquad (A7c)$$

Considering the next-order multiple-scales PDEs to determine the unknown averages $\bar{x}_1$, $\bar{y}_1$, and $\bar{z}_1$ and again averaging over $T$, we simply obtain

$$\frac{d\bar{x}_1}{dt} = 0, \quad \frac{d\bar{y}_1}{dt} = 0, \quad \frac{d\bar{z}_1}{dt} = 0, \qquad (A8)$$

noting that the forcing of the original system of differential equations contributes only at $O(1)$ as the swimmer orientation is prescribed, in contrast to the analysis of Sec. IV. Hence, the second-order correction to the average swimmer dynamics in three dimensions is simply a constant, so that drift away from this solution appears only in the higher-order terms. These constants are determined by enforcing the initial conditions of $x_1(0, 0) = y_1(0, 0) = z_1(0, 0) = 0$ in Eq. (A6), giving

$$\bar{x}_1 = \langle I_x \rangle_{t=0}, \quad \bar{y}_1 = \langle I_y \rangle_{t=0}, \quad \bar{z}_1 = \langle I_z \rangle_{t=0}. \qquad (A9)$$

## APPENDIX B: AN EXAMPLE MICROSWIMMER

Though established in *Leishmania mexicana* [31], the efficacy of Jeffery's orbit approximations for microswimmers in shear flow is far from being universally established. Here we numerically evidence the validity of utilizing a Jeffery's orbit model for the study of human spermatozoa, adopting the morphology of Gallagher and Smith [47] with the model beat pattern derived from the work of Dresdner and Katz [3], as shown in Fig. 6(a). Fixing a dimensionless shear rate of $\gamma = 1$, we modify the regularized-Stokeslet implementation of Gallagher and Smith [47] to accommodate the background shear flow, simulating the motion of a beating spermatozoon as it yaws due to the high-frequency flagellar beat. For full details of the methodology and swimmer morphology, we refer the interested reader to the publication of Gallagher and Smith [47] and the accompanying open-source implementation. In Fig. 6(a) we illustrate the spermatozoan beating gait by plotting snapshots of its flagellum relative to the swimmer body, the latter illustrated by the nodes of a computational mesh and with axis scales shown equal. Here we have made use of a standard sinusoidal traveling wave for the beating gait, given in terms of a flagellar arclength parameter $s \ge 0$ and swimmer-fixed reference coordinates $(\hat{x}, \hat{y})$ as

$$\hat{x}(s, T) = s, \quad \hat{y}(s, T) = (\alpha + \beta s)[\sin(ks - T) - \sin(-T)] \qquad (B1)$$

for wavelength $k = 2\pi$ and constants $\alpha = 0.1087$ and $\beta = 0.0543$ [3]. From simulation of the swimmer with this gait, temporarily in the absence of background flow and with reference to the notation of the main text, we compute the quantities $z_p = 0.9430 - 2 \times 10^{-4}i$ and $z_A = 0.78151 - 4 \times 10^{-4}i$, so that this gait approximately satisfies the conditions of Sec. II B and Sec. IV. Hence, upon assuming the applicability of the simple models explored in this work, this typical beating
<figure><img src="image_9.png" /></figure>

FIG. 6. Orbits of a model spermatozoon in shear flow. (a) A computational representation of the model spermatozoon and its beat, shown at multiple points over the period in a swimmer-fixed reference frame. (b) Plotted for $\omega \in \{0.5, 1, 20\}$, we showcase the evolution of the swimmer orientation $\phi$ as a function of time in a shear flow with $\gamma = 1$. We observe a remarkable agreement between the solutions, which lie approximately within an envelope of the temporal moving average, except the lowest frequency case, which is expected from our analysis. Fitting a Jeffery's orbit to the highest frequency case, finding $B_e \approx 0.93573$, we plot this solution as a blue dashed curve, demonstrating the remarkable validity of a Jeffery's orbit model. Inset: We plot the derived yawing function $G(T)$, from which we compute $z_p = 0.9430 - 2 \times 10^{-4}i$ and $z_A = 0.78151 - 4 \times 10^{-4}i$, so that both $\theta_p$ and $\theta_A$ are approximately zero and the slow-time motion is unbiased. Further, noting the sinusoidal form of $G$, this yawing function satisfies the conditions of Sec. IV, so that the long-time accumulation of error occurs at $O(1/\omega^2)$ or higher.

pattern would lead to negligible directional bias and is free from accumulating error until at least $O(1/\omega^2)$.

In Fig. 6(b) we plot the evolution of the sperm orientation, defined as the angle $\phi$ between a body-fixed reference frame and a fixed laboratory frame, with the origin of the body frame being at the center of the sperm head, for multiple values of the angular beat frequency $\omega$. The angle $\phi$ is defined such that $\langle G \rangle = 0$. In all but the lowest-frequency case, and concordant with the presented multiple-scales analysis, the angular evolution of the swimmer can be seen to be well approximated by a fitted Jeffery's orbit, shown as a blue dashed curve with $B_e \approx 0.93573$, with the trajectories lying within an envelope of the fitted solution. Further, as the beat frequency of spermatozoa often lies towards the upper range of values considered in this plot, such a Jeffery's orbit approximation is expected to be well suited for modeling the orientation of a beating spermatozoon in a shear flow, though further theoretical justification is required in order to establish this result in more generality.

## APPENDIX C: ANTIPERIODIC SWIMMERS AND ANTIPERIODIC YAW

Here, we detail how swimmers with particular symmetry give rise to antiperiodic yawing in the absence of external influences, focusing on the case where swimmer orientation can be described by a single angle expressed only as a function of the fast time $T$. Consider a swimmer whose geometrical configuration is antiperiodic, by which we mean that the configuration at time $T + \pi$ is a precise reflection (in a fixed plane) of its configuration at time $T$ when viewed in a swimmer-fixed reference frame, such as that shown in Fig. 6(a). Consider also an angle $\phi$ that measures the orientation of the swimmer-fixed frame with respect to a fixed laboratory frame, noting that there are many possible definitions of $\phi$, of which that from Fig. 1(b) is an example, and $\phi \mapsto \phi + D$ is another for any constant $D$. The rate of change of $\phi$ with time, written as $d\phi/dT$, is easily seen to be antiperiodic, noting the boundary-value nature of the Stokes swimming problem and the assumed antiperiodicity of the boundary conditions and swimmer geometry. This gives rise to a family of
derived yawing functions

$$G(T) = \underbrace{\int_{0}^{T} \frac{d\phi}{dT} ds}_{G_{1}(T)} + C, \qquad (C1)$$

where this freedom in choosing the constant of integration C is equivalent to shifting $\phi$ by a constant, capturing the many possible definitions of the swimmer orientation. Exploiting the antiperiodicity of the derivative, we can write the yawing functions as

$$G(T) = \begin{cases} G_1(T) + C, & T \in [0, \pi), \\ G_1(\pi) + C - G_1(T - \pi), & T \in [\pi, 2\pi), \end{cases} \qquad (C2)$$

recalling the $2\pi$-periodicity of G and noting that this satisfies $dG/dT = d\phi/dT$ for all T. From this, it is easy to see that we may ensure that G is antiperiodic by imposing the simple condition $G_1(\pi) + 2C = 0$, so that $G(T + \pi) + G(T) = 0$. Hence, for a swimmer with antiperiodic geometry over the course of its swimming gait, we may readily define its orientation, via C, such that its yawing function is antiperiodic. Indeed, this definition is equivalent to choosing $\phi$ such that $\langle G \rangle = 0$, as can be seen directly by computing the average of $G(T)$ as written in Eq. (C2), finding that this average vanishes if and only if $G_1(\pi) + 2C = 0$.

Of note, the above analysis extends to the antiderivatives of any antiperiodic function, with the antiderivative of an antiperiodic function being antiperiodic if and only if the average of the antiderivative is zero over a full period, such as that shown in Fig. 3(b).

## APPENDIX D: NUMERICALLY EXPLORING THE EFFECTS OF FREQUENCY ON ERROR

To support the multiple-scales analysis for the planar Jeffery's orbit in the presence of sinusoidal yawing, we numerically solve the perturbed system of Eq. (18c) and compare the solution to the leading-order multiple-scales approximation of Eq. (26). Computation was performed using the in-built ordinary differential equation solver ode45 in MATLAB [40] with error tolerances of $10^{-9}$. In Fig. 7, we plot the base ten logarithm of the infinity norm error between the full solution $\phi$ and the leading-order multiple-scales solution $\phi_0$, denoted by $\log_{10} \|\phi - \phi_0\|_{\infty}$, against $\log_{10} \omega$, comparing the solutions over ten periods of the unperturbed Jeffery's orbit. Immediately evident is the approximate $1/\omega$ scaling of the error with the frequency of oscillation, in line with the analysis of Sec. III and confirming the predicted accuracy of the multiple-scales approximation over this timescale. In particular, this order of error reflects the difference between the initial conditions of the leading-order and higher-order approximation, which scale with $1/\omega$.

## APPENDIX E: DETAILS OF HIGHER-ORDER DRIFT

In the analysis of Sec. IV, we remarked that agreement between the average leading-order dynamics of Eq. (28) and the higher-order averaged dynamics of Eq. (36) required that $\tilde{z}_p - z_p = o(1/\omega)$ and $\tilde{z}_A - z_A = o(1/\omega)$. Given the leading-order agreement, it remains to show that the subleading terms of $\tilde{z}_p$ and $\tilde{z}_A$ vanish, subject to constraints on the yawing function G. Expanding out the exponential integrands of the averages that define these complex numbers in Eq. (37) and retaining only the relevant terms, we have

$$\tilde{z}_p - z_p = -\frac{1}{\omega} \langle [I_{\phi} - \langle I_{\phi} \rangle] \sin G - i [I_{\phi} - \langle I_{\phi} \rangle] \cos G \rangle + o(1/\omega), \qquad (E1a)$$

$$\tilde{z}_A - z_A = -\frac{2}{\omega} \langle [I_{\phi} - \langle I_{\phi} \rangle] \sin 2G - i [I_{\phi} - \langle I_{\phi} \rangle] \cos 2G \rangle + o(1/\omega). \qquad (E1b)$$
<figure><img src="image_10.png" /><figcaption>FIG. 7. Error in the multiple-scales solution for planar Jeffery's orbits with sinusoidal yawing. For a range of values of $\omega$, we numerically solve the perturbed system (18c), comparing the solution with that of the multiple-scales analysis of Eq. (26). We plot, as black crosses, the base ten logarithm of the infinity norm error in the orientation as a function of $\log_{10} \omega$, clearly observing the predicted $1/\omega$ behavior as $\omega \to \infty$. Here the dynamics are simulated over ten periods of the unperturbed Jeffery's orbit, taking $B = 0.5$, $\gamma = 1$, and $A = 0.1$. A linear fit is shown as a gray, dashed line.</figcaption></figure>

Hence, we must compute the averages of terms that involve the quantity $I_{\phi} - \langle I_{\phi} \rangle$ multiplied by trigonometric functions of $G$ and $2G$. To do so, it is helpful to recall the definition of $I_{\phi}$ from Eq. (32), which can be expanded as

$$I_{\phi}(t, T) = \frac{\gamma B}{2} [\cos 2\bar{\phi}_{0}(t)I_{c}(T) - \sin 2\bar{\phi}_{0}(t)I_{s}(T)], \quad (E2a)$$

where

$$I_{c}(T) := \int_{0}^{T} [\langle \cos 2G \rangle - \cos 2G(s)] ds, \quad I_{s}(T) := \int_{0}^{T} [\langle \sin 2G \rangle - \sin 2G(s)] ds, \quad (E2b)$$

utilizing the general expression for the leading-order angular solution $\phi_0$ of Eq. (22). Thus, $I_{\phi}$ is a linear combination of $I_c$ and $I_s$, with coefficients that are independent of $T$, so that $I_{\phi} - \langle I_{\phi} \rangle$ is a linear combination of $I_c - \langle I_{c} \rangle$ and $I_s - \langle I_{s} \rangle$. Hence, the terms of Eq. (E1) will vanish if we can demonstrate that the following eight terms each have zero average over a period in $T$:

$$[I_c - \langle I_c \rangle] \sin G, \quad [I_c - \langle I_c \rangle] \cos G, \quad [I_s - \langle I_s \rangle] \sin G, \quad [I_s - \langle I_s \rangle] \cos G, \quad (E3a)$$

$$[I_c - \langle I_c \rangle] \sin 2G, \quad [I_c - \langle I_c \rangle] \cos 2G, \quad [I_s - \langle I_s \rangle] \sin 2G, \quad [I_s - \langle I_s \rangle] \cos 2G. \quad (E3b)$$

Here, the first row of terms pertains to the translational motion, while the second row concerns the angular evolution. We will show that $[I_c - \langle I_c \rangle] \cos 2G$ and $[I_s - \langle I_s \rangle] \sin 2G$ have identically zero averages in all cases, while the remaining terms have vanishing average due to parity arguments, following the imposition of two appropriate symmetry conditions on $G$.

We begin with the two terms that do not require additional conditions on $G$. Instead of $[I_c - \langle I_c \rangle] \cos 2G$, consider the modified expression $[I_c - \langle I_c \rangle][\langle \cos 2G \rangle - \cos 2G]$, which has the same average up to sign as the original by virtue of $\langle I_c - \langle I_c \rangle \rangle = 0$. Recalling the definition of $I_c$ and
defining $F(T) = I_c - \langle I_c \rangle$, so that $F'(T) = \langle \cos 2G \rangle - \cos 2G(T)$, we see that we can write

$$ \langle [I_c - \langle I_c \rangle] [\langle \cos 2G \rangle - \cos 2G] \rangle = \frac{1}{2\pi} \int_0^{2\pi} FF' dT = \frac{1}{4\pi} [F^2]_0^{2\pi} = 0 \quad (E4) $$

as $F$ is periodic, a property inherited from $G$. Hence, the average of $[I_c - \langle I_c \rangle] \cos 2G$ is zero. Analogous reasoning applies to $[I_s - \langle I_s \rangle] \sin 2G$, so that it also has zero mean, with no conditions on $G$ necessary.

Next, we impose an antiperiodicity condition on $G$; precisely the condition that was demonstrated to remove the leading-order orientational biases in the translational and angular dynamics in Sec. II B and Sec. III, respectively. This antiperiodicity immediately means that $\sin G$ and $\sin 2G$ are also antiperiodic, while $\cos G$ and $\cos 2G$ are periodic with period $\pi$. From this, we will now seek to establish the antiperiodicity of $I_s - \langle I_s \rangle$ and the $\pi$-periodicity of $I_c - \langle I_c \rangle$, which will lead to the vanishing of all but two of the remaining terms in Eq. (E3) by parity.

In fact, the $\pi$-periodicity of $I_c$, and therefore $I_c - \langle I_c \rangle$, follows immediately from its integrand $\langle \cos 2G \rangle - \cos 2G$ being $\pi$-periodic and having zero mean over a period in $T$. The $\pi$-antiperiodicity of $I_s - \langle I_s \rangle$ follows from the concluding remarks of Appendix C, with its integrand $\langle \sin 2G \rangle - \sin 2G$ being $\pi$-antiperiodic, and having zero mean over a $2\pi$-period in $T$.

The above parity arguments, which assume only that $G$ is antiperiodic, are sufficient to render the average of all of the remaining terms of Eq. (E3) zero, with the exception of the translational terms $[I_c - \langle I_c \rangle] \cos G$ and $[I_s - \langle I_s \rangle] \sin G$. Hence, we can conclude that the antiperiodicity of $G$ is sufficient to ensure that the average leading-order angular dynamics are equal to the higher-order averaged angular dynamics, with any deviation occurring at $O(1/\omega^2)$ or higher.

An additional condition is required in order to draw the same conclusion of the translational dynamics. Here, once again seeking to make use of a parity argument, we impose the additional constraint that $G$ is even about some $T_e^*$, which we note is satisfied by sinusoidal yawing. This implies that both $\cos G$ and $\sin G$ are even about the same $T_e^*$, while also ensuring that $I_c - \langle I_c \rangle$ and $I_s - \langle I_s \rangle$ are odd about $T_e^*$. These latter conclusions require some exposition. Consider $I_c - \langle I_c \rangle$, and write

$$ I_c(T_e^* + T) - \langle I_c \rangle = \left( \int_0^{T_e^*} + \int_{T_e^*}^{T_e^*+T} \right) [\langle \cos 2G \rangle - \cos 2G(s)] ds - \langle I_c \rangle. \quad (E5) $$

The second integral in the above is odd in $T$, with both $\langle \cos 2G \rangle$ and $\cos 2G$ being even about $T_e^*$. Hence, averaging over $T$, we have

$$ 0 = \int_0^{T_e^*} [\langle \cos 2G \rangle - \cos 2G(s)] ds - \langle I_c \rangle. \quad (E6) $$

Inserting this relation into Eq. (E5) gives

$$ I_c(T_e^* + T) - \langle I_c \rangle = \int_{T_e^*}^{T_e^*+T} [\langle \cos 2G \rangle - \cos 2G(s)] ds, \quad (E7) $$

which we have already noted is odd in $T$. Hence, $I_c - \langle I_c \rangle$ is odd about $T_e^*$. An entirely analogous analysis yields that $I_s - \langle I_s \rangle$ is odd about $T_e^*$.

Hence, by parity, the remaining terms of Eq. (E3) have vanishing average, and we can conclude that, if $G$ is both antiperiodic and even, the average leading-order dynamics of Eq. (28) are equal to the higher-order averaged dynamics of Eq. (36), with corrections arising at $O(1/\omega^2)$ or higher.

<footnote>[1] M. Gazzola, M. Argentina, and L. Mahadevan, Scaling macroscopic aquatic locomotion, Nat. Phys. 10, 758 (2014).</footnote>
[2] C. Moreau, K. Ishimoto, E. A. Gaffney, and B. J. Walker, Control and controllability of microswimmers by a shearing flow, *R. Soc. Open Sci.* **8**, 211141 (2021).

[3] R. D. Dresdner and D. F. Katz, Relationships of mammalian sperm motility and morphology to hydrodynamic aspects of cell function, *Biol. Reprod.* **25**, 920 (1981).

[4] D. J. Smith, E. A. Gaffney, H. Gadêlha, N. Kapur, and J. C. Kirkman-Brown, Bend propagation in the flagella of migrating human sperm, and its modulation by viscosity, *Cell Motil. Cytoskeleton* **66**, 220 (2009).

[5] T. Hyakutake, K. Sato, and K. Sugita, Study of bovine sperm motility in shear-thinning viscoelastic fluids, *J. Biomech.* **88**, 130 (2019).

[6] M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gompper, and H. Rieger, Computational models for active matter, *Nat. Rev. Phys.* **2**, 181 (2020).

[7] O. Chepizhko and T. Franosch, Random motion of a circle microswimmer in a random environment, *New J. Phys.* **22**, 073022 (2020).

[8] M. Bär, R. Großmann, S. Heidenreich, and F. Peruani, Self-propelled rods: Insights and perspectives for active matter, *Annu. Rev. Condens. Matter Phys.* **11**, 441 (2020).

[9] C. Pozrikidis, *A Practical Guide to Boundary Element Methods with the Software Library BEMLIB* (CRC Press, Boca Raton, FL, 2002).

[10] D. J. Smith, E. A. Gaffney, J. R. Blake, and J. C. Kirkman-Brown, Human sperm accumulation near surfaces: A simulation study, *J. Fluid Mech.* **621**, 289 (2009).

[11] B. J. Walker, R. J. Wheeler, K. Ishimoto, and E. A. Gaffney, Boundary behaviours of *Leishmania mexicana*: A hydrodynamic simulation study, *J. Theor. Biol.* **462**, 311 (2019).

[12] H. Shum, E. A. Gaffney, and D. J. Smith, Modelling bacterial behaviour close to a no-slip plane boundary: The influence of bacterial geometry, *Proc. R. Soc. London A* **466**, 1725 (2010).

[13] M. W. Rostami and S. D. Olson, Fast algorithms for large dense matrices with applications to biofluids, *J. Comput. Phys.* **394**, 364 (2019).

[14] N. Taketoshi, T. Omori, and T. Ishikawa, Elasto-hydrodynamic interaction of two swimming spermatozoa, *Phys. Fluids* **32**, 101901 (2020).

[15] E. Lauga and T. R. Powers, The hydrodynamics of swimming microorganisms, *Rep. Prog. Phys.* **72**, 096601 (2009).

[16] K. Ishimoto and E. A. Gaffney, A study of spermatozoan swimming stability near a surface, *J. Theor. Biol.* **360**, 187 (2014).

[17] B. J. Walker, K. Ishimoto, and E. A. Gaffney, Pairwise hydrodynamic interactions of synchronized spermatozoa, *Phys. Rev. Fluids* **4**, 093101 (2019).

[18] K. Ishimoto and E. A. Gaffney, Hydrodynamic clustering of human sperm in viscoelastic fluids, *Sci. Rep.* **8**, 15600 (2018).

[19] A. Creppy, O. Praud, X. Druart, P. L. Kohnke, and F. Plouraboué, Turbulence of swarming sperm, *Phys. Rev. E* **92**, 032722 (2015).

[20] J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär, and R. E. Goldstein, Fluid Dynamics of Bacterial Turbulence, *Phys. Rev. Lett.* **110**, 228102 (2013).

[21] C.-K. Tung, C. Lin, B. Harvey, A. G. Fiore, F. Ardon, M. Wu, and S. S. Suarez, Fluid viscoelasticity promotes collective swimming of sperm, *Sci. Rep.* **7**, 3152 (2017).

[22] D. M. Woolley, Motility of spermatozoa at surfaces, *Reproduction* **126**, 259 (2003).

[23] D. M. Woolley, R. F. Crockett, W. D. I. Groom, and S. G. Revell, A study of synchronisation between the flagella of bull spermatozoa, with related observations, *J. Exp. Biol.* **212**, 2215 (2009).

[24] H. Shum and E. A. Gaffney, Hydrodynamic analysis of flagellated bacteria swimming near one and between two no-slip plane boundaries, *Phys. Rev. E* **91**, 033012 (2015).

[25] G. Rosser, R. E. Baker, J. P. Armitage, and A. G. Fletcher, Modelling and analysis of bacterial tracks suggest an active reorientation mechanism in *Rhodobacter sphaeroides*, *J. R. Soc., Interface* **11**, 20140320 (2014).

[26] D. M. Richards and R. G. Endres, Target shape dependence in a simple model of receptor-mediated endocytosis and phagocytosis, *Proc. Natl. Acad. Sci. USA* **113**, 6113 (2016).
[27] C. Moreau and K. Ishimoto, Driving a microswimmer with wall-induced flow, *Micromachines* **12**, 1025 (2021).

[28] J. Słomka, U. Alcolombri, E. Secchi, R. Stocker, and V. I. Fernandez, Encounter rates between bacteria and small sinking particles, *New J. Phys.* **22**, 043016 (2020).

[29] D. Saintillan and M. J. Shelley, Instabilities and Pattern Formation in Active Particle Suspensions: Kinetic Theory and Continuum Simulations, *Phys. Rev. Lett.* **100**, 178103 (2008).

[30] M.-C. Kim and C. Klapperich, A new method for simulating the motion of individual ellipsoidal bacteria in microfluidic devices, *Lab Chip* **10**, 2464 (2010).

[31] B. J. Walker, K. Ishimoto, R. J. Wheeler, and E. A. Gaffney, Response of monoflagellate pullers to a shearing flow: A simulation study of microswimmer guidance, *Phys. Rev. E* **98**, 063111 (2018).

[32] G. B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, *Proc. R. Soc. London A* **102**, 161 (1922).

[33] F. P. Bretherton, The motion of rigid particles in a shear flow at low Reynolds number, *J. Fluid Mech.* **14**, 284 (1962).

[34] H. Brenner, The Stokes resistance of an arbitrary particle III, *Chem. Eng. Sci.* **19**, 631 (1964).

[35] K. Ishimoto, Helicoidal particles and swimmers in a flow at low Reynolds number, *J. Fluid Mech.* **892**, A11 (2020).

[36] K. Ishimoto, Jeffery orbits for an object with discrete rotational symmetry, *Phys. Fluids* **32**, 081904 (2020).

[37] K. Ishimoto, H. Gadêlha, E. A. Gaffney, D. J. Smith, and J. Kirkman-Brown, Coarse-Graining the Fluid Flow Around a Human Sperm, *Phys. Rev. Lett.* **118**, 124501 (2017).

[38] C. M. Bender and S. A. Orszag, *Advanced Mathematical Methods for Scientists and Engineers I* (Springer, New York, 1999).

[39] U. W. Goodenough, Motile detergent-extracted cells of Tetrahymena and Chlamydomonas. *J. Cell Biol.* **96**, 1610 (1983).

[40] L. F. Shampine and M. W. Reichelt, The MATLAB ODE suite, *SIAM J. Sci. Comput.* **18**, 1 (1997).

[41] V. Singh, D. L. Koch, and A. D. Stroock, Rigid ring-shaped particles that align in simple shear flow, *J. Fluid Mech.* **722**, 121 (2013).

[42] N. S. Borker and D. L. Koch, Slender body theory for particles with non-circular cross-sections with application to particle dynamics in shear flows, *J. Fluid Mech.* **877**, 1098 (2019).

[43] P. Chen and Q. Zhang, Dynamical solutions for migration of chiral dna-type objects in shear flows, *Phys. Rev. E* **84**, 056309 (2011).

[44] J. Fries, J. Einarsson, and B. Mehlig, Angular dynamics of small crystals in viscous flow, *Phys. Rev. Fluids* **2**, 014302 (2017).

[45] I. R. Thorp and J. R. Lister, Motion of a non-axisymmetric particle in viscous shear flow, *J. Fluid Mech.* **872**, 532 (2019).

[46] See http://dx.doi.org/10.5287/bodleian:xveKBqVPm.

[47] M. T. Gallagher and D. J. Smith, Meshfree and efficient modeling of swimming cells, *Phys. Rev. Fluids* **3**, 053101 (2018).