2.1. Filament model

The filament has length $L$ and cross-sectional radius $a$ , and is clamped at its base to a no-slip planar surface. The filament is forced at its free end by a follower force, a compressive tangential force, as shown in figure 1(a). The filament centreline, $\boldsymbol{Y}(s,t),$ is parameterised by arc-length $s\in [0,L]$ and evolves over time $t\in [0,\infty )$ . Filament deformation is described by the orthonormal local basis $\{\boldsymbol{\hat {t}}(s,t),\boldsymbol{\hat {\mu }}(s,t),\boldsymbol{\hat {\nu }}(s,t)\}$ . For the filament dynamics, we use the framework outlined in Schöller et al. (Reference Schöller, Townsend, Westwood and Keaveny2021). Due to the small length scales and velocity scales associated with ciliary motion, the Reynolds number is small and so we neglect the inertia of the fluid and the filament. Consequently, force and moment balances along the filament length yield

(2.1) \begin{align} \frac {\partial \boldsymbol{\varLambda }}{\partial s}+\boldsymbol{f}^H &=\boldsymbol{0}, \end{align}

(2.2) \begin{align} \frac {\partial \boldsymbol{M}}{\partial s}+\boldsymbol{{\hat {t}}} \times \boldsymbol{\varLambda }+\boldsymbol{\tau }^H &=\boldsymbol{0}, \end{align}

where $\boldsymbol{f}^H(s,t)$ and $\boldsymbol{\tau }^H(s,t)$ are the forces and torques per unit length, respectively, exerted by the surrounding fluid and $\boldsymbol{\varLambda }(s,t)$ and $\boldsymbol{M}(s,t)$ are the internal force and moment, respectively, acting on the filament cross-section (see figure 1 b). The internal forces act as Lagrange multipliers enforcing the kinematic constraint:

(2.3) \begin{equation} \frac {\partial \boldsymbol{Y}}{\partial s} =\boldsymbol{\hat {t}}, \end{equation}

and the internal moments are given by the constitutive relation (Schöller et al. Reference Schöller, Townsend, Westwood and Keaveny2021):

(2.4) \begin{equation} \boldsymbol{M}(s,t) = K_B \left (\hat {\boldsymbol{t}} \boldsymbol{\cdot }\frac {\partial \hat {\boldsymbol{\nu }}}{\partial s} - \tilde {\kappa }_0\right ) \hat {\boldsymbol{\mu }} + \beta K_B\left (\hat {\boldsymbol{\mu }} \boldsymbol{\cdot }\frac {\partial \hat {\boldsymbol{t}}}{\partial s} \right ) \hat {\boldsymbol{\nu }}+ K_T \left ( \hat {\boldsymbol{\nu }} \boldsymbol{\cdot }\frac {\partial \hat {\boldsymbol{\mu }}}{\partial s}\right ) \hat {\boldsymbol{t}}, \end{equation}

where $\tilde {\kappa }_0 = \kappa _0/L$ defines the intrinsic curvature, such that the rest configuration of the filament is a planar arc with constant curvature $\kappa_{0}$ , as sketched in figure 1(c). Therefore, larger values of $\kappa _0$ correspond to higher-curvature rest configurations for the filament. Also defined in (2.4) are $K_B$ , the filament bending stiffness in the direction of intrinsic curvature, $\beta K_B$ , the bending stiffness in the perpendicular plane, and $K_T$ , the twisting stiffness. Choosing $\beta \gt 1$ leads to preferred bending about the $\hat{\boldsymbol{\mu}}$ -axis and $\beta \lt 1$ leads to preferred bending about the $\hat {\boldsymbol{\nu }}$ -axis. Correspondingly, setting $\beta \rightarrow \infty$ should restrict the filament dynamics to the $(\hat {\boldsymbol{\nu }},\hat {\boldsymbol{t}})$ plane entirely. Choosing $\beta = 1, \kappa _0 = 0$ corresponds to the isotropic, symmetric case considered in Clarke et al. (Reference Clarke, Hwang and Keaveny2024).

The base of the filament $(s=0)$ is clamped to a no-slip planar surface, fixing the position and local orientation of the base to be $\boldsymbol{Y}(0,t) = \boldsymbol{0}$ and $(\hat {\boldsymbol{t}}(0,t),\hat {\boldsymbol{\mu }}(0,t),\hat {\boldsymbol{\nu }}(0,t)) = (\hat {\boldsymbol{z}},\hat {\boldsymbol{y}},\hat {\boldsymbol{x}})$ . The distal end $(s=L)$ is moment-free, $\boldsymbol{M}(L,t)=\boldsymbol{0}$ , and is driven by the follower force:

(2.5) \begin{equation} \boldsymbol{\varLambda }(L,t) = -\frac {f K_B}{L^2} \hat {\boldsymbol{t}}(L,t), \end{equation}

