2.1. Elastohydrodynamic model

As shown in figure 1, we consider an inextensible elastic filament of length $L_*$ and bending stiffness $B_*$ clamped normally to a flat wall at one end and subjected to a tangential compressive force of magnitude $\mathcal {F}_*$ at the other. The filament is immersed in a background liquid of viscosity $\eta _*$ . We neglect inertia of the fluid and the filament and model the latter as a homogeneous Kirchhoff rod having a circular cross-section of radius $\kappa L_*$ (Landau et al. Reference Landau, Pitaevskii, Kosevich and Lifshitz2012), further assuming that it is sufficiently slender such that the hydrodynamics can be described by resistive force theory (Graham Reference Graham2018). Throughout the paper, a subscript asterisk indicates a dimensional quantity; a superscript asterisk denotes complex conjugation; and $r$ and $i$ subscripts indicate the real and imaginary parts of a complex-valued quantity, respectively.

Figure 1.

(a) Dimensional schematic. (b) Forces and moments acting on an infinitesimal filament segment.

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

(2.1a,b) \begin{equation} \frac {\partial {\boldsymbol {F}}_*}{\partial s_*}+{\boldsymbol {f}}_*^H={\mathbf {0}}, \quad \frac {\partial {\boldsymbol {M}}_*}{\partial s_*}+{\hat {\boldsymbol {t}}}\times {\boldsymbol {F}}_*={\mathbf {0}}, \end{equation}

in which

(2.2) \begin{equation} {\boldsymbol {f}}_*^H=-\frac {4\pi \eta _*}{\ln (1/\kappa )}\left (\mathsf {I}-\frac {1}{2}{\hat {\boldsymbol {t}}}{\hat {\boldsymbol {t}}}\right )\boldsymbol {\cdot } \frac {\partial {\boldsymbol {r}}_*}{\partial t_*} \end{equation}

is the hydrodynamic force per unit length acting on the filament (see figure 1 b). The local anisotropic drag law (2.2) constitutes a leading-order asymptotic approximation of the hydrodynamics in the slender-filament limit $\kappa \to 0$ ; under that approximation, known as resistive force theory, the interaction of the filament with itself and with the wall is neglected despite the relative error associated with that approximation being only logarithmically small, scaling as $1/\ln (1/\kappa )$ (Graham Reference Graham2018).

The above governing equations are supplemented by initial conditions (which we do not specify at this stage), boundary conditions and a constitutive relation for the internal moment. The boundary conditions at the clamped end of the filament are

(2.3a,b) \begin{equation} {\boldsymbol {r}}_*={\mathbf {0}}, \quad {\hat {\boldsymbol {t}}}=\boldsymbol {\hat {\imath }} \quad \text {at} \quad s_*=0, \end{equation}

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

(2.4a,b) \begin{equation} {\boldsymbol {F}}_*= -\mathcal {F}_*{\hat {\boldsymbol {t}}}, \quad {\boldsymbol {M}}_*={\mathbf {0}} \quad \text {at} \quad s_*=L_*, \end{equation}

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

(2.5) \begin{equation} {\boldsymbol {M}}_*=B_*{\hat {\boldsymbol {t}}}\times \frac {\partial {\hat {\boldsymbol {t}}}}{\partial s_*}, \end{equation}

representing the resistance of the filament to bending, which is cross-sectionally isotropic given the assumption of a circular cross-section. In general, the constitutive relation includes an additional term representing the resistance of the filament to twisting. In the present scenario, however, where there are no external moments acting in the tangential direction and the bending is cross-sectionally isotropic, the twist term in the constitutive relation vanishes identically; this is verified in Appendix A following Landau et al. (Reference Landau, Pitaevskii, Kosevich and Lifshitz2012). The absence of twist in this scenario was also remarked on by Ling et al. (Reference Ling, Guo and Kanso2018). In contrast, Clarke et al. (Reference Clarke, Hwang and Keaveny2024) include twist as their model does not involve a thin-flament approximation of the hydrodynamics, so that viscous moments act in the tangential direction.

2.2. Dimensionless formulation

It is convenient to adopt a dimensionless convention where lengths are normalised by $L_*$ , forces by $B_*/L_*^2$ , moments by $B_*/L_*$ and time by $4\pi \eta _*(L_*^4/B_*)/\ln (1/\kappa )$ . Dimensionless fields are denoted similarly to their dimensional counterparts, only with the subscript asterisks omitted. The problem thus consists of the partial differential equations

(2.6a-d) \begin{equation} \frac {\partial {\boldsymbol {r}}}{\partial s}={\hat {\boldsymbol {t}}}, \quad \frac {\partial {\boldsymbol {F}}}{\partial s}-\left (\mathsf {I}-\frac {1}{2}{\hat {\boldsymbol {t}}}{\hat {\boldsymbol {t}}}\right )\boldsymbol {\cdot } \frac {\partial {\boldsymbol {r}}}{\partial t} ={\mathbf {0}}, \quad \frac {\partial {\boldsymbol {M}}}{\partial s}+{\hat {\boldsymbol {t}}}\times {\boldsymbol {F}}={\mathbf {0}}, \quad {\boldsymbol {M}}={\hat {\boldsymbol {t}}}\times \frac {\partial {\hat {\boldsymbol {t}}}}{\partial s}; \end{equation}

the inextensibility constraint

(2.7) \begin{equation} {\hat {\boldsymbol {t}}}\boldsymbol {\cdot }{\hat {\boldsymbol {t}}} = 1; \end{equation}

and the boundary conditions

(2.8a,b) \begin{equation} {\boldsymbol {r}}={\mathbf {0}}, \quad {\hat {\boldsymbol {t}}}=\boldsymbol {\hat {\imath }} \quad \text {at} \quad s=0 \end{equation}

and

(2.9a,b) \begin{equation} {\boldsymbol {F}}= -\mathcal {F}{\hat {\boldsymbol {t}}}, \quad {\boldsymbol {M}}={\mathbf {0}} \quad \text {at} \quad s=1, \end{equation}

wherein the dimensionless parameter

(2.10) \begin{equation} \mathcal {F}=\frac {\mathcal {F}_*L_*^2}{B_*} \end{equation}

measures the magnitude of the follower force. The problem depends on this sole parameter, aside for the initial conditions which we do not specify at this stage.

For all $\mathcal {F}$ , the above problem possesses the steady solution

(2.11a-d) \begin{equation} {\boldsymbol {r}}=s\boldsymbol {\hat {\imath }}, \quad {\hat {\boldsymbol {t}}} = \boldsymbol {\hat {\imath }}, \quad {\boldsymbol {F}} =-\mathcal {F}\boldsymbol {\hat {\imath }}, \quad {\boldsymbol {M}}={\mathbf {0}}, \end{equation}

where the filament is undeformed (straight) and under uniform compression by the follower force. Our aim will be to study the nonlinear dynamics of the filament for $\mathcal {F}$ near the buckling threshold where this steady ‘base state’ first becomes unstable.

3.1. Eigenvalue problem

As a preliminary step, we must characterise the loss of stability of the base state. To this end, we assume small-magnitude perturbations of the form

(3.1a-d) \begin{equation} {\boldsymbol {r}}-s\boldsymbol {\hat {\imath }} = e^{\lambda t}{\boldsymbol{\tilde{r}}} + \text {c.c.} \quad {\hat {\boldsymbol {t}}}-\boldsymbol {\hat {\imath }}= e^{\lambda t}\tilde {\boldsymbol {t}} + \text {c.c.}, \quad {\boldsymbol {F}}+\boldsymbol {\hat {\imath }}\mathcal {F}= e^{\lambda t}\tilde {{\boldsymbol {F}}} + \text {c.c.}, \quad {\boldsymbol {M}}= e^{\lambda t}\tilde {{\boldsymbol {M}}} + \text {c.c.}, \end{equation}

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

(3.2a-d) \begin{equation} {\boldsymbol{\tilde{r}}}^{\prime}=\tilde {\boldsymbol {t}}, \quad \tilde {{\boldsymbol {F}}}^{\prime}-\lambda \left (\mathsf {I}-\frac {1}{2}\boldsymbol {\hat {\imath }}\boldsymbol {\hat {\imath }}\right )\boldsymbol {\cdot } {\boldsymbol{\tilde{r}}}={\mathbf {0}}, \quad \tilde {{\boldsymbol {M}}}^{\prime}+\mathcal {F}\boldsymbol {\hat {\imath }}\times \tilde {\boldsymbol {t}}+\boldsymbol {\hat {\imath }}\times \tilde {{\boldsymbol {F}}}={\mathbf {0}}, \quad \tilde {{\boldsymbol {M}}}=\boldsymbol {\hat {\imath }}\times \tilde {\boldsymbol {t}}^{\prime}, \end{equation}

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

(3.3) \begin{equation} \boldsymbol {\hat {\imath }}\boldsymbol {\cdot }\tilde {\boldsymbol {t}} =0; \end{equation}

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

(3.4a-d) \begin{equation} {\boldsymbol{\tilde{r}}}={\mathbf {0}}, \quad \tilde {\boldsymbol {t}}={\mathbf {0}} \quad \text {at} \quad s=0; \qquad \tilde {{\boldsymbol {F}}}= -\mathcal {F}\tilde {\boldsymbol {t}}, \quad \tilde {{\boldsymbol {M}}}={\mathbf {0}} \quad \text {at} \quad s=1. \end{equation}

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

(3.5) \begin{equation} \tilde {x}^{\prime\prime\prime\prime}+\mathcal {F}\tilde {x}^{\prime\prime}+\lambda \tilde {x}=0 \end{equation}

and the boundary conditions

(3.6a-d) \begin{equation} \tilde {x}=0, \quad \tilde {x}^{\prime}=0 \quad \text {at} \quad s=0; \qquad \tilde {x}^{\prime\prime}=0, \quad \tilde {x}^{\prime\prime\prime}=0 \quad \text {at} \quad s=1. \end{equation}

The eigenvalue problem consisting of (3.5) and (3.6) was first derived and solved by De Canio et al. (Reference De Canio, Lauga and Goldstein2017) under the assumption of planar deformations. While the general solution to (3.5) is readily expressed in closed form, calculating the eigenvalues $\lambda$ as a function of $\mathcal {F}$ must, in general, be done numerically. In agreement with De Canio et al. (Reference De Canio, Lauga and Goldstein2017), we thereby find that the base state is monotonically stable for $\mathcal {F}\lt \mathcal {F}_{o}\approx 20.051$ , oscillatory stable for $\mathcal {F}_0\lt \mathcal {F}\lt \mathcal {F}_c\approx 37.695$ and oscillatory unstable for $\mathcal {F}\gt \mathcal {F}_c$ ; at the instability threshold, $\mathcal {F}=\mathcal {F}_c$ , the complex-conjugate pair of neutral eigenvalues is $\lambda =\pm i\omega$ , wherein $\omega \approx 191.26$ (see figure 2 a). We denote by $\varphi (s)=\varphi _r(s)+i\varphi _i(s)$ the eigenfunction normalised to unity at $s=1$ that is associated with the threshold eigenvalue $i\omega$ (see figure 2 b); $\varphi ^*(s)$ is accordingly the similarly normalised eigenfunction associated with the conjugate eigenvalue $-i\omega$ .

