1. Introduction
The translocation of molecular motors attached to biological filaments, such as microtubule or actin filaments, can give rise to rich dynamics and drive flows, as in ciliary beating or cytoplasmic streaming (Shelley Reference Shelley2016; Lauga Reference Lauga2020; Stein et al. Reference Stein, De Canio, Lauga, Shelley and Goldstein2021). In particular, such active filaments may buckle under sufficiently strong compressive forces exerted by the motors. In that scenario, the tangential nature of the loading suggests an oscillatory buckling instability leading to spontaneous dynamics. This is in contrast to the classical Euler buckling of a passive rod under a fixed load, where a monotonic instability leads to steady deformation.
De Canio et al. (Reference De Canio, Lauga and Goldstein2017) introduced a fundamental model consisting of a slender elastic filament clamped normal to a wall at one end and subjected to a tangential (compressive) ‘follower’ force at the other, which is meant to phenomenologically represent the force exerted on the filament by a motor walking towards its tip. (They also considered a generalised model purporting to account for the opposite force applied by the motor on the fluid.) De Canio et al. (Reference De Canio, Lauga and Goldstein2017) modelled the filament as an inertialess Kirchhoff rod (Landau et al. Reference Landau, Pitaevskii, Kosevich and Lifshitz2012), assumed that the filament motion is restricted to a plane and approximated the fluid–structure interaction according to ‘resistive force theory’, a local anisotropic drag law valid at zero Reynolds number and to leading order in the slender-filament limit (Lauga Reference Lauga2020). By performing a linear stability analysis, De Canio et al. (Reference De Canio, Lauga and Goldstein2017) showed that their model filament undergoes symmetry breaking via a supercritical Hopf bifurcation; that is, the undeformed configuration of the filament becomes unstable at a critical follower-force magnitude, above which small perturbations grow in an oscillatory fashion. In that regime, De Canio et al. (Reference De Canio, Lauga and Goldstein2017) demonstrated via initial-value simulations the approach to a limit cycle at late times, corresponding to the filament exhibiting periodic ‘beating’ oscillations.
The planar-motion assumption was subsequently relaxed by Ling et al. (Reference Ling, Guo and Kanso2018), who considered a similar set-up to that of De Canio et al. (Reference De Canio, Lauga and Goldstein2017) but allowing for three-dimensional deformations and generalising from a localised follower force at the tip to arbitrarily distributed axial loads. Revisiting the case of a localised follower force, Ling et al. (Reference Ling, Guo and Kanso2018) strikingly demonstrated through initial-value simulations that in a follower-force interval above the stability threshold the filament generally exhibits (at late times) non-planar ‘whirling’ oscillations – periodic states for which a steadily rotating frame exists in which the deformed filament appears fixed; in that interval, the beating states observed by De Canio et al. (Reference De Canio, Lauga and Goldstein2017) could only be attained by choosing the initial conditions to be precisely planar. At higher values of the follower force, i.e. away from the instability threshold, the simulations of Ling et al. (Reference Ling, Guo and Kanso2018) uncovered a secondary transition to a regime where planar beating is, in fact, preferred. Distributing the load was shown to modify the nature of that secondary transition. Their simulations also uncovered additional regimes exhibiting yet right dynamics, including chaotic-like phenomena.
To investigate the onset of spontaneous beating and whirling, Ling et al. (Reference Ling, Guo and Kanso2018) carried out a linear stability analysis allowing for general three-dimensional perturbations from the undeformed filament configuration. Similarly to De Canio et al. (Reference De Canio, Lauga and Goldstein2017), however, they found that the neutral modes at the instability threshold describe planar-beating oscillations of the filament. On that basis, Ling et al. (Reference Ling, Guo and Kanso2018) concluded that‘in the linear regime near the straight equilibrium, the filament undergoes planar deformations’, thus leaving unexplained the emergence of whirling and its preference over planar beating.
Partially motivated by this gap in understanding, Clarke et al. (Reference Clarke, Hwang and Keaveny2024) have recently carried out a thorough computational study of a generalised follower force model – accounting for Stokes hydrodynamics fully rather than through resistive force theory. Upon revisiting the linear stability of the undeformed configuration, Clarke et al. (Reference Clarke, Hwang and Keaveny2024) made the following key observations. While the space of neutral linear modes at the instability threshold can indeed be spanned by planar-beating modes (similar in nature to those found by De Canio et al. (Reference De Canio, Lauga and Goldstein2017) and Ling et al. (Reference Ling, Guo and Kanso2018)), their superposition allowing for differences in phase and the plane of motion generally produces non-planar whirling-like modes. (The false conclusion of Ling et al. (Reference Ling, Guo and Kanso2018) that the filament is restricted to planar deformations in the linear regime may accordingly be attributed to a misinterpretation of their linear analysis.) In light of symmetry, the linear beating modes are identical up to a rotation about the axis of the undeformed configuration (as well as magnitude and phase differences). The onset of instability thus occurs via a degenerate double-Hopf bifurcation, where two identical pairs of complex-conjugate eigenvalues become unstable simultaneously.
Employing methods of computational dynamical systems, Clarke et al. (Reference Clarke, Hwang and Keaveny2024) went beyond a linear stability analysis of the undeformed configuration and initial-value simulations – by computing and tracing the solution branches of periodic and quasi-periodic states, as well as determining the stability of periodic states via Floquet analysis. Their results reveal that two families of periodic states of the nonlinear problem, planar-beating states and whirling states, simultaneously emerge at the double-Hopf bifurcation. In an interval of the follower force above the threshold, the whirling states are stable while the planar-beating states are unstable under out-of-plane perturbations. Richer dynamical phenomena are uncovered at higher values of the follower force; in particular, the secondary transition to a regime of stable beating was shown to occur via a branch of quasi-periodic solutions connecting the whirling and beating solution branches.
Complementary to the primarily numerical studies mentioned above, we here employ the method of multiple scales to carry out a weakly nonlinear analysis of the follower force model near its instability threshold. The main outcome is a reduced-order model in the form of a nonlinear ‘amplitude equation’ governing the slow dynamics of small-magnitude oscillatory perturbations near the threshold. This amplitude equation is of higher dimension than the standard Stuart–Landau equation corresponding to a conventional (non-degenerate) Hopf instability (Drazin & Reid Reference Drazin and Reid2004). We analyse and illustrate our reduced-order model in order to derive new insights into the onset of spontaneous dynamics in the follower force model. We note that nonlinear amplitude equations have been employed to study various other scenarios of elastic filaments buckling in a viscous fluid, such as the onset of planar buckling of a clamped filament in a stagnation flow (Guglielmini et al. Reference Guglielmini, Kushwaha, Shaqfeh and Stone2012) or the early-time buckling of a highly flexible filament in a compressional flow (Chakrabarti et al. Reference Chakrabarti, Liu, LaGrone, Cortez, Fauci, du Roure, Saintillan and Lindner2020).
We shall consider the follower force model in its most elementary form, though allowing for non-planar deformations which is clearly crucial in light of the above discussion. We accordingly follow De Canio et al. (Reference De Canio, Lauga and Goldstein2017) and Ling et al. (Reference Ling, Guo and Kanso2018) by modelling the hydrodynamics via resistive force theory; follow De Canio et al. (Reference De Canio, Lauga and Goldstein2017) and Clarke et al. (Reference Clarke, Hwang and Keaveny2024) by assuming a follower force localised at the tip; and follow Ling et al. (Reference Ling, Guo and Kanso2018) and Clarke et al. (Reference Clarke, Hwang and Keaveny2024) by disregarding the biologically motivated entrainment flow included by De Canio et al. (Reference De Canio, Lauga and Goldstein2017). The preceding computational studies suggest that these modelling choices do not qualitatively influence the dynamics of a single active filament in the vicinity of the instability threshold.
The paper proceeds as follows. In § 2, we formulate the problem. In § 3, we revisit the linear stability of the undeformed configuration and the general linear approximation at the instability threshold. In § 4, we present, analyse and illustrate the reduced-order model, as well as validate the theory against numerical data provided by Eric E. Keaveny. The details of the weakly nonlinear analysis leading to the amplitude equation are retrospectively presented in § 5. We conclude in § 6 by discussing several future directions.
2. Problem formulation
2.1. Elastohydrodynamic model
As shown in figure 1, we consider an inextensible elastic filament of length
$L_*$
and bending stiffness
$B_*$
clamped normally to a flat wall at one end and subjected to a tangential compressive force of magnitude
$\mathcal {F}_*$
at the other. The filament is immersed in a background liquid of viscosity
$\eta _*$
. We neglect inertia of the fluid and the filament and model the latter as a homogeneous Kirchhoff rod having a circular cross-section of radius
$\kappa L_*$
(Landau et al. Reference Landau, Pitaevskii, Kosevich and Lifshitz2012), further assuming that it is sufficiently slender such that the hydrodynamics can be described by resistive force theory (Graham Reference Graham2018). Throughout the paper, a subscript asterisk indicates a dimensional quantity; a superscript asterisk denotes complex conjugation; and
$r$
and
$i$
subscripts indicate the real and imaginary parts of a complex-valued quantity, respectively.