where $f$ is the non-dimensional follower force strength. Through this non-dimensionalisation, the length is incorporated in the non-dimensional parameter, $f$ . The impact of filament aspect ratio, $L/a$ , in the symmetric and isotropic case ( $\kappa _0 = 0$ and $\beta = 1$ ) was investigated in Clarke et al. (Reference Clarke, Hwang and Keaveny2024). We note that an alternative definition of the non-dimensional follower force can be reached by balancing with the elastic force in the $\hat {\boldsymbol{\nu }}$ direction. This would lead to a rescaled non-dimensional follower force, $\tilde {f} = \beta f$ .

We discretise the filament into $N$ segments of length $\Delta L$ , such that segment $i$ has centre $\boldsymbol{Y}_i$ and an orientation defined by the local frame, $\{\boldsymbol{\hat {t}}_i,\boldsymbol{\hat {\mu }}_i,\boldsymbol{\hat {\nu }}_i\}$ for $i=1,\ldots ,N$ . We discretise (2.1), (2.2) and (2.3), using central differencing and rearrange to obtain

(2.6) \begin{align} \boldsymbol{F}_i^{C}+\boldsymbol{F}_i^H &=\boldsymbol{0}, \end{align}

(2.7) \begin{align} \boldsymbol{T}_i^{E}+\boldsymbol{T}_i^{C}+\boldsymbol{T}_i^H &=\boldsymbol{0} , \end{align}

for segment $i=1,\ldots ,N$ , and

(2.8) \begin{equation} \boldsymbol{Y}_{i+1}-\boldsymbol{Y}_i-\frac {\Delta L}{2}\left (\boldsymbol{\hat {t}}_i+\boldsymbol{\hat {t}}_{i+1}\right ) =\boldsymbol{0}, \end{equation}

for $i=1,\ldots ,N-1$ . Here $\boldsymbol{T}_i^{E} = \boldsymbol{M}_{i+1/2} - \boldsymbol{M}_{i-1/2}$ is the elastic torque, $\boldsymbol{F}_i^{C} = \boldsymbol{\varLambda }_{i+1/2} - \boldsymbol{\varLambda }_{i-1/2}$ and $\boldsymbol{T}_i^C =(\Delta L/2) \hat {\boldsymbol{t}}_i \times (\boldsymbol{\varLambda }_{i+1/2} + \boldsymbol{\varLambda }_{i-1/2})$ are the constraint forces and torques, respectively, and $\boldsymbol{F}_i^H = \Delta L \boldsymbol{f}_i^H$ and $\boldsymbol{T}_i^H = \Delta L \boldsymbol{\tau }_i^H$ are the hydrodynamic force and torque, respectively. The internal forces, $\boldsymbol{\varLambda }_{i+1/2},$ act as Lagrange multipliers enforcing the constraint (2.8). The internal moments, $\boldsymbol{M}_{i+1/2}$ , are defined by a discretised version of the constitutive law given in (2.4).

As we neglect fluid inertia due to the low Reynolds number, the fluid velocity is governed by the Stokes equations. Therefore, the translational and rotational velocities of the segments are linearly related to the hydrodynamic forces and torques on the segments through the mobility matrix, $\mathcal{M}$ . Explicitly, this can be expressed via the mobility relation:

(2.9) \begin{equation} \left (\begin{array}{l} \boldsymbol{V} \\ \boldsymbol{\varOmega }\end{array}\right )= \mathcal{M} \left (\begin{array}{l}-\boldsymbol{F}^{H} \\ -\boldsymbol{T}^{H}\end{array}\right ), \end{equation}

where $\mathcal{M}$ is the $6N \times 6N$ mobility matrix, $\boldsymbol{V}$ and $\boldsymbol{\varOmega }$ are $6N \times 1$ vectors containing the translational and angular velocities of all segments, respectively, and $\boldsymbol{F}^H$ and $\boldsymbol{T}^H$ are $6N\times 1$ vectors containing the hydrodynamic forces and torques of all segments, respectively, which are defined through the force and moment balances (2.6) and (2.7). In this work, we use the Rotne–Prager–Yamakawa (RPY) tensor (Wajnryb et al. Reference Wajnryb, Mizerski, Zuk and Szymczak2013), adapted to include the no-slip condition at $z=0$ (Swan & Brady Reference Swan and Brady2007), to define the mobility matrix, $\mathcal{M}$ . Here, each segment is modelled as a particle with hydrodynamic radius equal to the cross-sectional radius of the filament, $a$ . Various hydrodynamic models have been used in this context. The follower force model appears to give rise to the same dynamics in the forcing ranges we consider, regardless of the details of the hydrodynamic model (Ling et al. Reference Ling, Guo and Kanso2018). For a comparison between RPY and non-local slender body theory, we direct the reader to Maxian, Mogilner & Donev (Reference Maxian, Mogilner and Donev2021).