Figure 2.

(a) Complex growth rate $\lambda =\lambda _r+i\lambda _i$ as a function of $\mathcal {F}$ for the most unstable eigenfunctions of the eigenvalue problem consisting of (3.5) and (3.6). (b) The threshold eigenfunction $\varphi (s)=\varphi _r(s)+i\varphi _i(s)$ corresponding to the imaginary growth rate $\lambda =i\omega$ .

3.2. Double-Hopf bifurcation

The eigenfunctions at the instability threshold describe periodic beating oscillations of the filament in the plane through the origin and normal to $\boldsymbol {\hat {\imath }}\times {\hat {\boldsymbol {e}}}_1$ , of arbitrary magnitude and phase. Since $\tilde {y}(s)$ satisfies the same eigenvalue problem as $\tilde {x}(s)$ , there is clearly an independent pair of eigenfunctions representing similar beating oscillations rotated by $\pi /2$ about the undeformed configurations, i.e. in the plane through the origin and normal to $\boldsymbol {\hat {\imath }}\times {\hat {\boldsymbol {e}}}_2$ . As first pointed out by Clarke et al. (Reference Clarke, Hwang and Keaveny2024), this means that the filament loses stability via a double-Hopf bifurcation, where two pairs of complex-conjugate eigenvalues become unstable simultaneously; since the multiplicity here is induced by symmetry, the pairs are identical and have associated with them a full set of eigenfunctions. A key observation of Clarke et al. (Reference Clarke, Hwang and Keaveny2024) is that the linear superposition of beating modes allowing for differences in the phase and plane of motion generally produces non-planar whirling oscillations of the filament.

3.3. Linear approximation at the threshold

It will be convenient for us to employ complex vectors, whose real and imaginary parts are conventional 3D vectors. Let ${\boldsymbol {a}}={\boldsymbol {a}}_r+i{\boldsymbol {a}}_i$ and ${\boldsymbol {b}}={\boldsymbol {b}}_r+i{\boldsymbol {b}}_i$ be two such vectors. The dot and cross products are generalised from the case of real vectors as ${\boldsymbol {a}}\boldsymbol {\cdot }{\boldsymbol {b}}={\boldsymbol {a}}_r\boldsymbol {\cdot }{\boldsymbol {b}}_r-{\boldsymbol {a}}_i\boldsymbol {\cdot }{\boldsymbol {b}}_i+i({\boldsymbol {a}}_r\boldsymbol {\cdot }{\boldsymbol {b}}_i+{\boldsymbol {a}}_i\boldsymbol {\cdot }{\boldsymbol {b}}_r)$ and ${\boldsymbol {a}}\times {\boldsymbol {b}}={\boldsymbol {a}}_r\times {\boldsymbol {b}}_r-{\boldsymbol {a}}_i\times {\boldsymbol {b}}_i+i({\boldsymbol {a}}_r\times {\boldsymbol {b}}_i+{\boldsymbol {a}}_i\times {\boldsymbol {b}}_r)$ . An appropriate inner product is $\langle {\boldsymbol {a}},{\boldsymbol {b}}\rangle ={\boldsymbol {a}}^*\boldsymbol {\cdot }{\boldsymbol {b}}$ , which differs from the dot product. The magnitude of a complex vector $\boldsymbol {a}$ is accordingly defined as $|{\boldsymbol {a}}|=\sqrt {{\boldsymbol {a}}^*\boldsymbol {\cdot }{\boldsymbol {a}}}=\sqrt {|{\boldsymbol {a}}_r|^2+|{\boldsymbol {a}}_i|^2}$ .

Linear whirling oscillations formed by the superposition of planar-beating modes, can be represented by a complex-vector amplitude $\tilde{\boldsymbol{{A}}}$ parallel to the wall. Indeed, noting from (3.2) that an eigenfunction ${\boldsymbol{\tilde{r}}}={\hat {\boldsymbol {e}}}_1\varphi$ is associated with perturbations $\tilde {\boldsymbol {t}}={\hat {\boldsymbol {e}}}_1\varphi^{\prime}$ , $\tilde {{\boldsymbol {F}}}=-{\hat {\boldsymbol {e}}}_1(\varphi ^{\prime\prime\prime}+\mathcal {F}_c\varphi^{\prime})$ and $\tilde {{\boldsymbol {M}}}=\boldsymbol {\hat {\imath }}\times {\hat {\boldsymbol {e}}}_1\varphi ^{\prime\prime}$ , a general linear approximation at the threshold $\mathcal {F}=\mathcal {F}_c$ can be expressed as (cf. (3.1))

(3.7a) \begin{align} {\boldsymbol {r}}-s\boldsymbol {\hat {\imath }} & = e^{i\omega t}\,\tilde {\!{\boldsymbol {A}}}\varphi + \text {c.c.}, \end{align}

(3.7b) \begin{align} {\hat {\boldsymbol {t}}}-\boldsymbol {\hat {\imath }} & = e^{i\omega t}\,\tilde {\!{\boldsymbol {A}}}\varphi^{\prime} + \text {c.c.}, \end{align}

(3.7c) \begin{align} {\boldsymbol {F}}+\boldsymbol {\hat {\imath }}\mathcal {F}_c & = -e^{i\omega t}\,\tilde {\!{\boldsymbol {A}}}(\varphi ^{\prime\prime\prime}+\mathcal {F}_c\varphi^{\prime}) + \text {c.c.},\end{align}

(3.7d) \begin{align} {\boldsymbol {M}} & = e^{i\omega t}\boldsymbol {\hat {\imath }}\times \,\tilde {\!{\boldsymbol {A}}}\varphi ^{\prime\prime} + \text {c.c.},\end{align}

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

(3.8) \begin{equation} {\boldsymbol {r}}(1,t) -\boldsymbol {\hat {\imath }} = \,\tilde {\!{\boldsymbol {A}}}e^{i\omega t}+ \text {c.c}=2\,\tilde {\!{\boldsymbol {A}}}_r\cos \omega t-2\,\tilde {\!{\boldsymbol {A}}}_i\sin \omega t, \end{equation}

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

The properties of a general tip ellipse can be derived following Lindell (Reference Lindell1983). The lengths of the semi-major and semi-minor axes can be expressed as

(3.9) \begin{equation} \sqrt {|\,\tilde {\!{\boldsymbol {A}}}|^2+|\,\tilde {\!{\boldsymbol {A}}}\times \,\tilde {\!{\boldsymbol {A}}}^*|}\pm \sqrt {|\,\tilde {\!{\boldsymbol {A}}}|^2-|\,\tilde {\!{\boldsymbol {A}}}\times \,\tilde {\!{\boldsymbol {A}}}^*|}, \end{equation}

respectively, giving the ellipse area as $2\pi |\,\tilde {\!{\boldsymbol {A}}}\times \,\tilde {\!{\boldsymbol {A}}}^*|$ . The above-mentioned condition for a circular orbit, $\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}=0$ , can be inferred from (3.9) by noting the identity $|\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}|^2=|\,\tilde {\!{\boldsymbol {A}}}|^4-|\,\tilde {\!{\boldsymbol {A}}}\times \,\tilde {\!{\boldsymbol {A}}}^*|^2$ . For non-circular ellipses, the real and imaginary parts of $\,\tilde {\!{\boldsymbol {A}}}/(\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}})^{1/2}$ point in the directions of the major and minor axes, respectively (both the forward and backward axis directions, given the multiplicity of the square-root function). These formulae can alternatively be expressed in terms of the real vectors $\,\tilde {\!{\boldsymbol {A}}}_r$ and $\,\tilde {\!{\boldsymbol {A}}}_i$ , using the identities $\,\tilde {\!{\boldsymbol {A}}}\times \,\tilde {\!{\boldsymbol {A}}}^*=-2i\,\tilde {\!{\boldsymbol {A}}}_r\times \,\tilde {\!{\boldsymbol {A}}}_i$ and $\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}=|\,\tilde {\!{\boldsymbol {A}}}_r|^2-|\,\tilde {\!{\boldsymbol {A}}}_i|^2+2i\,\tilde {\!{\boldsymbol {A}}}_r\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}_i$ . Hence, the ellipse is a line if $\,\tilde {\!{\boldsymbol {A}}}_r\times \,\tilde {\!{\boldsymbol {A}}}_i={\mathbf {0}}$ ; it is a circle if both $\,\tilde {\!{\boldsymbol {A}}}_r\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}_i=0$ and $|\,\tilde {\!{\boldsymbol {A}}}_r|=|\,\tilde {\!{\boldsymbol {A}}}_i|$ .

The trajectory $c\,\tilde {\!{\boldsymbol {A}}}e^{i\omega t}+\text {c.c.}$ , wherein $c$ is a complex scalar, is identical to the elliptical tip trajectory (3.8) up to a dilation $|c|$ and phase $\angle c$ . It thus follows from (3.7a ) that, in fact, all points along the filament centreline carry out geometrically similar elliptical trajectories, obtained from the tip trajectory by a dilation $|\varphi (s)|$ and phase shift $\angle \varphi (s)$ . From (3.7 b–d), the tangent unit vector, force and moment perturbations along the filament can be analogously constructed by appropriate scale and phase modulations of the elliptical tip trajectory (as well as a $\pi /2$ rotation in the case of the moment perturbation).

It can be confirmed from the above formulae that in the special case of circular whirling, $\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}=0$ , any given point along the filament performs a steady circular motion parallel to the wall. Furthermore, in that case the geometry of the deformed filament and the force and moment distributions along it appear fixed in a frame of reference rotating about the undeformed configuration at a steady angular frequency $\omega$ .

4.1. Separation of time scales near the instability threshold

According to the linear theory, at the instability threshold the filament approaches periodic beating/whirling oscillations of the form (3.7), with the complex-vector amplitude $\,\tilde {\!{\boldsymbol {A}}}$ , assumed small in magnitude, fixed by the initial conditions. In fact, the fate of small-magnitude oscillations at the instability threshold is actually governed by weak nonlinear effects which become important at late times – invalidating the linear theory! Specifically, nonlinear terms formed of cubic products of the linear approximation (or its derivatives) may resonantly excite the neutral linear modes, thereby producing perturbations growing like $t|\,\tilde {\!{\boldsymbol {A}}}|^3$ . Over $O(1/|\,\tilde {\!{\boldsymbol {A}}}|^{2})$ times – long compared with the natural period $2\pi /\omega$ – such nonlinear perturbations would become comparable to the linear approximation and may accordingly influence the dynamics.

On the other hand, the linear theory predicts that small-magnitude perturbations exponentially decay for $\mathcal {F}\lt \mathcal {F}_c$ , or exponentially grow, in general, for $\mathcal {F}\gt \mathcal {F}_c$ . As $\mathcal {F}\to \mathcal {F}_c$ , the growth rate $\lambda _r$ of the most unstable linear modes vanishes like $\mathcal {F}-\mathcal {F}_c$ (see figure 2 a), suggesting that following a transient (during which more stable modes decay) the filament exhibits oscillations of the form (3.7) – modulated over long $O(1/|\mathcal {F}-\mathcal {F}_c|)$ times. We accordingly identify the distinguished scaling $|\,\tilde {\!{\boldsymbol {A}}}|=\text {ord}(|\mathcal {F}-\mathcal {F}_c|^{1/2})$ – conventional to symmetry-breaking bifurcations – such that linear and nonlinear effects both have a leading-order effect. For relatively larger amplitudes, or at the threshold, nonlinear effects dominate the dynamics; for relatively smaller amplitudes, linear growth or decay is instead dominant.