Figure 1. (a) Dimensional schematic. (b) Forces and moments acting on an infinitesimal filament segment.
We represent the filament by its centreline
${\boldsymbol {r}}_*(s_*,t_*)$
, where
$s_*$
is the arc length measured from the clamped end and
$t_*$
denotes the time. Since the filament is inextensible, the tangent vector
${\hat {\boldsymbol {t}}}=\partial {\boldsymbol {r}}_*/\partial s_*$
is constrained to be a unit vector. Let
${\boldsymbol {F}}_*(s_*,t_*)$
and
${\boldsymbol {M}}_*(s_*,t_*)$
be the cross-sectionally averaged internal (elastic) force and moment, respectively; with inertia neglected, they satisfy the equilibrium balances

in which

is the hydrodynamic force per unit length acting on the filament (see figure 1
b). The local anisotropic drag law (2.2) constitutes a leading-order asymptotic approximation of the hydrodynamics in the slender-filament limit
$\kappa \to 0$
; under that approximation, known as resistive force theory, the interaction of the filament with itself and with the wall is neglected despite the relative error associated with that approximation being only logarithmically small, scaling as
$1/\ln (1/\kappa )$
(Graham Reference Graham2018).
The above governing equations are supplemented by initial conditions (which we do not specify at this stage), boundary conditions and a constitutive relation for the internal moment. The boundary conditions at the clamped end of the filament are

where
$\boldsymbol {\hat {\imath }}$
is the unit vector normal to the wall pointing into the fluid and we choose the reference point for the position vector
${\boldsymbol {r}}_*$
to coincide with the clamping point. The boundary conditions at the other end of the filament are

the first prescribing the follower force and the second stating that there is no external moment acting at the tip. Lastly, the constitutive law for the internal moment can be written as (Landau et al. Reference Landau, Pitaevskii, Kosevich and Lifshitz2012)

representing the resistance of the filament to bending, which is cross-sectionally isotropic given the assumption of a circular cross-section. In general, the constitutive relation includes an additional term representing the resistance of the filament to twisting. In the present scenario, however, where there are no external moments acting in the tangential direction and the bending is cross-sectionally isotropic, the twist term in the constitutive relation vanishes identically; this is verified in Appendix A following Landau et al. (Reference Landau, Pitaevskii, Kosevich and Lifshitz2012). The absence of twist in this scenario was also remarked on by Ling et al. (Reference Ling, Guo and Kanso2018). In contrast, Clarke et al. (Reference Clarke, Hwang and Keaveny2024) include twist as their model does not involve a thin-flament approximation of the hydrodynamics, so that viscous moments act in the tangential direction.
2.2. Dimensionless formulation
It is convenient to adopt a dimensionless convention where lengths are normalised by
$L_*$
, forces by
$B_*/L_*^2$
, moments by
$B_*/L_*$
and time by
$4\pi \eta _*(L_*^4/B_*)/\ln (1/\kappa )$
. Dimensionless fields are denoted similarly to their dimensional counterparts, only with the subscript asterisks omitted. The problem thus consists of the partial differential equations

the inextensibility constraint

and the boundary conditions

and

wherein the dimensionless parameter

measures the magnitude of the follower force. The problem depends on this sole parameter, aside for the initial conditions which we do not specify at this stage.
For all
$\mathcal {F}$
, the above problem possesses the steady solution

where the filament is undeformed (straight) and under uniform compression by the follower force. Our aim will be to study the nonlinear dynamics of the filament for
$\mathcal {F}$
near the buckling threshold where this steady ‘base state’ first becomes unstable.
3. Linear theory
3.1. Eigenvalue problem
As a preliminary step, we must characterise the loss of stability of the base state. To this end, we assume small-magnitude perturbations of the form

where
$\lambda$
is a complex growth rate, the tilde-decorated vector fields are complex functions of
$s$
and
$\text {c.c.}$
stands for complex conjugate. Substituting (3.1) into (2.6) yields, upon linearisation, the set of ordinary differential equations

where henceforth a prime denotes an ordinary derivative with respect to
$s$
(e.g.
${\boldsymbol{\tilde{r}}}^{\prime}=\mathrm {d}{\boldsymbol{\tilde{r}}}/\mathrm {d}s$
). Furthermore, linearisation of the constraint (2.7) gives

and linearisation of the boundary conditions (2.8) and (2.9) gives

It is readily seen that the tilde-decorated perturbations are all perpendicular to
$\boldsymbol {\hat {\imath }}$
. In particular, let
${\boldsymbol{\tilde{r}}}={\hat {\boldsymbol {e}}}_1\tilde {x}+{\hat {\boldsymbol {e}}}_2\tilde {y}$
, wherein
$\tilde {x}(s)$
and
$\tilde {y}(s)$
are complex-valued functions and
${\hat {\boldsymbol {e}}}_1$
and
${\hat {\boldsymbol {e}}}_2$
are unit vectors such that
$\{{\hat {\boldsymbol {e}}}_1,{\hat {\boldsymbol {e}}}_2,\boldsymbol {\hat {\imath }}\}$
constitutes a right-handed orthogonal system. Substituting this Cartesian decomposition into (3.2) and (3.4), we find that
$\tilde {x}$
and
$\tilde {y}$
are decoupled, separately satisfying the same eigenvalue problem. In terms of
$\tilde {x}$
, say, that problem consists of the ordinary differential equation

and the boundary conditions


Figure 2. (a) Complex growth rate
$\lambda =\lambda _r+i\lambda _i$
as a function of
$\mathcal {F}$
for the most unstable eigenfunctions of the eigenvalue problem consisting of (3.5) and (3.6). (b) The threshold eigenfunction
$\varphi (s)=\varphi _r(s)+i\varphi _i(s)$
corresponding to the imaginary growth rate
$\lambda =i\omega$
.
The eigenvalue problem consisting of (3.5) and (3.6) was first derived and solved by De Canio et al. (Reference De Canio, Lauga and Goldstein2017) under the assumption of planar deformations. While the general solution to (3.5) is readily expressed in closed form, calculating the eigenvalues
$\lambda$
as a function of
$\mathcal {F}$
must, in general, be done numerically. In agreement with De Canio et al. (Reference De Canio, Lauga and Goldstein2017), we thereby find that the base state is monotonically stable for
$\mathcal {F}\lt \mathcal {F}_{o}\approx 20.051$
, oscillatory stable for
$\mathcal {F}_0\lt \mathcal {F}\lt \mathcal {F}_c\approx 37.695$
and oscillatory unstable for
$\mathcal {F}\gt \mathcal {F}_c$
; at the instability threshold,
$\mathcal {F}=\mathcal {F}_c$
, the complex-conjugate pair of neutral eigenvalues is
$\lambda =\pm i\omega$
, wherein
$\omega \approx 191.26$
(see figure 2
a). We denote by
$\varphi (s)=\varphi _r(s)+i\varphi _i(s)$
the eigenfunction normalised to unity at
$s=1$
that is associated with the threshold eigenvalue
$i\omega$
(see figure 2
b);
$\varphi ^*(s)$
is accordingly the similarly normalised eigenfunction associated with the conjugate eigenvalue
$-i\omega$
.
3.2. Double-Hopf bifurcation
The eigenfunctions at the instability threshold describe periodic beating oscillations of the filament in the plane through the origin and normal to
$\boldsymbol {\hat {\imath }}\times {\hat {\boldsymbol {e}}}_1$
, of arbitrary magnitude and phase. Since
$\tilde {y}(s)$
satisfies the same eigenvalue problem as
$\tilde {x}(s)$
, there is clearly an independent pair of eigenfunctions representing similar beating oscillations rotated by
$\pi /2$
about the undeformed configurations, i.e. in the plane through the origin and normal to
$\boldsymbol {\hat {\imath }}\times {\hat {\boldsymbol {e}}}_2$
. As first pointed out by Clarke et al. (Reference Clarke, Hwang and Keaveny2024), this means that the filament loses stability via a double-Hopf bifurcation, where two pairs of complex-conjugate eigenvalues become unstable simultaneously; since the multiplicity here is induced by symmetry, the pairs are identical and have associated with them a full set of eigenfunctions. A key observation of Clarke et al. (Reference Clarke, Hwang and Keaveny2024) is that the linear superposition of beating modes allowing for differences in the phase and plane of motion generally produces non-planar whirling oscillations of the filament.
3.3. Linear approximation at the threshold
It will be convenient for us to employ complex vectors, whose real and imaginary parts are conventional 3D vectors. Let
${\boldsymbol {a}}={\boldsymbol {a}}_r+i{\boldsymbol {a}}_i$
and
${\boldsymbol {b}}={\boldsymbol {b}}_r+i{\boldsymbol {b}}_i$
be two such vectors. The dot and cross products are generalised from the case of real vectors as
${\boldsymbol {a}}\boldsymbol {\cdot }{\boldsymbol {b}}={\boldsymbol {a}}_r\boldsymbol {\cdot }{\boldsymbol {b}}_r-{\boldsymbol {a}}_i\boldsymbol {\cdot }{\boldsymbol {b}}_i+i({\boldsymbol {a}}_r\boldsymbol {\cdot }{\boldsymbol {b}}_i+{\boldsymbol {a}}_i\boldsymbol {\cdot }{\boldsymbol {b}}_r)$
and
${\boldsymbol {a}}\times {\boldsymbol {b}}={\boldsymbol {a}}_r\times {\boldsymbol {b}}_r-{\boldsymbol {a}}_i\times {\boldsymbol {b}}_i+i({\boldsymbol {a}}_r\times {\boldsymbol {b}}_i+{\boldsymbol {a}}_i\times {\boldsymbol {b}}_r)$
. An appropriate inner product is
$\langle {\boldsymbol {a}},{\boldsymbol {b}}\rangle ={\boldsymbol {a}}^*\boldsymbol {\cdot }{\boldsymbol {b}}$
, which differs from the dot product. The magnitude of a complex vector
$\boldsymbol {a}$
is accordingly defined as
$|{\boldsymbol {a}}|=\sqrt {{\boldsymbol {a}}^*\boldsymbol {\cdot }{\boldsymbol {a}}}=\sqrt {|{\boldsymbol {a}}_r|^2+|{\boldsymbol {a}}_i|^2}$
.
Linear whirling oscillations formed by the superposition of planar-beating modes, can be represented by a complex-vector amplitude
$\tilde{\boldsymbol{A}}$
parallel to the wall. Indeed, noting from (3.2) that an eigenfunction
${\boldsymbol{\tilde{r}}}={\hat {\boldsymbol {e}}}_1\varphi$
is associated with perturbations
$\tilde {\boldsymbol {t}}={\hat {\boldsymbol {e}}}_1\varphi^{\prime}$
,
$\tilde {{\boldsymbol {F}}}=-{\hat {\boldsymbol {e}}}_1(\varphi ^{\prime\prime\prime}+\mathcal {F}_c\varphi^{\prime})$
and
$\tilde {{\boldsymbol {M}}}=\boldsymbol {\hat {\imath }}\times {\hat {\boldsymbol {e}}}_1\varphi ^{\prime\prime}$
, a general linear approximation at the threshold
$\mathcal {F}=\mathcal {F}_c$
can be expressed as (cf. (3.1))