After solving for the translational and rotational velocities of each segment, we integrate the system forwards in time. To describe the evolution of segment orientations, we utilise unit quaternions, $\boldsymbol{q} = (q_0,q_1,q_2,q_3)$ with $||\boldsymbol{q}||=1$ (Schöller et al. Reference Schöller, Townsend, Westwood and Keaveny2021; Hanson Reference Hanson2006). Segment $i$ is defined to have a corresponding unit quaternion, $\boldsymbol{q}_i$ , which maps the standard basis, $(\hat {\boldsymbol{x}},\hat {\boldsymbol{y}},\hat {\boldsymbol{z}})$ , to the segment’s local basis at time $t,$ $(\hat {\boldsymbol{t}}_i(t)\hat {\boldsymbol{\mu }}_i(t)\hat {\boldsymbol{\nu }}_i(t))$ , through a rotation matrix, $\boldsymbol{R}(\boldsymbol{q})$ (see Schöller et al. (Reference Schöller, Townsend, Westwood and Keaveny2021) for further details). The evolution of each quaternion is described by

(2.10) \begin{equation} \frac {{\textrm d}\boldsymbol{q}_i}{{\textrm d} t} = \frac {1}{2}\left ( 0, \boldsymbol{\varOmega }_i \right ) \bullet \boldsymbol{q}_i, \end{equation}

where $\bullet$ is the quaternion product, defined for two quaternions $\boldsymbol{p} = (p_0, \tilde {\boldsymbol{p}})$ and $\boldsymbol{q} = (q_0, \tilde {\boldsymbol{q}})$ as $\boldsymbol{p} \bullet \boldsymbol{q} = (p_0 q_0 - \tilde {\boldsymbol{p}}\boldsymbol{\cdot }\tilde {\boldsymbol{q}}, \; p_0\tilde {\boldsymbol{q}} + q_0 \tilde {\boldsymbol{p}} + \tilde {\boldsymbol{p}} \times \tilde {\boldsymbol{q}})$ .

Pairing (2.10) with an equation describing the evolution of segment positions, namely

(2.11) \begin{equation} \frac {{\textrm d}\boldsymbol{Y}_i}{{\textrm d} t} = \boldsymbol{V}_i, \end{equation}

we can discretise our system temporally using a second-order, geometric backward differential formula scheme, as described in Schöller et al. (Reference Schöller, Townsend, Westwood and Keaveny2021). These discretised equations, as well as the constraint given by (2.8), define a nonlinear system which we can solve iteratively using Broyden’s method (Broyden Reference Broyden1965).

In this work, the filament is discretised into $N=20$ segments of length $\Delta L = 2.2a$ , such that the filament has length $L = N\Delta L = 44a.$ This separation has been chosen to avoid overlap between neighbouring segments when the filament bends. The twist stiffness is taken to be equal to the bending stiffness in the $\hat {\boldsymbol{\mu }}$ direction, $K_B=K_T$ . In the isotropic case, where $\beta =1$ and $\kappa _0=0$ , the magnitude of the twist has little effect on the dynamics (Clarke et al. Reference Clarke, Hwang and Keaveny2024). We investigate the effect of varying three non-dimensional parameters: the intrinsic curvature, $\kappa _0$ , as in (2.4); the bending anisotropy, quantified by $\beta ,$ as in (2.4); and the non-dimensional follower force, $f,$ as defined in (2.5). We focus on ranges of $\kappa _0$ and $\beta$ that follow values explored in other studies (Rallabandi et al. Reference Rallabandi, Wang and Potomkin2022; Chakrabarti & Saintillan Reference Chakrabarti and Saintillan2019b , Reference Chakrabarti and Saintillana ; Sartori et al. Reference Sartori, Geyer, Scholich, Jülicher and Howard2016; Chakrabarti et al. Reference Chakrabarti, Fürthauer and Shelley2022; Cass & Bloomfield-Gadêlha Reference Cass and Bloomfield-Gadêlha2023; Lindemann & Lesich Reference Lindemann and Lesich2016; Wang et al. Reference Wang, Gsell, D’Ortona and Favier2023), and a range of $f$ which allows us to access and explore the whirling and beating states. All initial-value problems (IVPs) use the initial condition of a vertically upright filament, with random perturbations to the force on each segment. All time scales are non-dimensionalised by the filament’s relaxation time, $\tau$ , the characteristic time scale associated with a bent filament returning to its undeformed configuration, which we obtain numerically as described in Clarke et al. (Reference Clarke, Hwang and Keaveny2024).

2.2. Bifurcation and stability analyses