The time-scale separation at and near the instability threshold can be exploited towards systematically deriving a reduced-order nonlinear model governing the dynamics of small-magnitude oscillations in this regime. This is carried out in § 5 by means of a weakly nonlinear analysis of the governing equations employing the method of multiple scales; while conceptually standard, the derivation is presented fully so as to highlight some technical details peculiar to the present set-up and, more importantly, to facilitate future extensions to the theory as discussed in § 6. In the present section, we simply quote the reduced-order model and subsequently analyse and illustrate it in various scenarios.

4.2. Amplitude equation

The reduced-order model consists of the linear approximation (3.7), now understood to hold in a vicinity of the instability threshold and with $\,\tilde {\!{\boldsymbol {A}}}$ evolving with time according to the amplitude equation

(4.1) \begin{equation} \frac {\mathrm {d}\,\tilde {\!{\boldsymbol {A}}}}{\mathrm {d}t}=\alpha \,\tilde {\!{\boldsymbol {A}}}^*\,\tilde {\!{\boldsymbol {A}}}\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}+\beta \,\tilde {\!{\boldsymbol {A}}}{\,\tilde {\!{\boldsymbol {A}}}}^*\boldsymbol {\cdot }\,\tilde {\!{\boldsymbol {A}}}+(\mathcal {F}-\mathcal {F}_c)\gamma \,\tilde {\!{\boldsymbol {A}}}, \end{equation}

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

(4.2a–c) \begin{equation} \alpha \approx -451.038 -i327.599,\quad \beta \approx -517.974+i353.952, \quad \gamma \approx 7.34480+i5.34285. \end{equation}

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

We prefer to work with rescaled, order-unity quantities. Thus, consistently with the weakly nonlinear analysis in § 5, we define

(4.3) \begin{equation} \mathcal {F}=\mathcal {F}_c+\varepsilon \chi , \end{equation}

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

(4.4) \begin{equation} \,\tilde {\!{\boldsymbol {A}}}(t)=\varepsilon ^{1/2}{\boldsymbol {A}}(T), \end{equation}

wherein

(4.5) \begin{equation} T=\varepsilon t \end{equation}

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

(4.6) \begin{equation} \frac {\mathrm {d}{{\boldsymbol {A}}}}{\mathrm {d}T}=\alpha {{\boldsymbol {A}}}^*{{\boldsymbol {A}}}\boldsymbol {\cdot }{{\boldsymbol {A}}}+\beta {{\boldsymbol {A}}}{{{\boldsymbol {A}}}}^*\boldsymbol {\cdot }{{\boldsymbol {A}}}+\chi \gamma {{\boldsymbol {A}}}. \end{equation}

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

It is useful to note the following three relations:

(4.7a) \begin{align} \frac {\mathrm {d}}{\mathrm {d}T}({\boldsymbol {A}}^*\times {\boldsymbol {A}}) & =2(\beta _r|{\boldsymbol {A}}|^2+\chi \gamma _r){\boldsymbol {A}}^*\times {\boldsymbol {A}},\end{align}

(4.7b) \begin{align} \frac {\mathrm {d}}{\mathrm {d}T}({\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}})& = 2\left [(\alpha +\beta )|{\boldsymbol {A}}|^2+\chi \gamma \right ]{\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}},\end{align}

(4.7c) \begin{align} \frac {\mathrm {d}}{\mathrm {d}T}|{\boldsymbol {A}}|^2 & = 2(\alpha _r+\beta _r)|{\boldsymbol {A}}|^4-2\alpha _r|{\boldsymbol {A}}^*\times {\boldsymbol {A}}|^2+2\chi \gamma _r|{\boldsymbol {A}}|^2.\end{align}

They readily follow from (4.6): we obtain (4.7a ) by subtracting the cross product of ${\boldsymbol {A}}^*$ and (4.6) from the conjugate of that product; we obtain (4.7b ) by considering the dot product of $\boldsymbol {A}$ and (4.6), noting the identity $2\boldsymbol{{A}} \boldsymbol{\cdot} \mathrm{d}\boldsymbol{{A}}/\mathrm{d} T = \mathrm{d}(\boldsymbol{{A}} \boldsymbol{\cdot} \boldsymbol{{A}})/\mathrm{d} T$ ; and we obtain (4.7c ) by adding the dot product of ${\boldsymbol {A}}^*$ and (4.6) to the conjugate of that product. Equation (4.7a ) shows that if the fast-scale tip orbit is initially a line (corresponding to planar motion of the filament), then it remains so at all later slow times. Similarly, (4.7b ) shows that if the fast-scale tip orbit is initially circular then it remains so at all later slow times.

If the filament is confined to a plane (by the initial conditions), the complex-vector amplitude can be written ${\boldsymbol {A}}(T)={\hat {\boldsymbol {e}}} A(T)$ , where $A(T)$ is a complex-scalar amplitude and $\hat {\boldsymbol {e}}$ a unit vector parallel to the wall; the filament in that case performs planar-beating oscillations with $A(T)$ determining their magnitude and phase. In that restricted scenario, the amplitude equation (4.6) reduces to

(4.8) \begin{equation} \frac {\mathrm {d}A}{\mathrm {d}T}=(\alpha +\beta )|A|^2A+\chi \gamma A, \end{equation}

which has the form of the Stuart–Landau amplitude equation corresponding to a conventional (non-degenerate) Hopf instability (Drazin & Reid Reference Drazin and Reid2004).

We proceed in § 4.3 to find solutions of the amplitude equation (4.6) representing steady and periodic states; we also compare these solutions against numerical data provided by Eric E. Keaveny. In §§ 4.4–4.6, we analyse the stability of these states within the framework of the weakly nonlinear theory. In § 4.7, we present simulations of (4.6).

4.3. Periodic states

We consider ‘quasi-steady’ solutions of the amplitude equation (4.6) in the form

(4.9) \begin{equation} {\boldsymbol {A}}=e^{i\nu T}{\boldsymbol{\bar{A}}}, \end{equation}

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

(4.10) \begin{equation} \alpha {\boldsymbol{\bar{A}}}^*{\boldsymbol{\bar{A}}}\boldsymbol {\cdot }{\boldsymbol{\bar{A}}}+\beta {\boldsymbol{\bar{A}}}{\boldsymbol{\bar{A}}}^*\boldsymbol {\cdot } {\boldsymbol{\bar{A}}}+(\chi \gamma -i\nu ){\boldsymbol{\bar{A}}}={\mathbf {0}}. \end{equation}

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

4.3.1. Planar-beating states

In the planar-beating case, we can write

(4.11) \begin{equation} {\boldsymbol{\bar{A}}}={ {\bar{A}}}{\hat {\boldsymbol {e}}}, \end{equation}

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

(4.12) \begin{equation} { {\bar{A}}}=e^{i\vartheta }|{{\bar{A}}}|, \end{equation}

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

(4.13) \begin{equation} (\alpha +\beta )|{{\bar{A}}}|^2+\chi \gamma -i\nu =0, \end{equation}

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

(4.14a,b) \begin{equation} |{{\bar{A}}}|^2=-\frac {\chi \gamma _r}{\alpha _r+\beta _r}, \quad \nu =\chi \gamma _i-\chi \gamma _r\frac {\alpha _i+\beta _i}{\alpha _r+\beta _r}. \end{equation}

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

(4.15) \begin{equation} {\boldsymbol {r}}(1,t) -\boldsymbol {\hat {\imath }} = 2|{{\bar{A}}}|(\mathcal {F}-\mathcal {F}_c)^{1/2} \cos [(\omega +|\mathcal {F}-\mathcal {F}_c| \nu ) t+\vartheta ]{\hat {\boldsymbol {e}}}, \end{equation}

where the rescaled tip displacement $2|{{\bar{A}}}|\doteq 0.17412$ and angular-frequency correction $\nu \doteq 5.5426$ are calculated using (4.2) and (4.14).

4.3.2. Circular-whirling states

In the circular-whirling case, ${\boldsymbol{\bar{A}}}\boldsymbol {\cdot }{\boldsymbol{\bar{A}}}=0$ , the real vectors ${\boldsymbol {A}}_r$ and ${\boldsymbol {A}}_i$ are orthogonal and of equal magnitude. We accordingly write

(4.16) \begin{equation} {\boldsymbol{\bar{A}}}=\frac {1}{\sqrt {2}}|{\boldsymbol{\bar{A}}}|({\hat {\boldsymbol {e}}} \pm i\boldsymbol {\hat {\imath }}\times {\hat {\boldsymbol {e}}}), \end{equation}

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

(4.17) \begin{equation} \beta |{\boldsymbol{\bar{A}}}|^2+\gamma \chi -i\nu =0, \end{equation}

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

(4.18a,b) \begin{equation} |{\boldsymbol{\bar{A}}}|^2 = -\frac {\chi \gamma _r}{\beta _r}, \quad \nu = \chi \gamma _i-\chi \gamma _r\frac {\beta _i}{\beta _r}. \end{equation}

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

(4.19) \begin{align} {\boldsymbol {r}}(1,t) -\boldsymbol {\hat {\imath }} & = \sqrt {2}|{\boldsymbol{\bar{A}}}|(\mathcal {F}-\mathcal {F}_c)^{1/2}\left \{{\hat {\boldsymbol {e}}}\cos \left [(\omega +|\mathcal {F}-\mathcal {F}_c| \nu ) t\right ] \right . \nonumber\\ & \quad \left . \mp \boldsymbol {\hat {\imath }}\times {\hat {\boldsymbol {e}}} \sin \left [(\omega +|\mathcal {F}-\mathcal {F}_c|\nu )t\right ] \right \}, \end{align}

where the rescaled tip radius $\sqrt {2}|{\boldsymbol{\bar{A}}}|\doteq 0.16840$ and angular-frequency correction $\nu \doteq 10.362$ are calculated using (4.2) and (4.18).

4.3.3. Bifurcation diagram and comparison with numerical solutions of the full problem

In figure 3, we plot maximal tip displacement versus dimensionless follower force for the steady and periodic states predicted by the weakly nonlinear theory. As evident in (4.15) and (4.19), the maximal tip displacement $\propto \sqrt {\mathcal {F}-\mathcal {F}_c}$ for both the planar-beating and circular-whirling periodic states, the constant of proportionality being slightly larger for the planar-beating states. (We note that this square-root scaling near the threshold was hypothesised as well as numerically demonstrated by Clarke et al. (Reference Clarke, Hwang and Keaveny2024).) The stability characteristics indicated in the bifurcation diagram are derived and discussed in the following subsections.

Figure 3.