where we have ignored decaying stable modes. Given the normalisation
$\varphi (1)=1$
, the corresponding tip trajectory is

which traces a general ellipse parallel to the wall and centred about the nominal tip position, the vector
$\mathrm {Im}(\,\tilde {\!{\boldsymbol {A}}}\times \,\tilde {\!{\boldsymbol {A}}}^*)$
indicating the direction of motion according to the right-hand rule. It is easy to see that for
$\,\tilde {\!{\boldsymbol {A}}}\times \,\tilde {\!{\boldsymbol {A}}}^*={\mathbf {0}}$
, the ellipse degenerates to a line (parallel to
$\,\tilde {\!{\boldsymbol {A}}}_r$
and
$\,\tilde {\!{\boldsymbol {A}}}_i$
), of length
$4|\,\tilde {\!{\boldsymbol {A}}}|$
. For
$\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}=0$
, the ellipse constitutes a circle of radius
$\sqrt {2}|\,\tilde {\!{\boldsymbol {A}}}|$
.
The properties of a general tip ellipse can be derived following Lindell (Reference Lindell1983). The lengths of the semi-major and semi-minor axes can be expressed as

respectively, giving the ellipse area as
$2\pi |\,\tilde {\!{\boldsymbol {A}}}\times \,\tilde {\!{\boldsymbol {A}}}^*|$
. The above-mentioned condition for a circular orbit,
$\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}=0$
, can be inferred from (3.9) by noting the identity
$|\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}|^2=|\,\tilde {\!{\boldsymbol {A}}}|^4-|\,\tilde {\!{\boldsymbol {A}}}\times \,\tilde {\!{\boldsymbol {A}}}^*|^2$
. For non-circular ellipses, the real and imaginary parts of
$\,\tilde {\!{\boldsymbol {A}}}/(\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}})^{1/2}$
point in the directions of the major and minor axes, respectively (both the forward and backward axis directions, given the multiplicity of the square-root function). These formulae can alternatively be expressed in terms of the real vectors
$\,\tilde {\!{\boldsymbol {A}}}_r$
and
$\,\tilde {\!{\boldsymbol {A}}}_i$
, using the identities
$\,\tilde {\!{\boldsymbol {A}}}\times \,\tilde {\!{\boldsymbol {A}}}^*=-2i\,\tilde {\!{\boldsymbol {A}}}_r\times \,\tilde {\!{\boldsymbol {A}}}_i$
and
$\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}=|\,\tilde {\!{\boldsymbol {A}}}_r|^2-|\,\tilde {\!{\boldsymbol {A}}}_i|^2+2i\,\tilde {\!{\boldsymbol {A}}}_r\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}_i$
. Hence, the ellipse is a line if
$\,\tilde {\!{\boldsymbol {A}}}_r\times \,\tilde {\!{\boldsymbol {A}}}_i={\mathbf {0}}$
; it is a circle if both
$\,\tilde {\!{\boldsymbol {A}}}_r\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}_i=0$
and
$|\,\tilde {\!{\boldsymbol {A}}}_r|=|\,\tilde {\!{\boldsymbol {A}}}_i|$
.
The trajectory
$c\,\tilde {\!{\boldsymbol {A}}}e^{i\omega t}+\text {c.c.}$
, wherein
$c$
is a complex scalar, is identical to the elliptical tip trajectory (3.8) up to a dilation
$|c|$
and phase
$\angle c$
. It thus follows from (3.7a
) that, in fact, all points along the filament centreline carry out geometrically similar elliptical trajectories, obtained from the tip trajectory by a dilation
$|\varphi (s)|$
and phase shift
$\angle \varphi (s)$
. From (3.7
b–d), the tangent unit vector, force and moment perturbations along the filament can be analogously constructed by appropriate scale and phase modulations of the elliptical tip trajectory (as well as a
$\pi /2$
rotation in the case of the moment perturbation).
It can be confirmed from the above formulae that in the special case of circular whirling,
$\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}=0$
, any given point along the filament performs a steady circular motion parallel to the wall. Furthermore, in that case the geometry of the deformed filament and the force and moment distributions along it appear fixed in a frame of reference rotating about the undeformed configuration at a steady angular frequency
$\omega$
.
4. Weakly nonlinear theory
4.1. Separation of time scales near the instability threshold
According to the linear theory, at the instability threshold the filament approaches periodic beating/whirling oscillations of the form (3.7), with the complex-vector amplitude
$\,\tilde {\!{\boldsymbol {A}}}$
, assumed small in magnitude, fixed by the initial conditions. In fact, the fate of small-magnitude oscillations at the instability threshold is actually governed by weak nonlinear effects which become important at late times – invalidating the linear theory! Specifically, nonlinear terms formed of cubic products of the linear approximation (or its derivatives) may resonantly excite the neutral linear modes, thereby producing perturbations growing like
$t|\,\tilde {\!{\boldsymbol {A}}}|^3$
. Over
$O(1/|\,\tilde {\!{\boldsymbol {A}}}|^{2})$
times – long compared with the natural period
$2\pi /\omega$
– such nonlinear perturbations would become comparable to the linear approximation and may accordingly influence the dynamics.
On the other hand, the linear theory predicts that small-magnitude perturbations exponentially decay for
$\mathcal {F}\lt \mathcal {F}_c$
, or exponentially grow, in general, for
$\mathcal {F}\gt \mathcal {F}_c$
. As
$\mathcal {F}\to \mathcal {F}_c$
, the growth rate
$\lambda _r$
of the most unstable linear modes vanishes like
$\mathcal {F}-\mathcal {F}_c$
(see figure 2
a), suggesting that following a transient (during which more stable modes decay) the filament exhibits oscillations of the form (3.7) – modulated over long
$O(1/|\mathcal {F}-\mathcal {F}_c|)$
times. We accordingly identify the distinguished scaling
$|\,\tilde {\!{\boldsymbol {A}}}|=\text {ord}(|\mathcal {F}-\mathcal {F}_c|^{1/2})$
– conventional to symmetry-breaking bifurcations – such that linear and nonlinear effects both have a leading-order effect. For relatively larger amplitudes, or at the threshold, nonlinear effects dominate the dynamics; for relatively smaller amplitudes, linear growth or decay is instead dominant.
The time-scale separation at and near the instability threshold can be exploited towards systematically deriving a reduced-order nonlinear model governing the dynamics of small-magnitude oscillations in this regime. This is carried out in § 5 by means of a weakly nonlinear analysis of the governing equations employing the method of multiple scales; while conceptually standard, the derivation is presented fully so as to highlight some technical details peculiar to the present set-up and, more importantly, to facilitate future extensions to the theory as discussed in § 6. In the present section, we simply quote the reduced-order model and subsequently analyse and illustrate it in various scenarios.
4.2. Amplitude equation
The reduced-order model consists of the linear approximation (3.7), now understood to hold in a vicinity of the instability threshold and with
$\,\tilde {\!{\boldsymbol {A}}}$
evolving with time according to the amplitude equation