The filament model offers a method to generate numerical solutions given an initial configuration. In order to understand the bifurcations between different states that arise from the model, we need to first obtain these states and then assess their stability. We use a Jacobian-free Newton–Kyrlov subspace method (JFNK) to find the stable or unstable steady-state and time-periodic solutions, as described in Viswanath (Reference Viswanath2007), Willis (Reference Willis2019) and Clarke et al. (Reference Clarke, Hwang and Keaveny2024). The JFNK seeks solutions which minimise the error between the initial state and the state after being advanced for a period of time. Therefore, we can only use this to identify steady and time-periodic solutions – not, for instance, quasiperiodic or chaotic dynamics. The system is solved using Newton’s method, where the generalised minimal residual (GMRES) method is employed in each Newton iteration to solve for the update. Once we have a solution for a single set of parameters, we can then use continuation to track the solution branch for different values of $f, \kappa _0$ or $\beta$ , including values of these parameters for which the solution may be unstable.

After finding the steady-state or time-periodic base state, we can analyse its stability numerically. This is achieved by using the Arnoldi method to establish the eigenvalues of the system’s Jacobian, for steady-state solutions, or the Floquet matrix exponent, if the state is time-periodic. Further details of this method can be found in Clarke et al. (Reference Clarke, Hwang and Keaveny2024).

It may be that two solutions are stable for the same values of parameters. If this is the case, we expect for there to be an unstable state which lies on the boundary of the two stable solutions in the state space. We can identify the unstable solution by using the bisection algorithm outlined in Skufca, Yorke & Eckhardt (Reference Skufca, Yorke and Eckhardt2006), which iteratively bisects between initial conditions leading to either stable state. While this method allows us to identify an unstable solution lying on the boundary surface (separatrix) between the stable states in the state space, we only observe this transiently, before the filament converges to one of the stable states. Furthermore, if we expect multiple unstable solutions to exist for a set of parameters, bisection can only compute a solution stable (or attractive) along the directions tangential to the boundary surface between the two stable states.

3.1. Planar dynamics

We begin by considering an isotropic filament with finite intrinsic curvature, i.e. $\beta = 1,$ and $ \kappa _0 \neq 0$ , with its motion restricted to the plane in which it is deformed by its intrinsic curvature. The corresponding phase diagram of the planar filament states in $\kappa _0{-}f$ space, where the plane is aligned with the filament’s intrinsic curvature, is obtained by running one IVP for each set of parameters. Each IVP has an initial condition of a vertically upright filament with random, small perturbations to the forcing on a subset of the segments. We summarise the results in figure 2(a). For this study we consider $\kappa _0 \in [0,3\pi /4]$ . This range of values is appropriate as it includes values used in the literature in numerical studies (Wang et al. Reference Wang, Gsell, D’Ortona and Favier2023), in particular to recover dynamics similar to those observed experimentally (Chakrabarti & Saintillan Reference Chakrabarti and Saintillan2019a ; Sartori et al. Reference Sartori, Geyer, Scholich, Jülicher and Howard2016). These IVPs reveal that we observe two types of behaviour: steady states below a critical value of the forcing, or self-sustained beating oscillations when this critical value is exceeded, matching observations in Wang et al. (Reference Wang, Gsell, D’Ortona and Favier2023).

For $\kappa _0=0$ , the steady state is a vertically upright filament. However, with increasing $\kappa _0 \neq 0,$ the steady state becomes non-trivial. For $f=0,$ it is given by the prescribed intrinsic curvature, such that

(3.1) \begin{equation} \boldsymbol{t}(s) = (\sin (\kappa _0 s/L), 0, \cos (\kappa _0 s/L)), \;\;\; s\in [0,L], \end{equation}

as shown in figure 2(b). When $f\gt 0,$ this steady state changes; the tip of the filament is pushed by the follower force, until a new steady state is reached in which the elastic forces balance the follower force. This gives a non-trivial steady state, as shown for $\kappa _0=3\pi /4$ in figure 2(c).

Figure 2.

(a) Results from 2-D numerical simulations of a filament with intrinsic curvature $(\kappa _0 \geqslant 0,$ $ \beta = 1)$ . The black line is drawn to indicate the approximate boundaries of the solution regions. For $f\lt f^* \approx 35.3$ , filaments remain in their steady non-trivial state which varies with (b) $\kappa _0$ (shown for $f=0,$ $\kappa _0 \in \{0,\pi /4,\pi /2,\pi /4 \}$ ) and (c) $f$ (shown for $\kappa _0 = 3\pi /4,$ $ f \in \{0,5,10,20\}$ ). For $f\gt f^*$ filaments undergo asymmetric oscillations for $\kappa _0\gt 0$ , as shown for (d) $\kappa _0=3\pi /4$ .