(a) Bifurcation diagram showing maximal dip displacement (normalised by filament length) against the dimensionless follower-force magnitude $\mathcal {F}$ , for the (i) undeformed steady state; (ii) planar-beating periodic states; and (iii) circular-whirling periodic states. The curves depict the predictions of the weakly nonlinear theory (see § 4): the undeformed state is stable preceding the bifurcation (solid black line) and unstable following it (dashed black line); the planar-beating states are longitudinally stable, up to phase shifts, but transversally unstable (dash-dotted red curve); and the circular-whirling states are stable up to phase shifts (solid blue curve). The symbols depict numerical data provided by Eric E. Keaveny (see § 4.3.3). (b) Zoomed-in comparison between the theory and numerical data depicted on a log–log scale.

With the purpose of validating the present weakly nonlinear theory, Eric E. Keaveny has conducted simulations based on the numerical methodology of Schoeller et al. (Reference Schoeller, Townsend, Westwood and Keaveny2021), adopting the same formulation as Clarke et al. (Reference Clarke, Hwang and Keaveny2024) only with the full Stokes hydrodynamics included in their study simplified to resistive force theory – consistently with the present formulation (cf. (2.2)). Figure 3 shows remarkable agreement, in terms of the maximal tip displacements, between these simulations and the theoretical predictions.

4.4. Linear stability of the undeformed state

The stability characteristics of the undeformed state, as well as the planar-beating and circular-whirling periodic states, are readily determined within the framework of the weakly nonlinear theory. We begin by revisiting the stability of the undeformed state, which has already been analysed exactly in § 3. To this end, we consider small perturbations of the complex-vector amplitude:

(4.20) \begin{equation} {\boldsymbol {A}}(T)={\boldsymbol {a}}(T), \quad |{\boldsymbol {a}}|\ll 1. \end{equation}

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

(4.21) \begin{equation} \frac {\mathrm {d}{\boldsymbol {a}}}{\mathrm {d}T}=\chi \gamma {\boldsymbol {a}}, \end{equation}

which possesses the general solution

(4.22) \begin{equation} {\boldsymbol {a}}(T)=e^{\chi \gamma T}{\boldsymbol {a}}(0). \end{equation}

In this linear approximation, the elliptical tip orbit at slow time $T$ is obtained from the initial orbit defined by ${\boldsymbol {a}}(0)$ via a dilation by $\exp (\chi \gamma _r T)$ . Since $\gamma _r\gt 0$ , perturbations grow for $\chi \gt 0$ and decay for $\chi \lt 0$ . The case $\chi =0$ , where nonlinear terms determine the stability of the undeformed state, will be considered in § 4.6.

It is instructive to repeat the above stability analysis, this time starting from the real dynamical system implied by the linearised amplitude equation (4.21):

(4.23) \begin{equation} \frac {\mathrm {d}}{\mathrm {d}T}\left (\begin{array}{c}{\boldsymbol {a}}_r \\ {\boldsymbol {a}}_i \end{array}\right ) = \chi \left (\begin{array}{c@{\quad}c} \gamma _r & -\gamma _i \\ \gamma _i & \gamma _r\end{array}\right )\left (\begin{array}{c}{\boldsymbol {a}}_r \\ {\boldsymbol {a}}_i \end{array}\right ). \end{equation}

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

(4.24) \begin{equation} \left (\begin{array}{c}{\boldsymbol {a}}_r \\{\boldsymbol {a}}_i \end{array}\right )={\rm e}^{\sigma T}\left (\begin{array}{c}\tilde {{\boldsymbol {a}}}_r \\ \tilde {{\boldsymbol {a}}}_i \end{array}\right )+\text {c.c.}, \end{equation}

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

(4.25) \begin{equation} \left (\begin{array}{c@{\quad}c} \chi \gamma _r-\sigma & -\chi \gamma _i \\\chi \gamma _i & \chi \gamma _r-\sigma \end{array}\right )\left (\begin{array}{c}\tilde {{\boldsymbol {a}}}_r \\ \tilde {{\boldsymbol {a}}}_i\end{array}\right )=\left (\begin{array}{c}{\mathbf {0}} \\ {\mathbf {0}}\end{array}\right ). \end{equation}

The solution space is spanned by the eigenvectors

(4.26a,b) \begin{equation} \left (\begin{array}{c}\tilde {{\boldsymbol {a}}}_r \\ \tilde {{\boldsymbol {a}}}_i \end{array}\right )^{({1,2})}=\left (\begin{array}{c}{\hat {\boldsymbol {e}}}_{1} \\ \mp i{\hat {\boldsymbol {e}}}_{1} \end{array}\right ), \quad \left (\begin{array}{c}\tilde {{\boldsymbol {a}}}_r \\ \tilde {{\boldsymbol {a}}}_i \end{array}\right )^{({3,4})}=\left (\begin{array}{c}{\hat {\boldsymbol {e}}}_{2} \\ \mp {\hat {\boldsymbol {e}}}_{2} \end{array}\right ), \end{equation}

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

(4.27a,b) \begin{equation} \sigma ^{({1,3})}=\chi \gamma , \quad \sigma ^{({2,4})}=\chi \gamma ^*. \end{equation}

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

4.5. Linear stability of the periodic states

We next consider the stability of the planar-beating and circular-whirling states. In these cases, it is convenient to write the perturbation as

(4.28) \begin{equation} {\rm e}^{-i\nu T}{\boldsymbol {A}}(T) = {\boldsymbol{\bar{A}}}+{\boldsymbol {a}}(T), \quad |{\boldsymbol {a}}|\ll 1, \end{equation}

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

(4.29) \begin{equation} \frac {\mathrm {d}{\boldsymbol {a}}}{\mathrm {d}T}=2\alpha {\boldsymbol{\bar{A}}}^*{\boldsymbol{\bar{A}}}\boldsymbol {\cdot } {\boldsymbol {a}}+\alpha {\boldsymbol {a}}^*{\boldsymbol{\bar{A}}}\boldsymbol {\cdot } {\boldsymbol{\bar{A}}}+\beta {\boldsymbol{\bar{A}}}{\boldsymbol{\bar{A}}}\boldsymbol {\cdot } {\boldsymbol {a}}^*+\beta {\boldsymbol{\bar{A}}}{\boldsymbol {a}}\boldsymbol {\cdot }{\boldsymbol{\bar{A}}}^*+(\gamma -i\nu +\beta |{\boldsymbol{\bar{A}}}|^2){\boldsymbol {a}}. \end{equation}

4.5.1. Linear stability of the planar-beating states

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

(4.30a,b) \begin{equation} \frac {\mathrm {d}a_{\parallel }}{\mathrm {d}T} = (\alpha +\beta )|{{\bar{A}}}|^2\left (a_{\parallel }+{\rm e}^{2i\vartheta }a_{\parallel }^*\right ),\quad \frac {\mathrm {d}a_{\perp }}{\mathrm {d}T}=\alpha |{{\bar{A}}}|^2\left ({\rm e}^{2i\vartheta }a_{\perp }^*-a_{\perp }\right ), \end{equation}

wherein $|{{\bar{A}}}|$ is provided by (4.14a ).

We first consider parallel perturbations. Separating (4.30a ) into its real and imaginary parts, we find the two-dimensional real system

(4.31) \begin{equation} \dfrac {\mathrm {d}}{\mathrm {d}T}\left (\begin{array}{c}{a_{\parallel }}_r\\{a_{\parallel }}_i\end{array}\right )=|{{\bar{A}}}|^2\left (\begin{array}{c@{\quad}c}\zeta _r(1+\cos 2\vartheta )-\zeta _i\sin 2\vartheta & \zeta _r\sin 2\vartheta -\zeta _i(1-\cos 2\vartheta ) \\\zeta _r\sin 2\vartheta +\zeta _i(1+\cos 2\vartheta )&\zeta _r(1-\cos 2\vartheta )+\zeta _i\sin 2\vartheta \end{array}\right )\left (\begin{array}{c}{a_{\parallel }}_r\\{a_{\parallel }}_i\end{array}\right ), \end{equation}

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

(4.32a,b) \begin{equation} \sigma _{\parallel }^{({1})} = 0, \quad \sigma _{\parallel }^{({2})}=-2\gamma _r, \end{equation}

where we have used (4.14a ) for $|{{\bar{A}}}|$ . The zero eigenvalue $\sigma _{\parallel }^{({1})}$ reflects the arbitrariness of the phase $\vartheta$ . Since $\gamma _r\gt 0$ , the eigenvalue $\sigma _{\parallel }^{({2})}$ has a negative real part. Thus, planar beating is stable to parallel perturbations apart from the expected neutrality with respect to phase.

We next consider perpendicular perturbations. Separating (4.30b ) into its real and imaginary parts, we find the two-dimensional real system

(4.33) \begin{align} & \dfrac {\mathrm {d}}{\mathrm {d}T}\left (\begin{array}{c}{a_{\perp }}_r\\ {a_{\perp }}_i\end{array}\right )\nonumber\\ & =|{{\bar{A}}}|^2\left (\begin{array}{c@{\quad}c}\alpha _r(\cos 2\vartheta -1)-\alpha _i\sin 2\vartheta & \alpha _r\sin 2\vartheta +\alpha _i(1+\cos 2\vartheta ) \\\alpha _r\sin 2\vartheta +\alpha _i(\cos 2\vartheta -1) & -\alpha _r(1+\cos 2\vartheta )+\alpha _i\sin 2\vartheta \end{array}\right )\left (\begin{array}{c}{a_{\perp }}_r\\ {a_{\perp }}_i\end{array}\right ). \end{align}

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

(4.34a,b) \begin{equation} \sigma _{\perp }^{({1})} = 0, \quad \sigma _{\perp }^{({2})}=\frac {2\alpha _r\gamma _r}{\alpha _r+\beta _r}, \end{equation}

where we have used (4.14a ) for $|{{\bar{A}}}|$ . The zero eigenvalue $\sigma _{\perp }^{({1})}$ reflects the arbitrariness of the direction $\hat {\boldsymbol {e}}$ . Since $\gamma _r\gt 0$ , whereas $\alpha _r,\beta _r\lt 0$ , the eigenvalue $\sigma _{\perp }^{({2})}$ has a positive real part. Thus, planar beating is unstable under perpendicular perturbations – and therefore unstable under general perturbations.

4.5.2. Stability of the circling-whirling states

For the circular-whirling states (§ 4.3.2), the form of ${\boldsymbol{\bar{A}}}$ is given by (4.16), with the plus/minus sign determining the direction of motion and the unit vector $\hat {\boldsymbol {e}}$ representing a phase; without loss of generality, we only consider the stability of circular whirling in the clockwise direction, corresponding to choosing the plus sign. Writing ${\boldsymbol {a}}=a_1{\hat {\boldsymbol {e}}}+a_2\boldsymbol {\hat {\imath }}\times {\hat {\boldsymbol {e}}}$ , the linearised amplitude equation (4.29) gives, upon eliminating $\nu$ using the steady relation (4.17), the coupled pair of equations

(4.35a) \begin{eqnarray} \frac {\mathrm {d}a_1}{\mathrm {d}T}=\alpha |{\boldsymbol{\bar{A}}}|^2(a_1+ia_2)+\frac {1}{2}\beta |{\boldsymbol{\bar{A}}}|^2(a_1+a_1^*+ia_2^*-ia_2), \end{eqnarray}