in which the complex coefficients
$\alpha$
,
$\beta$
and
$\gamma$
are provided in Appendix B as quadratures involving the eigenfunction
$\varphi(s)$
, yielding the numerical values

In accordance with the above scaling remarks, this reduced-order model holds in the limit of small-magnitude perturbations
$|\,\tilde {\!{\boldsymbol {A}}}|\ll 1$
, with
$\mathcal {F}-\mathcal {F}_c=O(|\,\tilde {\!{\boldsymbol {A}}}|^2)$
. Furthermore, we note that the first two (nonlinear) terms on the right-hand side of (4.1) imply the long time scale
$1/|\,\tilde {\!{\boldsymbol {A}}}|^2$
, whereas the third (linear) term implies the long time scale
$1/|\mathcal {F}-\mathcal {F}_c|$
.
We prefer to work with rescaled, order-unity quantities. Thus, consistently with the weakly nonlinear analysis in § 5, we define

with
$0\lt \varepsilon \ll 1$
a small positive parameter and
$\chi$
real, and introduce the rescalings

wherein

is a slow time coordinate. The amplitude equation (4.1) then reads as

While only the product
$\varepsilon \chi$
is meaningful, the added freedom is convenient for jointly scaling the cases where the follower-force magnitude is at or near the instability threshold; in the latter case, the freedom can be removed by setting
$\chi =\mathrm {sgn}(\mathcal {F}-\mathcal {F}_c)$
.
It is useful to note the following three relations:



They readily follow from (4.6): we obtain (4.7a
) by subtracting the cross product of
${\boldsymbol {A}}^*$
and (4.6) from the conjugate of that product; we obtain (4.7b
) by considering the dot product of
$\boldsymbol {A}$
and (4.6), noting the identity
$2\boldsymbol{A} \boldsymbol{\cdot} \mathrm{d}\boldsymbol{A}/\mathrm{d} T = \mathrm{d}(\boldsymbol{A} \boldsymbol{\cdot} \boldsymbol{A})/\mathrm{d} T$
; and we obtain (4.7c
) by adding the dot product of
${\boldsymbol {A}}^*$
and (4.6) to the conjugate of that product. Equation (4.7a
) shows that if the fast-scale tip orbit is initially a line (corresponding to planar motion of the filament), then it remains so at all later slow times. Similarly, (4.7b
) shows that if the fast-scale tip orbit is initially circular then it remains so at all later slow times.
If the filament is confined to a plane (by the initial conditions), the complex-vector amplitude can be written
${\boldsymbol {A}}(T)={\hat {\boldsymbol {e}}} A(T)$
, where
$A(T)$
is a complex-scalar amplitude and
$\hat {\boldsymbol {e}}$
a unit vector parallel to the wall; the filament in that case performs planar-beating oscillations with
$A(T)$
determining their magnitude and phase. In that restricted scenario, the amplitude equation (4.6) reduces to

which has the form of the Stuart–Landau amplitude equation corresponding to a conventional (non-degenerate) Hopf instability (Drazin & Reid Reference Drazin and Reid2004).
We proceed in § 4.3 to find solutions of the amplitude equation (4.6) representing steady and periodic states; we also compare these solutions against numerical data provided by Eric E. Keaveny. In §§ 4.4–4.6, we analyse the stability of these states within the framework of the weakly nonlinear theory. In § 4.7, we present simulations of (4.6).
4.3. Periodic states
We consider ‘quasi-steady’ solutions of the amplitude equation (4.6) in the form

wherein
${\boldsymbol{\bar{A}}}$
is a time-independent complex vector parallel to the wall and
$\nu$
a real scalar. Such solutions represent periodic states of the filament at a corrected angular frequency
$\omega +\nu \chi ^{-1}(\mathcal {F}-\mathcal {F}_c)$
. From (4.6), the reduced amplitude
${\boldsymbol{\bar{A}}}$
satisfies

The trivial solution
${\boldsymbol{\bar{A}}}={\mathbf {0}}$
corresponds to the undeformed state of the filament. By crossing (4.10) with
${\boldsymbol{\bar{A}}}$
, we find that non-trivial quasi-steady solutions must satisfy either
${\boldsymbol{\bar{A}}}\times {\boldsymbol{\bar{A}}}^*={\mathbf {0}}$
, corresponding to planar beating, or
${\boldsymbol{\bar{A}}}\boldsymbol {\cdot }{\boldsymbol{\bar{A}}}=0$
, corresponding to circular whirling.
4.3.1. Planar-beating states
In the planar-beating case, we can write

wherein
$\hat {\boldsymbol {e}}$
is a unit vector parallel to the wall and
$\bar{A}$
is a non-vanishing complex scalar. Representing the latter in polar form,

wherein
$\vartheta$
represents a real phase, we find from (4.10) that the magnitude
$|{{\bar{A}}}|$
satisfies

whereas the direction
$\hat {\boldsymbol {e}}$
and phase
$\vartheta$
are arbitrary. Considering the real and imaginary parts of (4.13), we readily find

Given the coefficient values provided by (4.2), the right-hand side of (4.14a
) implies that the planar-beating states exist only in the case
$\chi \gt 0$
, namely when the follower-force magnitude is above the instability threshold. Setting
$\chi =1$
, without loss of generality, the corresponding tip trajectories follow from (3.8), using (4.4), (4.9), (4.11) and (4.12), as

where the rescaled tip displacement
$2|{{\bar{A}}}|\doteq 0.17412$
and angular-frequency correction
$\nu \doteq 5.5426$
are calculated using (4.2) and (4.14).
4.3.2. Circular-whirling states
In the circular-whirling case,
${\boldsymbol{\bar{A}}}\boldsymbol {\cdot }{\boldsymbol{\bar{A}}}=0$
, the real vectors
${\boldsymbol {A}}_r$
and
${\boldsymbol {A}}_i$
are orthogonal and of equal magnitude. We accordingly write

where the magnitude
$|{\boldsymbol{\bar{A}}}|\ne 0$
; the plus (minus) sign indicates clockwise (counter clockwise) motion (viewed from above the wall); and
${\hat {\boldsymbol {e}}}$
is a unit vector parallel to the wall which represents a phase. From (4.10), the magnitude
$|{\boldsymbol{\bar{A}}}|\ne 0$
satisfies

whereas the direction of motion and the phase unit vector
$\hat {\boldsymbol {e}}$
are arbitrary. Considering the real and imaginary parts of (4.17), we readily find

Similarly to the planar-beating states, we see from the right-hand side of (4.18a
) and the coefficient values provided by (4.2) that the circular-whirling states exist only in the case
$\chi \gt 0$
, similarly to the planar-beating states. Setting
$\chi =1$
, without loss of generality, the corresponding tip trajectories follow from (3.8), using (4.4), (4.9) and (4.16), as

where the rescaled tip radius
$\sqrt {2}|{\boldsymbol{\bar{A}}}|\doteq 0.16840$
and angular-frequency correction
$\nu \doteq 10.362$
are calculated using (4.2) and (4.18).
4.3.3. Bifurcation diagram and comparison with numerical solutions of the full problem
In figure 3, we plot maximal tip displacement versus dimensionless follower force for the steady and periodic states predicted by the weakly nonlinear theory. As evident in (4.15) and (4.19), the maximal tip displacement
$\propto \sqrt {\mathcal {F}-\mathcal {F}_c}$
for both the planar-beating and circular-whirling periodic states, the constant of proportionality being slightly larger for the planar-beating states. (We note that this square-root scaling near the threshold was hypothesised as well as numerically demonstrated by Clarke et al. (Reference Clarke, Hwang and Keaveny2024).) The stability characteristics indicated in the bifurcation diagram are derived and discussed in the following subsections.