When the force exceeds a critical value, the filament buckles and we observe planar oscillations. To find this critical value for each $\kappa _0,$ we track the steady state using the JFNK and use stability analysis. For $\kappa _0=0,$ this reveals that the buckling occurs at $f = f^* \approx 35.3,$ through a Hopf bifurcation. For $\kappa _0\gt 0$ , the bifurcation occurs through the same mechanism, but we observe a small increase in the critical forcing of approximately $0.5\,\%$ from $\kappa _0=0$ to $\kappa _0=3\pi /4.$ For $\kappa _0 = 0,$ numerical simulations after the bifurcation show that the planar beating oscillations are symmetric. However, this symmetry is broken by intrinsic curvature; for $\kappa _0 \gt 0,$ simulations reveal asymmetric beating patterns as shown in figure 2(d), and observed in Wang et al. (Reference Wang, Gsell, D’Ortona and Favier2023). We plot the frequencies of these solutions for various $\kappa _0$ in Appendix C. Close to the bifurcation this asymmetry is strongest, but further away from the bifurcation, where the follower force is strong, the effect of intrinsic curvature is less pronounced and the beating becomes more symmetric.

3.2. Non-planar dynamics: initial-value problems

When the filament is free to deform in any direction, we observe dramatic changes in filament dynamics as we vary $\kappa _0$ and $f$ . We summarise the filament states in figure 3(a). For the smallest values of the forcing, we observe the same steady states as in the planar case. However, we now find that the filament exhibits fully 3-D dynamics after buckling and the critical value of the forcing has a dependence on $\kappa _0$ .

Figure 3.

(a) Results of numerical simulations of a filament with intrinsic curvature in three dimensions $(\kappa _0 \geqslant 0, \beta = 1)$ . The black lines are drawn to indicate the approximate boundaries of the solution regions. (b) Snapshots of filaments over one period undergoing whirling dynamics for three values of intrinsic curvature close to ( $f=40$ ) and further from ( $f=100$ ) the bifurcation. We see that asymmetries in whirling for $\kappa _0\gt 0$ are more pronounced close to the bifurcation, and for higher $\kappa _0.$ (c) Snapshots of filaments over one period for $f=200,$ increasing $\kappa _0$ from $0$ to $3\pi /4$ . The dynamics changes from planar beating for $\kappa _0 = 0$ (left) to P1 for $\kappa _0 \gt 0$ (right).

For $\kappa _0 = 0$ (Ling et al. Reference Ling, Guo and Kanso2018; Clarke et al. Reference Clarke, Hwang and Keaveny2024), we observe the symmetric whirling state immediately after the buckling; a rigid-body rotation of the centreline, with the filament tip tracing a circle over one period. Increasing the forcing further, whirling is replaced by a quasiperiodic solution which we term QP1 – an elliptical whirling which rotates unidirectionally. Solution QP1 only exists for a small region of $f,$ before being replaced with planar beating. Planar beating is also observed in the two-dimensional (2-D) simulations, but is now symmetric (as $\kappa _0=0$ ) and can occur in any plane (due to the rotational symmetry of the problem).