(4.35b) \begin{eqnarray} \frac {\mathrm {d}a_2}{\mathrm {d}T}=\alpha |{\boldsymbol{\bar{A}}}|^2(a_2-ia_1)+\frac {1}{2}\beta |{\boldsymbol{\bar{A}}}|^2(a_2-a_2^*+ia_1+ia_1^*), \end{eqnarray}

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

(4.36) \begin{equation} \frac {\mathrm {d}}{\mathrm {d}T}\left (\begin{array}{c}{a_1}_r \\{a_1}_i \\{a_2}_r \\{a_2}_i\end{array}\right )=|{\boldsymbol{\bar{A}}}|^2\left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \alpha _r+\beta _r & -\alpha _i & -\alpha _i & \beta _r-\alpha _r \\\alpha _i+\beta _i & \alpha _r & \alpha _r & \beta _i-\alpha _i \\ \alpha _i-\beta _i & \alpha _r & \alpha _r & -(\alpha _i+\beta _i) \\ \beta _r-\alpha _r & \alpha _i & \alpha _i & \alpha _r+\beta _r\end{array}\right )\left (\begin{array}{c}{a_1}_r \\{a_1}_i \\{a_2}_r \\{a_2}_i\end{array}\right ). \end{equation}

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

(4.37) \begin{equation} \sigma (\sigma -2|{\boldsymbol{\bar{A}}}|^2\beta _r)\left (\sigma ^2-4\alpha _r|{\boldsymbol{\bar{A}}}|^2\sigma +4|{\boldsymbol{\bar{A}}}|^4|\alpha |^2\right )=0, \end{equation}

whose solutions yield the growth-rate eigenvalues

(4.38a-d) \begin{equation} \sigma ^{({1})}=0, \quad \sigma ^{({2})}=-2\gamma _r, \quad \sigma ^{({3})}= - 2\frac {\gamma _r}{\beta _r}\alpha , \quad \sigma ^{({4})}= - 2\frac {\gamma _r}{\beta _r}\alpha ^*, \end{equation}

where we have used (4.18a ) for $|{\boldsymbol{\bar{A}}}|$ . The zero eigenvalue $\sigma ^{({1})}$ reflects the arbitrariness of the phase unit vector $\hat {\boldsymbol {e}}$ . Since $\gamma _r\gt 0$ , whereas $\alpha _r,\beta _r\lt 0$ , the remaining eigenvalues have negative real parts. Thus, the circular-whirling states are stable apart from the expected neutrality with respect to phase.

The predicted stability characteristics of the different states are indicated in the bifurcation diagram (figure 3 a). The stability of the circular-whirling states and the instability of the planar-beating states (under general three-dimensional perturbations) are in accordance with the Floquet analyses numerically performed by Clarke et al. (Reference Clarke, Hwang and Keaveny2024). As first established by Clarke et al. (Reference Clarke, Hwang and Keaveny2024), these stability characteristics serve to rationalise the selection of whirling over beating in initial-value simulations not restricted to planar deformations (Ling et al. Reference Ling, Guo and Kanso2018).

4.6. Nonlinear stability of the undeformed state

It remains to determine the nonlinear stability of the undeformed state at the linear-instability threshold. Setting $\chi =0$ , the amplitude equation (4.6) reduces to

(4.39) \begin{equation} \frac {\mathrm {d}{\boldsymbol {A}}}{\mathrm {d}T}=\alpha {\boldsymbol {A}}^*{\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}}+\beta {\boldsymbol {A}}{\boldsymbol {A}}^*\boldsymbol {\cdot }{\boldsymbol {A}}. \end{equation}

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

(4.40a,b) \begin{equation} \frac {\mathrm {d}}{\mathrm {d}T}|{\boldsymbol {A}}|=(\alpha _r+\beta _r)|{\boldsymbol {A}}|^3-\alpha _r|{\boldsymbol {B}}|^2|{\boldsymbol {A}}|^{-1}, \quad \frac {\mathrm {d}}{\mathrm {d}T}|{\boldsymbol {B}}|=2\beta _r|{\boldsymbol {A}}|^2|{\boldsymbol {B}}|, \end{equation}

where we denote ${\boldsymbol {B}}={\boldsymbol {A}}^*\times {\boldsymbol {A}}$ .

In the special case of planar perturbations, where $\boldsymbol {B}$ vanishes, (4.40a ) reduces to

(4.41) \begin{equation} \frac {\mathrm {d}|{\boldsymbol {A}}|}{\mathrm {d}T}=(\alpha _r+\beta _r)|{\boldsymbol {A}}|^3. \end{equation}

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

(4.42) \begin{equation} \left |{\boldsymbol {A}}(T)\right |= \frac {|{\boldsymbol {A}}(0)|}{\sqrt {1-2|{\boldsymbol {A}}(0)|^{2}(\alpha _r+\beta _r)T}}\sim \frac {1}{\sqrt {-2(\alpha _r+\beta _r)T}}\quad \text {as} \quad T\to \infty , \end{equation}

revealing algebraic decay with the $-1/2$ power of time.

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

(4.43a,b) \begin{equation} |{\boldsymbol {A}}(T)|\sim \frac {1}{\sqrt {-2\beta _rT}}, \quad |{\boldsymbol {B}}(T)|\sim -\frac {1}{2\beta _rT} \quad \text {at} \quad T\to \infty . \end{equation}

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

(4.44) \begin{equation} \sqrt {|{\boldsymbol {A}}|^2+|{\boldsymbol {B}}|}\pm \sqrt {|{\boldsymbol {A}}|^2-|{\boldsymbol {B}}|} = \frac {1}{\sqrt {-\beta _rT}} + o(T^{-1/2}) \quad \text {as} \quad T\to \infty . \end{equation}

Thus, at late times the filament exhibits circular whirling on the fast scale, with the tip radius decaying algebraically with the $-1/2$ power of time.

4.7. Numerical illustrations

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

Figure 4.

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

Figure 5.

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

Figure 6.

Algebraic decay of perturbations at the linear-instability threshold $\mathcal {F}=\mathcal {F}_c$ (see § 4.6). The plot compares $|{\boldsymbol {A}}(T)|$ in the planar scenario ${\boldsymbol {A}}(0)={\hat {\boldsymbol {e}}}_x$ (thick solid curve) and non-planar scenario ${\boldsymbol {A}}(0)=({\hat {\boldsymbol {e}}}_x+i\sqrt {3}{\hat {\boldsymbol {e}}}_y)/2$ (thick dash-dotted curve) with the power-law asymptotes (4.42) and (4.43a ), where $\{{\hat {\boldsymbol {e}}}_x,{\hat {\boldsymbol {e}}}_y,\boldsymbol {\hat {\imath }}\}$ is a right-handed system of unit vectors. For the non-planar scenario, we also depict the magnitude $|{\boldsymbol {B}}|=|{\boldsymbol {A}}^*\times {\boldsymbol {A}}|$ (thick dashed curve) and its power-law asymptote (4.43b ), and in the inset a top view (looking towards the wall, with ${\hat {\boldsymbol {e}}}_x$ pointing to the right) of the tip trajectory scaled by $\varepsilon ^{1/2}$ , with $\varepsilon =0.5$ .

In figure 4, we consider the case where the initial fast-time oscillation is planar. In light of (4.7a ), the motion of the filament remains planar at all times. We observe the multiple-scales evolution of the dynamics towards the periodic planar-beating states (see §§ 4.3.1 and 4.5.1), for $\mathcal {F}\gt \mathcal {F}_c$ , or the undeformed configuration, for $\mathcal {F}\lt \mathcal {F}_c$ , with faster transients for larger $|\mathcal {F}-\mathcal {F}_c|$ . In figure 5, we consider an example where the initial fast-time oscillation is nearly planar and $\mathcal {F}\gt \mathcal {F}_c$ . We observe the slow-time expansion and rotation of the fast-scale elliptical orbit of the tip as it approaches the circular orbit corresponding to the periodic circular-whirling states (see §§ 4.3.2 and 4.5.2). Lastly, in figure 6, we demonstrate the algebraic decay of both planar and non-planar perturbations at the instability threshold, $\mathcal {F}=\mathcal {F}_c$ (see § 4.6).

5.1. Multiple-scale expansions near the instability threshold

In this section, we derive the amplitude equation (4.6), which together with the general representation (3.7) for the linear approximation at the instability threshold constitutes the weakly nonlinear theory that we have presented, analysed and illustrated in § 4. The derivation consists of a weakly nonlinear analysis in the multiple-scales near-threshold regime identified in § 4.1.

We write $\mathcal {F}=\mathcal {F}_c+\varepsilon \chi$ for the dimensionless follower force, as in (4.3), and consider the long-time dynamics of small-magnitude oscillations for $0\lt \varepsilon \ll 1$ , holding $\chi$ fixed. Following the method of multiple scales (Hinch Reference Hinch1991), we introduce the two-scale extension $\underline {{\boldsymbol {r}}}(s,\tau ,T)$ of the centreline position vector ${{\boldsymbol {r}}}(s,t)$ , wherein the ‘fast time’ $\tau$ and ‘slow time’ $T$ are treated as independent coordinates, such that $\underline {{\boldsymbol {r}}}(s,\tau ,T)={\boldsymbol {r}}(s,t)$ on the ‘physical diagonal’:

(5.1) \begin{equation} \tau =t \quad \text {and} \quad T=\varepsilon t; \end{equation}

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

(5.2) \begin{equation} \frac {\partial }{\partial t}\Rightarrow \frac {\partial }{\partial \tau }+\varepsilon \frac {\partial }{\partial T}. \end{equation}

We posit the ‘weakly nonlinear expansion’

(5.3) \begin{equation} \underline {{\boldsymbol {r}}}(s,\tau ,T)={{\boldsymbol {r}}}_0(s)+\varepsilon ^{1/2}{{\boldsymbol {r}}}_{1/2}(s,\tau ,T)+\varepsilon {{\boldsymbol {r}}}_1(s,\tau ,T)+\varepsilon ^{3/2}{\boldsymbol {r}}_{3/2}(s,\tau ,T)+\cdots , \end{equation}

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

(5.4) \begin{equation} {\boldsymbol {r}}_0=s\boldsymbol {\hat {\imath }}, \quad {\hat {\boldsymbol {t}}}_0 = \boldsymbol {\hat {\imath }}, \quad {\boldsymbol {F}}_0=-\mathcal {F}_c\boldsymbol {\hat {\imath }}, \quad {\boldsymbol {M}}_0={\mathbf {0}}. \end{equation}

A key idea underlying the method is to exploit the added freedom associated with the extension to two time scales in order to ensure that the weakly nonlinear expansions remain asymptotically ordered. We note that the scaling of time by $1/\varepsilon$ and the expansion in half-powers of $\varepsilon$ are suggested by the scaling arguments given in § 4.1.

5.2. The $O(\varepsilon ^{1/2})$ problem

At $O(\varepsilon ^{1/2})$ , the governing partial differential equations (2.6), together with the time-derivative transformation (5.2), give

(5.5a) \begin{equation} \frac {\partial {\boldsymbol {r}}_{1/2}}{\partial s}-\boldsymbol {t}_{1/2}={\mathbf {0}}, \end{equation}