Figure 3. (a) Bifurcation diagram showing maximal dip displacement (normalised by filament length) against the dimensionless follower-force magnitude
$\mathcal {F}$
, for the (i) undeformed steady state; (ii) planar-beating periodic states; and (iii) circular-whirling periodic states. The curves depict the predictions of the weakly nonlinear theory (see § 4): the undeformed state is stable preceding the bifurcation (solid black line) and unstable following it (dashed black line); the planar-beating states are longitudinally stable, up to phase shifts, but transversally unstable (dash-dotted red curve); and the circular-whirling states are stable up to phase shifts (solid blue curve). The symbols depict numerical data provided by Eric E. Keaveny (see § 4.3.3). (b) Zoomed-in comparison between the theory and numerical data depicted on a log–log scale.
With the purpose of validating the present weakly nonlinear theory, Eric E. Keaveny has conducted simulations based on the numerical methodology of Schoeller et al. (Reference Schoeller, Townsend, Westwood and Keaveny2021), adopting the same formulation as Clarke et al. (Reference Clarke, Hwang and Keaveny2024) only with the full Stokes hydrodynamics included in their study simplified to resistive force theory – consistently with the present formulation (cf. (2.2)). Figure 3 shows remarkable agreement, in terms of the maximal tip displacements, between these simulations and the theoretical predictions.
4.4. Linear stability of the undeformed state
The stability characteristics of the undeformed state, as well as the planar-beating and circular-whirling periodic states, are readily determined within the framework of the weakly nonlinear theory. We begin by revisiting the stability of the undeformed state, which has already been analysed exactly in § 3. To this end, we consider small perturbations of the complex-vector amplitude:

Linearisation of the amplitude equation (4.6) in the case
${\boldsymbol{\bar{A}}}={\mathbf {0}}$
gives

which possesses the general solution

In this linear approximation, the elliptical tip orbit at slow time
$T$
is obtained from the initial orbit defined by
${\boldsymbol {a}}(0)$
via a dilation by
$\exp (\chi \gamma _r T)$
. Since
$\gamma _r\gt 0$
, perturbations grow for
$\chi \gt 0$
and decay for
$\chi \lt 0$
. The case
$\chi =0$
, where nonlinear terms determine the stability of the undeformed state, will be considered in § 4.6.
It is instructive to repeat the above stability analysis, this time starting from the real dynamical system implied by the linearised amplitude equation (4.21):

As
${\boldsymbol {a}}_r$
and
${\boldsymbol {a}}_i$
are real vectors parallel to the wall, this system is four-dimensional. By substituting into (4.23) the modal ansatz

wherein
$\sigma$
represents a rescaled growth rate, we find the homogeneous system

The solution space is spanned by the eigenvectors

${\hat {\boldsymbol {e}}}_1$
and
${\hat {\boldsymbol {e}}}_2$
being a pair of orthogonal unit vectors parallel to the wall, with corresponding growth-rate eigenvalues

In accordance with the stability analysis in § 3, we find two identical pairs of complex-conjugate growth rates whose real part is positive for
$\chi \gt 0$
and negative for
$\chi \lt 0$
. The eigenvector pairs (4.27a
) and (4.27b
) correspond to planar beating along the
${\hat {\boldsymbol {e}}}_1$
and
${\hat {\boldsymbol {e}}}_2$
directions, respectively. The linearised whirling motion (4.22) is readily retrieved by general superposition of solutions in the form (4.24).
4.5. Linear stability of the periodic states
We next consider the stability of the planar-beating and circular-whirling states. In these cases, it is convenient to write the perturbation as

wherein the steady complex-vector amplitude
${\boldsymbol{\bar{A}}}$
and angular-frequency correction
$\nu$
are defined as in (4.9). Since these states exist only in the case
$\chi \gt 0$
, we set
$\chi =1$
, without loss of generality, for the remainder of this subsection. Linearisation of the amplitude equation (4.6) then gives, upon eliminating the quasi-steady relation (4.10),

4.5.1. Linear stability of the planar-beating states
For the planar-beating states (§ 4.3.1), the form of
${\boldsymbol{\bar{A}}}$
is given by (4.11) and (4.12), which feature the unit vector
$\hat {\boldsymbol {e}}$
(indicating the plane of beating) and phase
$\vartheta$
. Decomposing the perturbation into components parallel and perpendicular to the beating,
${\boldsymbol {a}}=a_{\parallel }{\hat {\boldsymbol {e}}}+a_{\perp }\boldsymbol {\hat {\imath }}\times {\hat {\boldsymbol {e}}}$
, the linearised amplitude equation (4.29) yields, upon eliminating
$\nu$
using the steady relation (4.13), the uncoupled pair of equations

wherein
$|{{\bar{A}}}|$
is provided by (4.14a
).
We first consider parallel perturbations. Separating (4.30a ) into its real and imaginary parts, we find the two-dimensional real system

wherein
$\zeta =\alpha +\beta$
. Seeking solutions in the form
$\{{a_\parallel }_r,{a_\parallel }_i\}=\{\tilde {a}{{}_{\parallel }}_r,\tilde {a}{{}_{\parallel }}_i\}\exp (\sigma _{\parallel }T) + \text {c.c.}$
, we find the growth-rate eigenvalues

where we have used (4.14a
) for
$|{{\bar{A}}}|$
. The zero eigenvalue
$\sigma _{\parallel }^{({1})}$
reflects the arbitrariness of the phase
$\vartheta$
. Since
$\gamma _r\gt 0$
, the eigenvalue
$\sigma _{\parallel }^{({2})}$
has a negative real part. Thus, planar beating is stable to parallel perturbations apart from the expected neutrality with respect to phase.
We next consider perpendicular perturbations. Separating (4.30b ) into its real and imaginary parts, we find the two-dimensional real system

Seeking solutions in the form
$\{{a_\perp }_r,{a_\perp }_i\}=\{\tilde {a}{{}_{\perp }}_r,\tilde {a}{{}_{\perp }}_i\}\exp (\sigma _{\perp }T) + \text {c.c.}$
, we readily find the growth-rate eigenvalues

where we have used (4.14a
) for
$|{{\bar{A}}}|$
. The zero eigenvalue
$\sigma _{\perp }^{({1})}$
reflects the arbitrariness of the direction
$\hat {\boldsymbol {e}}$
. Since
$\gamma _r\gt 0$
, whereas
$\alpha _r,\beta _r\lt 0$
, the eigenvalue
$\sigma _{\perp }^{({2})}$
has a positive real part. Thus, planar beating is unstable under perpendicular perturbations – and therefore unstable under general perturbations.
4.5.2. Stability of the circling-whirling states
For the circular-whirling states (§ 4.3.2), the form of
${\boldsymbol{\bar{A}}}$
is given by (4.16), with the plus/minus sign determining the direction of motion and the unit vector
$\hat {\boldsymbol {e}}$
representing a phase; without loss of generality, we only consider the stability of circular whirling in the clockwise direction, corresponding to choosing the plus sign. Writing
${\boldsymbol {a}}=a_1{\hat {\boldsymbol {e}}}+a_2\boldsymbol {\hat {\imath }}\times {\hat {\boldsymbol {e}}}$
, the linearised amplitude equation (4.29) gives, upon eliminating
$\nu$
using the steady relation (4.17), the coupled pair of equations


wherein
$|{\boldsymbol {A}}|$
is provided by (4.18a
). Decomposing (4.35) into their real and imaginary parts, we find the four-dimensional real system

Seeking solutions in the form
$\{{a_1}_r,{a_1}_i,{a_2}_r,{a_2}_i\}=\{\tilde {a}{{}_{1}}_r,\tilde {a}{{}_{1}}_i,\tilde {a}{{}_{2}}_r,\tilde {a}{{}_{2}}_i\}\exp (\sigma T) + \text {c.c.}$
, we find the dispersion relation

whose solutions yield the growth-rate eigenvalues

where we have used (4.18a
) for
$|{\boldsymbol{\bar{A}}}|$
. The zero eigenvalue
$\sigma ^{({1})}$
reflects the arbitrariness of the phase unit vector
$\hat {\boldsymbol {e}}$
. Since
$\gamma _r\gt 0$
, whereas
$\alpha _r,\beta _r\lt 0$
, the remaining eigenvalues have negative real parts. Thus, the circular-whirling states are stable apart from the expected neutrality with respect to phase.
The predicted stability characteristics of the different states are indicated in the bifurcation diagram (figure 3 a). The stability of the circular-whirling states and the instability of the planar-beating states (under general three-dimensional perturbations) are in accordance with the Floquet analyses numerically performed by Clarke et al. (Reference Clarke, Hwang and Keaveny2024). As first established by Clarke et al. (Reference Clarke, Hwang and Keaveny2024), these stability characteristics serve to rationalise the selection of whirling over beating in initial-value simulations not restricted to planar deformations (Ling et al. Reference Ling, Guo and Kanso2018).
4.6. Nonlinear stability of the undeformed state
It remains to determine the nonlinear stability of the undeformed state at the linear-instability threshold. Setting
$\chi =0$
, the amplitude equation (4.6) reduces to

