1. Introduction
Living cells generate fluid motion in and around them to facilitate fundamental biological processes such as motility, sensing and feeding. These fluid flows are essential to overcome the limitations of molecular diffusion for mass transport, and further, they must be effective in a reversible, viscosity-dominated environment where the effects of inertia are negligible (Lighthill Reference Lighthill1976; Lauga Reference Lauga2020). One way of creating these flows is through cargo-carrying motor proteins that move along biological filaments – so-called active filaments (Mallik et al. Reference Mallik, Carter, Lex, King and Gross2004; Cho & Vale Reference Cho and Vale2012). Within the cell, motor proteins move along microtubules, generating effective surface velocities by pushing the surrounding fluid as they move. Interactions between neighbouring active filaments enable coordinated fluid flows within the cell, which are critical, for example, to the proper development of fruit fly eggs (Stein et al. Reference Stein, De, G., E., M. and Goldstein2021; Htet & Lauga Reference Htet and Lauga2023; Dutta et al. Reference Dutta, Farhadifar, Lu, Kabacaoğlu, Blackwell, Stein, Lakonishok, Gelfand, Shvartsman and Shelley2024), and to the control and positioning of metaphase spindles during cell division (Wu et al. Reference Wu, Kabacaoğlu, Nazockdast, Chang, Shelley and Needleman2024; Young et al. Reference Young, Herrera, Zhang, Farhadifar and Shelley2025). Active filaments also comprise the axoneme, the internal structure of cilia, where motor proteins generate the internal shear that produces cilium beating. Cilia are responsible for fluid motion outside cells in processes that include the locomotion of microscopic organisms and cells (Berner Reference Berner1993; Gaffney et al. Reference Gaffney, Gadêlha, Smith, Blake and Kirkman-Brown2011), and the pumping of fluids in the vital organs of more complex organisms (Sleigh, Blake & Liron Reference Sleigh, Blake and Liron1988; Lyons, Saridogan & Djahanbakhch Reference Lyons, Saridogan and Djahanbakhch2006; Faubel et al. Reference Faubel, Westendorf, Bodenschatz and Eichele2016).
Experiments suggest that the piconewton-scale active forces exerted by the motor proteins are sufficient to deform the biofilament along which they move (Svoboda & Block Reference Svoboda and Block1994). As a result, active filament motion arises from the fluid–structure interactions driven by the stresses associated with motor protein activity. Since the Reynolds number associated with filament motion is low, models for active filaments can take advantage of theoretical and computational methods for slender, flexible bodies in Stokes flow. These often couple either Euler–Bernoulli or Kirchhoff rod theories with hydrodynamic models based on local resistive force theory (Gray & Hancock Reference Gray and Hancock1955; Lighthill Reference Lighthill1976) and non-local slender-body theory (Johnson Reference Johnson1980; Tornberg & Shelley Reference Tornberg and Shelley2004; Nazockdast et al. Reference Nazockdast, Rahimian, Zorin and Shelley2017), or instead methods based around regularised Stokeslets, the Rotne–Prager–Yamakawa (RPY) tensor, or immersed-boundary and force-coupling methods (Lim Reference Lim2010; Olson, Lim & Cortez Reference Olson, Lim and Cortez2013; Delmotte et al. Reference Delmotte, Climent and Plouraboué2015a ; Hall-McNair et al. Reference Hall-McNair, Montenegro-Johnson, Gadêlha, Smith and Gallagher2019; Jabbarzadeh & Fu Reference Jabbarzadeh and Fu2020; Walker, Ishimoto & Gaffney Reference Walker, Ishimoto and Gaffney2020; Maxian et al. Reference Maxian, Mogilner and Donev2021; Schoeller et al. Reference Schoeller, Townsend, Westwood and Keaveny2021; Maxian et al. Reference Maxian, Sprinkle, Peskin and Donev2022; Fuchter & Bloomfield-Gadêlha Reference Fuchter and Bloomfield-Gadêlha2023) that resolve hydrodynamic interactions between different elements of the filament. While these approaches provide descriptions of the passive elastic and viscous forces that are present, the active forces that represent motor activity and drive filament motion must also be incorporated.
A first approach to incorporating motor activity is to include a compressive, tangential load (Beck Reference Beck1952; Herrmann & Bungay Reference Herrmann and Bungay1964; Langthjem & Sugiyama Reference Langthjem and Sugiyama2000), the so-called follower force, to represent the forces that the motors exert on the filament. Using resistive force theory, De Canio et al. (Reference De Canio, Lauga and Goldstein2017) considered a single follower force at the free end of a filament that is clamped at its other end, and analysed its planar dynamics as the magnitude of the follower force increases. They showed that at a critical load, the filament buckles and undergoes a supercritical Hopf bifurcation leading to a time-periodic beating state. Allowing for fully three-dimensional (3-D) dynamics, Ling, Guo & Kanso (Reference Ling, Guo and Kanso2018) instead observed a whirling state just after buckling and planar beating arises at higher values of the follower force. They showed also how these states can change depending on where the follower forces are placed along the filament. To better understand the differences between the planar and fully-3-D results, Clarke, Hwang & Keaveny (Reference Clarke, Hwang and Keaveny2024) performed a computational dynamical systems analysis of the tip-driven model, and found that at buckling, both the planar beating and whirling solutions emerge through a double Hopf bifurcation, with whirling being stable. This scenario was later confirmed by Schnitzer (Reference Schnitzer2025), who performed a weakly nonlinear analysis around the bifurcation. Additionally, Clarke et al. (Reference Clarke, Hwang and Keaveny2024) found that a stable quasi-periodic state provides the transition from whirling to beating that occurs at higher follower forces. At even higher follower force values, beating is found to be unstable, and various quasi-periodic and chaotic states are observed instead (Ling et al. Reference Ling, Guo and Kanso2018; Krishnamurthy & Prakash Reference Krishnamurthy and Prakash2023; Clarke et al. Reference Clarke, Hwang and Keaveny2024).
In addition to follower forces acting on the filament, the motor proteins will exert the opposite force on the fluid, resulting in entrainment, or a surface flow, in the vicinity of the filament. The additional effect of the surface flow has been highlighted in modelling the collective dynamics of microtubules involved in cytoplasmic streaming (Stein et al. Reference Stein, De, G., E., M. and Goldstein2021) and is essential to the onset of large-scale flows. Despite this, less is known about single-filament dynamics when a surface flow is included. When it has been incorporated into the follower force model, surface flows have been shown to counter the compressive effects of the follower force and shift buckling to higher follower force values (De Canio et al. Reference De Canio, Lauga and Goldstein2017). In exploring a closely-related active colloidal filament model, where a spherical squirmer is at the tip of an otherwise passive filament, Laskar & Adhikari (Reference Laskar and Adhikari2017) found several dynamic states, such as beating and whirling, qualitatively similar to those that arise in the tip-driven, follower force model.
In this paper, we present a comparative dynamical systems study to identify filament states and their bifurcations for models that include or ignore surface flows. We also consider cases where activity is localised to the filament’s free end, or distributed along the filament length. Our computations reveal that filament states and bifurcations in all models are qualitatively similar at low actuation, though the inclusion of surface flows shifts bifurcations to higher activity values, consistent with previous studies. At higher actuation, however, we find that it is the distribution of activity that plays a key role in the transitions and resulting states. Actuation distributed uniformly along the filament length leads to a single, time-dependent state that remains stable even at very high values of actuation. When activity is concentrated at the free end, however, the filament undergoes a complex series of bifurcations culminating in chaotic behaviour. We attribute these differences in dynamic states at high actuation to the internal stress profiles that arise in the different models.
The paper is structured as follows. In § 2, we present the different active filament models. In this section, we also describe the numerical methods that we use to solve the resulting equations, as well as the computational dynamical systems tools that we employ to obtain different states, assess their stability, and analyse quasi-periodic solutions. In § 3, we present the main results of this study; we discuss and compare the bifurcation diagrams and dynamic states emerging from the active filament models, dividing results between the low and high actuation regimes. Finally, in § 4 we present the conclusions of our study.
2. Methods
2.1. Active filament models
2.1.1. Force and moment balances
We begin by presenting the modelling framework that we use to describe filament dynamics when driven by either follower forces or surface flows. In the follower force model, the filament is subject to a distribution of compressive tangential forces, and entrainment is ignored. Our implementation follows directly from Schoeller et al. (Reference Schoeller, Townsend, Westwood and Keaveny2021) and Clarke et al. (Reference Clarke, Hwang and Keaveny2024), where only a tip-driven follower force was considered. In the surface flow model, the effects of both compressive forces and entrainment are incorporated by introducing an effective surface flow based on the squirmer model (Lighthill Reference Lighthill1952; Blake Reference Blake1971), similar to the model considered by Laskar & Adhikari (Reference Laskar and Adhikari2017).
For both the follower force and surface flow models, we consider an inextensible filament whose base is held fixed at the origin and aligned with
$\hat {\boldsymbol{e}}_3$
, as illustrated in figure 1. The basis vectors
$\hat {\boldsymbol{e}}_1$
,
$\hat {\boldsymbol{e}}_2$
,
$\hat {\boldsymbol{e}}_3$
provide the
$x$
-,
$y$
- and
$z$
-directions, respectively, in the lab frame. The filament has length
$L$
and cross-sectional radius
$a$
, and bending and twisting rigidities
$K_B$
and
$K_T$
, respectively. In this work, we take
$K_T=K_B$
. Based on previous results for the tip-driven follower force model (Clarke et al. Reference Clarke, Hwang and Keaveny2024), we anticipate the dynamics and bifurcations for this particular value to be representative of those for
$K_T/K_B \geqslant 1$
. This corresponds to a large portion of the biologically relevant range, as we discuss in Appendix A. The filament is submerged in an unbounded fluid of viscosity
$\eta$
. The filament’s centreline position is denoted by
$\boldsymbol{Y}(s,t)$
, where
$s\in [0,L]$
is the arc length, and
$t\in [0,\infty )$
is time. We describe filament deformation using a local orthonormal basis
$\left \{\hat {\boldsymbol{t}}(s,t),\hat {\boldsymbol{\mu }}(s,t),\hat {\boldsymbol{\nu }}(s,t)\right \}$
at each material point along the centreline. Here,
$\hat {\boldsymbol{\mu }}(s,t)$
and
$\hat {\boldsymbol{\nu }}(s,t)$
span the filament cross-section at
$s$
, while
$\hat {\boldsymbol{t}}(s,t)$
is constrained to be the unit tangent vector through
\begin{equation} \hat {\boldsymbol{t}}(s,t)=\dfrac {\partial \boldsymbol{Y}(s,t)}{\partial s}. \end{equation}
(a) An illustration of an active filament with cross-sectional radius
$a$
. The filament is clamped at its base in an unbounded domain and has a free end. Here,
$\hat {\boldsymbol{e}}_1$
,
$\hat {\boldsymbol{e}}_2$
and
$\hat {\boldsymbol{e}}_3$
are the unit vectors in the
$x$
-,
$y$
- and
$z$
-directions. (b,c,d,e) Diagrams describing the actuation models that we consider: (b) tip-driven follower force, (c) distributed follower force, (d) tip-driven surface flow, and (e) distributed surface flow. In the follower force models, tangential forces are applied to the filament (black arrows). For the surface flow models, tangential filament velocities and associated surface flows (blue arrows) are specified.
In the over-damped limit where inertia is negligible, the equations of motion are provided by the force and moment balances:
\begin{equation} \begin{aligned} \dfrac {\partial \boldsymbol{\varLambda }}{\partial s}+\boldsymbol{f}^H+\boldsymbol{f}^A&=0,\\ \dfrac {\partial \boldsymbol{M}}{\partial s}+\hat {\boldsymbol{t}}\times \boldsymbol{\varLambda }+\boldsymbol{\tau }^H&=0, \end{aligned} \end{equation}
where
$\boldsymbol{f}$
and
$\boldsymbol{\tau }$
denote force and torque (per unit length) acting on the filament, and superscripts
$(\cdot)^{H}$
and
$(\cdot)^{A}$
indicate the hydrodynamic and active contributions. The internal force and internal moment on the filament cross-section are denoted by
$\boldsymbol{\varLambda }(s,t)$
and
$\boldsymbol{M}(s,t)$
, respectively. The internal moment is given by the constitutive law (Landau & Lifshitz Reference Landau and Lifshitz1986)
\begin{equation} \boldsymbol{M}(s,t)=K_B\left (\hat {\boldsymbol{t}}\times \dfrac {\partial \hat {\boldsymbol{t}}}{\partial s}\right )+K_T\left (\hat {\boldsymbol{\nu }} \boldsymbol{\cdot }\dfrac {\partial \hat {\boldsymbol{\mu }}}{\partial s}\right ), \end{equation}
while
$\boldsymbol{\varLambda }(s,t)$
enforces the constraint (2.1).
We describe below how we impose the active and hydrodynamic forces and torques after discretising the filament into segments. For the hydrodynamics, we model the low Reynolds number mobility matrices for the segments using solutions to the Stokes equations for spherical particles with radius equal to the cross-sectional radius
$a$
. In order for our results to connect with the previous literature – e.g. Laskar et al. (Reference Laskar, Singh, Ghose, Jayaraman, Kumar and Adhikari2013), Chelakkot et al. (Reference Chelakkot, Gopinath, Mahadevan and Hagan2014), Bayly & Dutcher (Reference Bayly and Dutcher2016), Laskar & Adhikari (Reference Laskar and Adhikari2017), De Canio et al. (Reference De Canio, Lauga and Goldstein2017) and Krishnamurthy & Prakash (Reference Krishnamurthy and Prakash2023) – we utilise mobility relations for an unbounded fluid rather than those that incorporate a nearby no-slip surface. This choice allows us to isolate and compare the effects of differences in actuation on resulting filament motions.
2.1.2. Numerical discretisation
The filament is discretised into
$N$
segments of length
$\Delta L=2.2 a$
such that
$L=N\,\Delta L$
(see figure 2). In all cases, the filaments have
$N=20$
segments, yielding filament aspect ratio
$a/L=0.0227$
. We discretise (2.2) using a second-order central difference method, and after multiplying the resulting equations by
$\Delta L$
, we obtain the discrete force and moment balances
\begin{equation} \begin{aligned} \boldsymbol{F}_n^C+\boldsymbol{F}_n^H+\boldsymbol{F}_n^A & =\boldsymbol{0}, \\ \boldsymbol{T}_n^E+\boldsymbol{T}_n^C+\boldsymbol{T}_n^H & =\boldsymbol{0}, \end{aligned} \end{equation}
for each segment
$n$
, where the superscripts
$E$
,
$C$
,
$H$
and
$A$
denote elastic, constraint, hydrodynamic and active forces and torques, respectively. We have
\begin{align} \boldsymbol{F}_n^C=&\boldsymbol{\varLambda }_{n+1/2}-\boldsymbol{\varLambda }_{n-1/2}, \end{align}
\begin{align} \boldsymbol{T}_n^C=&\dfrac {\Delta L}{2}\hat {\boldsymbol{t}}_n\times \left (\boldsymbol{\varLambda }_{n+1/2}+\boldsymbol{\varLambda }_{n-1/2}\right ), \end{align}
\begin{align} \boldsymbol{T}_n^E=&\boldsymbol{M}_{n+1/2}-\boldsymbol{M}_{n-1/2}, \end{align}
A discretised portion of the filament, showing segments of length
$\Delta L$
. The hydrodynamic mobilities of the segments are approximated by those of spherical particles of radius
$a$
.
with half-indices representing the function values at the segment ends. The hydrodynamic force on segment
$n$
is
$\boldsymbol{F}_n^H=\boldsymbol{f}_n^H\Delta L$
, and the torque is
$\boldsymbol{T}_n^H=\boldsymbol{\tau }_n^H\Delta L$
. The hydrodynamic forces and torques and the active force
$\boldsymbol{F}_n^A=\boldsymbol{f}_n^A\Delta L$
differ in the follower force and surface flow models, as we explain below.
2.1.3. Follower force model
In the follower force model, filament motion is driven by the active forces
$\boldsymbol{F}^A_n$
on each
$n$
. We will consider two different distributions of follower force. For the tip-driven follower force model, the active force is applied only at the tip and is zero otherwise, namely,
$\boldsymbol{F}^A_N=-F^A\hat {\boldsymbol{t}}_N$
and
$\boldsymbol{F}^A_n=\boldsymbol{0}$
,
$n\neq N$
. In the distributed follower force model, we apply the same active force to each segment, such that
$\boldsymbol{F}^A_n=-F^A\hat {\boldsymbol{t}}_n$
for all
$n$
.
Once the non-hydrodynamic forces are known, the resulting filament motion can be determined. Due to the low Reynolds number associated with biofilament motion, the hydrodynamic forces and torques on the filament segments will be linearly related to the translational velocity
$\boldsymbol{V}_{\!n}$
and angular velocity
$\boldsymbol{\varOmega }_n$
on each segment
$n$
through
\begin{equation} \begin{bmatrix} \boldsymbol{V}\\ \boldsymbol{\varOmega } \end{bmatrix}=\boldsymbol{\mathcal{M}}\,\begin{bmatrix} -\boldsymbol{F}^H\\ -\boldsymbol{T}^H \end{bmatrix}, \end{equation}
where
$\boldsymbol{V} = [\boldsymbol{V}_1^{\textrm{T}},\boldsymbol{V}_2^{\textrm{T}},\ldots ,\boldsymbol{V}_N^{\textrm{T}}]^{\textrm{T}}\in \mathbb{R}^{3N\times 1}$
and
$\boldsymbol{\varOmega }= [\boldsymbol{\varOmega }_1^{\textrm{T}},\boldsymbol{\varOmega }_2^{\textrm{T}}, \ldots ,\boldsymbol{\varOmega }_N^{\textrm{T}} ]^{\textrm{T}}\in \mathbb{R}^{3N\times 1}$
, and similarly
$\boldsymbol{F}^H,\,\boldsymbol{T}^H\in \mathbb{R}^{3N\times 1}$
are the vectors containing all components of hydrodynamic force and torque, respectively, on all filament segments. Appearing in (2.8) is the configuration-dependent mobility matrix
$\boldsymbol{\mathcal{M}}\in \mathbb{R}^{6N\times 6N}$
relating the forces and torques on the segments and the resulting motion. Following previous work (Schoeller et al. Reference Schoeller, Townsend, Westwood and Keaveny2021; Clarke et al. Reference Clarke, Hwang and Keaveny2024), we utilise the pairwise RPY tensor (Wajnryb et al. Reference Wajnryb, Mizerski, Zuk and Szymczak2013) with the hydrodynamic radius set to
$a$
to provide the entries of the mobility matrix and hydrodynamic interactions between the segments.
Finally, we must ensure that the resulting filament motion satisfies the boundary conditions at the clamped and free ends. In our discretisation, the clamped end condition is enforced by additional Lagrange multipliers to ensure
$\boldsymbol{Y}_1(t)=\boldsymbol{0}$
and
$\hat {\boldsymbol{t}}_1(t)=\hat {\boldsymbol{e}}_3$
. At the free end, additional conditions are not needed as the follower force at the tip and moment-free condition are incorporated by setting
$\boldsymbol{F}^A_N=-F^A\hat {\boldsymbol{t}}_N$
with
$\boldsymbol{\varLambda }_{N+1/2}= \boldsymbol{0}$
and
$\boldsymbol{M}_{N+1/2}=\boldsymbol{0}$
, respectively.
2.1.4. The surface flow model
In the surface flow model where entrainment effects are included, we have
$\boldsymbol{F}^A_n=\boldsymbol{0}$
for all
$n$
, and introduce activity through a mobility relation based on the squirmer model (Lighthill Reference Lighthill1952; Blake Reference Blake1971). As a result of the surface flow, each segment
$n$
attempts to move in the direction
$-\hat {\boldsymbol{t}}_{n}$
at speed
$U = 2B_1/3$
, where
$B_1$
is the activity parameter that controls the surface flow strength. As a result, each active segment
$n$
generates degenerate quadrupolar (potential dipolar) flow
\begin{equation} \boldsymbol{u}_n(\boldsymbol{x})=\dfrac {1}{4\pi \eta r_{n}^3}\left (\mathbb{I}-3\,\hat {\boldsymbol{r}}_{n}\hat {\boldsymbol{r}}_{n}\right )\boldsymbol{H}_n, \end{equation}
where
$\boldsymbol{r}_n=\boldsymbol{x}-\boldsymbol{Y}_{\!n}$
,
$r_n=\| r_n\|$
,
$\hat {\boldsymbol{r}}_{n}=\boldsymbol{r}_n/r_n$
, and
$\boldsymbol{H}_n=(4/3)\pi \eta a^3\,B_1\hat {\boldsymbol{t}}_n$
is the degenerate quadrupole strength. Using this flow field, we can incorporate the effects of the surface flow into the segment mobility such that the resulting segment velocities are given by
\begin{equation} \begin{bmatrix} \boldsymbol{V}\\ \boldsymbol{\varOmega } \end{bmatrix}=\boldsymbol{\mathcal{M}}\,\begin{bmatrix} -\boldsymbol{F}^H\\ -\boldsymbol{T}^H \end{bmatrix}+\begin{bmatrix} \boldsymbol{\mathcal{S}}&\boldsymbol{0}\\\boldsymbol{0}&\boldsymbol{0} \end{bmatrix}\begin{bmatrix} \boldsymbol{H}\\ \boldsymbol{0} \end{bmatrix}, \end{equation}
where
$\boldsymbol{\mathcal{S}}\in \mathbb{R}^{3N\times 3N}$
is the squirming matrix, and
$\boldsymbol{H}=[\boldsymbol{H}_1^{\textrm{T}}, \boldsymbol{H}_2^{\textrm{T}}, \ldots , \boldsymbol{H}_N^{\textrm{T}}]^{\textrm{T}} \in \mathbb{R}^{3N\times 1}$
(see Appendix B for details). The off-diagonal components of
$\boldsymbol{\mathcal{S}}$
are obtained by evaluating (2.9), while the diagonal entries are of the form
$-1/(2\pi \eta a^3)$
and account for the self-generated velocity. By changing
$\boldsymbol{H}_n$
, we can control the distribution of activity along the filament. For the tip-driven surface flow model, we take
$\boldsymbol{H}_N=(4/3)\pi \eta a^3B_1\hat {\boldsymbol{t}}_N$
and
$\boldsymbol{H}_n = \boldsymbol{0}$
for
$n \neq N$
. When actuation is uniform along the filament length in the distributed surface flow model, we have
$\boldsymbol{H}_n=(4/3)\pi \eta a^3B_1 \hat {\boldsymbol{t}}_n$
for all
$n$
.
As in the follower force cases, we apply a clamped-end boundary condition at the base by enforcing
$\boldsymbol{Y}_1(t)=\boldsymbol{0}$
and
$\hat {\boldsymbol{t}}_1(t)=\hat {\boldsymbol{e}}_3$
. At the free end, we again have
$\boldsymbol{\varLambda }_{N+1/2}=\boldsymbol{0}$
and
$\boldsymbol{M}_{N+1/2}=\boldsymbol{0}$
.
We note that the surface flow model corresponds to the filament model used in Dutta et al. (Reference Dutta, Farhadifar, Lu, Kabacaoğlu, Blackwell, Stein, Lakonishok, Gelfand, Shvartsman and Shelley2024) to study cytoplasmic streaming. In their model, filament compression arises from an imposed follower force, but due to the opposite motor force on the fluid, the flow produced by the filament is only that due to the internal stresses. Equivalently, by dividing the force balance by the drag coefficient, one can instead interpret the follower force in their model as a self-generated velocity, and the internal stresses as the forces required to keep the filament elements from translating. This is the scenario described by the surface flow model, with the difference being that we also include higher-order force–quadrupole terms associated with the self-generated velocity.
2.1.5. Differential–algebraic system and time stepping
Once we have computed the segment translational and angular velocities, we can update their positions and orientations. For the orientations, we recognise that the local basis
$\{\hat {\boldsymbol{t}}_n,\hat {\boldsymbol{\mu }}_n,\hat {\boldsymbol{\nu }}_n\}$
that is linked to filament deformation is related to the basis in the lab frame
$\{\hat {\boldsymbol{e}}_1,\hat {\boldsymbol{e}}_2,\hat {\boldsymbol{e}}_3\}$
through a rotation. Thus we track the evolution of these vectors using unit quaternions
$\boldsymbol{q}=[q_0,\overline {\boldsymbol{q}}]=[q_0,q_1,q_2,q_3]\in \mathbb{R}^4$
, with
$\|\boldsymbol{q}\|^2=q_0^2+q_1^2+q_2^2+q_3^2=1$
. The quaternions map the standard basis to the local frame at each point along the filament via a matrix rotation
$\boldsymbol{R}(\boldsymbol{q}_{\!n})= [\hat {\boldsymbol{t}}_n,\hat {\boldsymbol{\mu }}_n,\hat {\boldsymbol{\nu }}_n ]$
, where
\begin{equation} \boldsymbol{R}(\boldsymbol{q})=\begin{bmatrix} 1-2q_2^2-2q_3^2&2(q_1q_2-q_0q_3)&2(q_1q_3+q_0q_2)\\ 2(q_1q_2+q_0q_3)&1-2q_1^2-2q_3^2&2(q_2q_3-q_0q_1)\\ 2(q_1q_3-q_0q_2)&2(q_2q_3+q_0q_1)&1-2q_1^2-2q_2^2 \end{bmatrix}. \end{equation}
To update the positions and orientations of the segments, we must integrate the following nonlinear differential–algebraic system:
\begin{align} \dfrac {{\textrm d}\boldsymbol{Y}_{\!n}}{\rm dt} = & \boldsymbol{V}_{\!n}, \end{align}
\begin{align} \dfrac {{\textrm d}\boldsymbol{q}_{\!n}}{{\textrm d}t}=&\dfrac {1}{2}[0,\boldsymbol{\varOmega }_n]\bullet \boldsymbol{q}_{\!n}, \end{align}
\begin{align} \boldsymbol{0}=&\boldsymbol{Y}_{n+1}-\boldsymbol{Y}_{\!n}-\dfrac {\Delta L}{2}\left (\hat {\boldsymbol{t}}_n+\hat {\boldsymbol{t}}_{n+1}\right ), \end{align}
for each segment
$n$
, where
$\boldsymbol{q}\bullet \boldsymbol{p}=[q_0,\overline {\boldsymbol{q}}]\bullet [p_0,\overline {\boldsymbol{p}}]=[q_0p_0-\overline {\boldsymbol{q}}\boldsymbol{\cdot }\overline {\boldsymbol{p}},p_0\overline {\boldsymbol{q}}+q_0\overline {\boldsymbol{p}}+\overline {\boldsymbol{q}}\times \overline {\boldsymbol{p}}]$
is the quaternion product, and the algebraic equations that must be satisfied are the discrete form of (2.1). We employ a second-order geometric backward differentiation scheme to discretise the equations in time. We solve the resulting nonlinear system of equations using Broyden’s method (Broyden Reference Broyden1965) with tolerance
$10^{-12}$
. A detailed study of this numerical technique is presented in Schoeller et al. (Reference Schoeller, Townsend, Westwood and Keaveny2021).
To summarise, at each time step, we obtain
$\boldsymbol{M}$
from (2.3),
$\boldsymbol{F}^H$
and
$\boldsymbol{T}^H$
using (2.4),
$\boldsymbol{V}$
and
$\boldsymbol{\varOmega }$
from (2.8) or (2.10) depending on the model,
$\boldsymbol{Y}$
from (2.12), and
$\boldsymbol{q}$
from (2.13). We treat
$\boldsymbol{\varLambda }$
as a Lagrange multiplier for this differential–algebraic system, so that at each time step, the inextensibility constraint (2.14) is satisfied. In our simulations, the typical time step size is
$\Delta t=T/400$
, where
$T$
is the period of one cycle for periodic states, or based on the strongest frequency in non-periodic states.
2.1.6. Non-dimensional parameters
In the rest of this paper, we use the filament length
$l^*=L$
as the length scale, relaxation time
$t^*={\eta \,L^4}/{K_B}$
as the time scale, and
$F^*={K_B}/{L^2}$
as the force scale. Accordingly, we introduce
\begin{equation} f=\dfrac {F^A{\kern-1pt}L^2}{K_B},\quad b_1=\dfrac {\eta L^3B_1}{K_B} \end{equation}
as the non-dimensional actuation parameters for follower force models and surface flow models, respectively.
2.2. Dynamical systems analysis
In this subsection, we briefly describe the computational tools that we employ to perform a dynamical systems analysis of the filament models. In order to use these tools, we must first introduce an appropriate state vector. Unfortunately, the unit quaternions are not suitable since the sum of two unit quaternions is not itself a unit quaternion. In our case, it is therefore convenient to instead employ the effective Lie algebra element (Clarke et al. Reference Clarke, Hwang and Keaveny2024)
${\boldsymbol{U}}_n\in \mathbb{R}^3$
for each segment
$n$
to describe the rotation of the frame
$\{\hat {\boldsymbol{t}}_n,\hat {\boldsymbol{\mu }}_n,\hat {\boldsymbol{\nu }}_n\}$
relative to a state where the filament is straight and upright such that
$\hat {\boldsymbol{t}}_n=\hat {\boldsymbol{e}}_3$
,
$\hat {\boldsymbol{\mu }}_n=\hat {\boldsymbol{e}}_1$
and
$\hat {\boldsymbol{\nu }}_n=\hat {\boldsymbol{e}}_2$
for all
$n$
. In simple terms, an effective Lie algebra element is a 3-D generalisation of the tangent angle formulation that can be employed in two dimensions. Since the quaternion corresponding to the upright state is
\begin{equation} \boldsymbol{q}_{\textit{eq}}=\begin{bmatrix} \dfrac {1}{\sqrt {2}},0,-\dfrac {1}{\sqrt {2}},0 \end{bmatrix}, \end{equation}
with a conjugate
$\boldsymbol{q}^*_{\textit{eq}}=\begin{bmatrix} q_0,-\overline {\boldsymbol{q}} \end{bmatrix}$
, and the rotation of segment
$n$
relative to this state is given by the quaternion product
$\boldsymbol{p}=\boldsymbol{q}_{\!n}\bullet \boldsymbol{q}^*_{\textit{eq}}=\begin{bmatrix} p_0,\overline {\boldsymbol{p}} \end{bmatrix}$
, the effective Lie algebra element associated with a quaternion
$\boldsymbol{q}_{\!n}$
is
\begin{equation} {\boldsymbol{U}}_n = 2\arccos {\left (p_0\right )}\,\hat {\boldsymbol{p}}, \end{equation}
where
$\hat {\boldsymbol{p}}={\overline {\boldsymbol{p}}}/{\|\overline {\boldsymbol{p}}\|}$
. The state vector is then given by
${\boldsymbol{U}} = [{\boldsymbol{U}}_1^{\textrm{T}},{\boldsymbol{U}}_2^{\textrm{T}},\ldots ,{\boldsymbol{U}}_N^{\textrm{T}}]^{\textrm{T}}\in \mathbb{R}^{3N\times 1}$
, with
${\boldsymbol{U}}_1 = 0$
to incorporate the clamped-end condition. We refer the reader to Appendix C for more about properties of the effective Lie algebra element.
2.2.1. Jacobian-free Newton–Krylov method
The filament models display several time-periodic solutions after the steady, upright state becomes unstable. To find these time-periodic solutions, we first shift our thinking and consider the filament simulation as a means of computing the flow
$\boldsymbol{\phi }^t({\boldsymbol{U}})$
that maps solutions
${\boldsymbol{U}}(0)$
to
${\boldsymbol{U}}(t)$
. We then seek effective Lie algebra elements that satisfy
\begin{equation} \boldsymbol{F}({\boldsymbol{U}},T)=\boldsymbol{\phi }^{\textrm{T}}({\boldsymbol{U}})-{\boldsymbol{U}}=\boldsymbol{0} \end{equation}
for some non-zero, but unknown period
$T$
. We obtain solutions to (2.18) using a Newton–Krylov method with a GMRES-hook step, described in Viswanath (Reference Viswanath2007, Reference Viswanath2009). Our implementation closely follows that of Willis (Reference Willis2019). We continue the iterations until
$\| \boldsymbol{F}({\boldsymbol{U}},T)\|\lt 10^{-10}$
, with
$\| {\cdot }\|$
as the 2-norm operator.
2.2.2. Linear stability analysis
Along with finding time-periodic solutions, we would like to assess their stability, as well as the stability of the steady, upright state. The time evolution of the state vector can be expressed as
\begin{equation} \dfrac {{\textrm d} {\boldsymbol{U}}}{{\textrm d} t}=\boldsymbol{f}({\boldsymbol{U}};t). \end{equation}
In linear stability analysis, we are interested in the evolution of small perturbations
$\delta {\boldsymbol{U}}$
around a (steady or time-periodic) base state
${\boldsymbol{U}}_0$
. Growing perturbations indicate unstable base states, whereas dampened perturbations indicate linearly stable base states. We linearise (2.19) around
${\boldsymbol{U}}_0$
, and obtain the linearised dynamics for the perturbation as
\begin{equation} \dfrac {{\textrm d}\delta {\boldsymbol{U}}}{{\textrm d} t}=\boldsymbol{A}\delta {\boldsymbol{U}}, \end{equation}
where
$\boldsymbol{A}=\left . {\partial \boldsymbol{f}}/{\partial {\boldsymbol{U}}}\right |_{{\boldsymbol{U}}_0}$
.
Equation (2.20) is linear and has the principal fundamental matrix solution
$\boldsymbol{\varPhi }(t;{\boldsymbol{U}}_0(t))$
such that the general solution is
$\delta {\boldsymbol{U}}(t)=\boldsymbol{\varPhi }(t;{\boldsymbol{U}}_0(t))\,\delta {\boldsymbol{U}}_0$
. For a steady base state with
${\textrm d}{\boldsymbol{U}}_0/{\textrm d}t=\boldsymbol{0}$
,
\begin{equation} \boldsymbol{\varPhi }(t;{\boldsymbol{U}}_0(t))={\rm e}^{t\boldsymbol{A}}. \end{equation}
If
${\boldsymbol{U}}_0$
is instead time-periodic such that
${\boldsymbol{U}}_0(t)={\boldsymbol{U}}_0(t+T)$
, then Floquet’s theorem states that there is a Floquet normal form of the fundamental solution matrix,
\begin{equation} \boldsymbol{\varPhi }(t;{\boldsymbol{U}}_0)=\boldsymbol{P}(t)\,{\rm e}^{t\boldsymbol{B}}, \end{equation}
with a
$T$
-periodic vector matrix
$\boldsymbol{P}(t)$
, and a matrix operator
$\boldsymbol{B}$
(Perko Reference Perko2008). Notice that under the action of linearised Poincaré maps
${\boldsymbol{U}}(t)\mapsto {\boldsymbol{U}}(t+T)$
, a periodic state is steady, and the dynamics of a perturbation after one period is governed by the monodromy matrix
${\rm e}^{T\boldsymbol{B}}=\boldsymbol{\varPhi }^{-1}(0;{\boldsymbol{U}}_0)\,\boldsymbol{\varPhi }(T;{\boldsymbol{U}}_0)$
. This is equivalent to the fundamental matrix solution evaluated at
$t = T$
. For periodic solutions, the matrix
$\boldsymbol{B}$
is analogous to the matrix
$\boldsymbol{A}$
in the steady case.
The stability of a base state depends on the eigenvalues
$\lambda _i$
of
$\boldsymbol{A}$
or
$\boldsymbol{B}$
. A solution is linearly stable if all of its eigenvalues have non-positive real parts. We calculate the eigenvalues
$\mu _i$
of
$\boldsymbol{\varPhi }(T;{\boldsymbol{U}}_0)$
, where
$T$
is the period for periodic states and a sufficiently small time for steady base states. We compute the eigenvalues of
$\boldsymbol{\varPhi }(T;{\boldsymbol{U}}_0)$
using Arnoldi iteration, which does not require explicit computation of
$\boldsymbol{\varPhi }(T;{\boldsymbol{U}}_0)$
, but requires only its action on a set of vectors. At each Arnoldi iteration, an upper Hessenberg matrix is calculated, whose eigenvalues approximate those of
$\boldsymbol{\varPhi }(T;{\boldsymbol{U}}_0)$
. We record the first 13 eigenvalues as a vector at each Arnoldi iteration, and continue until the norm of the difference between consecutive iterations is lower than
$10^{-8}$
.
The eigenvalues
$\lambda _i$
are related to the eigenvalues of
$\boldsymbol{\varPhi }(T;{\boldsymbol{U}}_0)$
via
\begin{equation} \lambda _i=\dfrac {1}{T}\ln {(\mu _i)}. \end{equation}
2.2.3. Spectral proper orthogonal decomposition
Along with the upright steady-state and time-periodic solutions, the filament models can also admit more complex solutions, such as quasi-periodic states. To study these states, we extract their dominant motions at each frequency
$\omega$
using spectral proper orthogonal decomposition (SPOD). We utilise the implementation developed in Towne, Schmidt & Colonius (Reference Towne, Schmidt and Colonius2018), and summarise the procedure in our context below.
Given effective Lie algebra elements
${\boldsymbol{U}}(s,t)$
,
$\boldsymbol{V}(s,t)$
and the inner product
$\langle {\boldsymbol{U}},\boldsymbol{V}\rangle _{s,t}=\int_{-\infty }^\infty \int_0^L\boldsymbol{V}^*(s,t)\,{\boldsymbol{U}}(s,t)\,{\textrm d}s\,{\textrm d}t$
, with the superscript
$(\cdot)^*$
denoting the Hermitian transpose, SPOD modes are obtained by solving the optimisation problem
\begin{equation} \max _\phi \dfrac {E\{\lvert \langle {\boldsymbol{U}}(s,t),\boldsymbol{\phi }(s,t)\rangle _{s,t}\rvert ^2\}}{\langle \boldsymbol{\phi }(s,t),\boldsymbol{\phi }(s,t)\rangle _{s,t}}, \end{equation}
where
$E\{{\cdot }\}$
denotes ensemble average, for some function
$\boldsymbol{\phi }(s,t)$
that approximates
${\boldsymbol{U}}(s,t)$
on average. With a Fourier transform
$\boldsymbol{\phi }(s,t)=\int _{-\infty }^\infty \boldsymbol{\psi }(s,\omega )\,{\rm e}^{{\rm i}2\pi \omega t}\, {\textrm d}\omega$
, the optimisation problem in (2.24) can be formulated at each frequency
$\omega$
. The solution to the optimisation problem is obtained by solving the eigenvalue problem
\begin{equation} \int_0^L\boldsymbol{S}(s,s^\prime ,\omega ) \,\boldsymbol{\psi }(s^\prime ,\omega )\,{\textrm d}s^\prime =\lambda (\omega )\,\boldsymbol{\psi }(s,\omega ), \end{equation}
where
$\boldsymbol{S}(s,s^\prime ,\omega )$
is the cross-spectral density tensor, defined as the Fourier transform of the two-point space–time correlation tensor. For further details, the reader may refer to Towne et al. (Reference Towne, Schmidt and Colonius2018).
Given time series data of the effective Lie algebra element, we block the data into multiple realisations using Welch’s method, and take the Fourier transform with a Hamming window to reduce spectral leakage due to non-periodicity of each block. By ensemble averaging these Fourier realisations, the cross-spectral density matrix is obtained for each frequency, and the SPOD modes are calculated by performing an eigenvalue decomposition of the covariance operator in the frequency domain. The SPOD analysis based on these methods is performed using the software from Towne et al. (Reference Towne, Schmidt and Colonius2018), available at https://github.com/SpectralPOD/spod_matlab.
2.2.4. Bisection method
As we discuss in § 3, we find a range of actuation for which filament dynamics exhibits bistability. Here, the two periodic states of whirling and beating are connected through an unstable quasi-periodic state QP1. To compute QP1, we use the bisection algorithm from Skufca, Yorke & Eckhardt (Reference Skufca, Yorke and Eckhardt2006) and Schneider & Eckhardt (Reference Schneider and Eckhardt2009).
Given two initial conditions
${\boldsymbol{U}}_{{w}}$
and
${\boldsymbol{U}}_{{b}}$
that converge with time to whirling and beating, respectively, we compute the bisected initial condition as
\begin{equation} {\boldsymbol{U}}_{\textit{QP}}=\dfrac {{\boldsymbol{U}}_{{w}}+{\boldsymbol{U}}_{{b}}}{2}. \end{equation}
With this initial condition, we evolve the system to long times, and check if the filament dynamics is attracted to either of the periodic states. Once the solution reaches (criteria provided in § C.4) either of the states, we update
${\boldsymbol{U}}_{{w}}$
or
${\boldsymbol{U}}_{{b}}$
accordingly, and perform a bisection again according to (2.26). We repeat this process until
$\|{\boldsymbol{U}}_{{w}}-{\boldsymbol{U}}_{{b}}\|\lt 10^{-8}$
.
3. Results and discussions
We investigate the state space of the four different active filament models as we vary the actuation. Follower force models use the dimensionless follower force magnitude
$f$
as the actuation parameter, while surface flow models use the dimensionless surface flow magnitude
$b_1$
. These two parameters are directly related. A segment with follower force magnitude
$f$
has active force
$F^A=f{\!K}_{\!B}/\!L^2$
, while a surface flow segment with
$b_1$
experiences drag
$F_d=4\pi ab_1K_B/L^3$
. Equating these two forces yields that the follower force and the surface flow are related through
$b_1\approx 3.5f$
.
In exploring filament dynamics, we divide our presentation of the results into two actuation ranges: low-to-moderate and high. Primary and secondary filament instabilities occur within the low-to-moderate actuation regime, where previous studies have revealed similar filament dynamics (De Canio et al. Reference De Canio, Lauga and Goldstein2017; Laskar & Adhikari Reference Laskar and Adhikari2017; Ling et al. Reference Ling, Guo and Kanso2018; Clarke et al. Reference Clarke, Hwang and Keaveny2024). At high actuation, these states are unstable, and we observe greater differences in the dynamics given by the models.
3.1. Dynamic states at low-to-moderate actuation
We begin by classifying the dynamic states up to moderate actuation, and expand on previous results by directly comparing the different models. As stated earlier, we consider two actuation models, namely follower force (FF) and surface flow (SF), and two actuation profiles, namely concentrated at the tip (Tip) and uniformly distributed along the length of the filament (Dis). In figure 3, we present the bifurcation diagram for each model up to moderate actuation, which shows the possible filament states and transitions between these states. We indicate linearly stable states with solid lines and unstable states with dashed lines. We use the dimensionless mean bending energy
\begin{equation} \overline{E_b}=\dfrac {1}{T}\int_0^{\textrm{T}}\int_0^L\dfrac {1}{2}\dfrac {\partial ^2\boldsymbol{Y}(s,t)}{\partial s^2} \boldsymbol{\cdot }\dfrac {\partial ^2\boldsymbol{Y}(s,t)}{\partial s^2}\,{\textrm d}s\,{\textrm d}t \end{equation}
to distinguish different states. By comparing figures 3(a,c) with figures 3(b,d), we see that the bending energy is overall lower for the tip-driven cases than for the distributed cases. This implies that models with a distributed actuation profile produce filament states with higher curvature.
Bifurcation diagrams up to moderate actuation for (a) tip-driven follower force, (b) distributed follower force, (c) tip-driven surface flow, and (d) distributed surface flow models. Stable states are marked with solid lines, whereas unstable ones use dashed lines. We observe similar bifurcation patterns and states in low-to-moderate actuation levels. We plot the non-trivial filament states in figure 4.
Apart from small quantitative differences, the filament models exhibit similar states and dynamics in the low-to-moderate actuation regime. In all active filament models, we observe four common states: upright, whirling, beating and QP1. We plot these states (except upright) in figure 4. The increasing direction of the time is indicated by darker filament colours. The pathline of the filament tip is a closed curve for periodic states, whereas for QP1 it lies on a torus. To aid the viewer, we also plot the shadow of the filament at its base. Upright is a steady state (fixed point); whirling and beating are time-periodic; and QP1 is a quasi-periodic state.
Non-trivial states observed in all active filament models up to moderate actuation: (a) whirling, (b) beating, (c) QP1. All states are given at
$b_1=284$
for the distributed surface flow model (Dis SF). The filament’s colour gets darker as time passes. The filament’s shadow is given in grey, and its tip pathline is given in blue. We provide a video description of the states in supplementary movie 1.
From the bifurcation diagram in figure 3, we see that regardless of the type and distribution of actuation, we observe a double Hopf bifurcation when the upright state becomes unstable, creating whirling and beating (see § 3.2 for further details). Later, whirling undergoes a Hopf bifurcation, and beating undergoes a pitchfork bifurcation; QP1 connects to whirling at one end, and beating at the other (see § 3.3) through these bifurcations.
3.2. Primary instability of the upright filament
The upright filament is the trivial state and corresponds to a straight filament aligned with
$\hat {\boldsymbol{e}}_3$
. Actuation produces internal stresses, and like beams under compression, when these stresses exceed a threshold, the filament buckles. While the critical actuation value is model-dependent, this instability is observed in each model. We compute the eigenvalues of the linearised operator around the upright state, and in figure 5 plot the first two eigenvalues for each model as functions of the actuation parameter. The real and imaginary parts of the eigenvalues are coloured red and blue, respectively. The vertical dashed line indicates when the real part of the leading eigenvalue is zero. The leading eigenvalues and their dependence on the actuation parameter are qualitatively similar in all cases. We note, however, that the critical value of actuation is higher in the surface flow models. Consistent with De Canio et al. (Reference De Canio, Lauga and Goldstein2017), this indicates that the hydrodynamic entrainment provides a mechanism that delays the bifurcation.
Leading two eigenvalues of the upright state for (a) tip-driven follower force, (b) distributed follower force, (c) tip-driven surface flow, and (d) distributed surface flow models. Real and imaginary parts of the eigenvalues are shown in red and blue, respectively. The double Hopf bifurcation is marked with a vertical dashed line.
In all cases, as the actuation increases, we observe that identical sets of two real eigenvalue branches (only one set is shown in figure 5) merge. This point marks the transition of the upright filament from a stable node to a stable focus, i.e. the two sets of two real eigenvalues transition to two sets of a complex conjugate eigenvalue pair. Due to the symmetry of the upright state and no preferred directionality in the
$xy$
-plane, the leading eigenvalues are identical complex conjugate pairs, with one eigenmode corresponding to bending in the
$xz$
-plane, and the other to bending in the
$yz$
-plane.
These two sets of complex conjugate eigenmodes eventually become unstable at the same actuation value. As described in Clarke et al. (Reference Clarke, Hwang and Keaveny2024), when this occurs, two periodic states, whirling and beating, emerge through a double Hopf bifurcation, with whirling being stable, and beating being unstable. In all models, whirling involves a filament with time-independent curvature rotating about the
$z$
-axis at a constant speed (figure 4
a), while beating is a planar state where the tip traces an infinity-symbol-shaped path (figure 4
b). These states have been observed in the tip-driven follower force model (Chelakkot et al. Reference Chelakkot, Gopinath, Mahadevan and Hagan2014; De Canio et al. Reference De Canio, Lauga and Goldstein2017; Clarke et al. Reference Clarke, Hwang and Keaveny2024) and the tip-driven active colloid model (Laskar et al. Reference Laskar, Singh, Ghose, Jayaraman, Kumar and Adhikari2013; Laskar & Adhikari Reference Laskar and Adhikari2017).
As stated above, due to the rotational symmetry of the upright state, there are two decoupled orthonormal eigenmodes at the instability threshold, which individually describe planar beating oscillations. Without loss of generality, we can define these two orthonormal eigenmodes using the effective Lie algebra element: one for the
$xz$
-plane,
$\boldsymbol{\nu }_x(s)=[0,v(s),0]^{\textrm{T}}$
, and one for the
$yz$
-plane,
$\boldsymbol{\nu }_y(s)=[v(s),0,0]^{\textrm{T}}$
, where
$v(s)$
describes the filament shape. For surface flow models, periodic solutions close to the bifurcation take the form
\begin{equation} {\boldsymbol{U}}(s,t)\approx \left (A_x\,\boldsymbol{\nu }_x(s)+A_y\,\boldsymbol{\nu }_y(s)\right )\sqrt {b_1-b_1^*}\,{\rm e}^{{\rm i}\omega ^*t}+\text{c.c.}, \end{equation}
where
$b_1^*$
and
$\omega ^*$
are the surface flow magnitude and the frequency at the bifurcation, respectively. The amplitudes
$A_x, A_y\in \mathbb{C}$
are determined from the initial value
${\boldsymbol{U}}(s,0)$
of the periodic state. The form of the solution is the same for the follower force model with
$b_1$
replaced by
$f$
(see Clarke et al. Reference Clarke, Hwang and Keaveny2024).
For the motion to remain planar, the eigenmodes must remain in phase, so
$\operatorname {Im} (\langle A_x,A_y\rangle )=0$
for beating. For whirling,
$A_y = {\rm i}A_x$
, meaning that one mode lags behind the other by a
$\pi /2$
phase. These results are derived using the connection between filament states and effective Lie algebra elements in Appendix C.
In figure 6, we plot the leading eigenmode at the bifurcation in the
$xz$
-plane at different phases. The curved filament shape in each model is associated with its buckling; however, there are some qualitative differences between the modes for the different models.
Leading eigenmodes at the double Hopf bifurcation for (a) tip-driven follower force, (b) distributed follower force, (c) tip-driven surface flow, and (d) distributed surface flow models. The plots show the eigenmode on the
$xz$
-plane at different phases. A linear combination of the mode with its counterpart on the
$yz$
-plane creates a beating state, while a linear combination with a
$\pi /2$
phase lag creates a whirling state, as depicted in figure 4.
We investigate further the buckling instability by considering the force balance in the upright state. Linear stability analysis using resistive force theory (De Canio et al. Reference De Canio, Lauga and Goldstein2017; Ling et al. Reference Ling, Guo and Kanso2018; Schnitzer Reference Schnitzer2025) reveals that the internal stress distribution appears in the linear operator, and hence contributes to the eigenvalues of the primary instability. Motivated by this, we determine the internal stress for each model, with the aim of explaining the differences in the mode shapes in figure 6.
In follower force models, the hydrodynamic force and torque are zero due to the absence of induced surface flow and the filament being at rest. Thus the internal force results solely from balancing the prescribed active force, given by
$\boldsymbol{\varLambda }(s)=\boldsymbol{\varLambda }(L)+\int _s^L\boldsymbol{f}^A(s')\,{\textrm d}s'$
with the appropriate boundary condition at
$\boldsymbol{\varLambda }(L)$
. For the tip-driven case, we take
$\boldsymbol{f}^A(s)=\boldsymbol{0}$
and
$\boldsymbol{\varLambda }(L)=-F^A\,\hat {\boldsymbol{e}}_3$
, where
$F^A$
is the prescribed dimensional follower force magnitude. We note that this is equivalent to our numerical framework where instead activity is incorporated through the body force rather than the boundary condition; see Ling et al. (Reference Ling, Guo and Kanso2018). Consequently, the internal force distribution is
$\boldsymbol{\varLambda }(s)=-F^A\,\hat {\boldsymbol{e}}_3$
. In the distributed case, there is a uniform force distribution
$\boldsymbol{f}^A(s)=-F^A/\Delta L\,\hat {\boldsymbol{e}}_3$
. Hence
$\boldsymbol{\varLambda }(s)=-F^A(L-s)/\Delta L\,\hat {\boldsymbol{e}}_3$
.
In the surface flow model, however, the hydrodynamics plays a role as the internal stress is obtained by solving a resistance problem. For the surface flow model, the filament motion is given by
\begin{equation} \begin{bmatrix} \boldsymbol{V}-\boldsymbol{V}^{\textrm{sf}}\\ \boldsymbol{\varOmega } \end{bmatrix}=\boldsymbol{\mathcal{M}}\begin{bmatrix} -\boldsymbol{F}^H\\ -\boldsymbol{T}^H \end{bmatrix}, \end{equation}
where
$\boldsymbol{V}^{\textrm{sf}}=\boldsymbol{\mathcal{S}}\boldsymbol{H}$
is the surface velocity. In the upright state,
$\boldsymbol{V}=\boldsymbol{0}$
and
$\boldsymbol{\varOmega }=\boldsymbol{0}$
, yet non-zero hydrodynamic forces arise due to surface flow activity. As a result, internal stresses are given by the resistance problem
\begin{equation} \begin{bmatrix} \boldsymbol{F}^C\\\boldsymbol{T}^C \end{bmatrix}=-\boldsymbol{\mathcal{M}}^{-1}\begin{bmatrix} \boldsymbol{V}^{\textrm{sf}}\\\boldsymbol{0} \end{bmatrix}. \end{equation}
After solving this linear system for the tip-driven and distributed surface velocities, the internal stress can be determined by relating
$\boldsymbol{F}^C$
to
$\boldsymbol{\varLambda }$
using (2.5).
In figure 7, we plot the internal stress along the filament at buckling for the different models. For the distributed models, the stress is accumulated at the base and decays monotonically as one moves towards the free end. For the tip-driven follower force model, the stress is constant along the length. For the tip-driven surface flow case, however, after first increasing only gradually with arc length, the stress exhibits a peak at the free end. In this case, there is a single active segment that attempts to move towards the filament base. As a result, that segment and those that are nearby (due to the non-local nature of the hydrodynamics) require a larger force in order not to move. Once sufficiently far from this segment at the free end, the stress profile attains a near constant value.
Force distribution along the filament for different active filament models at the double Hopf bifurcation. Follower force models are plotted in red, while surface flow models are plotted in blue. Tip-driven models are given as dashed lines, whereas distributed models are plotted as solid lines.
As a result of the differences in these stress distributions, the eigenmodes for tip-driven models have a wider tip span (see figures 6 a,c), whereas distributed models show narrower tips and wider bases (see figures 6 b,d). Figure 7 also shows that the buckling threshold is generally higher for surface flow models. The surface flow, in opposing direction of compression, allows the filament to support larger loads, thus delaying the bifurcation.
3.3. Secondary instabilities and quasi-periodic transition
The double Hopf bifurcation creates two limit cycles: whirling and beating. In this subsection, we examine the instabilities associated with these periodic states. Near the double Hopf bifurcation, whirling is stable, while beating is unstable. However, as shown in figure 3, beating becomes stable and whirling becomes unstable at moderate actuation. We focus on understanding the changes in the stability of these periodic states, and the origin of this transition.
Using the Jacobian-free Newton–Krylov algorithm described in § 2.2.1, we track each state as the actuation increases. Increased actuation leads to higher internal stresses, filament deformation and bending energy, as shown in figure 3(a) for the tip-driven follower force model. For whirling, the radius of the circular path of each point along the filament increases, while for beating, the filament sweeps out a greater area. We also observe that periodic states tend to have higher frequencies as the actuation increases.
We perform a linear stability (Floquet) analysis for each state along these branches. The leading eigenvalues from the linear stability analysis of the whirling state are shown in figure 8. The leading eigenvalue is complex, meaning that whirling becomes unstable through a Hopf bifurcation. We plot the leading eigenvector of the monodromy matrix at different phases for each filament model in figure 9. These eigenvectors represent the instability dynamics, with the time-periodicity of whirling removed. Comparing the mode shapes for the tip-driven and distributed models (figures 9(a,b) and 9(c,d)), we observe that the unstable mode shape narrows towards the filament tip for distributed models, while it is widest at the tip for tip-driven models.
Leading eigenvalues of the whirling state for (a) tip-driven follower force, (b) distributed follower force, (c) tip-driven surface flow, and (d) distributed surface flow models. Real and imaginary parts of the leading eigenvalues are shown in red and blue, respectively. The Hopf bifurcation is marked with a vertical dashed line.
Leading eigenmodes at the Hopf bifurcation of whirling for (a) tip-driven follower force, (b) distributed follower force, (c) tip-driven surface flow, and (d) distributed surface flow models. Different filament colours correspond to different phases of the mode.
Since whirling can be viewed as a steady state in a rotating frame of reference, compression and shear forces are time-invariant in that frame. Figure 10 shows the force distribution for each active filament model. Compressive forces are transferred along the filament and balanced by the constraint force at the base. Shear is balanced by elastic forces that arise due to filament curvature, along with the viscous forces as a result of filament motion. We find that shear is much weaker than compression along the filament. As for the steady-state profile, we observe that distributed models have an almost linearly decaying tangential stress distribution (see figure 10). Due to its suspected presence in the linearised operator, this distribution of tangential stress produces a bottom-heavy mode shape, as seen in figure 9(d) for the distributed surface flow model. In contrast, the tip-driven models exhibit an increase in the tangential stress with arc length. Again, the surface flow model exhibits an additional peak at the free end. This form of the distribution results in a mode shape with a wide tip (figure 9 c).
Compression distribution along the length of a filament for different active filament models at the Hopf bifurcation. Follower force models are plotted in red, while surface flow models are plotted in blue. Tip-driven models are given as dashed lines, whereas distributed models are plotted as solid lines. Shear forces are much weaker than compression, and do not significantly contribute to the dynamics.
We also assess the stability of planar beating, which appears in all models as an unstable solution after the filament buckles. Figure 11 shows the leading eigenvalues from the stability analysis for the beating state. In all cases, the leading eigenvalue is purely real, so beating becomes stable through a pitchfork bifurcation. For our stability calculations, we use base states that beat in the
$xz$
-plane. Additionally, as the resulting eigenmodes will be dependent on the phase of the base state, we set
$t=0$
in our stability calculations to correspond to the instance when the filament tip attains its maximum amplitude in the
$-x$
-direction. The unstable modes at the bifurcation for each model are shown in figure 12. We see that the eigenmodes are planar, but in the
$yz$
-plane, orthogonal to the base state. This indicates that as actuation decreases, planar beating becomes unstable through out-of-plane tipping.
Leading eigenvalues of the beating state for (a) tip-driven follower force, (b) distributed follower force, (c) tip-driven surface flow, and (d) distributed surface flow models. Real and imaginary parts of the leading eigenvalues are shown in red and blue, respectively. The pitchfork bifurcation is marked with a vertical dashed line.
Leading eigenmodes at the pitchfork bifurcation for (a) tip-driven follower force, (b) distributed follower force, (c) tip-driven surface flow, and (d) distributed surface flow models. The base beating state for each model is at the same phase (when the tip of the filament is at its leftmost position) to remove phase dependency of the Floquet modes. The instability mechanism is out-of-plane tipping, as the most unstable mode is a planar filament state lying orthogonal to the beating plane.
These two bifurcations transition filament dynamics from a range where only whirling is stable to one where only beating is stable. Bridging these two bifurcations is the quasi-periodic state QP1; see figure 3. In this study, we find that while QP1 is a generic feature of active filament models, its stability appears to be sensitive to the details of the model. In the tip-driven surface flow model, whirling and beating undergo supercritical bifurcations, with a stable QP1 emerging from the Hopf bifurcation and annihilating at the pitchfork bifurcation. In contrast, all other models feature subcritical bifurcations, leading to bistability between whirling and beating, and an unstable QP1. In these cases, we obtain the unstable QP1 state using the bisection algorithm. From figure 3, we observe that distributed models have a wider range of bistability. In Clarke et al. (Reference Clarke, Hwang and Keaveny2024), the tip-driven follower force filament anchored to a no-slip surface was instead found to exhibit a stable QP1 state, indicating that this change in hydrodynamics is sufficient to change the stability of QP1. In addition, it has been shown that even in an unbounded fluid, increasing the filament aspect ratio by adding more segments in the discrete filament model is also able to stabilise QP1 (Clarke Reference Clarke2025).
We track QP1 throughout its existence as actuation increases, plotting its shape over time for the tip-driven surface flow model in figure 13. We see that QP1 carries the characteristics of both whirling and beating (see supplementary movie 1). Close to the Hopf bifurcation, at
$b_1=741.1$
, QP1 is a wobbly whirling state (see figure 13
a), with the tip following a near-circular path. As the actuation increases, beating characteristics become stronger. At
$b_1=744.9$
, near the pitchfork bifurcation (see figure 13
c), QP1 beats in a plane that is slowly whirling about the
$z$
-axis. Its tip pathline is similar to that of a beating. For stable QP1 states, the strength of the whirling mode decreases with actuation, while that of the beating mode increases. Conversely, for unstable QP1 states, the strength of the whirling mode increases with actuation, while that of the beating mode decreases in the bistable regime.
The QP1 state throughout its existence for the tip-driven surface flow model: (a) close to Hopf bifurcation at
$b_1=741.1$
, (b) at an intermediate value
$b_1=743.0$
, (c) close to pitchfork bifurcation at
$b_1=744.9$
. The filament’s colour gets darker in the increasing direction of time. Filament shadow is given in grey, and its tip pathline is given in blue.
Examining the tip pathline of the QP1 state in figure 14, we see that it lies on a torus, which gives a qualitative indication of the quasi-periodicity of this solution. To examine this state quantitatively, we extract dominant frequencies and modes of QP1 using the SPOD approach described in § 2.2.3 (Towne et al. Reference Towne, Schmidt and Colonius2018). The SPOD spectrum (which is the leading eigenvalue given by SPOD at each frequency) of QP1 for the tip-driven surface flow model with
$b_1=741.1$
is shown in figure 15(a). Here, QP1 has two incommensurate dominant frequencies, consistent with its quasi-periodic nature. These frequencies,
$40.8\,\eta L^4/K_B$
and
$29.4\,\eta L^4/K_B$
, are close to those of the whirling and beating states,
$40.8\,\eta L^4/K_B$
and
$35.3\,\eta L^4/K_B$
, respectively, at the same actuation. The remaining peaks in the frequency spectrum reflect nonlinear interactions of the two dominant frequencies.
The tip pathline of the stable QP1 state at
$b_1=741.1$
for the tip-driven surface flow model. The filament tip follows a path on a 2-torus.
(a) Dimensionless SPOD frequency spectrum of the stable QP1 at
$b_1=741.1$
(close to the Hopf bifurcation) for the tip-driven surface flow model. The SPOD mode shapes associated with (b) the first and (c) the second dominant frequency. The filament shapes in different phases of a mode are plotted in different colours. Filament shadow is given in grey. The dimensionless frequencies of whirling and beating states at
$b_1=741.1$
are 40.8 and 35.3, respectively.
Figures 15(b,c) show the two leading SPOD modes of QP1. At this actuation, which is close to the Hopf bifurcation where whirling becomes unstable, the leading SPOD mode (figure 15
b) resembles a whirling state, which can be considered as the primary motion of QP1 at
$b_1=741.1$
. Figure 15(c) shows the second mode, which is the leading mode associated with the second dominant frequency. We notice that this mode is planar and has a shorter wavelength than the eigenmode at the double Hopf bifurcation (figure 6
c). This suggests that the second SPOD mode corresponds to a secondary buckling instability, which triggers the transition from whirling to beating. As the surface flow magnitude increases, the compression along the filament exceeds a critical load, leading to buckling with a wavenumber higher than that observed in the double Hopf bifurcation. This is consistent with the Floquet analysis, which also yielded higher wavenumber modes (see figure 12).
3.4. Instabilities and dynamics at high actuation
Up to moderate actuation, we observe that all filament models exhibit the same dynamic states. As shown in §§ 3.1–3.3, each model undergoes similar bifurcations that ultimately lead to a stable beating state. At higher actuations, however, beating becomes unstable, and we encounter instabilities that yield different states for the different models.
The states encountered at high actuation for the different models are summarised in the bifurcation diagrams shown in figure 16. In each, the stable beating state is plotted in red. The leading eigenvalues from the stability analysis of the beating state are shown in figure 17. Each model undergoes distinct bifurcations, and reaches a different state at the highest actuation values. These states are shown as green curves in figure 16: a chaotic state for tip-driven models, a periodic bent-beating state for the distributed follower force model, and a quasi-periodic state QP2 for the distributed surface flow model. We now describe the states exhibited by each model.
Bifurcation diagrams for (a) tip-driven follower force, (b) distributed follower force, (c) tip-driven surface flow, and (d) distributed surface flow models at high actuation. Green full circles mark bifurcations. At the lower limit of each subplot, a stable beating state is plotted. The final state exhibited by each model (which is model-dependent) is plotted in green.
Leading eigenvalues of the beating state for (a) tip-driven follower force, (b) distributed follower force, (c) tip-driven surface flow, and (d) distributed surface flow models at high actuation. Real and imaginary parts of the leading eigenvalues are shown in red and blue, respectively. The bifurcations are marked with vertical dashed lines.
3.4.1. Transition to chaotic states for tip-driven models
For the tip-driven models, we observe a myriad of states and transitions that emerge at high actuation, consistent with previous studies on the tip-driven follower force model (Ling et al. Reference Ling, Guo and Kanso2018; Clarke et al. Reference Clarke, Hwang and Keaveny2024). In the discussion that follows, we focus exclusively on the tip-driven surface flow model, thoroughly mapping its instabilities and bifurcations, as this case has not been reported previously.
To begin, we explore the dynamics when filament motion is restricted to a plane. Here, we find three states connected by two saddle-node bifurcations, shown in figure 18(a). The three coexisting beating states plotted in figure 19 are beating, inter-beating and buckled-beating. Beating and buckled-beating are bistable, and the unstable inter-beating solution forms the branch that joins these two states.
Close views of the transition regime of the bifurcation diagram in figure 16(c) for the tip-driven surface flow model. Green full circles mark bifurcations. (a) In two dimensions, we observe a saddle-node transition between three beating states: beating, inter-beating and buckled-beating. (b) In addition to these, we observe two more bifurcations in the 3-D model where bent-beating and QP3 emerge. After this transition, the state spaces of 2-D and 3-D dynamics are the same. At a higher surface flow value, we observe a transition to a planar chaotic state. We plot these states in figures 19 and 20.
(a) Beating, (b) inter-beating and (c) buckled-beating states coexisting in the transition regime at
$b_1=2220$
, obtained using the tip-driven surface flow model. All three of the states are planar, hence exist in both two and three dimensions. The tip pathline of the filament is plotted in blue. The filament colour gets darker as time increases. We provide a video corresponding to this figure in supplementary movie 2.
In our discussion, we set time such that a beat cycle starts and ends (
$t=0$
and
$t=T$
) when the filament tip attains its maximum amplitude in the
$-x$
-direction; the filament tip is on the
$z$
-axis at
$t=T/4$
and
$t=3T/4$
, and attains its maximum amplitude in the
$x$
-direction at
$t=T/2$
. In the beating state, the filament remains curved throughout the beat cycle (see figure 19
a). In inter-beating and buckled-beating, the filament is nearly upright at
$t=T/4$
and
$t=3T/4$
. This is especially noticeable in buckled-beating (figure 19
c). At these instances, the nearly straight filament experiences compression that triggers dynamics reminiscent of a higher wavenumber buckling event. This is reflected in the downward motion of the filament tip. Nonlinear effects follow this buckling transient, and the filament continues along the rest of its beat path. Comparing figures 19(b) and 19(c), we observe that buckling occurs along a shorter section of the filament in inter-beating, and begins after the tip crosses
$x=0$
. Tracking the unstable inter-beating branch shows that the buckled section of the filament decreases in length with increasing actuation.
When we allow for fully-3-D filament motion, more complex bifurcations (figure 18 b) are observed. Evaluating the stability of these two-dimensional (2-D) states in three dimensions, we find that beating becomes unstable before its annihilation, and buckled-beating emerges as an unstable state. Furthermore, we discover two new states: bent-beating and QP3, as visualised in figures 20(a,b).
Filament states at high actuation: (a) bent-beating, (b) QP3, (c) chaotic beating. The plotted states are obtained using the tip-driven surface flow model at
$b_1=2220$
,
$b_1=2270$
and
$b_1=4020$
. The tip pathline of the filament is plotted in blue. The filament colour gets darker in the increasing direction of time. We provide a video description of this figure in supplementary movie 2.
First, we observe a supercritical pitchfork bifurcation at
$b_1=2106$
, and the creation of bent-beating. Bent-beating is similar to beating, but the tip of the filament curves out of the beating plane (see figure 20
a). We find that the leading unstable mode is orthogonal to the plane of the beating state, leading to 3-D dynamics. The bent-beating branch (orange in figure 18
b) is stable throughout its existence. By increasing the surface flow magnitude, we see that the out-of-plane tip deflection peaks at
$b_1=2172$
then declines, with bent-beating eventually merging back into beating at
$b_1=2242$
.
A second unstable mode of beating emerges as the surface flow magnitude increases and produces a supercritical pitchfork bifurcation at
$b_1=2234$
. As this occurs before the annihilation of bent-beating, beating becomes unstable to two modes, and continues to be unstable throughout. Then QP3 emerges from this bifurcation (grey in figure 18
b), and is destroyed at
$b_1=2456$
in yet another supercritical pitchfork bifurcation with the buckled-beating state. Plotted in figure 20(b), QP3, albeit quasi-periodic, is different from QP1 and QP2 (discussed in § 3.4.2): QP3 exhibits dynamics that combine out-of-plane tipping, as seen with bent-beating, with the periodic buckling encountered in the 2-D periodic states in figure 19. The buckling mechanism is weak around
$b_1=2234$
(near beating), and becomes stronger as the actuation increases towards
$b_1=2456$
(near the buckled-beating state). Tipping–buckling coupling causes sudden changes in the beating plane, seen in the top view in figure 20(b).
We note that the bifurcation patterns are simpler for the tip-driven follower force model. When beating becomes unstable at
$f=358$
, we observe only the bent-beating state, compared to the many states observed in the tip-driven surface flow model (see figures 16
a,c). The leading eigenvalue is real in figure 17(a) (as it is in the case of the tip-driven surface flow model), and due to the symmetry breaking in the out-of-plane direction, we observe a pitchfork bifurcation. We suspect that the greater variety of states exhibited by the surface flow model is related to the stress peak at the filament tip that arises in this model, and presumably produces shorter-wavelength buckling events.
The final state observed in tip-driven models at the highest values of actuation that we consider is a chaotic state (green in figures 16 a,c). This state arises as multiple buckled-beating eigenmodes become unstable over a small actuation range. It exhibits high-wavenumber instabilities near the filament tip (see supplementary movie 2) and exists in both two and three dimensions (figure 20 c).
3.4.2. Dynamics of distributed actuation models at high actuation
We find that when the actuation profile is uniform along the filament, the bifurcation patterns at high actuation are very different from those of tip-driven models. In both of the distributed actuation models, we find that a single attracting state emerges after beating becomes unstable (see figures 16 b,d). As in § 3.4.1, we start by classifying the dynamics for the surface flow model.
Figure 17(d) shows the leading eigenvalue (a complex conjugate pair) at high actuation, revealing a supercritical Hopf bifurcation for the distributed surface flow model. We identify a new quasi-periodic state, QP2. Figure 21 tracks QP2 along its branch as the actuation increases. Near
$b_1=634$
, QP2 exhibits in-plane beating combined with an out-of-plane motion, evident from the nearly planar shadow of the filament with some out-of-plane components (see figure 21
a). Since this instability is absent in two dimensions, we conclude that the out-of-plane instability is the main mechanism driving the creation of QP2.
The QP2 state with increasing actuation for the distributed surface flow model, at (a)
$b_1=709.9$
, (b)
$b_1=2366$
and (c)
$b_1=4732$
. The filament’s colour gets darker as time passes. Filament shadow is given in grey. We provide a video description of this figure in supplementary movie 3.
As actuation increases, the shadow of the filament changes from planar to circular, indicating a balance between the primary modes of QP2. The filament forms a rotating helix-like shape (figure 21 c). Presumably as a result of the surface flow, the filament is inclined to move towards the base and generate further compression, resulting in a more shortened shape overall (compare figures 21(a) and 21(c)).
A notable feature of QP2 is the handedness of the filament coils. Coils form at the base and advect towards the tip. Given the handedness of the helix-like shape and the overall rotation of the filament, it is inclined to move towards its base. This further increases compression at the base, and triggers buckling. This creates kinks that also travel towards the tip and lead to the reversal of the coiling direction. Due to reversal, the handedness directions above and below a kink are opposite to each other.
We track the angular position
$\theta$
of the filament tip extracted from its
$x$
- and
$y$
-coordinates, and from it compute the angular speed
$\dot {\theta }$
. The angular speed of the filament tip in figure 22(a) shows reversals linked to changes in coiling direction at the base. We mark four instances, and plot the filament shape at those instances together with the rotation directions in figures 22(b,c). The time stamps are in sequential order and are separated by five tip reversals. In the first instance (blue circle), the filament has a positive angular speed with left-handed coils without kinks. In the second instance, the angular speed of the tip is negative while the base rotates anticlockwise. The filament has left-handed coils near the base, and right-handed coils close to the free end, with a visible kink in between. The third instance returns to fully left-handed coils (see supplementary movie 3(c)). These transient reversals of handedness are unique to QP2 and the distributed surface flow model. Increasing actuation further reveals that QP2 is robust, remaining the only stable state that we find. We also track the unstable beating state, and find no additional instabilities at higher actuation.
(a) Angular speed of the filament tip about the
$z$
-axis, for the QP2 state at
$b_1=3310$
for the distributed squirming model. The filament shapes in the four instances marked are plotted from (b) top view and (c) side view. The rotation directions of the filament are indicated by black arrows.
Although both distributed actuation models exhibit a single stable state at high actuation, the characteristics of the transition are different. The distributed follower force model undergoes a supercritical pitchfork bifurcation (see figure 17 b), leading to bent-beating (see figure 20(a), discussed in § 3.4.1). Unlike in the tip-driven follower force model, bent-beating is robust to increasing actuation. These attracting states in both distributed models suggest that the uniform actuation profile stabilises the filament against high-wavenumber instabilities. This is, presumably, due to the peak stress being near the base rather than the free end.
4. Conclusions
In this paper, we carried out a bifurcation and stability analysis of various active filament models, discovering new states and identifying transition mechanisms. We considered a single active filament clamped at its base in an unbounded domain driven by two actuation models (follower force and surface flow) and two actuation profiles (tip-driven and distributed). We described the instability mechanisms that lead to new states and transitions.
Our computational analysis provides physical insight into active filament dynamics for these different actuation models. Examining the filament dynamics starting from zero actuation, we found that all models exhibited similar initial states and transitions. This indicates that the qualitative features of the dynamics and their bifurcations are independent of actuation type or profile, with the mechanism of instability in each case being buckling in the presence of rotational symmetry. The motion of the surrounding fluid in the surface flow case had a stabilising effect, delaying the primary and secondary instabilities. Distributed actuation yielded a stress profile that has its maximum value at the filament base, leading to unstable mode shapes that are narrower at the tip. Bistability between whirling and beating was observed in all models except the tip-driven surface flow model. The existence of QP1 is a generic feature; however, its stability is sensitive to model details. Using SPOD, we found that the transition regime involving whirling, beating and QP1 is linked to secondary buckling, with compression playing the primary role in the bifurcations.
A more surprising outcome of our study occurs at high actuation, where we find that the filament dynamics differs due to actuation profile. Tip-driven models exhibit different periodic or chaotic solutions, while distributed models exhibit states robust to increasing actuation. In particular, we found an interesting coiled state (QP2) in the distributed surface flow model that featured a rotating helix-like shape that changes handedness. Physically, we suspect that these differences in the states exhibited by the models are linked to the stress distributions that we explored in more detail at lower actuation. There, the distributed models had a stress profile that decreased towards the free end. Tip-driven models had the opposite profile, with the stress increasing towards the free end. We suspect that this profile triggers high-wavenumber instabilities at high actuation. We note that the torque supplied by the internal stresses after buckling (i.e.
$\hat {\boldsymbol{t}}\times \boldsymbol{\varLambda }$
in (2.2)) in tip-driven models is generally expected to be larger at the tip than that in distributed models. Additionally, given the clamped boundary condition at the base, the tangent vector rotation is also less restricted at the filament tip. Thus the larger internal stress observed around the tip in tip-driven models, coupled with the ease with which the filament tip can rotate, likely results in large filament deformations, ultimately producing chaotic beating dynamics.
The filament models that we investigated are biologically relevant for understanding microtubule-generated flows. Microtubules have bending rigidity approximately 22 pN
$\unicode{x03BC}$
m
$^2$
, and can exceed 40
$\unicode{x03BC}$
m in length (Gittes et al. Reference Gittes, Mickey, Nettleton and Howard1993). Their motion is driven by dynein motors exerting forces up to 6 pN, which can induce buckling and drive dynamic states (Shingyoji et al. Reference Shingyoji, Higuchi, Yoshimura, Katayama and Yanagida1998; Ling et al. Reference Ling, Guo and Kanso2018). These correspond to follower force magnitudes up to
$f=440$
, or equivalently surface flow magnitudes up to approximately
$b_1=1530$
in our dimensionless units. These values fall within the parameter range studied. While the distributed surface flow model most faithfully reflects motor protein activity on flexible filaments involved in cytoplasmic streaming, the states that we observe are exhibited by many different active biofilaments. One well-known example is ciliary beating (Gibbons Reference Gibbons1981), which can be found throughout the natural world, including in the human body in organs such as the fallopian tubes (Xiao et al. Reference Xiao, Coppeta, Rogers, Isenberg, Zhu, Olalekan and McKinnon2017), the lungs (Smith et al. Reference Smith, Gaffney and Blake2008b
) and the brain (Worthington & Cathcart Reference Worthington and Cathcart1966; Guldbrandsen et al. Reference Guldbrandsen, Vethe, Farag, Oveland, Garberg, Berle, Myhr, Opsahl, Barsnes and Berven2014; Zhang et al. Reference Zhang, Guo, Zou, Yang, Zhang, Ji, Shao, Sun and Wang2015). While cilia are typically found to undergo beating that is largely planar (though asymmetric), nodal cilia (Smith et al. Reference Smith, Blake and Gaffney2008a
) involved in embryonic development instead exhibit whirling. Another perhaps less well-known example is filopodia – actin-rich, filamentous tubular structures that cells use to probe their surroundings. Motion of these structures is driven by actin polymerisation, as well as activity of the motor protein myosin along the actin. Myosin activity in particular drives filopodia whirling, as well as coiling and buckling dynamics (Leijnse et al. Reference Leijnse2022) reminiscent of the QP2 state found with the distributed surface flow model at high actuation. While the details regarding the underlying forcing mechanisms for cilia and filopodia differ from the models explored here, the emergence of similar dynamic states points to some overall generic features of these systems that may provide a common avenue for instabilities to arise.
While our study has focused on the filament dynamics, it is important to note that the resulting flow fields for the follower force and surface flow models differ significantly, even in cases where the filament dynamics appear similar. For example, in the upright state, the follower force is balanced by internal stresses, resulting in no net fluid flow. In contrast, the surface flow model generates considerable fluid motion. After buckling, filament motion in the follower force model will produce a fluid flow; however, the flow generated by the surface flow model remains considerably stronger. As an example, figure 23 shows the period-averaged flow magnitude of the beating state at actuation values where beating annihilates:
$f=396.5$
and
$f=109.7$
for the follower force models, and
$b_1=3691$
and
$b_1=634.1$
for the surface flow models. For tip-driven models, these states correspond to the onset of chaotic dynamics. All active filament models exhibit a
$1/r$
decay in the far field, with the distributed surface flow model producing the strongest flow.
The period-averaged flow magnitude of the beating state for different active filament models. The chosen actuation values are
$f=396.5$
and
$f=109.7$
for follower force models, and
$b_1=3691$
and
$b_1=634.1$
for surface flow models, which are at the bifurcation point (marked in figure 17) leading to the annihilation of the beating state. Regardless of the model, the flow decays as
$1/r$
.
An important scenario that we have not explored here is the filament being attached to a no-slip surface. We can perform a comparison for the tip-driven follower force model between the results from this study and those from Clarke et al. (Reference Clarke, Hwang and Keaveny2024), where a no-slip surface is included. In doing so, we see that at low actuation, the initial double Hopf bifurcation is largely unchanged. There is a modification of the transition between whirling and beating. While the no-slip surface allows for a stable QP1 branch, and without the surface it is unstable, in both cases the QP1 branch exists for only a very small range of follower force values. At higher forcing when beating becomes unstable, we find that the series of bifurcations and states reported here (critical value of the forcing and emergence of bent-beating followed by chaotic behaviour) corresponds very closely to those exhibited by filaments with a much higher aspect ratio (
$L/a = 220$
) when the no-slip surface is included. In Clarke et al. (Reference Clarke, Hwang and Keaveny2024), the aspect ratio was increased by increasing the length of the filament while keeping the cross-sectional radius fixed. This results in the forced tip being located further away from the surface, reducing its hydrodynamic effect. For lower-aspect-ratio filaments, the bent-beating state did not appear after beating became unstable, therefore we see that the nearby surface has an effect on the nature of this transition. For the surface flow model, especially for the distributed case, we expect the no-slip surface to have a stronger effect as it would alter substantially the flows generated by the internal stresses, even at low actuation before the filament buckles. This could, for example, reduce the stabilising effect of entrainment that produced a delay in the initial buckling. We expect the study of the hydrodynamic effects of a no-slip surface on the bifurcations to form an important aspect of our future research.
Finally, although the differences in the flow fields for the different models largely do not affect the qualitative dynamics and bifurcations at low actuation, this may not be the case when multiple filaments interact. Indeed, when multiple filaments are present and fluid entrainment is incorporated as in Dutta et al. (Reference Dutta, Farhadifar, Lu, Kabacaoğlu, Blackwell, Stein, Lakonishok, Gelfand, Shvartsman and Shelley2024), for example, hydrodynamic interactions are able to alter the bifurcations that are observed, and most notably change the initial double Hopf bifurcation to a pitchfork bifurcation where the filaments are found to bend in the same direction (Stein et al. Reference Stein, De, G., E., M. and Goldstein2021; Dutta et al. Reference Dutta, Farhadifar, Lu, Kabacaoğlu, Blackwell, Stein, Lakonishok, Gelfand, Shvartsman and Shelley2024). This is not observed with the follower force model where there is no entrainment and only coordinated beating or whirling has been seen to emerge (Westwood & Keaveny Reference Westwood and Keaveny2021; Su & Keaveny Reference Su and Keaveny2024). When the double Hopf bifurcation is suppressed, the fluid–structure interactions of fluid entraining filaments have been shown to contribute to the emergence of large-scale streaming flows that could aid transport within cells (Dutta et al. Reference Dutta, Farhadifar, Lu, Kabacaoğlu, Blackwell, Stein, Lakonishok, Gelfand, Shvartsman and Shelley2024). Our future work will focus on examining these and related bifurcations that occur in biofilament clusters, as well as bifurcations and instabilities in models driven by active shearing (Oriola, Gadêlha & Casademunt Reference Oriola, Gadêlha and Casademunt2017) and incorporating structural features such as anisotropy (Clarke, Hwang & Keaveny Reference Clarke, Hwang and Keaveny2025).
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2026.11259.
Funding
I.R.O. gratefully acknowledges funding by an Imperial College Department of Mathematics Roth Scholarship. B.C. gratefully acknowledges funding from an EPSRC scholarship (grant no. EP/W523872/1).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Effect of twist
In this work, we have set the the twist modulus equal to the bending modulus such that
$K_T/K_B = 1$
. Using values for the geometry and material parameters of microtubules from the literature, Clarke et al. (Reference Clarke, Hwang and Keaveny2024) estimated the range of this ratio to span several orders of magnitude,
$0.002 \lt K_T/K_B \lt 350$
. Based on this, Clarke et al. (Reference Clarke, Hwang and Keaveny2024) performed simulations for filaments driven by a tip-driven follower force with
$K_T/K_B = 0.01, 1, 100$
. Filament dynamics and the resulting states remained qualitatively similar, except for the case
$K_T/K_B = 0.01$
, where some differences in the whirling-to-beating transition were observed. Assuming this trend to also be exhibited by the other models that we explore in this paper, we expect the dynamics to be largely unaffected by choice of
$K_T/K_B$
provided that
$K_T/K_B \geqslant 1$
, which corresponds to a large portion of the relevant range.
Appendix B. Derivation of the slip velocity correction
B.1. Single force-free squirmer
The squirmer is a spherical swimming body that imposes instantaneous velocity boundary conditions to the surrounding fluid due to effective surface deformations (Lighthill Reference Lighthill1952; Blake Reference Blake1971). A force-free, neutral squirmer swimming in the
$\hat {\boldsymbol{e}}$
-direction applies the force quadruple (Ishikawa, Simmonds & Pedley Reference Ishikawa, Simmonds and Pedley2006; Delmotte et al. Reference Delmotte, Keaveny, Plouraboué and Climent2015b
)
$\boldsymbol{H}=-(4/3)\pi \eta a^3B_1\hat {\boldsymbol{e}}$
on the fluid, and generates a Stokes flow that satisfies
\begin{align} \boldsymbol{0}& =- \boldsymbol{\nabla }p+\eta\, {\nabla} ^2\boldsymbol{u}+\boldsymbol{H}\, {\nabla} ^2\boldsymbol{\delta }\left (\boldsymbol{r}\right ), \end{align}
\begin{align} \boldsymbol{0}& = \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}, \end{align}
where
$\boldsymbol{r}$
is the position relative to the squirmer centre,
$r=\|\boldsymbol{r}\|$
, and
$\hat {\boldsymbol{r}}=\boldsymbol{r}/r$
. The solution to this system reads
\begin{equation} \boldsymbol{u}(\boldsymbol{r})=\dfrac {1}{4\pi \eta r^3}\left (\mathbb{I}-3\,\hat {\boldsymbol{r}}\hat {\boldsymbol{r}}\right )\boldsymbol{H}. \end{equation}
We obtain the squirmer-induced velocity
$\boldsymbol{u}^\prime$
of a nearby particle with hydrodynamic radius
$a$
using Faxén’s relations (Kim & Karrila Reference Kim and Karrila2013):
\begin{equation} \boldsymbol{u}^\prime =\left (\mathbb{I}+\dfrac {a^2}{6} {\nabla} ^2\right )\boldsymbol{u}(\boldsymbol{r}). \end{equation}
Since in this case the contribution of the Laplacian term is zero, the induced velocity is simply
$\boldsymbol{u}(\boldsymbol{r})$
and given by (B4). On the other hand, the self-induced velocity of a squirmer is
$(2/3)B_1\hat {\boldsymbol{e}}$
.
B.2. Hydrodynamically coupled squirmers
For an
$N$
-squirmer system, we have
\begin{equation} \boldsymbol{0}=- \boldsymbol{\nabla }p+\eta\, {\nabla} ^2\boldsymbol{u}+\sum \limits _i\boldsymbol{H}_i\, {\nabla} ^2\boldsymbol{\delta }\left (\boldsymbol{x}-\boldsymbol{Y}_i\right ), \end{equation}
where
$\boldsymbol{H}_i=(4/3)\pi \eta a^3B_1\hat {\boldsymbol{t}}_i$
. Here, we choose the swimming direction
$\hat {\boldsymbol{e}} =-\hat {\boldsymbol{t}}$
such that each squirmer tries to move towards the base of the filament. Then the flow induced on segment
$j$
by segment
$i$
(
$i\neq j$
) is
\begin{equation} \boldsymbol{u}^\prime \left (\boldsymbol{r}_{ij}=\boldsymbol{Y}_i-\boldsymbol{Y}_j\right )=\dfrac {1}{4\pi \eta r_{ij}^3}\left (\mathbb{I}-3\,\hat {\boldsymbol{r}}_{ij}\hat {\boldsymbol{r}}_{ij}\right )\boldsymbol{H}_j=\boldsymbol{S}_{ij}\boldsymbol{H}_j \end{equation}
for a squirming interaction matrix
$\boldsymbol{S}_{ij}\in \mathbb{R}^{3\times 3}$
. We remark that
$i,j$
are not matrix indices, but denote segment indices. Similarly, self-interactions can be expressed as
\begin{equation} \boldsymbol{u}^\prime \left (\boldsymbol{r}_{ii}\right )=-\dfrac {1}{2\pi \eta a^3}\mathbb{I}\,\boldsymbol{H}_i=\boldsymbol{S}_{ii}\boldsymbol{H}_i. \end{equation}
Using (B6) and (B7), we can group all
$\boldsymbol{H}_i$
under a general swimming vector
$\boldsymbol{H}= [\boldsymbol{H}_1^{\textrm{T}},\ldots ,\boldsymbol{H}_N^{\textrm{T}} ]^{\textrm{T}}\in \mathbb{R}^{3N\times 1}$
and a general swimming matrix
$\boldsymbol{\mathcal{S}}\in \mathbb{R}^{3N\times 3N}$
with block entries
$\boldsymbol{S}_{ij}$
. Then the translational and angular velocities due to the segment surface flows can be expressed as
\begin{equation} \begin{bmatrix} \boldsymbol{V}^\prime \\ \boldsymbol{\varOmega }^\prime \end{bmatrix}=\begin{bmatrix} \boldsymbol{\mathcal{S}}&\boldsymbol{0}\\\boldsymbol{0}&\boldsymbol{0} \end{bmatrix}\begin{bmatrix} \boldsymbol{H} \\\boldsymbol{0} \end{bmatrix}. \end{equation}
Appendix C. The effective Lie algebra element
It is sufficient to track the tangent angle along the filament to describe the 2-D filament dynamics. However, in three dimensions, a different approach is required to describe the evolution of a local frame at each point along the filament centreline. In our numerical framework, following Clarke et al. (Reference Clarke, Hwang and Keaveny2024), we use unit quaternions to describe local rotations of the coordinate frame, and effective Lie algebra elements for the state space representation of the dynamics.
Given the equilibrium quaternion, which describes the upright filament, i.e.
\begin{equation} \boldsymbol{q}_{\textit{eq}}=\begin{bmatrix} \dfrac {1}{\sqrt {2}},0,-\dfrac {1}{\sqrt {2}},0 \end{bmatrix},\quad \text{with}\quad \boldsymbol{q}^*_{\textit{eq}}=\begin{bmatrix} \dfrac {1}{\sqrt {2}},0,\dfrac {1}{\sqrt {2}},0 \end{bmatrix}, \end{equation}
the effective Lie algebra element associated with a quaternion
$\boldsymbol{q}$
is defined as
\begin{equation} {\boldsymbol{U}} = 2\arccos {\left (p_0\right )}\,\hat {\boldsymbol{p}}, \end{equation}
where
$\boldsymbol{q}\bullet \boldsymbol{q}^*_{\textit{eq}}=\begin{bmatrix} p_0,\overline {\boldsymbol{p}} \end{bmatrix}$
and
$\hat {\boldsymbol{p}}= {\overline {\boldsymbol{p}}}/{\|\overline {\boldsymbol{p}}\|}$
.
C.1. Connection to the position vector
Given an effective Lie algebra element
${\boldsymbol{U}}(s_k)$
associated with the
$k$
th filament segment, the corresponding quaternion is given by
\begin{equation} \boldsymbol{q}_k=\exp {\left ({\boldsymbol{U}}(s_k)\right )}\bullet \boldsymbol{q}_{\textit{eq}}, \end{equation}
where
$\exp { ({\boldsymbol{U}} )}=\begin{bmatrix} \cos ({ {\|{\boldsymbol{U}}\|}/{2}}),\sin ( {{\|{\boldsymbol{U}}\|}/{2}})({{\boldsymbol{U}}}/{\|{\boldsymbol{U}}\|}) \end{bmatrix}$
. Then the centre position of the
$k$
th segment is found recursively as
\begin{equation} \boldsymbol{X}_k=\boldsymbol{X}_{k-1}+\dfrac {\Delta L}{2}\left (\boldsymbol{R}\left (\boldsymbol{q}_k\right ) \hat {\boldsymbol{e}}_1+\boldsymbol{R}\left (\boldsymbol{q}_{k-1}\right ) \hat {\boldsymbol{e}}_1\right ), \end{equation}
where the action of the rotation matrix on an arbitrary vector
$\boldsymbol{v}$
is defined as
\begin{equation} \begin{bmatrix} 0,\boldsymbol{R}\left (\boldsymbol{q}\right )\boldsymbol{v} \end{bmatrix}=\boldsymbol{q}\bullet \begin{bmatrix} 0,\boldsymbol{v} \end{bmatrix}\bullet \boldsymbol{q}^*. \end{equation}
Then the position of a segment centre is
\begin{equation} \begin{aligned} {\boldsymbol{X}}_k&={\boldsymbol{X}}_1+\dfrac {\Delta L}{2}\left (\boldsymbol{R}\left (\boldsymbol{q}_1\right )+\boldsymbol{R}\left (\boldsymbol{q}_k\right )+2\sum _2^{k-1}\boldsymbol{R}\left (\boldsymbol{q}_i\right )\right )\hat {\boldsymbol{e}}_1\\ &={\boldsymbol{X}}_1+\dfrac {\Delta L}{2}\textbf {vec}\left (\boldsymbol{q}_1\bullet \begin{bmatrix}0,\hat {\boldsymbol{e}}_1\end{bmatrix}\bullet \boldsymbol{q}^*_1+\boldsymbol{q}_k\bullet \begin{bmatrix}0,\hat {\boldsymbol{e}}_1\end{bmatrix}\bullet \boldsymbol{q}^*_k+2\sum _2^{k-1}\boldsymbol{q}_i\bullet \begin{bmatrix}0,\hat {\boldsymbol{e}}_1\end{bmatrix}\bullet \boldsymbol{q}^*_i\right ), \end{aligned} \end{equation}
where the
$\textbf {vec}$
operator takes the vector part of the quaternion. Notice that each quaternion multiplication term corresponds to the tangent vector of the respective filament segment. We investigate some special cases of the effective Lie algebra element to show its connection to the tangent vectors, and hence the position and geometry of the filament.
C.2. Special cases of effective Lie algebra element
C.2.1. Case I
If the effective Lie algebra element of a segment is of the form
\begin{equation} {\boldsymbol{U}}(s,t)=\begin{bmatrix}U_1(s,t),0,0\end{bmatrix}, \end{equation}
then the associated tangent vector of a segment is
\begin{equation} \begin{aligned} \hat {\boldsymbol{t}}(s,t)&=0\,\hat {\boldsymbol{e}}_1-\sin \|{\boldsymbol{U}}(s,t)\|\,\hat {\boldsymbol{e}}_2+\cos \|{\boldsymbol{U}}(s,t)\|\,\hat {\boldsymbol{e}}_3\\ &=0\,\hat {\boldsymbol{e}}_1-\sin \lvert U_1(s,t)\rvert \,\hat {\boldsymbol{e}}_2+\cos \lvert U_1(s,t)\rvert \,\hat {\boldsymbol{e}}_3. \end{aligned} \end{equation}
This form of effective Lie algebra element relates to planar motion in the
$yz$
-plane.
C.2.2. Case II
If the effective Lie algebra element of a segment is of the form
\begin{equation} {\boldsymbol{U}}(s,t)=\begin{bmatrix}0,U_2(s,t),0\end{bmatrix}, \end{equation}
then the associated tangent vector of a segment is
\begin{equation} \begin{aligned} \hat {\boldsymbol{t}}(s,t)&=\sin \|{\boldsymbol{U}}(s,t)\|\,\hat {\boldsymbol{e}}_1+0\,\hat {\boldsymbol{e}}_2+\cos \|{\boldsymbol{U}}(s,t)\|\,\hat {\boldsymbol{e}}_3\\ &=\sin \lvert U_2(s,t)\rvert \,\hat {\boldsymbol{e}}_1+0\,\hat {\boldsymbol{e}}_2+\cos \lvert U_2(s,t)\rvert \,\hat {\boldsymbol{e}}_3. \end{aligned} \end{equation}
This form of effective Lie algebra element relates to planar motion in the
$xz$
-plane.
C.2.3. Case III
A natural generalisation of the previous cases is their sum, i.e. an effective Lie algebra element of a segment is of the form
\begin{equation} {\boldsymbol{U}}(s,t)=\begin{bmatrix}U_1(s,t),U_2(s,t),0\end{bmatrix}. \end{equation}
The associated tangent vector of a segment is
\begin{equation} \hat {\boldsymbol{t}}(s,t)=\dfrac {U_2(s,t)}{\|{\boldsymbol{U}}(s,t)\|}\sin \|{\boldsymbol{U}}(s,t)\|\,\hat {\boldsymbol{e}}_1-\dfrac {U_1(s,t)}{\|{\boldsymbol{U}}(s,t)\|}\sin \|{\boldsymbol{U}}(s,t)\|\,\hat {\boldsymbol{e}}_2+\cos \|{\boldsymbol{U}}(s,t)\|\,\hat {\boldsymbol{e}}_3. \end{equation}
In the case where
$\|{\boldsymbol{U}}(s,t)\|$
is time-independent, such that
$U_1^2(s,t)+U_2^2(s,t)=c(s)$
for a function
$c(s)$
, the tangent vector follows a circular path restricted to a
$z=\text{constant}$
plane. This is a model for a whirling-like motion of the filament.
On the other hand, if the norm of the Lie algebra element changes but
${U_2}/{U_1}=\text{constant}$
(or
${U_1}/{U_2}=\text{constant}$
, whichever is well-defined), then the motion is restricted to a plane. The angle that the plane makes with the
$x$
-axis is
$\theta =-\arctan ({U_1}/{U_2})$
.
C.2.4. Case IV
The effective Lie algebra element in its most general form is
\begin{equation} {\boldsymbol{U}}(s,t)=\begin{bmatrix}U_1(s,t),U_2(s,t),U_3(s,t)\end{bmatrix}. \end{equation}
The resulting tangent vector is
\begin{equation} \begin{aligned} \hat {\boldsymbol{t}}(s,t)=&\left [\dfrac {U_1}{\|{\boldsymbol{U}}\|}\dfrac {U_3\left (1-\cos \|{\boldsymbol{U}}\|\right )}{\|{\boldsymbol{U}}\|}+\dfrac {U_2}{\|{\boldsymbol{U}}\|}\sin \|{\boldsymbol{U}}\|\right ]\,\hat {\boldsymbol{e}}_1\\&+\left [\dfrac {U_2}{\|{\boldsymbol{U}}\|}\dfrac {U_3\left (1-\cos \|{\boldsymbol{U}}\|\right )}{\|{\boldsymbol{U}}\|}-\dfrac {U_1}{\|{\boldsymbol{U}}\|}\sin \|{\boldsymbol{U}}\|\right ]\,\hat {\boldsymbol{e}}_2\\&+\left [\dfrac {U_3}{\|{\boldsymbol{U}}\|}\dfrac {U_3\left (1-\cos \|{\boldsymbol{U}}\|\right )}{\|{\boldsymbol{U}}\|}+\cos \|{\boldsymbol{U}}\|\right ]\,\hat {\boldsymbol{e}}_3. \end{aligned} \end{equation}
C.3. Connection to periodic states
As shown in § C.2, certain types of the effective Lie algebra element induce certain restrictions on the filament dynamics. Here, we generalise these results as two examples of filament states.
C.3.1. Decomposition of beating state
Beating is a planar state, which can be described as a time-periodic bending of the filament. Since it is essentially a 2-D state, a segment’s tangent vector can be expressed by only one time-dependent variable:
\begin{equation} \boldsymbol{\hat {{t}}}(s,t)=\sin {\left (\phi (s,t)\right )}\cos {\theta }\,\hat {\boldsymbol{e}}_1- \sin {\left (\phi (s,t)\right )}\sin {\theta }\,\hat {\boldsymbol{e}}_2+\cos {\left (\phi (s,t)\right )}\,\hat {\boldsymbol{e}}_3, \end{equation}
where
$\theta$
determines the beating plane, and
$\phi (s,t)$
is the beat function. From § C.2.3, we take
$U_1(s,t)=\phi (s,t)\sin \theta$
and
$U_2(s,t)=\phi (s,t)\cos \theta$
. Then the associated effective Lie algebra element is
\begin{equation} \begin{aligned} {\boldsymbol{U}}(s,t)&=\begin{bmatrix} \phi (s,t)\sin {\theta },\phi (s,t)\cos {\theta },0 \end{bmatrix}\\ &=\sin {\theta }\begin{bmatrix} \phi (s,t),0,0 \end{bmatrix}+\cos {\theta }\begin{bmatrix} 0,\phi (s,t),0 \end{bmatrix}. \end{aligned} \end{equation}
Notice that we restrict neither the filament shape nor its time-dependency. Here, we only assume (owing to periodicity) that
$\phi (s,t)=\phi (s,t+T)$
, where
$T$
is the beating period. This decomposition also demonstrates that the summation of two beating signals (with unit modulus) corresponds to rotations about the
$z$
-axis.
C.3.2. Decomposition of whirling state
The tangent vector of a whirling filament can be described without loss of generality as
\begin{equation} \hat {\boldsymbol{t}}(s,t)=A\,{\rm e}^{{\rm i}\omega t}\sin \psi (s)\,\hat {\boldsymbol{e}}_1-A\,{\rm e}^{{\rm i}\omega t}\sin \psi (s)\,\hat {\boldsymbol{e}}_2+\cos \psi (s)\,\hat {\boldsymbol{e}}_3+\text{c.c.} \end{equation}
for some arbitrary
$A\in \mathbb{C}$
with
$\lvert A\rvert =1$
. Here,
$\tan \psi (s)$
is the time-independent angle that the tangent vector makes from the meridional point of view. From § C.2.3, we take
$U_1^2(s,t)^2+U_2^2(s,t)^2=\psi ^2(s)$
, and obtain
$U_1(s,t)={\rm i}A\,{\rm e}^{{\rm i}\omega t}\,\psi (s)$
and
$U_2(s,t)=A\,{\rm e}^{{\rm i}\omega t}\,\psi (s)$
. Then any whirling Lie algebra element can be decomposed as
\begin{equation} \begin{aligned} {\boldsymbol{U}}(s,t)&=\begin{bmatrix} {\rm i}A\,\psi (s)\,{\rm e}^{{\rm i}\omega t},A\,\psi (s)\,{\rm e}^{{\rm i}\omega t},0 \end{bmatrix}+\text{c.c.}\\ &=\left ({\rm i}A\begin{bmatrix} \psi (s),0,0 \end{bmatrix}+A\begin{bmatrix} 0,\psi (s),0 \end{bmatrix}\right ) {\rm e}^{{\rm i}\omega t}+\text{c.c.} \end{aligned} \end{equation}
Each of these effective Lie algebra elements describes planar filament dynamics – one on the
$yz$
-plane, and another on the
$xz$
-plane. They are identical in shape, up to a
$\pi /2$
rotation about the
$z$
-axis. The periodic time dependence is separated from space; therefore, the filament shape is, in a sense, ‘frozen in time’, and whirling is a solid body rotation.
However, another interpretation of this result is that these effective Lie algebra elements individually describe a beating filament. The beating has the frequency
$\omega /(2\pi )$
, and the filament has a specific shape described by
$\phi (s)$
. The coefficients associated with each signal have the same magnitude, but are orthogonal to each other in the complex plane. Therefore, a
$\pi /2$
-delayed combination of two beating signals creates whirling. Furthermore, the sum of any two beating signals is not necessarily a whirling state, in that whirling requires time and space to be separated in the beating function.
C.4. Use in the bisection method
Given
${\boldsymbol{U}}_{{w}}(s,t)=[U_{1,{w}}(s,t),U_{2,{w}}(s,t),0]$
for whirling, we have
\begin{equation} U_{1,{w}}^2(s)+U_{2,{w}}^2(s)=c(s) \end{equation}
for some
$c(s)\in \mathbb{R}$
. Similarly, given
${\boldsymbol{U}}_{{b}}(s,t)=[U_{1,{b}}(s,t),U_{2,{b}}(s,t),0]$
for beating, we have
\begin{equation} \dfrac {U_{1,{b}}(s)}{U_{2,{b}}(s)}=c \end{equation}
for some constant
$c\in \mathbb{R}$
. After a sufficiently long time, a solution satisfies either (C19) or (C20), which we use as convergence criteria.
Open access