Increasing $\kappa _0 \gt 0,$ the initial buckling occurs at smaller values of the forcing. The whirling behaviour still exists but is no longer symmetric, as shown in figure 3(b) and supplementary movie 1 available at [https://doi.org/10.1017/jfm.2025.10787]. The frequencies of the whirling solutions are plotted in Appendix C. In fact, the asymmetry is strongest as we increase the intrinsic curvature, and choose forcing values closer to the initial bifurcation. Increasing the forcing has the effect of smoothing out the asymmetry, as the follower force dominates over the torques maintaining the intrinsic curvature.

Increasing $f$ further, we observe that whirling gives way to a new periodic solution for $\kappa _0 \gt 0$ , which we call P1. This consists of a bent-over beating, where the filament bends in the plane in which it would deform in the absence of the follower force, and beats in the plane perpendicular to it, as shown in figure 3(a) inset and supplementary movie 2. The intrinsic curvature breaks left–right symmetry, in such a way that the filament beats in the left half of the plane. In figure 3(c), we fix $f=200$ and increase $\kappa _0$ from $0$ to $3\pi /4$ to show how P1 continues from the planar beating solution. We see that increasing $\kappa _0$ increases the maximum perpendicular distance of the tip from the beat plane. We plot the frequencies of P1 for various $f$ and $\kappa _0$ in Appendix C. Varying $\kappa _0\gt 0$ , planar beating is no longer visible with IVPs alone – it appears that the 2-D asymmetric beating patterns observed in the planar case are unstable in the 3-D model.

For larger values of $\kappa _0$ , we note that the whirling state is also no longer visible using IVPs, regardless of the initial filament configuration. Furthermore, for large $\kappa _0$ and $f,$ we see a new periodic mode appearing, which we call P2. Mode P2 beats similarly to P1, but traces a distinct shape of 8 with its tip over one period (see figure 3 a, inset).

3.3. Bifurcation and stability analyses

To explore in more detail how filament dynamics changes with $\kappa _0,$ we track the steady and time-periodic states with JFNK and perform a stability analysis to establish bifurcations within the system. An overview of the resulting bifurcation diagrams for $f\in [0,200]$ is shown in figure 4, where we plot the non-dimensionalised mean bending energy, given by

(3.2) \begin{equation} \bar {E}_b = \frac {1}{T} \int _0^T \int _0^L \frac {1}{2} \kappa (s)^2 \;{\textrm d}s \;{\textrm d}t, \end{equation}

for each solution branch for various values of $\kappa _0$ . Here $\kappa (s)$ is the local curvature of the filament. While the bending energy is usually defined with respect to the intrinsic curvature, we instead use the standard bending energy described above for ease of comparison as we vary $\kappa _0$ . We now discuss the key changes in the bifurcation diagrams as we increase $\kappa _0$ . We provide a more detailed version of this discussion in Appendix B.1.

Figure 4.

Bifurcation diagrams showing the different solutions as we vary $f$ , and the corresponding mean bending energy of these solutions, $\bar {E}_b$ , for an isotropic filament ( $\beta = 1$ ) with intrinsic curvature for (a) $\kappa _0 = 0$ , (b) $\kappa _0 = \pi /4$ , (c) $\kappa _0 = \pi /2$ and (d) $\kappa _0 = 3\pi /4$ . Dashed (solid) lines refer to unstable (stable) solutions. Square markers are used to indicate the location of bifurcations as we vary $f$ . Inset, left, shows the bifurcation diagram near the initial buckling event and inset, right, shows the bifurcation diagram around the bistable regions, if these are present. As $\kappa _0$ increases we see the emergence of the P1 solution branch (b–d), from which whirling bifurcates (b,c, left inset). We observe bistability between P1 and whirling (b,c, right inset) and, for the largest values of $\kappa _0$ , whirling vanishes completely (d).

3.3.1. Isotropic case

We begin by recapping the isotropic case, where $\kappa _0 = 0$ , shown in figure 4(a). Focusing first on the initial buckling event, in Clarke et al. (Reference Clarke, Hwang and Keaveny2024) and Schnitzer (Reference Schnitzer2025) it was shown that two complex conjugate pairs of eigenmodes become unstable at $f=f^*$ (see figure 5 a). We hence identify this bifurcation as a double Hopf bifurcation.

Figure 5.

The real part of the four largest eigenvalues, $\lambda _i$ for $i=1,\ldots ,4$ , associated with the steady state for $\beta =1$ : (a) $\kappa _0 = 0$ and (b) $\kappa _0 = 3\pi /4$ . Analysing the corresponding eigenvectors at each bifurcation, obtained as described in the main text and displayed on the right, allows us to identify the branches corresponding to off-plane and in-plane beating. (c) The real part of the largest eigenvalue associated with the steady state for various values of $\kappa _0.$ We see that the eigenvalue becomes unstable for smaller $f$ as we increase $\kappa _0.$ (d) The real part of the largest eigenvalue associated with the P1 state for various values of $\kappa _0$ . For $\kappa _0 \gt \kappa _0^*$ , we see that the eigenvalue remains stable.

To understand the buckling bifurcation further, we visualise the eigenvectors. Defining each eigenvalue–eigenvector pair as $(\lambda _i, \boldsymbol{x}_i)$ , we note that the $\boldsymbol{x}_i$ are complex for $i=1,\ldots ,4$ . Therefore, to interpret their physical meaning, we take linear combinations of the eigenvectors $\boldsymbol{v}_1 = a(\boldsymbol{x}_1 + \boldsymbol{x}_2), \boldsymbol{v}_2 = a(\boldsymbol{x}_1 - \boldsymbol{x}_2)$ and similar for $\boldsymbol{v}_3$ and $\boldsymbol{v}_4$ using $\boldsymbol{x}_3$ and $\boldsymbol{x}_4$ , where $a$ is an arbitrary constant chosen to emphasise the structure of the eigenmode. As the magnitude of the state variable relates to angles of rotation of the local basis, taking larger values of $a$ results in more exaggerated filament deformation. We show these eigenvectors in figure 5 for $a=5.$ We observe that the first complex conjugate pair is associated with beating in the $(x,z)$ plane, and the second is associated with beating in the perpendicular plane, as shown by the eigenvectors in figure 5(a) (right). It is these modes that combine to give rise to the stable whirling solution observed after the bifurcation, as well as the unstable planar beating branch. For larger $f$ , whirling becomes unstable and gives rise to QP1. The QP1 branch exists for a small region of $f$ , and vanishes when planar beating stabilises for larger $f$ (see figure 4 a, inset). Further discussion on the bifurcation diagram in the isotropic case, including exploration of higher forcing values, can be found in Clarke et al. (Reference Clarke, Hwang and Keaveny2024).

3.3.4. Loss of whirling solution

In the previous sections we identify three key bifurcations: P1 becoming unstable, P1 subsequently restabilising and whirling becoming unstable. As we increase $\kappa _0$ , these three bifurcations occur within a smaller range of $f$ , until the intrinsic curvature reaches a critical value, and the bifurcations coalesce. For larger $\kappa _0,$ for instance $\kappa _0=3\pi /4$ as shown in figure 4(d), the whirling solution vanishes completely, and we only observe the stable P1 solution branch.

We can explore the disappearance of the whirling solution by performing Floquet analysis on the P1 solution for various values of $\kappa _0$ . The real part of the dominant eigenvalue for each $\kappa _0$ is shown in figure 5(d). We see that the two bifurcations, corresponding to P1 becoming unstable and then restabilising, get closer together as we increase $\kappa _0$ , until they merge at a critical value, $\kappa _0=\kappa _0^*$ . At this value, the whirling solution collapses onto the P1 branch, also terminating QP1 and QP1* in the process. Beyond this $\kappa _0,$ only P1 is stable in this region of $f.$ We can identify the loss of the whirling solution by again revisiting the initial buckling event. For $\kappa _0 \gt \kappa _0^*$ , it appears that the growth rate associated with the off-plane modes is so large compared with the growth rate of the in-plane modes that the whirling solution never exists. Instead, we only observe P1, matching observations from IVPs.

Figure 6.

(a) Results of numerical simulations of an anisotropic filament with no intrinsic curvature $(\kappa _0 = 0, \beta \geqslant 1)$ . The black lines are drawn to indicate the approximate boundaries of the solution regions. (b) Snapshots of filaments over one period undergoing whirling dynamics in the isotropic case, close to the bifurcation $(\beta = 1, f = 36)$ , an anisotropic case near the bifurcation $ (\beta =1.2, f=46)$ and an anisotropic case further from the bifurcation $(\beta = 1.2, f=100)$ . The three dynamics are indicated on the phase diagram in (a) using numbered circles.

4.1. Initial-value problems

In the previous section we investigated the effect of intrinsic curvature defined in the plane defined by $\hat {\boldsymbol{\nu }}$ and $\hat {\boldsymbol{t}}$ . Now, we set $\kappa _0 = 0$ and focus on the effect of introducing anisotropy through the bending stiffness by introducing the parameter $\beta$ , which describes the ratio of the bending stiffnesses about the $\hat {\boldsymbol{\nu }}$ and $\hat {\boldsymbol{\mu }}$ axes. Throughout this paper we focus on $\beta \gt 1,$ but the $\beta \lt 1$ case follows immediately by considering a rotation of the axes, and redefining the non-dimensional follower force as $\tilde {f} = \beta f$ (cf. (2.5)). In figure 6(a), we show the results of IVPs for $\beta \in [1,1.5]$ and $f\in [0,300].$ Regardless of $\beta$ and $f$ , the steady state is the vertically upright filament. For $\beta =1$ , this buckles into whirling, while for $\beta \gt 1$ the steady state buckles into planar beating in the $(x,z)$ plane. For moderate values of $\beta$ , beating becomes unstable as we increase the forcing, with whirling appearing in its place. The frequency of these whirling solutions for different $\beta$ and $f$ values can be found in Appendix C. Close to the bifurcation in which planar beating appears to become unstable, whirling is very similar to in-plane beating, with a small off-plane component. This causes the tip to trace an ellipse with a major axis that aligns with the $\hat {\boldsymbol{x}}$ axis, as shown in figure 6(b), middle, and supplementary movie 3. The large amplitude in the $\hat {\boldsymbol{x}}$ direction close to the bifurcation suggests that we can expect the whirling solution to branch from the planar beating solution branch, instead of the steady state. Further from the bifurcation, however, the off-plane component grows and whirling resembles the isotropic case more closely (see figure 6 b, right). For larger forcing values, whirling becomes unstable and planar beating in the $(x,z)$ plane reappears. As we increase $\beta$ , the region in which whirling is stable appears to decrease, until $\beta = 1.35$ , where we no longer observe whirling through IVPs. For these $\beta$ values, we only observe planar beating in the $(x,z)$ plane after buckling.

Figure 7.

(a) Bifurcation diagrams showing the different solutions as we vary $f$ , and the corresponding mean bending energy of these solutions, $\bar {E}_b$ , for an anisotropic filament with no intrinsic curvature ( $\kappa _0 = 0$ ) for $\beta = 1.2,\ 1.3$ and $1.4$ (top to bottom). Dashed (solid) lines refer to unstable (stable) solutions. Square markers are used to indicate the location of bifurcations as we vary $f$ . We see the separation of the in-plane and off-plane solution branches as $\beta$ increases, and observe bistability between in-plane beating and whirling (inset). For larger values of $\beta$ , the whirling solution branch vanishes completely (bottom). (b) The real part of the four largest eigenvalues associated with the steady state for $\beta = 1.5$ . As described in the main text, we visualise the eigenvectors, $\boldsymbol{x}_i$ for $i=1,\ldots ,4$ , by taking linear combinations of the complex conjugate pairs, for example $\boldsymbol{v}_1 = a(\boldsymbol{x}_1 + \boldsymbol{x}_2)$ , where $a$ is a constant. We plot the eigenvectors for $a=5$ to identify the branches corresponding to in-plane beating (bottom left) and off-plane beating (bottom right). (c) The real part of the second largest eigenvalue associated with off-plane beating for various values of $\beta .$ We associate this eigenvalue becoming stable with P4 collapsing onto the off-plane beating branch. (d) The real part of the largest eigenvalue associated with the in-plane beating state for various values of $\beta$ . For $\beta \gt \beta ^*$ , we see that the eigenvalue remains stable.

4.2. Bifurcation and stability analyses

Using JFNK to track the steady and periodic solutions and then performing a stability analysis on the subsequent branches allows us to generate the bifurcation diagrams shown in figure 7(a). We outline the key features of the bifurcation diagrams below, and provide further details in Appendix B.2.

Focusing first on the initial bifurcation, we see that the critical forcing value at which the steady state buckles remains at $f=f^*$ regardless of $\beta$ . At this point, the planar beating branch in the $(x,z)$ plane, which we will henceforth call ‘in-plane’ beating (as opposed to ‘off-plane’ beating, in the $(y,z)$ plane), is born through a Hopf bifurcation. We can also identify the unstable off-plane beating branch for each $\beta$ , which also appears to connect to the trivial steady state, bifurcating at $f= \beta f^*$ . These solutions are present regardless of $\beta$ , as shown in figure 7(a). These branches are linked to the initial buckling event; introducing $\beta \neq 1$ has the effect of separating the critical values at which the beating in either plane becomes unstable, similar to the effect of having $\kappa _0 \neq 0$ . This is a consequence of the non-dimensionalisation of the follower force, as given in (2.5). Accordingly, we expect eigenmodes to become unstable at $f=f^*$ (corresponding to buckling in the $(x,z)$ plane) and $f = \beta f^*$ (corresponding to buckling in the $(y,z)$ plane). Thus, the magnitude of $\beta$ dictates which plane we initially buckle into, and furthermore when successive buckling modes become activated. To verify the effect of $\beta$ on the initial buckling bifurcation, we plot the dominant eigenvalues for $\beta =1.5$ in figure 7(b). As expected, we find that the first bifurcation occurs at $f= f^*,$ with corresponding planar eigenmodes lying in the $(x,z)$ plane. Therefore, the filament is expected to undergo planar deformations in this plane, as is confirmed by the numerical simulations. The second mode, corresponding to buckling in the perpendicular plane, is found to occur at $f = \beta f^* \approx 53.0$ .

For moderate $\beta$ (see figure 7 a, top, middle) we find that, for a larger value of the forcing, in-plane beating becomes unstable and the whirling solution is born. Between the bifurcations in which in-plane beating becomes stable again and whirling becomes unstable, there is a region of bistability. Bisection in this region yields a single solution branch, which we call P4. As P4 is periodic, we can track the behaviour using JFNK. In doing so, we find that P4 connects to the off-plane beating branch, rather than the whirling branch, as shown in figure 7(a) (top, inset). Floquet analysis of the off-plane beating branch confirms this, showing that the second largest eigenvalue, $\lambda _2$ , becomes stable when P4 collapses onto the branch. By performing a Floquet analysis on the off-plane beating branch for various values of $\beta ,$ as shown in figure 7(c), we see that $\lambda _2$ becomes stable at higher $f$ as we increase $\beta$ , indicating that P4 joins the off-plane beating branch at higher values of $f$ as we increase $\beta$ . This trend continues until a critical value of $\beta ,$ beyond which P4 now connects to the whirling branch, after which both branches terminate. Consequently, for these values of $\beta ,$ the whirling branch only exists when it is stable (see figure 7 a, middle).

The three bifurcation points corresponding to in-plane beating becoming unstable, becoming stable again and the whirling branch disappearing all occur within smaller ranges of $f$ as we increase $\beta$ , until the three bifurcations collide and we see only one stable branch – in-plane beating, as shown in figure 7(a) (bottom) for $\beta = 1.4.$ Indeed, a stability analysis of the in-plane beating solution, shown in figure 7(d), confirms that the range of $f$ for which beating is unstable decreases with increasing $\beta ,$ up to a critical value of $\beta = \beta ^*$ where the two pitchfork bifurcations collide, which likely also coalesces with the P4–whirling bifurcation. Beyond this, beating is stable for all values of $f$ in our tested range and, as in the intrinsic curvature case, whirling ceases to exist entirely. Both the Floquet analysis and IVPs suggest that the whirling branch collapses onto the beating branch at $\beta =\beta ^*,$ resulting in the loss of P4 and whirling as distinct solutions. The disappearance of whirling can again be explained through the large separation in critical force values between in-plane and off-plane beating solutions bifurcating from the trivial state.