(In this case, the small parameter
$\varepsilon$
simply measures the magnitude of the perturbation, which scales with the
$1/2$
power of that parameter.) We wish to determine whether
$|{\boldsymbol {A}}|$
decays or grows as
$T\to \infty$
. To this end, we set
$\chi =0$
in (4.7a,c
) to find the pair of equations

where we denote
${\boldsymbol {B}}={\boldsymbol {A}}^*\times {\boldsymbol {A}}$
.
In the special case of planar perturbations, where
$\boldsymbol {B}$
vanishes, (4.40a
) reduces to

Since
$\alpha _r,\beta _r\lt 0$
, planar perturbations decay. Explicitly, integration of (4.41) gives

revealing algebraic decay with the
$-1/2$
power of time.
For initially non-planar perturbations we consider the pair of equations (4.40) as a nonlinear dynamical system for the magnitudes
$|{\boldsymbol {A}}|$
and
$|{\boldsymbol {B}}|$
. This system has no fixed points, and consideration of the phase space shows that all orbits in the first quadrant (
$|{\boldsymbol {A}}|,|{\boldsymbol {B}}|\gt 0$
) approach the origin. Positing a power-law behaviour for
$|{\boldsymbol {A}}|$
and
$ |{\boldsymbol {B}}|$
, we readily identify the late-time behaviours

Given these asymptotes, we further find using (3.9) that the semi-major and semi-minor axes of the fast-scale tip orbit behave as (up to an
$\varepsilon ^{1/2}$
scaling)

Thus, at late times the filament exhibits circular whirling on the fast scale, with the tip radius decaying algebraically with the
$-1/2$
power of time.
4.7. Numerical illustrations
The amplitude equation (4.6) is straightforward to integrate numerically starting from some initial amplitude
${\boldsymbol {A}}(0)$
. Given a solution
${\boldsymbol {A}}(T)$
, a function of the slow time
$T$
, the corresponding motion of the filament in three dimensions as well as the internal force and moment distributions can be calculated at any time
$t$
by making the substitution
$\,\tilde {\!{\boldsymbol {A}}}\Rightarrow \varepsilon ^{1/2}{\boldsymbol {A}}(T)$
in (3.7) (cf. (4.4)). For the purpose of illustrating the theory, we present in figures 4–6 sample solutions for several values of
$\mathcal {F}-\mathcal {F}_c=\varepsilon \chi$
and initial conditions.

Figure 4. Planar multiple-scale dynamics for the initial condition
${\boldsymbol {A}}(0)=0.02{\hat {\boldsymbol {e}}}_x$
and indicated values of
$\mathcal {F}-\mathcal {F}_c$
, where
${\hat {\boldsymbol {e}}}_x$
is a unit vector parallel to the wall. The plot depicts the tip displacement in the
${\hat {\boldsymbol {e}}}_x$
direction (solid curves), scaled by
$|\mathcal {F}-\mathcal {F}_c|^{1/2}$
, as a function of time, scaled by the natural filament period
$2\pi /\omega$
. The dashed curves depict the slow-time envelopes
$2|{\boldsymbol {A}}(T)|$
. The dash-dotted line marks the peak displacement corresponding to the planar-beating state.

Figure 5. Non-planar multiple-scale dynamics for the initial condition
${\boldsymbol {A}}(0)=0.03{\hat {\boldsymbol {e}}}_x+i0.005{\hat {\boldsymbol {e}}}_y$
and
$\mathcal {F}-\mathcal {F}_c=0.5$
, where
$\{{\hat {\boldsymbol {e}}}_x,{\hat {\boldsymbol {e}}}_y,\boldsymbol {\hat {\imath }}\}$
is a right-handed system of unit vectors. The panels depict a top view (looking towards the wall, with
${\hat {\boldsymbol {e}}}_x$
pointing to the right) showing the tip position (filled circle) and filament projection (thick solid line) at the indicated times, along with the instantaneous fast-scale tip orbit (dashed ellipse), tip trajectory starting from the previous time stamp (fading thin curves) and the radius corresponding to the circular-whirling states (dash-dotted circle); distances are scaled by
$\sqrt {\mathcal {F}-\mathcal {F}_c}$
.

Figure 6. Algebraic decay of perturbations at the linear-instability threshold
$\mathcal {F}=\mathcal {F}_c$
(see § 4.6). The plot compares
$|{\boldsymbol {A}}(T)|$
in the planar scenario
${\boldsymbol {A}}(0)={\hat {\boldsymbol {e}}}_x$
(thick solid curve) and non-planar scenario
${\boldsymbol {A}}(0)=({\hat {\boldsymbol {e}}}_x+i\sqrt {3}{\hat {\boldsymbol {e}}}_y)/2$
(thick dash-dotted curve) with the power-law asymptotes (4.42) and (4.43a
), where
$\{{\hat {\boldsymbol {e}}}_x,{\hat {\boldsymbol {e}}}_y,\boldsymbol {\hat {\imath }}\}$
is a right-handed system of unit vectors. For the non-planar scenario, we also depict the magnitude
$|{\boldsymbol {B}}|=|{\boldsymbol {A}}^*\times {\boldsymbol {A}}|$
(thick dashed curve) and its power-law asymptote (4.43b
), and in the inset a top view (looking towards the wall, with
${\hat {\boldsymbol {e}}}_x$
pointing to the right) of the tip trajectory scaled by
$\varepsilon ^{1/2}$
, with
$\varepsilon =0.5$
.
In figure 4, we consider the case where the initial fast-time oscillation is planar. In light of (4.7a
), the motion of the filament remains planar at all times. We observe the multiple-scales evolution of the dynamics towards the periodic planar-beating states (see §§ 4.3.1 and 4.5.1), for
$\mathcal {F}\gt \mathcal {F}_c$
, or the undeformed configuration, for
$\mathcal {F}\lt \mathcal {F}_c$
, with faster transients for larger
$|\mathcal {F}-\mathcal {F}_c|$
. In figure 5, we consider an example where the initial fast-time oscillation is nearly planar and
$\mathcal {F}\gt \mathcal {F}_c$
. We observe the slow-time expansion and rotation of the fast-scale elliptical orbit of the tip as it approaches the circular orbit corresponding to the periodic circular-whirling states (see §§ 4.3.2 and 4.5.2). Lastly, in figure 6, we demonstrate the algebraic decay of both planar and non-planar perturbations at the instability threshold,
$\mathcal {F}=\mathcal {F}_c$
(see § 4.6).
5. Weakly nonlinear analysis
5.1. Multiple-scale expansions near the instability threshold
In this section, we derive the amplitude equation (4.6), which together with the general representation (3.7) for the linear approximation at the instability threshold constitutes the weakly nonlinear theory that we have presented, analysed and illustrated in § 4. The derivation consists of a weakly nonlinear analysis in the multiple-scales near-threshold regime identified in § 4.1.
We write
$\mathcal {F}=\mathcal {F}_c+\varepsilon \chi$
for the dimensionless follower force, as in (4.3), and consider the long-time dynamics of small-magnitude oscillations for
$0\lt \varepsilon \ll 1$
, holding
$\chi$
fixed. Following the method of multiple scales (Hinch Reference Hinch1991), we introduce the two-scale extension
$\underline {{\boldsymbol {r}}}(s,\tau ,T)$
of the centreline position vector
${{\boldsymbol {r}}}(s,t)$
, wherein the ‘fast time’
$\tau$
and ‘slow time’
$T$
are treated as independent coordinates, such that
$\underline {{\boldsymbol {r}}}(s,\tau ,T)={\boldsymbol {r}}(s,t)$
on the ‘physical diagonal’:

we also introduce analogous extensions
$\underline {{\hat {\boldsymbol {t}}}}(s,\tau ,T)$
,
$\underline {{\boldsymbol {F}}}(s,\tau ,T)$
and
$\underline {{\boldsymbol {M}}}(s,\tau ,T)$
for the tangent unit vector
${\hat {\boldsymbol {t}}}(s,t)$
, internal force
${\boldsymbol {F}}(s,t)$
and internal moment
${\boldsymbol {M}}(s,t)$
, respectively. The extended fields satisfy the same problem as formulated in § 2, only with the time derivative in (2.6b
) transformed according to

We posit the ‘weakly nonlinear expansion’