(5.5b) \begin{equation} \frac {\partial {\boldsymbol {F}}_{1/2}}{\partial s}-\left (\mathsf {I}-\frac {1}{2}\boldsymbol {\hat {\imath }}\boldsymbol {\hat {\imath }}\right )\boldsymbol {\cdot } \frac {\partial {\boldsymbol {r}}_{1/2}}{\partial \tau }={\mathbf {0}}, \end{equation}

(5.5c) \begin{equation} \frac {\partial {\boldsymbol {M}}_{1/2}}{\partial s}+\mathcal {F}_c\boldsymbol {\hat {\imath }}\times \boldsymbol {t}_{1/2}+\boldsymbol {\hat {\imath }}\times {\boldsymbol {F}}_{1/2}={\mathbf {0}}, \end{equation}

(5.5d) \begin{equation} {\boldsymbol {M}}_{1/2}-\boldsymbol {\hat {\imath }}\times \frac {\partial \boldsymbol {t}_{1/2}}{\partial s}={\mathbf {0}}; \end{equation}

the inextensibility constraint (2.7) gives

(5.6) \begin{equation} \boldsymbol {\hat {\imath }}\boldsymbol {\cdot }\boldsymbol {t}_{1/2}=0; \end{equation}

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

(5.7a,b) \begin{equation} {\boldsymbol {r}}_{1/2}={\mathbf {0}}, \quad \boldsymbol {t}_{1/2}={\mathbf {0}} \quad \text {at} \quad s=0 \end{equation}

and

(5.8a,b) \begin{equation} {\boldsymbol {F}}_{1/2}+\mathcal {F}_c\boldsymbol {t}_{1/2}={\mathbf {0}}, \quad {\boldsymbol {M}}_{1/2}={\mathbf {0}} \quad \text {at} \quad s=1. \end{equation}

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

(5.9a) \begin{align} {\boldsymbol {r}}_{1/2} &= {\boldsymbol {A}}(T)\varphi (s){\rm e}^{i\omega \tau }+\text {c.c.},\end{align}

(5.9b) \begin{align} \boldsymbol {t}_{1/2} &= {\boldsymbol {A}}(T)\varphi ^{\prime}(s){\rm e}^{i\omega \tau }+\text {c.c.}, \end{align}

(5.9c) \begin{align} {\boldsymbol {F}}_{1/2}&= -{\boldsymbol {A}}(T)[\varphi ^{\prime\prime\prime}(s)+\mathcal {F}_c\varphi ^{\prime}(s)]{\rm e}^{i\omega \tau }+\text {c.c.},\end{align}

(5.9d) \begin{align}{\boldsymbol {M}}_{1/2}&= \boldsymbol {\hat {\imath }}\times {\boldsymbol {A}}(T)\varphi ^{\prime\prime}(s){\rm e}^{i\omega \tau }+\text {c.c.},\end{align}

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

5.3. The $O(\varepsilon )$ problem

At $O(\varepsilon )$ , the governing partial differential equations (2.6), together with the time-derivative transformation (5.2), give

(5.10a) \begin{align} \frac {\partial {\boldsymbol {r}}_1}{\partial s}-\boldsymbol {t}_1={\mathbf {0}}, \end{align}

(5.10b) \begin{align} \frac {\partial {\boldsymbol {F}}_1}{\partial s}-\left (\mathsf {I}-\frac {1}{2}\boldsymbol {\hat {\imath }}\boldsymbol {\hat {\imath }}\right )\boldsymbol {\cdot } \frac {\partial {\boldsymbol {r}}_1}{\partial \tau } = -\frac {1}{2}\left (\boldsymbol {t}_{1/2}\boldsymbol {\hat {\imath }}+\boldsymbol {\hat {\imath }}\boldsymbol {t}_{1/2}\right )\boldsymbol {\cdot }\frac {\partial {\boldsymbol {r}}_{1/2}}{\partial \tau }\\ \left \{ =-\boldsymbol {\hat {\imath }}\frac {i\omega }{2}\left ({\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}}\varphi ^{\prime}\varphi {\rm e}^{2i\omega \tau }-|{\boldsymbol {A}}|^2\varphi ^{\prime}\varphi ^*\right )+\text {c.c.}\right \}, \end{align}

(5.10c) \begin{align} \frac {\partial {\boldsymbol {M}}_1}{\partial s}+\mathcal {F}_c\boldsymbol {\hat {\imath }}\times \boldsymbol {t}_1+\boldsymbol {\hat {\imath }}\times {\boldsymbol {F}}_1 = -\boldsymbol {t}_{1/2}\times {\boldsymbol {F}}_{1/2} \quad \left \{ = {\boldsymbol {A}}\times {\boldsymbol {A}}^*\varphi ^{\prime}{\varphi}^{*\prime\prime\prime}+\text {c.c.} \right \}, \end{align}

(5.10d) \begin{align} {\boldsymbol {M}}_1-\boldsymbol {\hat {\imath }}\times \frac {\partial \boldsymbol {t}_1}{\partial s}=\boldsymbol {t}_{1/2}\times \frac {\partial \boldsymbol {t}_{1/2}}{\partial s} \quad \left \{={\boldsymbol {A}}\times {\boldsymbol {A}}^*\varphi ^{\prime}{\varphi}^{*\prime\prime} + \text {c.c.}\right \}; \end{align}

the inextensibility constraint (2.7) gives

(5.11) \begin{equation} \boldsymbol {\hat {\imath }}\boldsymbol {\cdot }\boldsymbol {t}_1 = -\frac {1}{2}\boldsymbol {t}_{1/2}\boldsymbol {\cdot }\boldsymbol {t}_{1/2} \quad \left \{= -\frac {1}{2}{\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}} {\varphi ^{\prime}}^2{\rm e}^{2i\omega \tau }-\frac {1}{2}|{\boldsymbol {A}}|^2|\varphi ^{\prime}|^2 + \text {c.c.}\right \}; \end{equation}

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

(5.12a,b) \begin{equation} {\boldsymbol {r}}_{1}={\mathbf {0}}, \quad \boldsymbol {t}_1={\mathbf {0}} \quad \text {at} \quad s=0 \end{equation}

and

(5.13a,b) \begin{equation} {\boldsymbol {F}}_1+\mathcal {F}_c\boldsymbol {t}_1=-\chi \boldsymbol {\hat {\imath }} \quad \left \{=-\frac {1}{2}\chi \boldsymbol {\hat {\imath }}+\text {c.c.}\right \}, \quad {\boldsymbol {M}}_1={\mathbf {0}} \quad \text {at} \quad s=1. \end{equation}

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

(5.14a) \begin{align} {\boldsymbol {r}}_1=-\boldsymbol {\hat {\imath }}\frac {1}{2}{\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}} {\rm e}^{2i\omega \tau }\int _0^s \{{\varphi ^{\prime}}(p)\}^2\,{\rm d}p-\boldsymbol {\hat {\imath }}\frac {1}{2}|{\boldsymbol {A}}|^2\int _0^s|\varphi ^{\prime}(p)|^2\,{\rm d}p + \text {c.c.}, \end{align}

(5.14b) \begin{align} \boldsymbol {t}_1=-\boldsymbol {\hat {\imath }}\frac {1}{2}{\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}} {\varphi ^{\prime}}^2{\rm e}^{2i\omega \tau }-\boldsymbol {\hat {\imath }}\frac {1}{2}|{\boldsymbol {A}}|^2|\varphi ^{\prime}|^2 + \text {c.c.},\end{align}

(5.14c) \begin{align} {\boldsymbol {F}}_1 & =\frac {1}{2}\boldsymbol {\hat {\imath }} {\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}} {\rm e}^{2i\omega \tau }\left [\mathcal {F}_c\{\varphi ^{\prime}(1)\}^2+\frac {i\omega }{2}\varphi ^2(1)-\frac {i\omega }{2}\varphi ^2(s)+i\omega \int _s^1{\rm d}p\,\int _0^p {\rm d}q\,\{\varphi ^{\prime}(q)\}^2\right ]\nonumber\\& \quad -\boldsymbol {\hat {\imath }}\frac {\chi }{2}+\boldsymbol {\hat {\imath }}\frac {\mathcal {F}_c}{2}|{\boldsymbol {A}}|^2|\varphi ^{\prime}(1)|^2-\boldsymbol {\hat {\imath }}\frac {i\omega }{2}|{\boldsymbol {A}}|^2\int _s^1\varphi ^{\prime}(p)\varphi ^*(p)\,{\rm d}p +\text {c.c.} \end{align}

(5.14d) \begin{align} {\boldsymbol {M}}_1={\boldsymbol {A}}\times {\boldsymbol {A}}^*\varphi ^{\prime}{\varphi}^{*\prime\prime} + \text {c.c.}, \end{align}

whose periodicity in $\tau$ confirms the regularity of the weakly nonlinear expansions to $O(\varepsilon )$ (cf. (5.3)). In (5.14), we have discarded all homogeneous solutions, including beating/whirling solutions of the form (5.9) as well as solutions decaying on the fast time scale – these do not effect the slow-time dynamics of the leading-order amplitude ${\boldsymbol {A}}(T)$ .

5.4. The $O(\varepsilon ^{3/2})$ problem

At $O(\varepsilon ^{3/2})$ , the governing partial differential equations (2.6), together with the time-derivative transformation (5.2), give

(5.15a) \begin{equation} \frac {\partial {\boldsymbol {r}}_{3/2}}{\partial s}-\boldsymbol {t}_{3/2}={\mathbf {0}}, \end{equation}

(5.15b) \begin{align} \frac {\partial {\boldsymbol {F}}_{3/2}}{\partial s} & -\left (\mathsf {I}-\frac {1}{2}\boldsymbol {\hat {\imath }}\boldsymbol {\hat {\imath }}\right )\boldsymbol {\cdot } \frac {\partial {\boldsymbol {r}}_{3/2}}{\partial \tau } = \left (\mathsf {I}-\frac {1}{2}\boldsymbol {\hat {\imath }}\boldsymbol {\hat {\imath }}\right )\boldsymbol {\cdot } \frac {\partial {\boldsymbol {r}}_{1/2}}{\partial T}-\frac {1}{2}\left (\boldsymbol {t}_{1/2}\boldsymbol {\hat {\imath }}+\boldsymbol {\hat {\imath }}\boldsymbol {t}_{1/2}\right )\boldsymbol {\cdot }\frac {\partial {\boldsymbol {r}}_{1}}{\partial \tau }\nonumber\\& \quad - \frac {1}{2}\left (\boldsymbol {t}_1\boldsymbol {\hat {\imath }}+\boldsymbol {\hat {\imath }}\boldsymbol {t}_1+\boldsymbol {t}_{1/2}\boldsymbol {t}_{1/2}\right )\boldsymbol {\cdot }\frac {\partial {\boldsymbol {r}}_{1/2}}{\partial \tau } \quad \left \{={\rm e}^{i\omega \tau }{\boldsymbol {f}}_{(2)}(s)+\text {c.c.} + \text {n.s.t.}\right \}, \end{align}

