1. Introduction
Bulk properties of particle suspensions depend on particle orientation, and particle--particle interactions can be neglected for sufficiently dilute suspensions (Leal & Hinch Reference Leal and Hinch1971; Saintillan & Shelley Reference Saintillan and Shelley2015). Hence, in many cases interaction with the local flow is a key factor in particle orientation. For small enough particles, viscous effects dominate and orientation is mainly forced by the local shear flow approximation. Thus, the rotational dynamics of single particles in viscous shear flows are of fundamental interest in fluid mechanics.
A classic result in fluid mechanics is that passive spheroids undergo non-uniform rotation in viscous shear flow. Their angular dynamics are governed by Jeffery’s equations (Jeffery Reference Jeffery1922; Taylor Reference Taylor1923), and their periodic but uneven rotations are called Jeffery’s orbits. The precise nature of the Jeffery’s orbit depends on the spheroid aspect ratio
$r \in (0, \infty )$
via the quantity
$(r^2 - 1)/(r^2 + 1)$
, with values of
$r$
further from unity causing more uneven dynamics.
More generally, Jeffery’s equations hold for a wider class of particles beyond spheroids, including passive axisymmetric objects (Bretherton Reference Bretherton1962; Brenner Reference Brenner1964), parameterised via a coefficient called the Bretherton parameter,
$B$
. The derived governing equations for axisymmetric objects are equivalent to Jeffery’s equations when equating
$B = (r^2 - 1)/(r^2 + 1)$
. Therefore, axisymmetric objects with
$B \in (-1,1)$
demonstrate angular dynamics in shear flow equivalent to those of a spheroid with an aspect ratio of
$r = \sqrt{(1+B)/(1-B)}$
.
Asymmetry of particles can induce fundamentally different behaviours to axisymmetric objects (Hinch & Leal Reference Hinch and Leal1979; Roggeveen & Stone Reference Roggeveen and Stone2022; Miara et al. Reference Miara, Vaquero-Stainer, Pihler-Puzović, Heil and Juel2024). For example, helicoidal objects are governed by modified versions of Jeffery’s equations, with extra terms characterised by two additional coefficients that account for chiral effects (Ishimoto Reference Ishimoto2020). In addition, the loss of axial symmetry caused by replacing spheroids with triaxial ellipsoids, and more generally to particles with two orthogonal planes of symmetry, can generate chaotic dynamics (Yarin, Gottlieb & Roisman Reference Yarin, Gottlieb and Roisman1997; Thorp & Lister Reference Thorp and Lister2019).
In the studies mentioned in the paragraph above, the particles are passive. However, particle activity makes interactions with fluid flow much more complicated (Lauga & Powers Reference Lauga and Powers2009; Wittkowski & Löwen Reference Wittkowski and Löwen2012; Elgeti, Winkler & Gompper Reference Elgeti, Winkler and Gompper2015; Saintillan Reference Saintillan2018; Junot et al. Reference Junot, Figueroa-Morales, Darnige, Lindner, Soto, Auradou and Clément2019). Recent work deriving the emergent behaviour of single active particles in shear flow has shown that self-propelled objects exhibiting fast-scale periodic motion can generate emergent slow angular dynamics in shear flow (Ishimoto Reference Ishimoto2023). For example, in two dimensions, the oscillatory yawing of ellipses (Walker et al. Reference Walker, Ishimoto, Gaffney, Moreau and Dalwadi2022), and, in three dimensions, the constant rotation of spheroids and helicoidal objects (Dalwadi et al. Reference Dalwadi, Moreau, Gaffney, Ishimoto and Walker2024a ,Reference Dalwadi, Moreau, Gaffney, Walker and Ishimoto b ). These studies use the method of multiple scales (Hinch Reference Hinch1991) to understand the nonlinear interaction between the fast self-propulsion and slow shear flow, and demonstrate that this interaction generates emergent angular dynamics equivalent to those of a passive particle. The calculated shapes of these equivalent effective particles depend on the type of fast motion and the original shapes. Notably, in these scenarios the effective shape maintains the hydrodynamic symmetries of the original particle. The method of multiple scales has also been used recently to understand the effective dynamics of particles in unsteady flow fields (Ma, Pujara & Thiffeault Reference Ma, Pujara and Thiffeault2022; Pujara & Thiffeault Reference Pujara and Thiffeault2023; Ventrella et al. Reference Ventrella, Pujara, Boffetta, Cencini, Thiffeault and De Lillo2023).
The specific type of activity we are interested in here is rapid yawing in three dimensions. Undulatory motion can be observed in many microswimmers, especially flagellates (Guasto, Rusconi & Stocker Reference Guasto, Rusconi and Stocker2012), for example,Chlamydomonas (Leptos et al. Reference Leptos, Chioccioli, Furlan, Pesci and Goldstein2023) and spermatozoon (Shaebani et al. Reference Shaebani, Wysocki, Winkler, Gompper and Rieger2020). The latter is particularly well studied due to the implications for understanding sperm motility, impacting our understanding of motility-based male fertility diagnostics, reproductive toxicology and basic sperm function (Gaffney et al. Reference Gaffney, Gadêlha, Smith, Blake and Kirkman-Brown2011; Walker et al. Reference Walker, Phuyal, Ishimoto, Tung and Gaffney2020).
For simplicity and analytic tractability, we neglect the complexities of how exactly the rapid motion arises at the microswimmer scale. With the goal of gaining insight into how rapid microscale motion interacts with a far-field shear flow, we consider a simple model of self-generated rapid yawing of a rigid spheroid in steady Stokes flow, with an accompanying self-generated translation. Some care must be taken in considering the limit of rapid yawing while using steady Stokes flow, since very fast oscillations can induce the inclusion of a time derivative via the unsteady Stokes equations (Clarke et al. Reference Clarke, Cox, Williams and Jensen2005). We justify our consideration of steady Stokes flow here by noting that the oscillatory Reynolds number
$\Omega L^2 / \nu$
(where
$\Omega$
is the frequency of rotation,
$L$
is the swimmer length and
$\nu$
is the kinematic viscosity of the fluid) is generally small for microswimmers, despite their fast self-generated motions (Lauga Reference Lauga2020). Hence, while our rapid yawing analysis is relevant for the typical parameter values of many microswimmers, we emphasize that it would formally break down for very large rotation rates (e.g.
$10^4$
Hz for a bacterium with a length scale of
$10 \, \mu$
m in water), when induced inertial terms would become important.
The emergent dynamics for yawing in three dimensions cannot be understood by simply combining results for yawing in two dimensions and constant rotation in three dimensions, for two main reasons. First, three-dimensional (3-D) orientation is governed by a system of three nonlinear equations, in comparison to just a single nonlinear equation in two dimensions. Second, yawing corresponds to a time-dependent rather than constant angular velocity, the latter being the case for constant rotation. This means that yawing in three dimensions is governed by a non-autonomous nonlinear 3-D dynamical system at leading order. Using a multiple scales analysis for systems, in this study we show that rapid yawing generates asymmetric emergent effects that are not present in the original particle. This emergent behaviour is fundamentally different to that arising from yawing in two dimensions (Walker et al. Reference Walker, Ishimoto, Gaffney, Moreau and Dalwadi2022) and constant rotation in three dimensions (Dalwadi et al. Reference Dalwadi, Moreau, Gaffney, Ishimoto and Walker2024a ,Reference Dalwadi, Moreau, Gaffney, Walker and Ishimoto b ). In particular, we demonstrate that the emergent asymmetry generated by rapid yawing can result in chaotic dynamics, which is not possible for passive spheroids nor for the emergent dynamics arising from rapid (constant) rotation.
We start in § 2 below by setting up the physical problem and equations of motion, including in § 2.1 a short summary of the main results we derive subsequently. We present our main analysis in § 3, where we derive the emergent rotational dynamics of the system, relegating some of the technical details to Appendices A and B. In § 4 we demonstrate that the asymmetry generated in the emergent equations can exhibit chaotic dynamics, which is a fundamentally different behaviour to that seen for passive spheroids. We then derive the emergent translational dynamics in § 5. Finally, we discuss our results and their wider implications in § 6.
2. Problem set-up
We consider the dynamics of a self-propelling rigid spheroid in a viscous (Stokes) fluid with an imposed far-field shear flow. We work in dimensionless quantities, scaling time with the inverse shear rate of the imposed far-field flow and space with the equatorial radius of the spheroid. The distance from the centre of the spheroid to its pole along the symmetry axis is
$r$
, and the spheroid self-generates a fast periodic yawing within a swimmer-fixed plane containing its symmetry axis. This yawing manifests through an unsteady angular velocity
$\boldsymbol {\Omega }(t)$
in a quiescent fluid, where
$t$
denotes time. The fast yawing means that the orientation of the swimmer varies rapidly in the laboratory frame. The spheroid also self-generates a translation
$\boldsymbol {V}(t)$
, periodic in a swimmer-fixed reference frame we define below, and with the same period as the yawing.
We define the spheroidal axis of symmetry via a swimmer-fixed axis of
$\hat {\boldsymbol {e}}_{1}$
, and define the self-generated yawing through its angular velocity
$\boldsymbol {\Omega }(t)$
, which is perpendicular to
$\hat {\boldsymbol {e}}_{1}$
. We then define
$\hat {\boldsymbol {e}}_{2}$
to be the direction of
$\boldsymbol {\Omega }$
, and write
\begin{align} \boldsymbol {\Omega }(t) = \Omega A \cos \left (\Omega t\right ) \hat {\boldsymbol {e}}_{2}, \end{align}
where
$\Omega \gg 1$
is the fast frequency of the yawing and
$A$
is the amplitude of the yawing. Finally, we define
$\hat {\boldsymbol {e}}_{3} = \hat {\boldsymbol {e}}_{1} \times \hat {\boldsymbol {e}}_{2}$
. We define the orthonormal basis of the laboratory frame to be
$\{\boldsymbol {e}_{1},\boldsymbol {e}_{2},\boldsymbol {e}_{3}\}$
, orientated in terms of the far-field shear flow that has velocity field
$\boldsymbol {u}(x,y,z) = y \boldsymbol {e}_{3}$
for coordinates
$(x,y,z)$
in the laboratory frame. These definitions are illustrated in figure 1. In this swimmer-fixed basis we also prescribe the translational velocity
\begin{align} \boldsymbol {V}(t) = \sum _{i = 1}^{3} V_i(t) \hat {\boldsymbol {e}}_{i}, \quad \text {with } V_i(t) = a_i + b_i \cos \left (\Omega t - \delta _i \right ), \end{align}
where, in the
$\hat {\boldsymbol {e}}_{i}$
direction,
$a_i$
is the average translational velocity,
$b_i$
is the amplitude of the translational velocity oscillation and
$\delta _i$
is the phase shift of this oscillation.
Schematic of (a) physical set-up and (b) Euler angle definitions, with laboratory (
$\boldsymbol {e}_{i}$
) and swimmer-fixed (
$\hat {\boldsymbol {e}}_{i}$
) frames denoted by black and green arrows, respectively. The Euler angle rotations in (b) occur in the order
$\phi$
,
$\theta$
,
$\psi$
. The swimmer self-generates a rapid yawing via the time-dependent angular velocity
$\boldsymbol {\Omega }(t) = \Omega A \cos (\Omega t) \hat {\boldsymbol {e}}_{2}$
(curved purple arrows), a time-dependent translational velocity
$\boldsymbol {V}(t)$
(straight purple arrow in panel a) and interacts with a far-field shear flow
$\boldsymbol {u} = y \boldsymbol {e}_{3}$
(blue arrows).
A key goal of our subsequent analysis is to understand the dynamics of the particle as it interacts with the far-field flow. To quantify these, we use Euler angles
$\theta \in [0, \pi )$
(pitch),
$\psi \text { mod } 2\pi$
(roll) and
$\phi \text { mod } 2\pi$
(yaw), illustrated in figure 1(b), formally defined via an
$xyx$
-Euler angle transformation
\begin{align} \hat {\boldsymbol {e}}(\theta , \psi , \phi ) = \boldsymbol {\textsf {M}}(\theta ,\psi ,\phi ) \boldsymbol {e}, \end{align}
where
$\hat {\boldsymbol {e}} := (\hat {\boldsymbol {e}}_{1}, \hat {\boldsymbol {e}}_{2}, \hat {\boldsymbol {e}}_{3})^{\intercal }$
is the swimmer-fixed basis, which depends on the orientation,
$\boldsymbol {e} := (\boldsymbol {e}_{1}, \boldsymbol {e}_{2}, \boldsymbol {e}_{3})^{\intercal }$
is the laboratory basis, and we define
\begin{align} \boldsymbol {\textsf {M}}(\theta ,\psi ,\phi ) := \left (\begin{array}{c|c|c} c_\theta & s_\phi s_\theta & - c_\phi s_\theta \\ s_\psi s_\theta & \hphantom {+}c_\phi c_\psi - s_\phi c_\theta s_\psi & \hphantom {+}s_\phi c_\psi + c_\phi c_\theta s_\psi \\ c_\psi s_\theta & - c_\phi s_\psi -s_\phi c_\theta c_\psi & -s_\phi s_\psi + c_\phi c_\theta c_\psi \end{array}\right ), \end{align}
using the shorthand notation
$s_{\theta } = \sin \theta$
, etc. Then, applying the model derivation of Dalwadi et al. (Reference Dalwadi, Moreau, Gaffney, Ishimoto and Walker2024a
,Reference Dalwadi, Moreau, Gaffney, Walker and Ishimoto
b
) to the motion (2.1), the resulting rotational dynamics for a spheroid in shear flow with self-induced yawing are
\begin{align} \frac {\mathrm {d} \theta }{\mathrm {d} t} &= \Omega A \cos \left (\Omega t\right ) \cos \psi + g_1(\theta ,\phi ; B), \end{align}
\begin{align} \frac {\mathrm {d} \psi }{\mathrm {d} t} &= - \Omega A \cos \left (\Omega t\right ) \dfrac {\cos \theta \sin \psi }{\sin \theta } + g_2(\theta ,\phi ; B), \end{align}
\begin{align} \frac {\mathrm {d} \phi }{\mathrm {d} t} &= \Omega A \cos \left (\Omega t\right ) \dfrac {\sin \psi }{\sin \theta } + g_3(\phi ; B), \end{align}
with arbitrary initial conditions. Here, the first terms on the right-hand sides represent the fast yawing motion, and the remaining functions
$g_i$
(
$i = 1, 2, 3$
, here and henceforth) represent the slow interaction with the far-field shear flow. The functions
$g_i$
that quantify this interaction are
\begin{align} g_1 = -\dfrac {B}{2} \cos \theta \sin \theta \sin 2 \phi , \qquad g_2 = \dfrac {B}{2} \cos \theta \cos 2 \phi , \qquad g_3 = \dfrac {1}{2} \left (1 - B \cos 2 \phi \right ), \end{align}
where
$B(r) = (r^2-1)/(r^2 + 1)$
is the Bretherton parameter (Bretherton Reference Bretherton1962), which takes values
$B \in (-1, 1)$
for a spheroid. If there were no yawing (
$A = 0$
), the system (2.4) would reduce exactly to Jeffery’s equations for the orientation of a passive spheroid in shear flow (Jeffery Reference Jeffery1922). We emphasize that the right-hand sides of Jeffery’s equations (2.5) are independent of the spheroid roll
$\psi$
, illustrating their axisymmetry.
The translational dynamics for the centre of mass of the spheroid (
$\boldsymbol {X} = X \boldsymbol {e}_{1} + Y \boldsymbol {e}_{2} + Z \boldsymbol {e}_{3}$
) in shear flow with self-induced translation are
\begin{align} \frac {\mathrm {d} \boldsymbol {X}}{\mathrm {d} t} = \boldsymbol {V}(t) + Y \boldsymbol {e}_{3}, \end{align}
with initial conditions
$\boldsymbol {X}(0) = \boldsymbol {0}$
, noting that we are free to prescribe the origin of the laboratory frame to be at the initial spheroid centre of mass without loss of generality. Note that
$\boldsymbol {V}(t)$
, defined in (2.2), is only a straightforward oscillation in the swimmer frame, and that (2.6) is strongly coupled to the angular dynamics via the evolution of the spheroid orientation through (2.4). However, the reverse is not true -- the angular dynamics (2.4) and (2.5) do not depend on the translational dynamics (2.6). As such, it will be helpful to first consider the angular dynamics in our analysis, then to investigate the translational dynamics.
The full dynamics of the non-autonomous nonlinear dynamical system (2.4)–(2.6) in the fast yawing limit
$\Omega \gg 1$
(black lines in figure 2) have two main effects that occur over distinct times cales. These are (a) yawing over a fast
$t = \mathit {O}(1/\Omega )$
time scale and (b) shear interaction over a slow
$t = \mathit {O}(1)$
time scale. In this study we investigate the emergent dynamics of (2.4)–(2.6) in the fast yawing limit where
$\Omega \gg 1$
(and all other parameters are of
$\mathit {O}(1)$
). This is a singular perturbation problem where the fast oscillatory effects are maintained over the slow shear time scale, i.e. the emergent dynamics cannot be obtained by simply ignoring the slow evolution due to the shear interaction (see red lines in figures 2 and 3). We therefore use the method of multiple scales to calculate the emergent effects.
Numerical solutions of the full rotational dynamics (2.4) and (2.5) (solid blue lines), compared with (a) ignoring the slow evolution, by setting
$g_i = 0$
in (2.4) (solid red lines) and (b) the asymptotic solutions, consisting of the leading-order solutions we derive in (3.5) and the emergent slow evolution equations we derive in (3.17), where the latter are solved numerically (dotted black lines). We use parameter values
$B = 0.9$
,
$A = 2$
and
$\Omega = 3$
with initial conditions
$(\theta ,\psi , \phi ) = (\pi /6, \pi /12, \pi /12)$
. We see that the emergent (asymptotic) dynamics we derive in the limit of large
$\Omega$
agree well with the full dynamics, even for moderate values of
$\Omega$
.
Numerical solutions of the full translational dynamics (2.6), which also depend on the solution to the full rotational dynamics (2.4)–(2.5) (solid blue lines), compared with (a) ignoring the slow evolution, by setting
$g_i = 0$
in (2.4) (solid red lines), and (b) the asymptotic solutions, from the emergent slow evolution equations we derive in (5.7), solved numerically (dotted black lines). We use the same parameter values as in figure 2, and additionally
$\boldsymbol {V}(t)$
is defined in (2.2), with
$(a_1, a_2, a_3) = (-0.2, 0.5, 0.2)$
,
$(b_1, b_2, b_3) = (0.2, 0.6, 0.5)$
,
$(\delta _1, \delta _2, \delta _3) = (\pi /2, \pi /4, -\pi /4)$
and initial conditions
$\boldsymbol {X}(0) = \boldsymbol {0}$
.
2.1. Summary of main results
We show that the rotational dynamics of rapidly yawing spheroids in three dimensions do not act as effective passive spheroids in general, in contrast to recent results for yawing in two dimensions (Walker et al. Reference Walker, Ishimoto, Gaffney, Moreau and Dalwadi2022) and constant rotation in three dimensions (Dalwadi et al. Reference Dalwadi, Moreau, Gaffney, Ishimoto and Walker2024a ). While we show that the orientation of active spheroids here can be related to an equivalent passive particle, this effective particle only has two planes of symmetry in general, thereby losing axial symmetry. This difference in the symmetries of the effective particle is particularly notable because passive particles with two planes of symmetry have been shown to demonstrate chaotic dynamics (Yarin et al. Reference Yarin, Gottlieb and Roisman1997; Thorp & Lister Reference Thorp and Lister2019). However, the classic (3-D) Jeffery equations for passive spheroids are integrable (Jeffery Reference Jeffery1922), constraining the dynamics to a two-dimensional (2-D) surface in phase space and, therefore, ruling out chaos (Thorp & Lister Reference Thorp and Lister2019).
Therefore, with the benefit of hindsight, it is helpful to record here the rotational dynamical equations for the appropriate class of asymmetric passive particles that will emerge; particles with two orthogonal planes of symmetry (Bretherton Reference Bretherton1962; Brenner Reference Brenner1964; Harris, Nawaz & Pittman Reference Harris, Nawaz and Pittman1979; Thorp & Lister Reference Thorp and Lister2019). These can be written as modified versions of Jeffery’s equations ((2.4) with
$A = 0$
) transformed via
\begin{align} g_i(\theta ,\phi ;B) \mapsto \tilde {g}_i(\theta ,\psi ,\phi ;B_1,B_2,B_3) = g_i\left (\theta ,\phi ;B_1\right ) + h_i(\theta ,\psi ,\phi ;B_1,B_2,B_3), \end{align}
where the
$h_i$
encode non-axisymmetric effects via their dependence on
$\psi$
, and are defined as
\begin{align} h_1(\theta ,\psi , \phi ; B_1, B_2) &= \dfrac {B_1 + B_2}{2} s_{\theta } s_{\psi }\left (c_\theta s_\psi s_{2 \phi } - c_\psi c_{2 \phi } \right ), \end{align}
\begin{align} h_2(\theta , \psi , \phi ; B_1, B_2, B_3) &= \dfrac {B_1 + B_2 + B_3}{2} c_{\theta } c_{\psi } \left (c_\theta s_{\psi } s_{2 \phi } - c_{\psi } c_{2 \phi } \right )\nonumber\\&\quad + \dfrac {B_3}{2} s_{\psi } \left ( c_\psi s_{2 \phi } + c_\theta s_{\psi } c_{2 \phi } \right ), \end{align}
\begin{align} h_3(\theta ,\psi , \phi ; B_1, B_2) &= \dfrac {B_1 + B_2}{2} \left ( c^2_{\psi } c_{2 \phi } - c_{\theta } s_{\psi } c_{\psi } s_{2\phi } \right ), \end{align}
again using the shorthand notation
$s_{\theta } = \sin \theta$
, etc. Importantly, the transformation (2.7a
) introduces three independent coefficients
$B_i$
in place of
$B$
. For an ellipsoid, the governing equations are the same as (2.7) (Hinch & Leal Reference Hinch and Leal1979), except that
$B_i$
are not independent. In fact, Jeffery (Reference Jeffery1922) showed that an ellipsoid with axes
$a$
,
$b$
,
$c$
has the explicit relationship
\begin{align} B_1 = \frac {c^2 - b^2}{c^2 + b^2}, \qquad B_2 = \frac {a^2 - c^2}{a^2 + c^2}, \qquad B_3 = \frac {b^2 - a^2}{b^2 + a^2}, \end{align}
from which it follows (Bretherton Reference Bretherton1962) that, for an ellipsoid,
$B_i$
must satisfy the relationship
\begin{align} B_1 B_2 B_3 + B_1 + B_2 + B_3 = 0. \end{align}
For an axisymmetric ellipsoid with
$a = b$
(i.e. a spheroid),
$B_1 = -B_2 = B$
and
$B_3 = 0$
, causing
$h_i= 0$
in (2.7) and removing all non-axisymmetric effects, as expected.
In our analysis, we find that the average orientation
$(\bar {\vartheta },\bar {\varPsi },\bar {\varphi })$
(defined appropriately below) evolves over the
$t = \mathit {O}(1)$
scale, with emergent dynamics governed by
\begin{align} \frac {\mathrm {d} \bar {\vartheta }}{\mathrm {d} t} &= g_1(\bar {\vartheta },\bar {\varphi }; \widehat {B}_1) + h_1(\bar {\vartheta },\bar {\varPsi },\bar {\varphi }; \widehat {B}_1,\widehat {B}_2), \end{align}
\begin{align} \frac {\mathrm {d} \bar {\varPsi }}{\mathrm {d} t} &= g_2(\bar {\vartheta },\bar {\varphi }; \widehat {B}_1) + h_2(\bar {\vartheta },\bar {\varPsi },\bar {\varphi }; \widehat {B}_1, \widehat {B}_2,\widehat {B}_3), \end{align}
\begin{align} \frac {\mathrm {d} \bar {\varphi }}{\mathrm {d} t} &= g_3(\bar {\varphi }; \widehat {B}_1) + h_3(\bar {\vartheta },\bar {\varPsi },\bar {\varphi }; \widehat {B}_1,\widehat {B}_2), \end{align}
with
$g_i$
and
$h_i$
defined in (2.5) and (2.7), but now where
$\widehat {B}_i$
are functions of
$A$
and
$B(r)$
, and do not generally satisfy (2.9). That is, we show that (i) a rapidly yawing active spheroid generates a loss of axial symmetry in its emergent dynamics (via the
$h_i$
terms), (ii) the emergent dynamics are equivalent to those of an effective passive particle with two planes of symmetry, and (iii) the effective shape is not an ellipsoid in general. The main part of our analysis involves deriving (2.10) from (2.4) and (2.5), including explicit representations of the average orientation
$(\bar {\vartheta },\bar {\varPsi },\bar {\varphi })$
and the functions
$\widehat {B}_i(A, B)$
. We also show that the emergent dynamics we derive can demonstrate chaos, and therefore, that these dynamics can be fundamentally different to those for any passive spheroid.
Finally, we also determine the emergent translational dynamics over the
$t = \mathit {O}(1)$
scale. Though we find that the spheroid translates with an effective velocity
$\widehat {\boldsymbol {V}}$
in addition to a straightforward shear contribution, this effective velocity is not simply the average of the physical velocity
$\boldsymbol {V}(t)$
in the swimmer-fixed frame. Importantly, we show that phase lags in the translational velocity oscillations can generate non-trivial contributions to the effective translational velocity.
3. Deriving the emergent rotational dynamics
3.1. Framework for the method of multiple scales
We start by analysing the emergent rotational dynamics (2.4) and (2.5) in the limit of rapid yawing, corresponding to
$\Omega \gg 1$
. We use the method of multiple scales (Hinch Reference Hinch1991) to derive equations for the emergent behaviour. We introduce the ‘fast’ time scale
$T = \mathit {O}(1)$
via
$T = \Omega t$
, and refer to the original time scale
$t$
as the ‘slow’ time scale. Using the method of multiple scales, we treat these two time scales as independent, transforming the time derivative as
\begin{align} \frac {\mathrm {d} }{\mathrm {d} t} \mapsto \Omega \frac {\partial }{\partial T} + \frac {\partial }{\partial t}. \end{align}
Under the time derivative mapping (3.1), the dynamical system for the swimmer orientation (2.4) is transformed to
\begin{align} \Omega \frac {\partial \theta }{\partial T} + \frac {\partial \theta }{\partial t} &= \Omega A \cos T \cos \psi + g_1(\theta ,\phi ;B), \end{align}
\begin{align} \Omega \frac {\partial \psi }{\partial T} + \frac {\partial \psi }{\partial t} &= - \Omega A \cos T \dfrac {\cos \theta \sin \psi }{\sin \theta } + g_2(\theta ,\phi ;B), \end{align}
\begin{align} \Omega \frac {\partial \phi }{\partial T} + \frac {\partial \phi }{\partial t} &= \Omega A \cos T \dfrac {\sin \psi }{\sin \theta } + g_3(\phi ;B). \end{align}
We expand each dependent variable as an asymptotic series in inverse powers of
$\Omega$
and as a function of both the fast and slow timescales, writing
\begin{align} y(T,t) \sim y_0(T,t) + \dfrac {1}{\Omega } y_1(T,t) \quad \text {as } \Omega \to \infty , \text { for } y \in \{\theta , \psi , \phi \}. \end{align}
3.2. Leading-order analysis
Using the asymptotic expansions (3.3) in the transformed governing equations (3.2), we obtain the leading-order (i.e.
$\mathit {O}(\Omega )$
) system
\begin{align} \frac {\partial \theta _0}{\partial T} = A \cos T \cos \psi _0, \qquad \frac {\partial \psi _0}{\partial T} = - A \cos T \dfrac {\cos \theta _0 \sin \psi _0}{\sin \theta _0}, \qquad \frac {\partial \phi _0}{\partial T} = A \cos T \dfrac {\sin \psi _0}{\sin \theta _0}. \end{align}
We derive an exact solution to the nonlinear, non-autonomous leading-order system (3.4) in Appendix A, obtaining
\begin{align} \cos \theta _0 &= \cos \bar {\vartheta } \cos f(T) - \sin \bar {\vartheta } \cos \bar {\varPsi } \sin f(T), \end{align}
\begin{align} \sin \theta _0 \sin \psi _0 &= \sin \bar {\vartheta } \sin \bar {\varPsi }, \end{align}
\begin{align} \tan (\phi _0 - \bar {\varphi }) &= \dfrac {\sin \bar {\varPsi } \sin f(T)}{\sin \bar {\vartheta } \cos f(T) + \cos \bar {\vartheta } \cos \bar {\varPsi } \sin f(T)}, \end{align}
defining
\begin{align} f(T;A) := A \sin T, \end{align}
and where
$\bar {\vartheta }(t)$
,
$\bar {\varPsi }(t)$
and
$\bar {\varphi }(t)$
are the three slow-time functions of integration that remain undetermined from our leading-order analysis. The goal of the next-order analysis in § 3.3 will be to derive the governing equations satisfied by
$\bar {\vartheta }$
,
$\bar {\varPsi }$
and
$\bar {\varphi }$
. The three degrees of freedom that arise from integrating (3.4) could be included as different combinations of
$\bar {\vartheta }$
,
$\bar {\varPsi }$
and
$\bar {\varphi }$
(for example, we could have used
$(\bar {\vartheta },\bar {\varPsi }) \mapsto (\bar {\vartheta },H)$
, where
$H(t) = \sin \bar {\vartheta } \sin \bar {\varPsi }$
); we choose the specific forms in (3.5) so that
$(\bar {\vartheta }, \bar {\varPsi }, \bar {\varphi })$
is associated with
$(\theta , \psi , \phi )$
and represents the average orientation direction of the spheroid over a single yawing oscillation. That is,
$\langle \hat {\boldsymbol {e}}_{i}(\theta ,\psi ,\phi ) \rangle \propto \hat {\boldsymbol {e}}_{i}(\bar {\vartheta },\bar {\varPsi },\bar {\varphi })$
, where we use the notation
$\langle \boldsymbol {\cdot } \rangle$
to denote the average of its argument over one fast-time oscillation, defined as
\begin{align} \left \lt y \right \gt = \dfrac {1}{2 \pi }\int _0^{2\pi } \! y \, \mathrm {d}T. \end{align}
3.3. Next-order system
Our remaining goal is to determine the governing equations satisfied by the slow-time functions
$\bar {\vartheta }(t)$
,
$\bar {\varPsi }(t)$
and
$\bar {\varphi }(t)$
. To do this, we must determine the solvability conditions required for the first-order correction (i.e.
$\mathit {O}(1)$
) terms in (3.2) after posing the asymptotic expansions (3.3). These
$\mathit {O}(1)$
terms are
\begin{align} \frac {\partial \theta _1}{\partial T} + \left (A \cos T\right ) \psi _1 \sin \psi _0 &= g_1(\theta _0,\phi _0) - \frac {\partial \theta _0}{\partial t}, \end{align}
\begin{align} \frac {\partial \psi _1}{\partial T} - \left (A \cos T\right ) \theta _1 \dfrac {\sin \psi _0}{\sin ^2 \theta _0} + \left (A \cos T\right ) \psi _1 \dfrac {\cos \theta _0 \cos \psi _0}{\sin \theta _0} &= g_2(\theta _0,\phi _0) - \frac {\partial \psi _0}{\partial t}, \end{align}
\begin{align} \frac {\partial \phi _1}{\partial T} + \left (A \cos T\right ) \theta _1 \dfrac {\cos \theta _0 \sin \psi _0}{\sin ^2 \theta _0} - \left (A \cos T\right ) \psi _1\dfrac {\cos \psi _0}{\sin \theta _0} &= g_3(\theta _0,\phi _0) - \frac {\partial \phi _0}{\partial t}. \end{align}
To derive the required solvability conditions, we use the method of multiple scales for systems; see, e.g., Dalwadi (Reference Dalwadi2014, pp. 127–128) or Dalwadi et al. (Reference Dalwadi, Chapman, Oliver and Waters2018, p. 22). Namely, the solvability condition for the system
$L \boldsymbol {X} = \boldsymbol {G}$
can be found by calculating the solution of
\begin{align} L^* \boldsymbol {Y} = \boldsymbol {0}, \end{align}
where
$L^*$
is the matrix differential-algebraic linear adjoint operator. The requisite solvability condition for (3.8) is then
\begin{align} \left \lt \boldsymbol {Y} \boldsymbol {\cdot } \boldsymbol {G} \right \gt = 0. \end{align}
Generally, each linearly independent solution of the homogeneous adjoint problem (3.9) will contribute one solvability condition. For (3.8), the homogeneous adjoint operator is the transpose of the matrix operator taking the adjoint of each element, and is therefore
\begin{align} L^* \!=\! \begin{pmatrix} \!\!\!- \partial _{T} & -A \cos T \sin \psi _0/\!\sin ^2 \theta _0 & A \cos T \cos \theta _0 \sin \psi _0 /\! \sin ^2 \theta _0 \\ A \cos T \sin \psi _0 & - \partial _{T} \!+\! A \cos T \cos \theta _0 \cos \psi _0/\sin \theta _0 & - A \cos T \cos \psi _0/\sin \theta _0 \\ 0 & 0 & - \partial _{T} \\ \!\!\end{pmatrix}\!. \end{align}
Hence, this problem is not self-adjoint, i.e.
$L \neq L^*$
.
Even though the system (3.9) with operator (3.11) is non-autonomous, we are able to solve it exactly by deducing appropriate nonlinear transformations. We present this analysis in Appendix B, where we calculate that (3.9), (3.11) has general periodic solution
\begin{align} \boldsymbol {Y} &= C_1 \begin{pmatrix} \cos \theta _0 \sin \psi _0 \\ \sin \theta _0 \cos \psi _0 \\ 0 \end{pmatrix} + C_2 \begin{pmatrix} \cos \psi _0 \\ - \sin \theta _0 \cos \theta _0 \sin \psi _0 \\ 0 \end{pmatrix} + C_3 \begin{pmatrix} 0 \\ \cos \theta _0 \\ 1 \end{pmatrix} \end{align}
for arbitrary constants
$C_1$
,
$C_2$
and
$C_3$
. The solution (3.12) can also be verified a posteriori by direct substitution into (3.9), (3.11) and applying (3.4).
Finally, we derive our required solvability conditions by substituting the adjoint solutions (3.12) into the general solvability condition (3.10), with
$\boldsymbol {G}$
defined as the vector right-hand side of (3.8), and setting the resulting coefficients of
$C_i$
to zero. This procedure yields the following three solvability conditions:
\begin{align} \left \lt \theta _{0 t} \cos \theta _0 \sin \psi _0 + \psi _{0 t} \sin \theta _0 \cos \psi _0 \right \gt &= \left \lt g_1 \cos \theta _0 \sin \psi _0 + g_2 \sin \theta _0 \cos \psi _0 \right \gt , \end{align}
\begin{align} \left \lt \theta _{0 t} \cos \psi _0 -\psi _{0 t} \sin \theta _0 \cos \theta _0 \sin \psi _0 \right \gt &= \left \lt g_1 \cos \psi _0 - g_2 \sin \theta _0 \cos \theta _0 \sin \psi _0 \right \gt , \end{align}
\begin{align} \left \lt \psi _{0t} \cos \theta _0 + \phi _{0t} \right \gt &= \left \lt g_2 \cos \theta _0 + g_3 \right \gt. \end{align}
Here the subscript
$t$
denotes partial differentiation with respect to
$t$
. To derive the emergent governing equations we seek, our remaining task is to evaluate the averages in (3.13) in terms of the slow-time functions
$\bar {\vartheta }$
,
$\bar {\varPsi }$
and
$\bar {\varphi }$
.
3.4. Evaluating the solvability conditions
Using our leading-order solutions (3.5), we evaluate the left-hand sides of (3.13) to obtain
\begin{align} \left \lt \theta _{0 t} \cos \theta _0 \sin \psi _0 + \psi _{0 t} \sin \theta _0 \cos \psi _0 \right \gt &= \bar {\vartheta }_t \cos \bar {\vartheta } \sin \bar {\varPsi } + \bar {\varPsi }_t \sin \bar {\vartheta } \cos \bar {\varPsi }, \end{align}
\begin{align} \left \lt \theta _{0 t} \cos \psi _0 -\psi _{0 t} \sin \theta _0 \cos \theta _0 \sin \psi _0 \right \gt &= \bar {\vartheta }_t \cos \bar {\varPsi } - \bar {\varPsi }_t \sin \bar {\vartheta } \cos \bar {\vartheta } \sin \bar {\varPsi }, \end{align}
\begin{align} \left \lt \psi _{0t} \cos \theta _0 + \phi _{0t} \right \gt &= \bar {\varPsi }_t \cos \bar {\vartheta } + \bar {\varphi }_t. \end{align}
After more algebra, we can also evaluate the right-hand sides of (3.13) to deduce that
\begin{align} &g_1 \cos \theta _0 \sin \psi _0 + g_2 \sin \theta _0 \cos \psi _0 \notag \\ &\qquad = -\dfrac {B}{2} \left \{s_{\bar {\vartheta }} s_{\bar {\varPsi }} \left [c_{\bar {\vartheta }}^2 c_f^2 - c_{\bar {\varPsi }}^2 (1 + c_{\bar {\vartheta }}^2)s_f^2 \right ] + c_{\bar {\vartheta }}^3 s_{2 \bar {\varPsi }} s_f c_f \right \} s_{2 \bar {\varphi }} \notag \\ &\qquad \quad -\dfrac {B}{2} \left \{ s_{\bar {\vartheta }} c_{\bar {\vartheta }} c_{\bar {\varPsi }} \left ( c_{2\bar {\varPsi }} s_f^2 - c_f^2\right ) + \left [s_{\bar {\vartheta }}^2 c_{\bar {\varPsi }}^2 - c_{\bar {\vartheta }}^2 c_{2 \bar {\varPsi }} \right ] s_f c_f \right \} c_{2 \bar {\varphi }} , \end{align}
\begin{align} &g_1 \cos \psi _0 - g_2 \sin \theta _0 \cos \theta _0 \sin \psi _0 \notag \\ &\qquad = \dfrac {B}{2} \left \{s_{\bar {\vartheta }} c_{\bar {\vartheta }} c_{\bar {\varPsi }} \left (c_{2 \bar {\varPsi }} s_f^2 - c_f^2 \right ) + \left [s_{\bar {\vartheta }}^2 c_{\bar {\varPsi }}^2 - c_{\bar {\vartheta }}^2 c_{2 \bar {\varPsi }} \right ] s_f c_f \right \} s_{2 \bar {\varphi }} \notag \\ &\qquad \quad -\dfrac {B}{2} \left \{s_{\bar {\vartheta }} s_{\bar {\varPsi }} \left [ c_{\bar {\vartheta }}^2 c_f^2 - c_{\bar {\varPsi }}^2 (1 + c_{\bar {\vartheta }}^2) s_f^2 \right ] + c_{\bar {\vartheta }}^3 s_{2 \bar {\varPsi }} s_f c_f \right \} c_{2 \bar {\varphi }}, \end{align}
\begin{align} &g_2 \cos \theta _0 + g_3 = \dfrac {1}{2} \left (1 - B \left \{ (c_{\bar {\vartheta }}^2 c_{\bar {\varPsi }}^2 - s_{\bar {\varPsi }}^2) s_f^2 + s_{\bar {\vartheta }}^2 c_f^2 + 2s_{\bar {\vartheta }} c_{\bar {\vartheta }} c_{\bar {\varPsi }} c_f s_f \right \}c_{2 \bar {\varphi }}\right ) \notag \\ &\qquad \qquad \qquad \, \, \, + B s_{\bar {\varPsi }} \left[ c_{\bar {\vartheta }} c_{\bar {\varPsi }} s_f^2 + s_{\bar{\vartheta}} c_f s_f \right] s_{2 \bar {\varphi }} , \end{align}
using the shorthand notation
$s_{\bar {\vartheta }} = \sin \bar {\vartheta }$
, etc., and where
$f$
is defined in (3.6). Then, noting that
$\langle c_f \rangle = J_0(A)$
and
$\langle s_f \rangle= 0$
, deduced via the integral representation of
$J_0(A)$
(the Bessel function of order zero) and parity arguments, we obtain
\begin{align} \left \lt s_f^2 \right \gt = \dfrac {1}{2}\left (1 - J_0(2 A) \right ), \qquad \left \lt c_f^2 \right \gt = \dfrac {1}{2}\left (1 + J_0(2 A) \right ), \qquad \left \lt s_f c_f \right \gt = 0. \end{align}
Taking the averages of (3.15) using (3.16) allows us to evaluate the right-hand sides of the solvability conditions (3.13). Then, combining this result with (3.14) for the left-hand sides of (3.13) and rearranging, we obtain
\begin{align} \frac {\mathrm {d} \bar {\vartheta }}{\mathrm {d} t} &= -\dfrac {B J_0(2A)}{2} s_{\bar {\vartheta }} c_{\bar {\vartheta }} s_{2\bar {\varphi }} -\dfrac {B (1 - J_0(2A))}{4} s_{\bar {\vartheta }} s_{\bar {\varPsi }} \left (c_{\bar {\vartheta }} s_{\bar {\varPsi }} s_{2 \bar {\varphi }} - c_{\bar {\varPsi }} c_{2 \bar {\varphi }} \right ), \end{align}
\begin{align} \frac {\mathrm {d} \bar {\varPsi }}{\mathrm {d} t} &= \dfrac {B J_0(2A)}{2} c_{\bar {\vartheta }} c_{2 \bar {\varphi }} + \dfrac {B (1 - J_0(2 A))}{4} s_{\bar {\varPsi }} \left ( c_{\bar {\varPsi }} s_{2 \bar {\varphi }} + c_{\bar {\vartheta }} s_{\bar {\varPsi }} c_{2 \bar {\varphi }} \right ), \end{align}
\begin{align} \frac {\mathrm {d} \bar {\varphi }}{\mathrm {d} t} &= \dfrac {1}{2} \left (1 - B J_0(2A) c_{2 \bar {\varphi }}\right ) - \dfrac {B (1 - J_0(2A))}{4} c_{\bar {\varPsi }} \left (c_{\bar {\varPsi }} c_{2 \bar {\varphi }} - c_{\bar {\vartheta }} s_{\bar {\varPsi }} s_{2 \bar {\varphi }} \right ), \end{align}
which are the key results of our analysis; the emergent slow evolution equations we have been seeking. Recombining the slow evolution equations (3.17) with the fast oscillations (3.5), we see the leading-order solutions agree very well with the full dynamics (blue and black lines in figure 2), even for values of
$\Omega$
as low as
$3$
. The agreement improves further as
$\Omega$
increases. Finally, we note that the emergent equations (3.17) have the functional form we claimed in (2.10), with analytically derived representations for the coefficients (illustrated in figure 4):
\begin{align} \widehat {B}_1 = B J_0(2 A), \qquad \widehat {B}_2 = -\dfrac {B}{2} \left (1 + J_0(2 A) \right ), \qquad \widehat {B}_3 = \dfrac {B}{2} \left (1 - J_0(2 A) \right ). \end{align}
The effective coefficients (3.18), obtained by comparing (3.17) with (2.10). The marked stars along the
$x$
axis are the values of
$A$
at which (6.1) is satisfied (i.e.
$J_0(2 A) = 0$
) and, therefore, when the effective shape is an ellipsoid. In fact, in these cases the effective shape is constrained further to a spheroid, whose aspect ratio is given in (6.2).
To summarise, we have demonstrated that the emergent rotational dynamics (3.17) gain an effective asymmetry, with effective coefficients derived in (3.18). This is most explicitly demonstrated through their dependence on
$\bar {\varPsi }$
in the terms proportional to
$B (1 - J_0(2A))$
on the right-hand sides of (3.17). As discussed in § 2.1, the emergent dynamics (3.17) are equivalent to those of a passive object with two orthogonal planes of symmetry. If the effective coefficients (3.18) were independent then it would immediately follow from the results of Yarin et al. (Reference Yarin, Gottlieb and Roisman1997) and Thorp & Lister (Reference Thorp and Lister2019) that chaotic behaviour was possible. Since the effective coefficients (3.18) are not independent, it is not immediately clear whether such behaviour is possible in the system we derive. In the next section we demonstrate that chaotic behaviour is possible in the system (3.17) and (3.18).
4. Chaotic behaviour in the emergent dynamics
To investigate the possibility of chaos in the system (3.17) and (3.18), we follow Hinch & Leal (Reference Hinch and Leal1979) and Yarin et al. (Reference Yarin, Gottlieb and Roisman1997) and Thorp & Lister (Reference Thorp and Lister2019) and define an appropriate Poincaré section. Specifically, we reduce the full 3-D continuous dynamics of (3.17) and (3.18) to a 2-D discrete dynamical system in
$(\bar {\varPsi },\bar {\vartheta })$
whenever
$\bar {\varphi } = n \pi$
for non-negative integers
$n$
. We can do this straightforwardly by solving (3.17a
) and (3.17b
) in terms of
$\bar {\varphi }$
, transforming the time derivatives via
\begin{align} \frac {\mathrm {d} }{\mathrm {d} t} \mapsto \frac {\mathrm {d} \bar {\varphi }}{\mathrm {d} t} \frac {\mathrm {d} }{\mathrm {d} \bar {\varphi }}, \end{align}
and using (3.17c) to evaluate
$\mathrm {d} \bar {\varphi }/ \mathrm {d} t \gt 0$
. The monotonic nature of
$\mathrm {d} \bar {\varphi }/ \mathrm {d} t$
follows from
$|B| \lt 1$
, and ensures that the transformation is well defined. For given system parameters
$A$
(amplitude of yawing),
$B = (r^2 - 1)/(r^2 + 1)$
(Bretherton parameter in terms of spheroid aspect ratio,
$r$
), and initial conditions
$(\bar {\varPsi }_0, \bar {\vartheta }_0) = (\bar {\varPsi }(0), \bar {\vartheta }(0))$
(setting
$\bar {\varphi }(0) = 0$
), this generates an iteration
$(\bar {\varPsi }_n,\bar {\vartheta }_n)$
.
In figures 5 and 6 we show Poincaré sections for
$A = 0$
(i.e. classic Jeffery orbits) and
$A = 0.25$
, respectively. These consist of a point shown for every iteration up to
$n=500$
for different initial conditions. Periodic orbits in
$(\bar {\varPsi }_n,\bar {\vartheta }_n)$
repeat exactly, and are represented by distinct points that repeat themselves after a finite number of iterations. Quasiperiodic orbits appear as one-dimensional (1-D) curves, i.e. orbits consist of points that densely cover a 1-D path but never exactly repeat themselves. These correspond to orbits with more than one periodic component, but with incommensurate frequencies. Finally, chaotic orbits appear as dense patchworks of points in a 2-D region of the Poincaré section.
Poincaré section for the classic Jeffery’s equations, which are equivalent to setting
$A = 0$
in our emergent equations (3.17) and (3.18). Using the Poincaré map outlined in the main text, we use
$A = 0$
,
$B = 0.99$
and iterate up to
$n = 500$
. The full 2-D phase space is obtained by exploiting reflectional symmetry across
$\bar {\vartheta } = \pi /2$
and translational symmetry in
$\bar {\varPsi } \mapsto \bar {\varPsi } + \pi$
. No chaos is possible for the classic Jeffery’s equations, as observed here.
Poincaré section for the emergent dynamics (3.17) and (3.18). Using the Poincaré map outlined in the main text, we use
$A = 0.25$
,
$B = 0.99$
and iterate up to
$n = 500$
. The full 2-D phase space is obtained by exploiting reflectional symmetry across
$\bar {\vartheta } = \pi /2$
and translational symmetry in
$\bar {\varPsi } \mapsto \bar {\varPsi } + \pi$
. Significant regions of chaos are observed in the upper third of the figure.
For the classic Jeffery’s equations (for passive spheroids), the orbits are constant in
$\bar {\vartheta }$
(figure 5). No chaos is possible in this case, which follows from the integrability of the classic 3-D Jeffery’s equations constraining the dynamics to a 2-D surface in phase space. The lack of chaotic dynamics for passive spheroids can be seen visually in figure 5. These features are well known for passive spheroids, and have been explored in detail as part of the more general analysis in Thorp & Lister (Reference Thorp and Lister2019).
In contrast to the classic Jeffery’s equations, we see that chaos is possible for the more general emergent equations we derive here (figure 6). Chaotic regions are plentiful in the upper third of the figure. In the lower two-thirds of the figure, the behaviour is mainly dominated by periodic and quasiperiodic orbits, and we note that these behaviours also appear within islands in the upper third. The quasiperiodic behaviours include orbits that sample all values of
$\bar {\varPsi } \in [- \pi /2, \pi /2]$
and orbits that only take a subset of values therein.
Finally, in figure 7 we also show Poincaré sections for different values of
$A$
, starting with the same initial condition. For lower values of
$A$
, the orbits are quasiperiodic. However, the orbits demonstrate chaos as
$A$
increases, before returning to a quasiperiodic orbit as
$A$
increases further.
Poincaré section for the emergent dynamics (3.17) and (3.18). Using the Poincaré map outlined in the main text, we use
$B = 0.99$
with initial conditions
$(\bar {\varPsi }_0, \bar {\vartheta }_0) = (0, 9\pi /20)$
(shown as a black asterisk) and iterate up to
$n = 1000$
. We use
$A \in \{0, 0.2, 0.25, 0.3, 0.5 \}$
, as described in the legend. For lower values of
$A$
, the orbits are quasiperiodic. However, the orbit is chaotic at
$A = 0.25$
and
$A = 0.3$
, before returning to quasiperiodicity for
$A = 0.5$
.
Hence, we have shown that chaotic behaviour is possible in the emergent dynamical system we have derived (3.17) and (3.18). This behaviour arises as a direct result of the asymmetry that is generated by the rapid yawing of the active spheroids we consider. This asymmetry does not arise for a rapid (constant) rotation of spheroids, for which the emergent dynamics are equivalent to effective passive spheroids. In this latter case, integrability of the system means that chaotic behaviour is not possible, as shown explicitly in Thorp & Lister (Reference Thorp and Lister2019).
5. Deriving the emergent translational dynamics
We now consider the emergent translational dynamics from the system (2.6) in the large-
$\Omega$
limit. The emergent translational dynamics we derive are significantly more straightforward than their rotational equivalents. In multiple scales form, the time derivative is again transformed via (3.1), and so (2.6) becomes
\begin{align} \Omega \frac {\partial \boldsymbol {X}}{\partial T} + \frac {\partial \boldsymbol {X}}{\partial t} = \boldsymbol {V}(T) + Y \boldsymbol {e}_{3}, \end{align}
emphasizing that
$\boldsymbol {X} = \boldsymbol {X}(T,t)$
in general, and explicitly writing the velocity
$\boldsymbol {V}$
from (2.2) in terms of the fast timescale as
\begin{align} \boldsymbol {V}(T) = \sum _{i = 1}^{3} V_i(T) \hat {\boldsymbol {e}}_{i}, \quad \text {where } V_i(T) := a_i + b_i \cos \left (T - \delta _i \right ). \end{align}
Expanding
$\boldsymbol {X} \sim \boldsymbol {X}_0 + (1/\Omega ) \boldsymbol {X}_1$
and substituting into (5.1), the
$\mathit {O}(\Omega )$
terms yield
\begin{align} \frac {\partial \boldsymbol {X}_0}{\partial T} = \boldsymbol {0}, \end{align}
which is trivially solved by
$\boldsymbol {X}_0 = \boldsymbol {X}_0(t)$
. At next order, the
$\mathit {O}(1)$
terms in (5.1) yield
\begin{align} \frac {\partial \boldsymbol {X}_1}{\partial T} + \frac {\mathrm {d} \boldsymbol {X}_0}{\mathrm {d} t} = \boldsymbol {V}(T) + Y_0 \boldsymbol {e}_{3}, \end{align}
where we note that the dependence of
$\boldsymbol {V}$
on
$\hat {\boldsymbol {e}}_{i}$
in (5.2) means that
$\boldsymbol {V}$
depends on the leading-order Euler angles. These are defined via the fast oscillations (3.5) and the slow evolution equations (3.17).
The emergent equations are obtained by averaging (5.4) over a
$2 \pi$
period in
$T$
. Performing this procedure, the fast-time derivative of
$\boldsymbol {X}_1$
vanishes (by design), and we are left with
\begin{align} \frac {\mathrm {d} \boldsymbol {X}_0}{\mathrm {d} t} = \left \lt \sum _{i = 1}^{3} \left [a_i + b_i \cos \left (T - \delta _i \right ) \right ] \hat {\boldsymbol {e}}_{i}(\theta _0, \psi _0, \phi _0) \right \gt + Y_0 \boldsymbol {e}_{3}, \end{align}
where we have replaced
$\boldsymbol {V}(T)$
in the average operator with its specific form via (5.2), to emphasize the dependence of
$\boldsymbol {V}(T)$
on the leading-order spheroid orientation. To evaluate this averaged velocity, we use the transformation between laboratory basis and swimmer-fixed basis (2.3) to write
$\hat {\boldsymbol {e}}_{i}$
explicitly in terms of
$(\theta _0, \psi _0, \phi _0)$
and the laboratory basis. Then, we combine this with our leading-order fast-time solutions (3.5) and directly calculate the following averages:
\begin{align} \left \lt \hat {\boldsymbol {e}}(\theta _0, \psi _0, \phi _0) \right \gt &= \boldsymbol {\textsf {D}}\, \tilde {\boldsymbol {e}}(\bar {\vartheta },\bar {\varPsi },\bar {\varphi }), \end{align}
\begin{align} \left \lt \hat {\boldsymbol {e}}(\theta _0, \psi _0, \phi _0) \sin T \right \gt &= \boldsymbol {\textsf {A}} \, \tilde {\boldsymbol {e}}(\bar {\vartheta },\bar {\varPsi },\bar {\varphi }), \end{align}
\begin{align} \left \lt \hat {\boldsymbol {e}}(\theta _0, \psi _0, \phi _0) \cos T \right \gt &= \boldsymbol {0}. \end{align}
Here
\begin{align} \boldsymbol {\textsf {D}} = \left (\begin{array}{c c c} J_0(A) & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & J_0(A) \end{array}\right ), \quad \boldsymbol {\textsf {A}} = \left (\begin{array}{c c c} 0 & 0 & J_1(A) \\ 0 & 0 & 0 \\ -J_1(A) & 0 & 0 \end{array}\right ), \quad \tilde {\boldsymbol {e}}(\bar {\vartheta },\bar {\varPsi },\bar {\varphi }) = \boldsymbol {\textsf {M}}(\bar {\vartheta },\bar {\varPsi },\bar {\varphi }) \boldsymbol {e}, \end{align}
noting the appearance of Bessel functions of order one in
$\boldsymbol {\textsf {A}}$
, along with the Bessel functions of order zero in
$\boldsymbol {\textsf {D}}$
that we have already seen appear in the rotational dynamics (though with a slightly different argument). Before we use these results to determine the average translational dynamics, it is instructive to note that
$\tilde {\boldsymbol {e}} = (\tilde {\boldsymbol {e}}_1, \tilde {\boldsymbol {e}}_2, \tilde {\boldsymbol {e}}_3)^{\intercal }$
is defined as the basis transformation
$\boldsymbol {\textsf {M}}$
(defined in (2.3)) applied to the laboratory basis, but evaluated using the slow evolution angles
$(\bar {\vartheta },\bar {\varPsi },\bar {\varphi })$
instead of the rapidly varying full Euler angles
$(\theta , \psi , \phi )$
. That is, (5.6a
) tells us that the average of the leading-order swimmer-fixed basis
$\langle \hat {\boldsymbol {e}} \rangle$
is the slow evolution angle basis, weighted by the diagonal matrix
$\boldsymbol {\textsf {D}}$
. We can therefore, in some sense, interpret
$(\bar {\vartheta },\bar {\varPsi },\bar {\varphi })$
as an ‘average orientation’ of the spheroid over the fast timescale.
Using the averaged results (5.6) in the solvability condition (5.5), we obtain the emergent governing equation for the translational dynamics
\begin{align} \frac {\mathrm {d} \boldsymbol {X}_0}{\mathrm {d} t} = \widehat {\boldsymbol {V}} + Y_0 \boldsymbol {e}_{3}, \end{align}
where the effective self-generated translational velocity
$\widehat {\boldsymbol {V}}$
is defined as
\begin{align} \widehat {\boldsymbol {V}} = \left [a_1 J_0(A) + b_3 J_1(A) \sin \delta _3 \right ] \tilde {\boldsymbol {e}}_1 + a_2 \tilde {\boldsymbol {e}}_2 + \left [a_3 J_0(A) - b_1 J_1(A) \sin \delta _1 \right ] \tilde {\boldsymbol {e}}_3. \end{align}
Importantly, we see that the effective velocity is constant in the averaged swimmer-fixed frame defined through
$\tilde {\boldsymbol {e}}(\bar {\vartheta },\bar {\varPsi },\bar {\varphi })$
. We discuss the implication of these results in § 6.
6. Discussion
We show that a rapidly yawing spheroidal swimmer interacting with a far-field shear flow generates non-axisymmetric emergent effects in its rotational dynamics, equivalent to those of a passive particle with two orthogonal planes of symmetry. With the caveat that the effective coefficients we derive in (3.18) are not independent, our emergent equations are equivalent to those that have been investigated in, e.g., Hinch & Leal (Reference Hinch and Leal1979), Thorp & Lister (Reference Thorp and Lister2019) and Yarin et al. (Reference Yarin, Gottlieb and Roisman1997). From these works, it is known that passive particles with this type of asymmetry can behave very differently to passive spheroids. We demonstrate that the emergent dynamics we derive here can exhibit chaotic behaviours, in stark contrast to passive spheroids.
Given that the effective coefficients we derive in (3.18) are not independent (notably,
$\widehat {B}_1 + \widehat {B}_2 + \widehat {B}_3 = 0$
), the effective shape described by our emergent equations is constrained to particles with two orthogonal planes of symmetry. A basic class of objects with this type of symmetry is ellipsoids. To check when our effective coefficients describe a passive ellipsoid, we use (2.9) to derive the requirement
\begin{align} 0 = \widehat {B}_1 \widehat {B}_2 \widehat {B}_3 + \widehat {B}_1 + \widehat {B}_2 + \widehat {B}_3 = -B^3 J_0(2 A) \left (1 - J_0^2(2 A) \right )/4. \end{align}
Therefore, the equivalent effective particle described by the emergent evolution equations (3.17) is not an ellipsoid in general, unless one of three specific conditions hold: (1)
$A = 0$
(i.e. no yawing), (2)
$B = 0$
(i.e.
$r = 1$
; the original spheroid is a sphere), or (3)
$J_0(2 A) = 0$
. There are infinitely many discrete values of
$A$
that generate the non-trivial case (3); for these scenarios we can use the relationship (2.8) between
$B_i$
and the ellipsoid axes to deduce that the effective passive shape becomes a spheroid with axes
$\widehat {a}$
,
$\widehat {b}$
,
$\widehat {b}$
, where
\begin{align} \widehat {a}/\widehat {b} = \sqrt {(r^2 + 3)/(3 r^2 + 1)}. \end{align}
Thus, in case (3) active prolate spheroids behave as passive oblate spheroids and vice versa. Notably, the effective aspect ratio (6.2) is the same as that which arises for a spheroid rapidly and uniformly rotating about an axis perpendicular to its symmetry axis (Dalwadi et al. Reference Dalwadi, Moreau, Gaffney, Ishimoto and Walker2024a
). This can be understood intuitively in the large-
$A$
limit, where
$J_0(2 A) \to 0$
, for which the large-amplitude yawing has a similar effect to uniform rotation. However, we emphasize that case (3) also occurs for the infinitely many finite roots of
$J_0(2 A) = 0$
.
The emergent loss of symmetry here is fundamentally different to the results of recent studies of different types of rapidly moving rigid bodies, e.g. yawing in two dimensions (Walker et al. Reference Walker, Ishimoto, Gaffney, Moreau and Dalwadi2022), and constant rotation in three dimensions (Dalwadi et al. Reference Dalwadi, Moreau, Gaffney, Ishimoto and Walker2024a ,Reference Dalwadi, Moreau, Gaffney, Walker and Ishimoto b ). While these studies do also show that their specific rapid motions in shear flow lead to emergent dynamics, the effective passive shapes generated all preserve the hydrodynamic symmetries of the original physical shapes. Moreover, the equivalence to an effective passive spheroid is generic for periodically shape-deforming swimmers in two dimensions (Gaffney et al. Reference Gaffney, Dalwadi, Moreau, Ishimoto and Walker2022). Our study shows that symmetry of the physical swimmer is not maintained for rapid yawing in three dimensions.
Given the fundamental difference between our 3-D results (3.17) and the generic 2-D behaviour (Gaffney et al. Reference Gaffney, Dalwadi, Moreau, Ishimoto and Walker2022), it is instructive to understand how our results collapse to the 2-D case. By constraining the swimmer pitch and yawing motions to the shear plane (
$\theta = \psi = \pi /2$
in the full dynamics (2.4)) and solely evolving the remaining equation for
$\phi$
, we reduce our set-up to the 2-D yawing problem considered in Walker et al. (Reference Walker, Ishimoto, Gaffney, Moreau and Dalwadi2022). This corresponds to fixing
$\bar {\vartheta } = \bar {\varPsi } = \pi /2$
in the slow variables, under which the remaining emergent equation for in-plane orientation
$\bar {\varphi }$
(3.17c) reduces significantly to
\begin{align} \frac {\mathrm {d} \bar {\varphi }}{\mathrm {d} t} = \dfrac {1}{2} \left (1 - B J_0(2A) \cos 2 \bar {\varphi } \right ). \end{align}
Therefore, restricting motion to the 2-D shear plane means that the active spheroid behaves as a passive spheroid with effective Bretherton parameter
$B J_0(2A)$
, in agreement with the 2-D results of Walker et al. (Reference Walker, Ishimoto, Gaffney, Moreau and Dalwadi2022). Specifically, the emergent asymmetry generated in the full 3-D emergent dynamics (3.17) vanishes in the constrained 2-D dynamics. Hence, we may conclude that the emergent asymmetry that arises is a 3-D effect generated by out-of-shear-plane interactions between the swimmer and the shear flow.
A natural question to ask is why do symmetries appear to be preserved in the emergent dynamics arising from some types of self-generated motion? Intuitively, it seems as though an important factor should be the average shape of a rapidly moving object, with any symmetries therein conserved in the emergent dynamics. For a spheroid rapidly rotating at a constant rate about an axis fixed in the swimmer frame, the average shape is axisymmetric and the emergent dynamics is equivalent to those of a passive axisymmetric object (Dalwadi et al. Reference Dalwadi, Moreau, Gaffney, Ishimoto and Walker2024a
). For the rapidly yawing spheroid considered here, the average shape is not axisymmetric in general. However, the average shape does have two planes of symmetry, and this symmetry is retained in the effective passive shape represented by the emergent dynamics we derive here. Curiously, however, the average shape does not tell the full story; it is not the shape of the equivalent passive object in general. This can be demonstrated by considering the aspect ratio (6.2) of the effective passive spheroid that arises when
$J_0(2A) = 0$
(including the large-
$A$
limit). In this scenario, the effective aspect ratio (6.2) is bounded within
$(1/\sqrt {3},\sqrt {3})$
, no matter how large or small the physical aspect ratio
$r$
, and is therefore different from the average shape in general. While this may seem surprising due to the linearity of the Stokes equations, the difference occurs because the Stokes equations are not linear in geometry. The general nature of the relationship between the average shape of the fast motion and any effective hydrodynamic shape therefore remains an open question.
Given the technical nature of our analysis, it is instructive to consider further the averages (3.13) required to determine the emergent equations. These averages are weighted non-trivially in a manner that is systematically determined by solving the (non-self) adjoint of the non-autonomous first-correction system (3.9), (3.11) in Appendix B. Requiring a technical analysis to determine the appropriate averages to take is not unusual in nonlinear multiple scales problems; since
$\langle ab \rangle \neq \langle a \rangle \langle b \rangle$
in general, intuitive arguments that do not properly account for nonlinearities may break down, and a key question in such problems is often which average one should take. In fact, the averages that arise from our analysis are more straightforward to interpret if we write the original functions
$g_i$
from (2.5) (which represent the slow interaction with the far-field shear flow) in terms of the angular velocity components of the spheroid in the swimmer-fixed frame
$(\hat {\omega }_1, \hat {\omega }_2, \hat {\omega }_3)$
:
\begin{align} g_1 = \hat {\omega }_2 \cos \psi - \hat {\omega }_3 \sin \psi , \qquad g_2 = \hat {\omega }_1 - g_3 \cos \theta , \qquad g_3 = \dfrac {\hat {\omega }_2 \sin \psi + \hat {\omega }_3 \cos \psi }{\sin \theta }. \end{align}
Substituting (6.4) into the right-hand sides of the averages (3.13), we see that the averages are linear combinations of the averaged quantities
\begin{align} \left \lt \hat {\omega }_1 \right \gt , \qquad \left \lt \hat {\omega }_2 \cos \theta - \hat {\omega }_3 \sin \theta \cos \psi \right \gt , \qquad \left \lt \hat {\omega }_2 \sin \theta \cos \psi + \hat {\omega }_3 \cos \theta \right \gt . \end{align}
Therefore, (6.5) provides a more physical interpretation of the averages we have systematically derived via our technical analysis; the appropriate averages are specific combinations of the angular velocity components of the spheroid. Specifically, the component along the symmetry axis (
$\hat {\omega }_1$
) is averaged without modification, but the components perpendicular to this axis (
$\hat {\omega }_2$
and
$\hat {\omega }_3$
) must be weighted in a manner that accounts for the plane of yawing.
We also derived the emergent translational dynamics (5.7) that arise from the combination of rapid yawing with self-generated translation of the spheroid centre of mass. We specifically consider self-generated motion that is oscillatory in a swimmer-fixed frame, with the same period as the yawing. Importantly, the effective translational velocity
$\widehat {\boldsymbol {V}}$
we derive in (5.7b
) is constant in the average orientation basis vectors
$\tilde {\boldsymbol {e}}_i(\bar {\vartheta },\bar {\varPsi },\bar {\varphi })$
fixed in the average swimmer frame. We note that the effective velocity in the direction normal to the yawing plane
$(\widehat {\boldsymbol {V}} \boldsymbol {\cdot } \tilde {\boldsymbol {e}}_2)$
is simply the average of the full translational velocity (5.2) in this direction. Therefore, the emergent translation in the direction normal to the yawing plane is essentially unaffected by yawing (as expected intuitively) and, being independent of
$b_2$
and
$\delta _2$
(the amplitude and phase shift of the oscillation in the direction normal to the yawing plane), also ignores any oscillation of the translational velocity in this direction. However, these properties do not carry over to the effective velocity in the yawing plane.
In the yawing plane, the emergent translational velocity has two key contributions. The first is due to the average of the full translational velocity in this plane, weighted by
$J_0(A)$
. Physically, this is due to the yawing causing the constant translation to generate a curved trajectory in physical space, reducing the magnitude of the net translation over a yawing period. The second is due to the interaction of the yawing motion with the oscillatory part of the full translational velocity in the yawing plane. Importantly, this contribution can only arise if the translational oscillations in the yawing plane are not in phase or antiphase with the yawing oscillation (i.e. it requires the phase shifts in the yawing plane
$\delta _1, \delta _3 \neq 0, \pi$
), as indicated by the dependence of the effective velocity
$\widehat {\boldsymbol {V}}$
on
$\sin \delta _1$
and
$\sin \delta _3$
(5.7b
). We also note that translational oscillations in one direction of the yawing plane generate orthogonal contributions within the yawing plane, and these contributions are weighted by
$J_1(A)$
. Physically, this is because yawing moves the swimmer-fixed basis in the yawing plane orthogonally within the laboratory frame, and so phase differences of translational oscillations in the swimmer frame can lead to orthogonal emergent contributions.
Since we have demonstrated that rapid yawing causes a spheroid to demonstrate rotational dynamics equivalent to those for an effective object with two planes of symmetry, and objects with this type of symmetry can have translation--rotation coupling, it is interesting to note that this coupling does not arise in the emergent translational dynamics here. That is, the rapid yawing by itself does not generate any emergent translation in (5.7). However, if the spheroid self-generates oscillatory self-directed motion then this can combine with the rapid yawing to generate orthogonal effective translation. We note that this effect is independent of the shear flow in the sense that it would still occur for a self-propelling spheroid in a fluid with no externally imposed flow, though it is implicitly affected by the shear flow via the interacting flow effect on the slow evolution of
$(\bar {\vartheta },\bar {\varPsi },\bar {\varphi })$
.
Finally, we note that our results are straightforward to generalize to several other scenarios. Since the translational results do not affect the leading-order rotational dynamics, it is straightforward to incorporate different types of self-translation into our results by calculating the resulting average in (5.5). For the rotational dynamics, the simplest generalization is that our results hold immediately for rapidly yawing general axisymmetric objects, now interpreting the Bretherton parameter as the measure of an effective physical aspect ratio (Bretherton Reference Bretherton1962; Brenner Reference Brenner1964). Our results can also be extended to consider general periodic yawing functions, essentially replacing
$\Omega A \cos (\Omega t )$
in (2.1) with a general periodic function
$\Omega f'(\Omega t)$
. In this case, all our analysis up to and including (3.15) still holds, and the corresponding versions of the slow evolution equations (3.17) can be obtained by simple evaluation of
$\langle s^2_f \rangle$
,
$\langle c^2_f \rangle$
and
$\langle s_f c_f \rangle$
in terms of the integrated function
$f(T)$
(imposing
$f(0) = 0$
, and where
$T = \Omega t$
).
For odd yawing functions, we have
$\langle s_f c_f \rangle = 0$
and the appropriate emergent slow evolution equations are (3.17), replacing
$J_0(2 A)$
with
$\langle e^{2 \mathrm {i} f} \rangle$
. If we additionally have
$\langle e^{2 \mathrm {i} f} \rangle = 0$
then we recover the non-trivial case (3) above; the effective shape again reduces to a spheroid with axes
$\widehat {a}$
,
$\widehat {b}$
,
$\widehat {b}$
and aspect ratio (6.2). In this scenario, the nonlinear transformations
\begin{align} s_{\bar {\vartheta }} s_{\bar {\varPsi }} \mapsto c_{\bar {\vartheta }}, \qquad s_{\bar {\vartheta }} c_{\bar {\varPsi }} \mapsto s_{\bar {\vartheta }} s_{\bar {\varPsi }}, \qquad c_{\bar {\vartheta }} s_{\bar {\varPsi }} s_{\bar {\varphi }} - c_{\bar {\varPsi }} c_{\bar {\varphi }} \mapsto s_{\bar {\vartheta }} s_{\bar {\varphi }} \end{align}
recover Jeffery’s equations directly. That is, they transform (3.17) into (2.10) with
$\widehat {B}_1 = - \widehat {B}_2 = -B/2$
and
$\widehat {B}_3 = 0$
. This observation provides a potential interpretation for the generation of asymmetry. Namely, considering
$e^{\mathrm {i} f} = \xi +\mathrm {i}\eta$
on the complex unit circle,
$\langle e^{2 \mathrm {i} f} \rangle = 0$
corresponds to
$\langle \xi ^2 \rangle = \langle \eta ^2 \rangle$
and
$\langle \xi \eta \rangle = 0$
, i.e. the mean square orientation having no preferred direction. Hence, emergent asymmetry can arise when there is some bias in the preferred mean square orientation of rapid motion.
Acknowledgements.
The author thanks E.A. Gaffney, K. Ishimoto, C. Moreau, and B.J. Walker for useful discussions, and the anonymous referees for their helpful comments.
Conflict of interests.
The author reports no conflict of interest.
Data accessibility.
The code used in this study is available at https://github.com/m-dalwadi/Yawing-spheroids.
Appendix A. Solving the leading-order system (3.4)
In this appendix we solve the nonlinear, non-autonomous leading-order system (3.4). We start by noting that the first two equations in (3.4) decouple from the third and exhibit fast-time conserved quantities. To see this, we divide the second equation by the first to obtain
\begin{align} \frac {\partial \psi _0}{\partial \theta _0} = - \cot \theta _0 \tan \psi _0. \end{align}
Then, we can integrate (A1) directly to deduce that
\begin{align} \sin \theta _0 \sin \psi _0 = \alpha (t), \end{align}
where
$\alpha (t)$
is a (slow-time) function of integration, i.e. a fast-time conserved quantity.
We then substitute (A2) into the third equation in (3.4) to obtain the governing equation
\begin{align} \frac {\partial \theta _0}{\partial T} = A \cos T \dfrac {\sqrt {\sin ^2 \theta _0 - \alpha ^2(t)}}{\sin \theta _0}, \end{align}
which can be rearranged to obtain
\begin{align} \int \dfrac {\sin \theta _0 \, \mathrm {d}\theta _0}{\sqrt {1 - \alpha ^2(t) - \cos ^2 \theta _0}} = \int A \cos T \, \mathrm {d}T = A \sin T + \beta (t), \end{align}
where
$\beta (t)$
is the second of three functions of integration. The integral on the left-hand side of (A4) can be calculated by direct substitution of
$\cos \theta _0 = \sqrt {1 - \alpha ^2} \cos u$
, yielding the solution
\begin{align} \cos \theta _0 = \sqrt {1 - \alpha ^2(t)} \cos (A \sin T + \beta (t)). \end{align}
Finally, we can solve the remaining leading-order equation by substituting (A2) and (A5) into the third equation of (3.4) to obtain
\begin{align} \frac {\partial \phi _0}{\partial T} = \dfrac {\alpha (t) A \cos T}{1 - (1 - \alpha ^2(t)) \cos ^2(A \sin T + \beta (t))}. \end{align}
We can integrate (A6) directly using the substitution
$\tan (A \sin T + \beta ) = \alpha \tan u$
, resulting in the solution
\begin{align} \alpha (t) \tan (\phi _0 - \gamma (t)) = \tan (A \sin T + \beta (t)), \end{align}
where
$\gamma (t)$
is the third and final function of integration from the leading-order solution.
The leading-order solutions (A2), (A5) and (A7) are the general solutions to the nonlinear, non-autonomous leading-order system (3.4). These are equivalent to (3.5) after appropriately redefining the functions of integration
$(\alpha (t), \beta (t), \gamma (t))$
into
$(\bar {\vartheta }(t), \bar {\varPsi }(t), \bar {\varphi }(t))$
. We redefine these fast-time conserved quantities so that the new functions of integration are equivalent to
$(\theta , \psi , \phi )$
in the limit
$A \to 0$
, in which case we expect to recover the original Jeffery’s equations ((2.4) with
$A = 0$
) with no fast-time variation. Specifically, we use the following transformations:
\begin{align} \alpha (t) &= \sin \bar {\vartheta }(t) \sin \bar {\varPsi }(t), \end{align}
\begin{align} \tan \beta (t) &= \tan \bar {\vartheta }(t) \cos \bar {\varPsi }(t), \end{align}
\begin{align} \tan \gamma (t) &= \dfrac {\cos \bar {\vartheta }(t) \tan \bar {\varPsi }(t) \tan \bar {\varphi }(t) - 1}{\cos \bar {\vartheta }(t) \tan \bar {\varPsi }(t) + \tan \bar {\varphi }(t)}. \end{align}
Appendix B. Solving the adjoint system (3.9), (3.11)
In this appendix we solve the linear, non-autonomous adjoint system (3.9), (3.11) for
$\boldsymbol {Y} = (P, Q, R)$
with periodic boundary conditions. This solution will allow us to generate the appropriate solvability conditions, and hence, emergent equations, via (3.10). We note that the solution method we present here is a modified version of that in Dalwadi et al. (Reference Dalwadi, Moreau, Gaffney, Ishimoto and Walker2024a
) for a similar (but autonomous) adjoint system.
Since we have a 3-D linear system, we seek three linearly independent solutions. We start by considering the bottom row of
$L^*$
in (3.11), which tells us that
$R = \text {constant}$
in all solutions. We can reduce our task to solving a 2-D system by setting this constant equal to zero in two of the linearly independent solutions, ensuring that the third solution will be able to generate our solution basis for
$R$
. For this third solution, we explicitly write
\begin{align} R = C_3, \end{align}
where
$C_3$
is an arbitrary constant. Substituting (B1) into (3.9), (3.11), we obtain
\begin{align} \frac {\mathrm {d} P}{\mathrm {d} T} &= \dfrac {A \cos T \sin \psi _0}{\sin ^2 \theta _0} \left [- Q + C_3 \cos \theta _0 \right ], \end{align}
\begin{align} \frac {\mathrm {d} Q}{\mathrm {d} T} &= A \cos T \left [P \sin \psi _0 + \dfrac {\cos \psi _0}{\sin \theta _0} \left (\cos \theta _0 Q - C_3 \right ) \right ]. \end{align}
Given that the right-hand side of (B2a
) does not depend on
$P$
, it is convenient to introduce
$Q = C_3 \cos \theta _0 + \tilde {Q}$
to reduce the complexity of the system. In this case, the time derivative generated by differentiating
$\cos \theta _0$
in (B2b
) cancels exactly with the inhomogeneous term on the right-hand side, and so (B2) becomes
\begin{align} \frac {\mathrm {d} P}{\mathrm {d} T} &= -\dfrac {A \cos T \sin \psi _0 }{\sin ^2 \theta _0} \tilde {Q}, \end{align}
\begin{align} \frac {\mathrm {d} \tilde {Q}}{\mathrm {d} T} &= A \cos T \left [P \sin \psi _0 + \dfrac {\cos \psi _0 \cos \theta _0}{\sin \theta _0} \tilde {Q} \right ]. \end{align}
That is, we have transformed (B2) into a homogeneous system, which has trivial solution
$(P, \tilde {Q}) = (0, 0)$
. Hence, we have generated one non-trivial solution to (3.9), (3.11), namely
\begin{align} \boldsymbol {Y} = C_3 (0, \cos \theta _0, 1)^{\intercal }. \end{align}
As noted above, the remaining two solutions can be determined by now setting
$C_3 = 0$
in (B2) to obtain the system
\begin{align} \frac {\mathrm {d} P}{\mathrm {d} T} &= -\dfrac {A \cos T \sin \psi _0}{\sin ^2 \theta _0}Q, \end{align}
\begin{align} \frac {\mathrm {d} Q}{\mathrm {d} T} &= A \cos T \left [P \sin \psi _0 + \dfrac {\cos \psi _0 \cos \theta _0}{\sin \theta _0} Q \right ]. \end{align}
To derive two linearly independent solutions to (B5), we seek nonlinear transformations to map the remaining adjoint problem (B5) into a homogeneous version of the first-correction system (3.8). The reason this is helpful is because the linear operator of (3.8) is a perturbed version of the (nonlinear) operator of the leading-order system (3.4). Hence, perturbed versions of the solutions we derived to the leading-order system (3.5) will solve the homogeneous version of the first-correction system (3.8). That is,
\begin{align} \theta _1 = \dfrac {1}{A \cos T}\frac {\partial \theta _0}{\partial T} = \cos \psi _0, \qquad \psi _1 = \dfrac {1}{A \cos T}\frac {\partial \psi _0}{\partial T} = - \dfrac {\cos \theta _0 \sin \psi _0}{\sin \theta _0} \end{align}
solve the homogeneous version of the first-correction system (3.8a ), (3.8b ).
To determine how to map (B5) into a homogeneous version of (3.8), we introduce the transformations
\begin{align} P(T) = \zeta _1(\theta _0) \tilde {P}(T), \quad Q(T) = \zeta _2(\theta _0) \tilde {Q}(T), \end{align}
where
$\zeta _i$
are the (as-of-yet) unknown nonlinear functions of
$\theta _0$
we seek. Substituting (B7) into (B5), we obtain
\begin{align} \frac {\mathrm {d} \tilde {P}}{\mathrm {d} T} &= -A \cos T \left [\dfrac {\zeta '_1}{\zeta _1} \cos \psi _0 \tilde {P} + \dfrac {\zeta _2}{\zeta _1} \dfrac { \sin \psi _0}{\sin ^2 \theta _0} \tilde {Q} \right ], \end{align}
\begin{align} \frac {\mathrm {d} \tilde {Q}}{\mathrm {d} T} &= A \cos T \left [\dfrac {\zeta _1}{\zeta _2} \sin \psi _0 \tilde {P} + \cos \psi _0 \left ( \dfrac { \cos \theta _0}{\sin \theta _0} - \dfrac {\zeta '_2}{\zeta _2}\right ) Q \right ]. \end{align}
Then, we compare (B8) to the homogeneous version of the first-correction system (3.8a
), (3.8b
). We see that
$(\tilde {P}, \tilde {Q}) = (\theta _1, \psi _1)$
when
\begin{align} \dfrac {\zeta '_1}{\zeta _1} = 0, \qquad \dfrac {\zeta '_2}{\zeta _2} = \dfrac {2 \cos \theta _0}{\sin \theta _0}, \qquad \dfrac {\zeta _2}{\zeta _1} = \sin ^2 \theta _0. \end{align}
That is, when
\begin{align} \zeta _1 = C_2, \qquad \zeta _2 = C_2 \sin ^2 \theta _0, \end{align}
where
$C_2$
is an arbitrary constant. Then substituting
$(\tilde {P}, \tilde {Q}) = (\theta _1, \psi _1)$
into (B7) and using the results (B6) and (B10) gives us the second linearly independent solution:
\begin{align} \boldsymbol {Y} = C_2 (\cos \psi _0, -\sin \theta _0 \cos \theta _0 \sin \psi _0, 0)^{\intercal }. \end{align}
The third and final linearly independent solution to the adjoint problem (3.9), (3.11) arises from comparing (B8) to the homogeneous version of the first-correction system (3.8a
), (3.8b
) and now looking for
$(\tilde {P}, \tilde {Q}) = (\psi _1, \theta _1)$
. This requires
\begin{align} \dfrac {\zeta '_1}{\zeta _1} = 0, \qquad \dfrac {\zeta '_2}{\zeta _2} = 2 \dfrac {\cos \theta _0}{\sin \theta _0}, \qquad \dfrac {\zeta _2}{\zeta _1} = \sin ^2 \theta _0, \end{align}
which we can solve straightforwardly to obtain
\begin{align} \zeta _1 = C_1 \sin \theta _0, \qquad \zeta _2 = -C_1 \sin \theta _0. \end{align}
Finally, substituting
$(\tilde {P}, \tilde {Q}) = (\psi _1, \theta _1)$
into (B7) and using the results (B6) and (B13) gives us the final linearly independent solution:
\begin{align} \boldsymbol {Y} = C_1 (\sin \theta _0 \cos \psi _0, \cos \theta _0 \sin \psi _0, 0)^{\intercal }. \end{align}
Together, the three linearly independent solutions (B4), (B11), (B14) to the adjoint problem (3.9), (3.11) give the general solution (3.12).
Open access