with analogous expansions defined for
$\underline {{\hat {\boldsymbol {t}}}}(s,\tau ,T)$
,
$\underline {{\boldsymbol {F}}}(s,\tau ,T)$
and
$\underline {{\boldsymbol {M}}}(s,\tau ,T)$
, the zeroth-order fields being the base state (2.11) evaluated at the threshold:

A key idea underlying the method is to exploit the added freedom associated with the extension to two time scales in order to ensure that the weakly nonlinear expansions remain asymptotically ordered. We note that the scaling of time by
$1/\varepsilon$
and the expansion in half-powers of
$\varepsilon$
are suggested by the scaling arguments given in § 4.1.
5.2.
The
$O(\varepsilon ^{1/2})$
problem
At
$O(\varepsilon ^{1/2})$
, the governing partial differential equations (2.6), together with the time-derivative transformation (5.2), give




the inextensibility constraint (2.7) gives

while the boundary conditions (2.8) and (2.9) give

and

The above problem is nothing but the full problem at the instability threshold linearised for small deformations (scaled by
$\varepsilon ^{1/2}$
). We may accordingly deduce the appropriate general solution from the linear theory of § 3:




wherein
${\boldsymbol {A}}(T)$
is a complex-vector amplitude parallel to the wall evolving on the slow time scale. As detailed in § 3, this solution form represents an arbitrary superposition of linearly neutral planar-beating modes, which generally gives rise to a three-dimensional whirling motion. We have not included in the general solution (5.9) potential contributions of any stable linear modes as these would exponentially decay on the fast time scale with no effect on the slow-time dynamics. We must proceed to higher order in
$\varepsilon$
to derive an evolution equation for
${\boldsymbol {A}}(T)$
.
5.3.
The
$O(\varepsilon )$
problem
At
$O(\varepsilon )$
, the governing partial differential equations (2.6), together with the time-derivative transformation (5.2), give




the inextensibility constraint (2.7) gives

while the boundary conditions (2.8) and (2.9) give

and

The homogeneous part of the above linear problem is the same as for the
$O(\varepsilon ^{1/2})$
problem. The inhomogeneous terms, calculated by substituting (5.9), are provided explicitly inside the curly brackets; they are normal to the wall and composed of contributions that are harmonic in the fast time
$\tau$
with angular frequencies
$0$
or
$2\omega$
. Accordingly, these forcing terms cannot resonantly excite homogeneous beating/whirling solutions of the form (5.9), which are parallel to the wall and harmonic in the fast time
$\tau$
with angular frequency
$\omega$
. Indeed, it is straightforward to derive the solutions




whose periodicity in
$\tau$
confirms the regularity of the weakly nonlinear expansions to
$O(\varepsilon )$
(cf. (5.3)). In (5.14), we have discarded all homogeneous solutions, including beating/whirling solutions of the form (5.9) as well as solutions decaying on the fast time scale – these do not effect the slow-time dynamics of the leading-order amplitude
${\boldsymbol {A}}(T)$
.
5.4.
The
$O(\varepsilon ^{3/2})$
problem
At
$O(\varepsilon ^{3/2})$
, the governing partial differential equations (2.6), together with the time-derivative transformation (5.2), give




the inextensibility constraint (2.7) gives

where the right-hand side vanishes since
$\boldsymbol {t}_1$
and
$\boldsymbol {t}_{1/2}$
are perpendicular; and the boundary conditions (2.8) and (2.9) give

and

The above linear problem has the same form as the problems encountered in the preceding orders, except that at the present order the inhomogeneous forcing terms include potentially resonant contributions – parallel to the wall and harmonic in the fast time
$\tau$
with the natural angular frequency
$\omega$
. To focus attention on these contributions, we express the forcing terms as shown inside the curly brackets. Therein, ‘
$\text {n.s.t.}$
’ stands for non-secular terms and represents all other contributions, while
${\boldsymbol {f}}_{(2)}$
,
${\boldsymbol {f}}_{(3)}$
and
${\boldsymbol {f}}_{(4)}$
are the ‘reduced forcing terms’




obtained by substituting (5.9) and (5.14); their numbering will become meaningful in the following subsection.
5.5. Solvability condition
5.5.1. Restricted linear problem
The
$\exp (i\omega \tau )$
Fourier component of the
$O(\varepsilon ^{3/2})$
problem, projected parallel to the wall, furnishes a ‘restricted’ linear problem in the form

where
$ \boldsymbol{\psi }(s) = [ \acute {{\boldsymbol {r}}}(s), \,\, \acute {\boldsymbol {t}}(s), \,\, \acute {{\boldsymbol {F}}}(s), \,\, \acute {{\boldsymbol {M}}}(s) ]^T$
denotes the unknown column-array field, whose elements are complex vector fields parallel to the wall;
$\mathsf{L}$
,
$\mathsf{S}$
and
$\mathsf{T}$
are the matrix-differential operators

and
$\boldsymbol{f}(s)=[{\boldsymbol {f}}_{(1)}(s),\,{\boldsymbol {f}}_{(2)}(s)\,{\boldsymbol {f}}_{(3)}(s),\,{\boldsymbol {f}}_{(4)}(s)]^T$
,
$\boldsymbol{g}=[{\boldsymbol {g}}_{(1)},\,{\boldsymbol {g}}_{(2)}]^T$
and
$\boldsymbol{h}=[{\boldsymbol {h}}_{(1)},\, {\boldsymbol {h}}_{(2)}]^T$
are column-array forcing terms, whose elements are vectors parallel to the wall (fields for
$\boldsymbol{f}(s)$
, constants for
$\boldsymbol{g}$
and
$\boldsymbol{h}$
). In the present scenario, we find from (5.15)–(5.18) that
${\boldsymbol {f}}_{(1)}(s)$
,
${\boldsymbol {h}}_{(2)}$
,
${\boldsymbol {g}}_{(1)}$
and
${\boldsymbol {g}}_{(2)}$
vanish, whereas
${\boldsymbol {f}}_{(2)}(s)$
,
${\boldsymbol {f}}_{(3)}(s)$
,
${\boldsymbol {f}}_{(4)}(s)$
and
${\boldsymbol {h}}_{(1)}$
are provided by (5.19). We know from the linear theory of § 3 that the restricted problem possesses non-trivial homogeneous solutions in the form
$ \boldsymbol{\psi } = [{\boldsymbol {a}}\varphi , \,\, {\boldsymbol {a}}\varphi ^{\prime}, \,\, -{\boldsymbol {a}}(\varphi ^{\prime\prime\prime}+\mathcal {F}_c\varphi ^{\prime}),\,\, \boldsymbol {\hat {\imath }}\times {\boldsymbol {a}}\varphi ^{\prime\prime} ]^T$
, with
$\boldsymbol {a}$
an arbitrary complex vector parallel to the wall. We therefore expect that solutions to the inhomogeneous restricted problem exist – whereby resonance is avoided in the
$O(\varepsilon ^{3/2})$
problem of § 5.4 – only under certain ‘solvability’ conditions on the forcing terms.
Below, we derive a necessary condition for existence of solutions to the restricted problem. (By alluding to the Fredholm alternative theorem for differential operators (Keener Reference Keener2018), it could be shown that the condition is also sufficient.) In deriving this solvability condition, we shall allow for the full form of the forcing terms in (5.20) in order to facilitate future generalisations of the theory as discussed in § 6.
5.5.2. Adjoint operators and homogeneous solutions
For any pair of column-vector fields
$ \boldsymbol{\psi }^{({1})}(s) = [ \acute {{\boldsymbol {r}}}^{({1})}(s), \,\, \acute {\boldsymbol {t}}^{({1})}(s), \,\, \acute {{\boldsymbol {F}}}^{({1})}(s), \,\, \acute {{\boldsymbol {M}}}^{({1})}(s) ]^T$
and
$\boldsymbol {\psi }^{({2})}(s) = [ \acute {{\boldsymbol {r}}}^{({2})}(s), \,\, \acute {\boldsymbol {t}}^{({2})}(s), \,\, \acute {{\boldsymbol {F}}}^{({2})}(s), \,\, \acute {{\boldsymbol {M}}}^{({2})}(s) ]^T$
, we define the inner product

Following Keener (Reference Keener2018), we seek adjoint operators
$\mathsf {L}^{\dagger }$
,
$\mathsf {S}^{\dagger }$
and
$\mathsf {T}^{\dagger }$
such that the factor

vanishes for any
$ \boldsymbol{\psi }^{({1})}$
satisfying the ‘direct boundary conditions’
$\mathsf {S} \boldsymbol{\psi }^{({1})}(0)= \textbf{0}$
and
$\mathsf {T} \boldsymbol{\psi }^{({1})}(1)= \textbf{0}$
, and
$ \boldsymbol{\psi }^{({2})}$
satisfying the ‘adjoint boundary conditions’
$\mathsf {S}^\dagger \boldsymbol{\psi }^{({2})}(0)= \textbf{0}$
and
$\mathsf {T}^\dagger \boldsymbol{\psi }^{({2})}(1)= \textbf{0}$
. To this end, we integrate by parts the product