(5.15c) \begin{align} & \frac {\partial {\boldsymbol {M}}_{3/2}}{\partial s}+\mathcal {F}_c\boldsymbol {\hat {\imath }}\times \boldsymbol {t}_{3/2}+\boldsymbol {\hat {\imath }}\times {\boldsymbol {F}}_{3/2} = -\boldsymbol {t}_{1/2}\times {\boldsymbol {F}}_{1}-\boldsymbol {t}_1\times {\boldsymbol {F}}_{1/2} \nonumber\\& \left \{={\rm e}^{i\omega \tau }{\boldsymbol {f}}_{(3)}(s)+\text {c.c.} + \text {n.s.t.}\right \}, \end{align}

(5.15d) \begin{align} {\boldsymbol {M}}_{3/2}-\boldsymbol {\hat {\imath }}\times \frac {\partial \boldsymbol {t}_{3/2}}{\partial s}=\boldsymbol {t}_{1/2}\times \frac {\partial \boldsymbol {t}_1}{\partial s}+\boldsymbol {t}_1\times \frac {\partial \boldsymbol {t}_{1/2}}{\partial s} \quad \left \{={\rm e}^{i\omega \tau }{\boldsymbol {f}}_{(4)}(s)+\text {c.c.} + \text {n.s.t.}\right \}; \end{align}

the inextensibility constraint (2.7) gives

(5.16) \begin{equation} \boldsymbol {t}_{3/2}\boldsymbol {\cdot }\boldsymbol {\hat {\imath }}=-\boldsymbol {t}_1\boldsymbol {\cdot }\boldsymbol {t}_{1/2}, \end{equation}

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

(5.17) \begin{equation} {\boldsymbol {r}}_{3/2}={\mathbf {0}}, \quad \boldsymbol {t}_{3/2}={\mathbf {0}} \quad \text {at} \quad s=0 \end{equation}

and

(5.18) \begin{equation} {\boldsymbol {F}}_{3/2}+\mathcal {F}_c\boldsymbol {t}_{3/2}=-\chi \boldsymbol {t}_{1/2} \quad \left \{={\rm e}^{i\omega \tau }{\boldsymbol {h}}_{(1)}+\text {c.c.} + \text {n.s.t.}\right \}, \quad {\boldsymbol {M}}_{3/2}={\mathbf {0}} \quad \text {at} \quad s=1. \end{equation}

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

(5.19a) \begin{align} {\boldsymbol {f}}_{(2)} & =\frac {\mathrm {d}{\boldsymbol {A}}}{\mathrm {d}T}\varphi +\frac {i\omega }{2}{\boldsymbol {A}}^*({\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}}){\varphi}^{*\prime}\int _0^s\{\varphi ^{\prime}(p)\}^2\,\mathrm {d}p\nonumber\\&\quad -\frac {i\omega }{2}\left [{\boldsymbol {A}}|{\boldsymbol {A}}|^2\varphi ^{\prime}({\varphi^{*\prime}}\varphi -\varphi ^{\prime}{\varphi ^*})+{\boldsymbol {A}}^*({\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}})\varphi |\varphi ^{\prime}|^2\right ], \end{align}

(5.19b) \begin{align} {\boldsymbol {f}}_{(3)} & = -\boldsymbol {\hat {\imath }}\times {\boldsymbol {A}}\varphi ^{\prime}\chi +\boldsymbol {\hat {\imath }}\nonumber\\ & \quad \times {\boldsymbol {A}}|{\boldsymbol {A}}|^2\left \{\mathcal {F}_c\varphi ^{\prime}(|\varphi ^{\prime}(1)|^2-|\varphi ^{\prime}|^2)-|\varphi ^{\prime}|^2\varphi ^{\prime\prime\prime}+\omega \varphi ^{\prime}\mathrm {Im}\int _s^1\varphi ^{\prime}(p)\varphi ^*(p)\,\mathrm {d}p\right \}\nonumber\\& \quad +\frac {1}{2}\boldsymbol {\hat {\imath }}\times {\boldsymbol {A}}^*({\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}})\left [-{\varphi ^{\prime}}^2(\varphi ^{\prime\prime\prime}+\mathcal {F}_c\varphi ^{\prime})^*+\mathcal {F}_c{\varphi^{*\prime}}\{\varphi ^{\prime}(1)\}^2+\frac {i\omega }{2}{\varphi^{*\prime}}\varphi ^2(1) \right .\nonumber\\ & \quad \left . -\frac {i\omega }{2}{\varphi^{*\prime}}\varphi ^2+i\omega {\varphi^{*\prime}}\int _s^1\mathrm {d}p\,\int _0^p \mathrm {d}q\,\{\varphi ^{\prime}(q)\}^2\right ], \end{align}

(5.19c) \begin{align} {\boldsymbol {f}}_{(4)}=\boldsymbol {\hat {\imath }}\times {\boldsymbol {A}} |{\boldsymbol {A}}|^2{\varphi ^{\prime}}^2{\varphi}^{*\prime\prime} + \boldsymbol {\hat {\imath }}\times {\boldsymbol {A}}^*({\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {A}})\left (\varphi ^{\prime\prime}|\varphi ^{\prime}|^2-\frac {1}{2}{\varphi ^{\prime}}^2{\varphi}^{*\prime\prime}\right ), \end{align}

(5.19d) \begin{align}{\boldsymbol {h}}_{(1)}=-\chi {\boldsymbol {A}}\varphi ^{\prime}(1), \end{align}

obtained by substituting (5.9) and (5.14); their numbering will become meaningful in the following subsection.

5.5. Solvability condition

5.5.1. Restricted linear problem

The $\exp (i\omega \tau )$ Fourier component of the $O(\varepsilon ^{3/2})$ problem, projected parallel to the wall, furnishes a ‘restricted’ linear problem in the form

(5.20a,b,c) \begin{equation} \mathsf {L}{\vec{\psi}} = {\vec{f}}, \quad \mathsf {S}{\vec{\psi}}(0)= {\vec{g}}, \quad \mathsf {T}{\vec{\psi}}(1)= {\vec{h}}, \end{equation}

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

(5.21a,b,c) \begin{equation} \mathsf {L} = \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\frac {{\rm d}}{{\rm d}s} & -1 & 0 & 0 \\[2pt]-i\omega & 0 & \frac {{\rm d}}{{\rm d}s} & 0 \\[2pt]0 & \mathcal {F}_c\boldsymbol {\hat {\imath }}\times & \boldsymbol {\hat {\imath }}\times & \frac {{\rm d}}{{\rm d}s} \\[2pt]0 & -\boldsymbol {\hat {\imath }}\times \frac {{\rm d}}{{\rm d}s} & 0 & 1\end{array}\right ), \quad \mathsf {S} = \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}1 & 0 & 0 & 0 \\0& 1 & 0 & 0 \end{array}\right ), \quad \mathsf {T} = \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & \mathcal {F}_c & 1 & 0 \\0& 0 & 0 & 1 \end{array}\right ); \end{equation}

and ${\vec{f}}(s)=[{\boldsymbol {f}}_{(1)}(s),\,{\boldsymbol {f}}_{(2)}(s)\,{\boldsymbol {f}}_{(3)}(s),\,{\boldsymbol {f}}_{(4)}(s)]^T$ , ${\vec{g}}=[{\boldsymbol {g}}_{(1)},\,{\boldsymbol {g}}_{(2)}]^T$ and ${\vec{h}}=[{\boldsymbol {h}}_{(1)},\, {\boldsymbol {h}}_{(2)}]^T$ are column-array forcing terms, whose elements are vectors parallel to the wall (fields for ${\vec{f}}(s)$ , constants for ${\vec{g}}$ and ${\vec{h}}$ ). In the present scenario, we find from (5.15)–(5.18) that ${\boldsymbol {f}}_{(1)}(s)$ , ${\boldsymbol {h}}_{(2)}$ , ${\boldsymbol {g}}_{(1)}$ and ${\boldsymbol {g}}_{(2)}$ vanish, whereas ${\boldsymbol {f}}_{(2)}(s)$ , ${\boldsymbol {f}}_{(3)}(s)$ , ${\boldsymbol {f}}_{(4)}(s)$ and ${\boldsymbol {h}}_{(1)}$ are provided by (5.19). We know from the linear theory of § 3 that the restricted problem possesses non-trivial homogeneous solutions in the form $ {\vec{\psi}} = [{\boldsymbol {a}}\varphi , \,\, {\boldsymbol {a}}\varphi ^{\prime}, \,\, -{\boldsymbol {a}}(\varphi ^{\prime\prime\prime}+\mathcal {F}_c\varphi ^{\prime}),\,\, \boldsymbol {\hat {\imath }}\times {\boldsymbol {a}}\varphi ^{\prime\prime} ]^T$ , with $\boldsymbol {a}$ an arbitrary complex vector parallel to the wall. We therefore expect that solutions to the inhomogeneous restricted problem exist – whereby resonance is avoided in the $O(\varepsilon ^{3/2})$ problem of § 5.4 – only under certain ‘solvability’ conditions on the forcing terms.

Below, we derive a necessary condition for existence of solutions to the restricted problem. (By alluding to the Fredholm alternative theorem for differential operators (Keener Reference Keener2018), it could be shown that the condition is also sufficient.) In deriving this solvability condition, we shall allow for the full form of the forcing terms in (5.20) in order to facilitate future generalisations of the theory as discussed in § 6.

5.5.2. Adjoint operators and homogeneous solutions

For any pair of column-vector fields $ {\vec{\psi}}^{({1})}(s) = [ \acute {{\boldsymbol {r}}}^{({1})}(s), \,\, \acute {\boldsymbol {t}}^{({1})}(s), \,\, \acute {{\boldsymbol {F}}}^{({1})}(s), \,\, \acute {{\boldsymbol {M}}}^{({1})}(s) ]^T$ and ${\vec{\psi}}^{({2})}(s) = [ \acute {{\boldsymbol {r}}}^{({2})}(s), \,\, \acute {\boldsymbol {t}}^{({2})}(s), \,\, \acute {{\boldsymbol {F}}}^{({2})}(s), \,\, \acute {{\boldsymbol {M}}}^{({2})}(s) ]^T$ , we define the inner product

(5.22) \begin{equation} \langle {\vec{\psi}}^{({1})}, {\vec{\psi}}^{({2})} \rangle = \int _0^1\left \{{{\boldsymbol {r}}^*}^{({1})}\boldsymbol {\cdot }{\boldsymbol {r}}^{({2})}+{\boldsymbol {t}^*}^{({1})}\boldsymbol {\cdot }\boldsymbol {t}^{({2})}+{{\boldsymbol {F}}^*}^{({1})}\boldsymbol {\cdot }{\boldsymbol {F}}^{({2})}+{{\boldsymbol {M}}^*}^{({1})}\boldsymbol {\cdot }{\boldsymbol {M}}^{({2})}\right \}\,\mathrm {d}s. \end{equation}

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

(5.23) \begin{equation} J( {\vec{\psi}}^{({1})}, {\vec{\psi}}^{({2})})=\langle \mathsf {L} {\vec{\psi}}^{({1})}, {\vec{\psi}}^{({2})} \rangle - \langle {\vec{\psi}}^{({1})}, \mathsf {L}^\dagger {\vec{\psi}}^{({2})}\rangle \end{equation}

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