and use the vector triple product to find

Thus, by inspection, the adjoint operators are

Associated with the adjoint operators is the adjoint homogeneous problem,

It is readily found to possess the non-trivial solutions

where
$\boldsymbol {b}$
is an arbitrary complex vector parallel to the wall and the adjoint eigenfunction
$\phi (s)$
satisfies the ordinary differential equation

and the boundary conditions

The unidimensional problem (5.29)–(5.30) defines the adjoint eigenfunction
$\phi (s)$
up to an arbitrary complex prefactor, which we conveniently choose such that
$\phi (0)=1$
; the eigenfunction can be readily expressed analytically in terms of the numerical values
$\mathcal {F}_c$
and
$\omega$
, or calculated numerically.
5.5.3. Solvability condition
In the definition (5.23) of the factor
$J( \boldsymbol{\psi }^{({1})}, \boldsymbol{\psi }^{({2})})$
, let
$ \boldsymbol{\psi }^{({1})}= \boldsymbol{\psi }$
be a solution to the inhomogeneous problem (5.20) and
$ \boldsymbol{\psi }^{({2})}= \boldsymbol{\zeta }$
any solution to the homogeneous adjoint problem (cf. (5.28)). Using (5.20a
) and (5.27a
), we find from (5.23) the relation

Since
$ \boldsymbol{\psi }$
satisfies inhomogeneous conditions, the factor
$J( \boldsymbol{\psi }, \boldsymbol{\zeta })$
does not necessarily vanish; rather, the boundary terms in (5.25) give, upon substituting (5.20b,c
) and (5.28)–(5.30),

From (5.20a ) we find, using the triple-vector-product identity,

Substituting (5.32) and (5.33) into (5.31), we find, given that
$\boldsymbol {b}$
is arbitrary,

This result constitutes a general solvability condition that could be applied to any inhomogeneous linear problem in the form (5.20). In the present scenario, where
${\boldsymbol {f}}_{(1)}(s)$
,
${\boldsymbol {h}}_{(2)}$
,
${\boldsymbol {g}}_{(1)}$
and
${\boldsymbol {g}}_{(2)}$
vanish, (5.34) reduces to

where we have also taken the complex conjugate. Substituting (5.19) for
${\boldsymbol {f}}_{(2)}(s)$
,
${\boldsymbol {f}}_{(3)}(s)$
,
${\boldsymbol {f}}_{(4)}(s)$
and
${\boldsymbol {h}}_{(1)}$
, we arrive at the amplitude equation (4.6) with the coefficients
$\alpha$
,
$\beta$
and
$\gamma$
defined by the quadratures presented in Appendix B.
6. Concluding remarks
We have developed a weakly nonlinear theory illuminating the onset of spontaneous beating and whirling in the follower force model. The biological inspiration for the follower force model suggests several extensions to the modelling, including cross-sectionally anisotropic bending; preferred curvature; axially distributed compressive forces; applied moments; non-local hydrodynamic corrections and entrainment of the flow by molecular motors; and interactions with ambient flows and neighbouring filaments. Sufficiently near the threshold, even weak perturbations such as the above may significantly influence the leading-order dynamics through resonant interactions with the linear beating modes. It would be relatively straightforward to generalise our weakly nonlinear analysis to such perturbative scenarios: as long as the homogeneous linearised operator remains unaffected, a generalised weakly nonlinear analysis could directly utilise the linear approximation at the threshold and the solvability condition developed herein. For non-weak perturbations that appreciably alter the near-threshold linearised problem, our approach remains applicable, though most steps of the analysis would need to be revisited.
Besides studying the onset of spontaneous dynamics under variations to the modelling, the present work could be extended in several other interesting directions. One is to carry out weakly nonlinear analyses near other critical points of the dynamics, e.g. near the bifurcations discovered by Clarke et al. (Reference Clarke, Hwang and Keaveny2024) where the quasi-periodic branch of solutions termed ‘QP1’ merges with the planar-beating and whirling states, modifying their stability. Another direction is to consider weak perturbation and interaction effects away from critical points; such perturbations can again have a leading-order influence over long times, e.g. manifested in slow phase variations or slow reorientation of the beating plane. Lastly, it may be of interest to consider the onset of spontaneous dynamics in the inertial version of the follower force model, where a finite-mass filament deforms in the absence of viscous effects. The inertial problem was widely considered in the elasticity literature (Beck Reference Beck1952; Langthjem & Sugiyama Reference Langthjem and Sugiyama2000), although mostly assuming planar deformations. A three-dimensional nonlinear analysis in the spirit of this work may accordingly prove illuminating.
Supplementary material
A supplementary material is available at https://doi.org/10.1017/jfm.2025.128.
Acknowledgements
The author is grateful to Dr Eric E. Keaveny of Imperial College London for introducing him to this problem and providing the numerical results presented in § 4.3.3.
Funding
The author acknowledges the support of the Leverhulme Trust through Research Project Grant RPG-2021-161.
Declaration of interests
The author reports no conflict of interest.
Appendix A. Absence of twist
When including resistance to twist, the constitutive equation (2.5) generalises to

where
$T_*$
is the twist stiffness and the unit-vector triplet
$\{\hat {\boldsymbol {\mu }},\hat {\boldsymbol {\nu }},{\hat {\boldsymbol {t}}}\}$
is a right-handed orthogonal material frame (Schoeller et al. Reference Schoeller, Townsend, Westwood and Keaveny2021). Following Landau et al. (Reference Landau, Pitaevskii, Kosevich and Lifshitz2012), we confirm below that the twist term in (A1) vanishes identically for the scenario considered in this paper. We also show how that condition could be used to calculate the material frame along the filament given its centreline, the latter independently determined by the bending problem formulated in § 2.
The dot product of the equilibrium relation (2.1b
) and
$\hat {\boldsymbol {t}}$
gives

Similarly, the dot product of the generalised constitutive relation (A1) and
$\partial {\hat {\boldsymbol {t}}}/\partial s_*$
gives

since
$\partial {\hat {\boldsymbol {t}}}/\partial s_*$
is perpendicular to
$\hat {\boldsymbol {t}}$
. Combining (A2) and (A3), we find

whereby integration together with the tip boundary condition (2.4b ) yields

It follows that the twist term in (A1) vanishes.
The last result implies the orthogonality relation

The remaining components of
$\partial \hat {\boldsymbol {\mu }}/\partial s_*$
can be found as follows. First, the constraint that
$\hat {\boldsymbol {\mu }}$
is a unit vector implies

Second, since
$\hat {\boldsymbol {t}}$
and
$\hat {\boldsymbol {\mu }}$
are orthogonal,

Combining (A6)–(A8), we find the equation

which is supplemented by the boundary condition

wherein
$\hat {\boldsymbol {e}}$
is an arbitrarily chosen unit vector parallel to the wall; since the filament is clamped at the wall,
$\hat {\boldsymbol {e}}$
is time-independent and can accordingly be identified with
$\hat {\boldsymbol {\mu }}$
in its undeformed configuration. Given
${\hat {\boldsymbol {t}}}(s_*,t_*)$
from the solution to the bending problem formulated in § 2, (A9) and (A10) could, in principle, be solved for
$\hat {\boldsymbol {\mu }}(s_*,t_*)$
. The remaining unit vector
$\hat {\boldsymbol {\nu }}(s_*,t_*)$
could then be obtained from
$\hat {\boldsymbol {\nu }}={\hat {\boldsymbol {t}}}\times \hat {\boldsymbol {\mu }}$
.
Appendix B. Coefficients appearing in the amplitude equation
The coefficients appearing in the amplitude equation (4.1), or (4.6), are defined as



wherein
$\xi = -\int _0^1{\phi^{*\prime}}\varphi \,\mathrm {d}s$
. Recall that
$\varphi (s)$
is the eigenfunction defined in § 3, along with the threshold follower-force and angular-frequency values
$\mathcal {F}_c$
and
$\varpi$
; and
$\phi (s)$
is the adjoint eigenfunction defined in § 5.5. Since all of these quantities are parameter-free, so are the coefficients
$\alpha$
,
$\beta$
and
$\gamma$
. See the supplementary material for a Mathematica notebook (Wolfram Research 2020) that calculates the eigenfunctions
$\varphi (s)$
and
$\phi (s)$
, and the critical values
$\mathcal {F}_c$
and
$\omega$
, and uses these to evaluate the above integrals giving the numerical values of the coefficients quoted in (4.2).