(5.24) \begin{align} & \langle \mathsf {L} {\vec{\psi}}^{({1})}, {\vec{\psi}}^{({2})} \rangle = \int _0^1 \left \{\left (\frac {{\rm d}{{{\acute {{\boldsymbol {r}}}}}}^{({1})}{}^*}{{\rm d}s}-{\acute {\boldsymbol {t}}}^{({1})}{}^*\right ) \boldsymbol {\cdot } \acute {{\boldsymbol {r}}}^{({2})}+\left (i\omega {\acute {{\boldsymbol {r}}}}^{({1})}{}^*+\frac {{\rm d}{\acute {{\boldsymbol {F}}}}^{({1})}{}^*}{{\rm d}s}\right )\boldsymbol {\cdot }\acute {\boldsymbol {t}}^{({2})} \right . \nonumber\\& \quad \left . +\left (\mathcal {F}_c\boldsymbol {\hat {\imath }}\times {\acute {\boldsymbol {t}}}^{({1})}{}^* + \boldsymbol {\hat {\imath }}\times {\acute {{\boldsymbol {F}}}}^{({1})}{}^*+\frac {{\rm d}{\acute {{\boldsymbol {M}}}}^{({1})}{}^*}{{\rm d}s}\right )\boldsymbol {\cdot }\acute {{\boldsymbol {F}}}^{({2})}+\left (-\boldsymbol {\hat {\imath }}\times \frac {{\rm d}{\acute {\boldsymbol {t}}}^{({1})}{}^*}{{\rm d}s}+{\acute {{\boldsymbol {M}}}}^{({1})}{}^*\right )\boldsymbol {\cdot }\acute {{\boldsymbol {M}}}^{({2})}\right \}\mathrm {ds} \end{align}

and use the vector triple product to find

(5.25) \begin{align} & \langle \mathsf {L} {\vec{\psi}}^{({1})}, {\vec{\psi}}^{({2})} \rangle = \left [\acute {{\boldsymbol {r}}}^{({2})}\boldsymbol {\cdot }\acute {{\boldsymbol {r}}}^{({1})}{}^*+\boldsymbol {\hat {\imath }}\times \acute {{\boldsymbol {M}}}^{({2})}\boldsymbol {\cdot } \acute {\boldsymbol {t}}^{({1})}{}^*+\acute {\boldsymbol {t}}^{({2})}\boldsymbol {\cdot }{\acute {{\boldsymbol {F}}}}^{({1})}{}^*+\acute {{\boldsymbol {F}}}^{({2})}\boldsymbol {\cdot }{\acute {{\boldsymbol {M}}}}^{({1})}{}^*\right ]_0^1 \nonumber\\& \quad +\int _0^{1}\left \{{\acute {{\boldsymbol {r}}}}^{({1})}{}^*\boldsymbol {\cdot }\left (-\frac {{\rm d}\acute {{\boldsymbol {r}}}^{({2})}}{{\rm d}s}+i\omega \acute {\boldsymbol {t}}^{({2})}\right )+{\acute {\boldsymbol {t}}}^{({1})}{}^*\boldsymbol {\cdot }\left (-\acute {{\boldsymbol {r}}}^{({2})}-\mathcal {F}_c\boldsymbol {\hat {\imath }}\times \acute {{\boldsymbol {F}}}^{({2})}-\boldsymbol {\hat {\imath }}\times \frac {{\rm d}\acute {{\boldsymbol {M}}}^{({2})}}{{\rm d}s}\right ) \right . \nonumber\\& \quad \left . +{\acute {{\boldsymbol {F}}}}^{({1})}{}^*\boldsymbol {\cdot }\left (-\frac {{\rm d}\acute {\boldsymbol {t}}^{({2})}}{{\rm d}s}-\boldsymbol {\hat {\imath }}\times \acute {{\boldsymbol {F}}}^{({2})}\right )+{\acute {{\boldsymbol {M}}}}^{({1})}{}^*\boldsymbol {\cdot }\left (-\frac {{\rm d}\acute {{\boldsymbol {F}}}^{({2})}}{{\rm d}s}+\acute {{\boldsymbol {M}}}^{({2})}\right )\right \}\,\mathrm {d}s. \end{align}

Thus, by inspection, the adjoint operators are

(5.26a,b,c) \begin{align} \mathsf {L}^\dagger & = \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}-\frac {{\rm d}}{{\rm d}s} & i\omega & 0 & 0 \\-1 & 0 & -\mathcal {F}_c\boldsymbol {\hat {\imath }}\times & -\boldsymbol {\hat {\imath }}\times \frac {{\rm d}}{{\rm d}s} \\0 & -\frac {{\rm d}}{{\rm d}s} & -\boldsymbol {\hat {\imath }}\times & 0 \\0 & 0 & -\frac {{\rm d}}{{\rm d}s} & 1\end{array}\right ), \quad \mathsf {S}^\dagger = \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & 1 & 0 & 0 \\0& 0 & 1 & 0 \end{array}\right ),\nonumber\\\mathsf {T}^\dagger & = \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}1 & 0 & 0 & 0 \\0& -\mathcal {F}_c & 0 & \boldsymbol {\hat {\imath }}\times \end{array}\right ). \end{align}

Associated with the adjoint operators is the adjoint homogeneous problem,

(5.27a,b,c) \begin{equation} \mathsf {L}^\dagger {\vec{\zeta}} = {\vec{0}}, \quad \mathsf {S}^\dagger {\vec{\zeta}}(0)= {\vec{0}}, \quad \mathsf {T}^\dagger {\vec{\zeta}}(1)= {\vec{0}}. \end{equation}

It is readily found to possess the non-trivial solutions

(5.28) \begin{equation} {\vec{\zeta}}=\left [{\boldsymbol {b}}\phi ,\,\, {\boldsymbol {b}}\frac {1}{i\omega }\phi ^{\prime},\,\, \boldsymbol {\hat {\imath }}\times {\boldsymbol {b}}\frac {1}{i\omega }\phi ^{\prime\prime}, \,\, \boldsymbol {\hat {\imath }}\times {\boldsymbol {b}}\frac {1}{i\omega }\phi ^{\prime\prime\prime}\right ]^T, \end{equation}

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

(5.29) \begin{equation} \phi ^{\prime\prime\prime\prime}+\mathcal {F}_c\phi ^{\prime\prime}-i\omega \phi =0 \end{equation}

and the boundary conditions

(5.30a,b,c,d) \begin{equation} \phi ^{\prime} =0, \quad \phi ^{\prime\prime} =0 \quad \text {at} \quad s=0; \quad \phi =0, \quad \phi ^{\prime\prime\prime}+\mathcal {F}_c\phi ^{\prime}=0 \quad \text {at} \quad s=1. \end{equation}

The unidimensional problem (5.29)–(5.30) defines the adjoint eigenfunction $\phi (s)$ up to an arbitrary complex prefactor, which we conveniently choose such that $\phi (0)=1$ ; the eigenfunction can be readily expressed analytically in terms of the numerical values $\mathcal {F}_c$ and $\omega$ , or calculated numerically.

5.5.3. Solvability condition

In the definition (5.23) of the factor $J( {\vec{\psi}}^{({1})}, {\vec{\psi}}^{({2})})$ , let $ {\vec{\psi}}^{({1})}= {\vec{\psi}}$ be a solution to the inhomogeneous problem (5.20) and $ {\vec{\psi}}^{({2})}= {\vec{\zeta}}$ any solution to the homogeneous adjoint problem (cf. (5.28)). Using (5.20a ) and (5.27a ), we find from (5.23) the relation

(5.31) \begin{equation} J( {\vec{\psi}}, {\vec{\zeta}})=\langle {\vec{f}}, {\vec{\zeta}}\rangle . \end{equation}

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

(5.32) \begin{equation} J( {\vec{\psi}}, {\vec{\zeta}})={\boldsymbol {b}}\boldsymbol {\cdot }\left (\frac {\phi ^{\prime}(1)}{i\omega }{\boldsymbol {h}}_{(1)}^*+\frac {\phi ^{\prime\prime}(1)}{i\omega }{\boldsymbol {h}}_{(2)}^*\times \boldsymbol {\hat {\imath }}-\phi (0){\boldsymbol {g}}_{(1)}^*+\frac {\phi ^{\prime\prime\prime}(0)}{i\omega }{\boldsymbol {g}}_{(2)}^*\right ). \end{equation}

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

(5.33) \begin{equation} \langle {\vec{f}}, {\vec{\zeta}}\rangle ={\boldsymbol {b}}\boldsymbol {\cdot }\int _0^1\left ({\boldsymbol {f}}_{(1)}^*\phi +{\boldsymbol {f}}_{(2)}^*\frac {\phi ^{\prime}}{i\omega }+{\boldsymbol {f}}_{(3)}^*\times \boldsymbol {\hat {\imath }}\frac {\phi ^{\prime\prime}}{i\omega }+{\boldsymbol {f}}_{(4)}^*\times \boldsymbol {\hat {\imath }}\frac {\phi ^{\prime\prime\prime}}{i\omega }\right )\,\mathrm {d} s. \end{equation}

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

(5.34) \begin{align} & \int _0^1\left ({\boldsymbol {f}}_{(1)}^*\phi +{\boldsymbol {f}}_{(2)}^*\frac {\phi ^{\prime}}{i\omega }+{\boldsymbol {f}}_{(3)}^*\times \boldsymbol {\hat {\imath }}\frac {\phi ^{\prime\prime}}{i\omega }+{\boldsymbol {f}}_{(4)}^*\times \boldsymbol {\hat {\imath }}\frac {\phi ^{\prime\prime\prime}}{i\omega }\right )\,\mathrm {d}s\nonumber\\& \qquad =\frac {\phi ^{\prime}(1)}{i\omega }{\boldsymbol {h}}_{(1)}^*+\frac {\phi ^{\prime\prime}(1)}{i\omega }{\boldsymbol {h}}_{(2)}^*\times \boldsymbol {\hat {\imath }}-\phi (0){\boldsymbol {g}}_{(1)}^*+\frac {\phi ^{\prime\prime\prime}(0)}{i\omega }{\boldsymbol {g}}_{(2)}^*. \end{align}

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

(5.35a–c) \begin{equation} \int _0^1\left ({\boldsymbol {f}}_{(2)}{\phi^{*\prime}}-\boldsymbol {\hat {\imath }}\times {\boldsymbol {f}}_{(3)}{\varphi}^{*\prime\prime}-\boldsymbol {\hat {\imath }}\times {\boldsymbol {f}}_{(4)}{\phi}^{*\prime\prime\prime}\right )\,{\rm d}s = {\boldsymbol {h}}_{(1)}{\phi^{*\prime}}(1), \end{equation}

where we have also taken the complex conjugate. Substituting (5.19) for ${\boldsymbol {f}}_{(2)}(s)$ , ${\boldsymbol {f}}_{(3)}(s)$ , ${\boldsymbol {f}}_{(4)}(s)$ and ${\boldsymbol {h}}_{(1)}$ , we arrive at the amplitude equation (4.6) with the coefficients $\alpha$ , $\beta$ and $\gamma$ defined by the quadratures presented in Appendix B.