2 Continuum model

Consider a nearly inextensible elastic filament of length $\ell$ and diameter $a$ rooted at the origin $O$ of a Cartesian coordinate system $(x,y,z)$ . Let $\{\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3}\}$ be the corresponding orthonormal basis. The filament is free to deform in the half-space fluid domain $z\geqslant 0$ , where $z=0$ corresponds to a no-slip solid wall (see figure 1). The centreline of the filament is denoted by the position vector $\boldsymbol{r}(s,t)$ , where $s$ and $t$ represent the arclength and time, respectively. The balance of forces and moments on a cross-section of the filament are given by Kirchhoff’s equations for an elastic rod (Audoly & Pomeau Reference Audoly and Pomeau2010)

Here, the prime $(\cdot )^{\prime }=\unicode[STIX]{x2202}(\cdot )/\unicode[STIX]{x2202}s$ denotes differentiation with respect to arclength $s$ , $\hat{\boldsymbol{t}}=\boldsymbol{r}^{\prime }/|\boldsymbol{r}^{\prime }|$ is the tangent unit vector along the filament centreline, $\boldsymbol{N}$ and $\boldsymbol{M}$ are the internal force and bending moment, respectively, and $-\boldsymbol{f}$ is the drag force per unit length exerted by the surrounding fluid on the filament ( $\boldsymbol{f}$ is the force per unit length exerted by the filament on the fluid). The Hookean constitutive relation between the bending moment $\boldsymbol{M}$ and the bending deformation (curvature) of the filament is given by $\boldsymbol{M}=B\,\hat{\boldsymbol{t}}\times \hat{\boldsymbol{t}}^{\prime }$ , where $B$ is the bending rigidity. The internal force $\boldsymbol{N}$ consists of a bending force and a constraint tension force that enforces the inextensibility condition. The constraint is satisfied in a weak form by considering an elastic filament with large tensile stiffness.

The filament is free at its tip $s=\ell$ and is actuated by an internal motor at its base $s=0$ that produces a torque $\boldsymbol{M}(0,t)=\boldsymbol{M}_{b}$ . The torque $\boldsymbol{M}_{b}$ is a configuration-dependent torque that switches direction when the base angle $\unicode[STIX]{x1D703}_{b}(t)$ of the filament, defined as $\unicode[STIX]{x1D703}_{b}=\arcsin (\boldsymbol{e}_{3}\times \hat{\boldsymbol{t}}(0,t)\boldsymbol{\cdot }\boldsymbol{e}_{2})$ , reaches predefined target orientations $\pm \unicode[STIX]{x1D6E9}$ . More specifically, we consider $\boldsymbol{M}_{b}=\unicode[STIX]{x1D6FC}M_{b}\boldsymbol{e}_{2}$ , where $M_{b}$ is a positive constant and $\unicode[STIX]{x1D6FC}\in \{-1,1\}$ is a state variable that defines the torque direction; $\unicode[STIX]{x1D6FC}$ changes from $1$ to $-1$ at $\unicode[STIX]{x1D703}_{b}=\unicode[STIX]{x1D6E9}$ and from $-1$ to $1$ at $\unicode[STIX]{x1D703}_{b}=-\unicode[STIX]{x1D6E9}$ (see figure 1). This torque model can be viewed as an extension to the geometric switch model for colloidal systems studied in Kotar et al. (Reference Kotar, Leoni, Bassetti, Lagomarsino and Cicuta2010) and Bruot et al. (Reference Bruot, Kotar, de Lillo, Lagomarsino and Cosentino2012), and as a simplified version of the ‘geometric clutch’ model proposed in Lindemann (Reference Lindemann1994). A target angle that acts as a geometric switch to drive elastic filaments along their entire length was used by Buchmann, Cortez & Fauci (Reference Buchmann, Cortez and Fauci2017) to model the power and recovery strokes of eukaryotic cilia.

The fluid motion is governed by the incompressible Stokes equation for zero-Reynolds-number flows,

Here, $p$ is the pressure field, $\unicode[STIX]{x1D707}$ is the fluid viscosity, $\boldsymbol{u}$ is the fluid velocity field, and $\boldsymbol{F}$ is the (Eulerian) force density exerted by the filament on the fluid. The last, $\boldsymbol{F}(\boldsymbol{x},t)$ , is related to the force per unit length $\boldsymbol{f}(s,t)$ as $\boldsymbol{F}(\boldsymbol{x},t)=\int _{s\in [0,\ell ]}\boldsymbol{f}(s,t)\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{r}(s,t))\,\text{d}s$ , where $\unicode[STIX]{x1D6FF}$ is the 3D Dirac $\unicode[STIX]{x1D6FF}$ -function and $\boldsymbol{x}$ is the position vector. These equations are subject to the no-slip condition $\boldsymbol{u}=0$ at the bounding wall $z=0$ . We take advantage of the small aspect ratio $a/\ell \ll 1$ of the filament to approximate the velocity at the filament boundary by the velocity along its centreline,

To fully determine the filament deformation $\boldsymbol{r}(s,t)$ given the moment $\boldsymbol{M}_{b}$ at the filament base, we need to solve the coupled fluid–filament system of equations (2.1)–(2.3). It is convenient for building an efficient numerical method to (i) write the moment equation in (2.1) in integral form and (ii) assume that the filament is quasi-inextensible (Teran, Fauci & Shelley Reference Teran, Fauci and Shelley2010; Chrispell, Fauci & Shelley Reference Chrispell, Fauci and Shelley2013; Olson & Fauci Reference Olson and Fauci2015). In particular, we integrate the moment equation in (2.1) from the filament free end at $\ell$ to any location $s$ along the filament, taking into account that $\boldsymbol{M}(\ell )=0$ and that along the filament $\boldsymbol{N}^{\prime }=\boldsymbol{f}$ . We get, after an integration by parts on the second term, that

We then write the force density $\boldsymbol{f}(s,t)$ applied by the filament on the surrounding fluid as $\boldsymbol{f}=\boldsymbol{f}^{\bot }+\boldsymbol{f}^{\Vert }$ . We assume that the force component $\boldsymbol{f}^{\Vert }=(\,\boldsymbol{f}\boldsymbol{\cdot }\hat{\boldsymbol{t}})\hat{\boldsymbol{t}}$ tangent to the filament’s centreline can be obtained explicitly by considering a large tensile stiffness  $K$ ,

thus ensuring that the filament’s length remains almost constant. We substitute (2.5) into (2.4) taking into account that $\boldsymbol{M}=B\,\boldsymbol{r}^{\prime }\times \boldsymbol{r}^{\prime \prime }$ to obtain an expression for $\boldsymbol{f}^{\bot }$ in terms of the position vector $\boldsymbol{r}(s)$ and its spatial derivatives. To this end, one gets both components of the force density $\boldsymbol{f}=\boldsymbol{f}^{\bot }+\boldsymbol{f}^{\Vert }$ in terms of the kinematic variables $\boldsymbol{r}$ and its derivatives. We substitute these expressions for $\boldsymbol{f}$ into (2.2) and we solve numerically subject to (2.3) to obtain the nonlinear dynamics of the filament, as discussed next.

Table 1.

Characteristic scales of the system.

To obtain non-dimensional counterparts to the equations of motion, we consider the dimensional scales associated with the fluid viscosity $\unicode[STIX]{x1D707}$ and cilium length $\ell$ . Because of the geometric switch model, the system does not have an intrinsic time scale. To remedy this, we consider the time scale $T=\ell ^{4}\unicode[STIX]{x1D707}/B$ arising from balancing the filament’s elasticity with the fluid viscosity. To this end, we consider the bending rigidity to be of the order $B=800~\text{pN}~\unicode[STIX]{x03BC}\text{m}^{2}$ , as reported in Xu et al. (Reference Xu, Wilson, Okamoto, Shao, Dutcher and Bayly2016) for wild-type Chlamydomonas flagella. A list of the dimensional parameters used to scale the equations of the motion is reported in table 1. Hereafter, all quantities are considered to be dimensionless unless otherwise stated.

3 Numerical method

We discretize the filament into a uniform chain of $N+1$ segments of length $a$ such that $\unicode[STIX]{x0394}s=\ell /N=a$ (figure 1 b). The segments are labelled from $n=0$ at the filament base to $n=N$ at its tip. The position vector is discretized by $\boldsymbol{r}_{n}=x_{n}\boldsymbol{e}_{1}+z_{n}\boldsymbol{e}_{3}$ and the local orientation $\unicode[STIX]{x1D703}_{n}$ of the tangent vector to segment $n$ is defined as the angle between the $z$ -axis and the vector $\unicode[STIX]{x0394}\boldsymbol{r}_{n}=(\boldsymbol{r}_{n+1}-\boldsymbol{r}_{n})$ . Equation (2.4) can be written in discrete form as

Here, $\boldsymbol{M}_{n}=-B[(\unicode[STIX]{x1D703}_{n}-\unicode[STIX]{x1D703}_{n-1})/\unicode[STIX]{x0394}s]\boldsymbol{e}_{2}$ for $1\leqslant n<N$ whereas $\boldsymbol{M}_{o}=\boldsymbol{M}_{b}$ at the no-slip wall $z=0$ .

We decompose the force $\boldsymbol{f}_{n}=\boldsymbol{f}_{n}^{\bot }+\boldsymbol{f}_{n}^{\Vert }$ exerted by segment $n$ on the surrounding fluid into two components, $\boldsymbol{f}_{n}^{\bot }$ and $\boldsymbol{f}_{n}^{\Vert }$ , that are perpendicular and parallel to $\unicode[STIX]{x0394}\boldsymbol{r}_{n-1}=(\boldsymbol{r}_{n}-\boldsymbol{r}_{n-1})$ . We substitute into (3.1) and rearrange the term containing $\boldsymbol{f}_{n}^{\bot }$ to the other side of the equation to get

Upon taking the cross-product with $\unicode[STIX]{x0394}\boldsymbol{r}_{n-1}/\Vert \unicode[STIX]{x0394}\boldsymbol{r}_{n-1}\Vert ^{2}$ , the above equation becomes

The parallel component $\boldsymbol{f}_{n}^{\Vert }$ is given by the discrete analogue to (2.5). Namely, for $1\leqslant n<N$ , one has

whereas for $n=N$ , one has

Using (3.3) and (3.4)–(3.5), $\boldsymbol{f}_{n}$ can be evaluated sequentially from the filament tip (by decreasing order of $n$ ) in terms of the filament kinematic variables $\boldsymbol{r}_{m}$ , $m\geqslant n-1$ .

Figure 2.

Filament dynamics for (a,c,e,g,i) $M_{b}=1$ and (b,d,f,h,j) $M_{b}=3$ . From top to bottom: snapshots of the filament deformations, configuration-dependent moment $M_{b}$ , angular velocity $\dot{\unicode[STIX]{x1D703}}_{b}$ , angle $\unicode[STIX]{x1D703}_{b}$ and phase $\unicode[STIX]{x1D719}$ at the base of the filament as functions of time. Note that the frequency of oscillations is not known a priori, only the switching angles $\unicode[STIX]{x1D6E9}=\pm 0.15\unicode[STIX]{x03C0}$ .

To solve (2.2)–(2.3), we use a one-dimensional distribution of regularized stokeslets along the centreline of the filament together with an ‘image’ distribution to impose the no-slip boundary conditions at the $z=0$ plane (Cortez & Varela Reference Cortez and Varela2015). The regularized stokeslets are placed at the centre $\boldsymbol{r}_{n}$ of each segment. The strength of the regularized stokeslet at $\boldsymbol{r}_{n}$ is equal to the discrete force $\boldsymbol{f}_{n}$ , and the fluid velocity generated by the filament at an arbitrary position $\boldsymbol{x}$ in the fluid domain is given by $\boldsymbol{u}(\boldsymbol{x})=\sum _{n=1}^{N}\boldsymbol{G}(\boldsymbol{x}-\boldsymbol{r}_{n})\boldsymbol{\cdot }\boldsymbol{f}_{n}$ , where $\boldsymbol{G}(\boldsymbol{x}-\boldsymbol{r}_{n})$ is the Green’s tensor for the regularized stokeslet near a wall (Ainley et al. Reference Ainley, Durkin, Embid, Boindala and Cortez2008). We substitute this expression for the fluid velocity into (2.3), recalling (3.3) and (3.4)–(3.5) to express $\boldsymbol{f}_{n}$ in terms of the filament position $\boldsymbol{r}_{m}$ . This yields a set of coupled equations for the filament dynamics that we evolve forwards in time using the forward Euler method. Initial conditions for this system are the configuration of the filament $\boldsymbol{r}(s,t)$ and the state variable $\unicode[STIX]{x1D6FC}=\pm 1$ .

In all numerical simulations, we fix the target angles at $\unicode[STIX]{x1D6E9}=0.15\unicode[STIX]{x03C0}$ , we use $N=20$ segments of length $a=1/20$ to discretize the filament, and we set the regularization parameter of the stokeslet to be equal to $a$ . The tensile stiffness $K=5000$ is a numerical parameter chosen to keep the length of the filament almost constant. The filaments are initialized in straight configuration normal to the wall. Here, the time step is $0.2\times 10^{-4}$ and the system is integrated to $T=20$ , which is sufficient to capture the long-time dynamics; see table 2 for a summary of all parameter values.

Figure 3.

Beating frequency $\unicode[STIX]{x1D714}$ and Sperm number $Sp$ of a single filament as a function of bending moment at the base $M_{b}$ .

Table 2.

Dimensionless parameters used in the simulations ( $\ell ,~\unicode[STIX]{x1D707}$ and $B$ are normalized to 1).

4 Deformation of a single filament

We examine the motion of a single filament actively driven at its base by $\boldsymbol{M}_{b}=\unicode[STIX]{x1D6FC}M_{b}\boldsymbol{e}_{2}$ ( $\unicode[STIX]{x1D6FC}=\pm 1$ ) that switches direction at $\unicode[STIX]{x1D6E9}=\pm 0.15\unicode[STIX]{x03C0}$ . Figure 2 shows the deformations and time evolution of the filament for two distinct values of the bending moment: (a,c,e,g,i) $M_{b}=1$ and (b,d,f,h,j) $M_{b}=3$ . In each case, the moment $\boldsymbol{M}_{b}$ is constant in magnitude between switches (figure 2 c,d). In response, the angular velocity $\dot{\unicode[STIX]{x1D703}}_{b}$ slows down after each switching event (figure 2 e,f). This decrease in $\dot{\unicode[STIX]{x1D703}}_{b}$ is induced by the internal restoring moment due to the filament elasticity, which acts in the opposite direction to $\boldsymbol{M}_{b}$ .

We introduce the phase variable $\unicode[STIX]{x1D719}\in [0,1]$ as a linear interpolation of $\unicode[STIX]{x1D703}_{b}$ between the two target angles $\pm \unicode[STIX]{x1D6E9}$ ,

By definition, the values of $\unicode[STIX]{x1D719}$ lie in $[0,0.5]$ when the base moment is positive ( $\unicode[STIX]{x1D6FC}=1$ ) and $\unicode[STIX]{x1D719}\in [0.5,1]$ when the base moment switches to negative ( $\unicode[STIX]{x1D6FC}=-1$ ). Thus, the phase variable $\unicode[STIX]{x1D719}$ is monotonic in time over one oscillation period (see figure 2 i,j) and, as a result, it can be viewed as a time reparametrization for examining the long-term periodic behaviour of the filament. We therefore label the snapshots in figure 2(a,b) by their phase $\unicode[STIX]{x1D719}$ .

Figure 4.

(a) Filament with tilted angles produce asymmetric beating patterns. Here $M_{b}=3$ , target angles are $\unicode[STIX]{x03C0}/4\pm \unicode[STIX]{x1D6E9}$ . The dashed line shows the average orientation $\unicode[STIX]{x03C0}/4$ . (b) The filament with $M_{b}=3$ recovers the same beating pattern when starting from a perturbed initial condition; the long-term beating pattern is shown in green for comparison. Filaments with lower base moment converge to the long-term beating pattern within fewer cycles (see supplementary movies 1–4 available at https://doi.org/10.1017/jfm.2017.816).

The frequency of switching in the drive is a property of the system that depends on the system’s parameters, including the magnitude of the base moment $M_{b}$ . Figure 2 indicates that for $M_{b}=1$ , the base angle $\unicode[STIX]{x1D703}_{b}$ takes longer time to reach $\pm \unicode[STIX]{x1D6E9}$ than for $M_{b}=3$ . To quantify the frequency $\unicode[STIX]{x1D714}$ of switching, we define it as the number of left-side switching events per unit time, averaged over late-time behaviour. Figure 3 depicts $\unicode[STIX]{x1D714}$ as a function of $M_{b}$ . The frequency $\unicode[STIX]{x1D714}$ increases monotonically with $M_{b}$ as well as the waviness of the filament. For $M_{b}$ larger than $1.7$ , the curvature along the filament $(0<s<\ell )$ changes sign at least once, indicating a clear travelling-wave pattern. The switching frequency can be mapped to a Sperm number $Sp=(\ell ^{4}\unicode[STIX]{x1D709}_{\bot }\unicode[STIX]{x1D714}/B)^{1/4}$ , defined as the ratio between the filament length and the elasto-viscous penetration length (Wiggins & Goldstein Reference Wiggins and Goldstein1998; Lagomarsino, Capuani & Lowe Reference Lagomarsino, Capuani and Lowe2003; Eloy & Lauga Reference Eloy and Lauga2012). Here, $\unicode[STIX]{x1D709}_{\bot }=4\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}/\!\ln (2\ell /a)$ is the local drag coefficient perpendicular to the filament direction (Lauga & Powers Reference Lauga and Powers2009). The Sperm number $Sp$ increases monotonically with $M_{b}$ , ranging from $1.6$ to $4.3$ for $M_{b}\in [0.5,3]$ . This range of Sperm numbers is consistent with those observed empirically in flagella and cilia. For example, the wild-type Chlamydomonas has a beating frequency of 60–80 Hz (Leptos et al. Reference Leptos, Wan, Polin, Tuval, Pesci and Goldstein2013), which yields a Sperm number $Sp\approx 3$ , given the characteristic parameters listed in table 1.

Although the Sperm number is comparable to that of cilia and flagella, the filament deformations deviate from those observed in nature in that the amplitude of the travelling wave decreases towards the tip. This is due to the fact that the model considers a driving moment at the base only, while, in many biological filaments, the moments are distributed along the whole filament. In the model, the beating pattern is related to the choice of the target angle $\unicode[STIX]{x1D6E9}$ ; larger $\unicode[STIX]{x1D6E9}$ produces beating patterns with high curvature. Further, the filament deformations are symmetric because the switching orientations $\pm \unicode[STIX]{x1D6E9}$ are equal and opposite. To break this symmetry, it suffices to tilt the angle about which the geometric switch is applied by setting $\unicode[STIX]{x1D6E9}_{\text{left}}=0.1\unicode[STIX]{x03C0}$ and $\unicode[STIX]{x1D6E9}_{\text{right}}=0.4\unicode[STIX]{x03C0}$ as shown in figure 4(a). Hereafter, we restrict our discussion to the symmetric case. Finally, we note that the long-term behaviour of the filament depends on $M_{b}$ but it is independent of the filament initial configuration, as illustrated in figure 4(b).

5 Synchronization of two filaments

We consider the behaviour of two hydrodynamically coupled, elastic filaments separated by a distance $d$ and subject to the same moment $M_{b}$ at their base. We set the separation distance to be large enough ( $d>0.4$ ) so that the two filaments do not intersect.

Figure 5.

Long-term dynamics of a pair of filaments with (a,c,e)  $M_{b}=1$ and (b,d,f) $M_{b}=3$ for $d=0.7$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.49$ . (a,b) Beating patterns at four different phases ( $\unicode[STIX]{x1D719}_{1}=0,0.25,0.5,0.75$ ). (c,d) Base angles $\unicode[STIX]{x1D703}_{b}$ as a function of time for $0<t<0.3$ (to highlight transient behaviour). (e,f) Phase difference $|\unicode[STIX]{x0394}\unicode[STIX]{x1D719}|$ for $0<t<2$ . See supplementary movies 5 and 6.

Figure 5 shows the long-term behaviour of two filaments that start in a nearly anti-phase configuration; the left filament, referred to as filament 1, is initially straight and moving to the right ( $\unicode[STIX]{x1D6FC}_{1}(0)=1$ ) whereas the right filament, or filament 2, is initially moving to the left ( $\unicode[STIX]{x1D6FC}_{2}(0)=-1$ ) such that the phase difference is equal to $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=\unicode[STIX]{x1D719}_{2}(0)-\unicode[STIX]{x1D719}_{1}(0)=0.49$ . Here, the state variables $\unicode[STIX]{x1D6FC}_{1}$ and $\unicode[STIX]{x1D6FC}_{2}$ and phase variables $\unicode[STIX]{x1D719}_{1}$ and $\unicode[STIX]{x1D719}_{2}$ for filaments 1 and 2 are defined as in (4.1). This initial configuration corresponds to a small perturbation away from the anti-phase configuration for which $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.5$ . The coupling between the two filaments is due to hydrodynamic interactions only.

The two filaments exhibit anti-phase synchronization for $M_{b}=1$ and $d=0.7$ as shown in figure 5(a,c,e), whereas for $M_{b}=3$ the two filaments depart from their anti-phase initial conditions and approach in-phase synchronization as shown in figure 5(b,d,f). These modes of synchronization are quantified in figure 5(e,f). In both cases, the shapes of the filaments show no significant difference compared to those exhibited by the single filaments.

To quantify the long-term synchronization mode between the two filaments, we adapt the synchronization order parameter $Q$ proposed in Kotar et al. (Reference Kotar, Leoni, Bassetti, Lagomarsino and Cicuta2010). Namely, we let

where $T$ is the total integration time, and $T^{\ast }$ is chosen to ensure that transient behaviour is excluded. By construction, one has $-1\leqslant Q\leqslant 1$ , where $Q=-1$ describes exactly in-phase motions while $Q=1$ corresponds to exactly anti-phase motions. In the simulations, we set $T^{\ast }=15$ and $T=20$ time units, respectively; the filaments are said to be in-phase if $Q\in [-1,-0.5]$ and anti-phase if $Q\in [0.5,1]$ .

Figure 6.

Synchronization order parameter $Q$ as a function of (a,b) bending moment $M_{b}$ and (c,d) separation distance $d$ . (a,c) In-phase initial conditions are shown in blue and anti-phase in red. (b,d) Mean $\langle Q\rangle$ (dotted line) and standard deviation $\text{SD}(Q)$ (grey error bars) corresponding to 20 random initial conditions; $Q$ values are shown as blue and red dots, with colour intensity proportional to the percentage of initial conditions resulting in these values.

Figure 7.

Beating frequency $\unicode[STIX]{x1D714}$ as a function of (a,b) bending moment $M_{b}$ and (c,d) separation distance $d$ . (a,c) In-phase initial conditions are shown in blue and anti-phase in red. (b,d) Mean $\langle \unicode[STIX]{x1D714}\rangle$ (dotted line) and standard deviation $\text{SD}(\unicode[STIX]{x1D714})$ (grey error bars) corresponding to 20 random initial conditions; $\unicode[STIX]{x1D714}$ values are shown as blue and red dots, with colour intensity proportional to the percentage of initial conditions resulting in these values. In all panels, the frequency of a single filament (dashed black line) is superimposed for comparison.

We fix the separation distance between the filaments at $d=0.7$ and investigate the effect of the bending moment $M_{b}$ on the long-term synchronization modes between the two filaments. We consider in-phase and anti-phase initial conditions $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.5$ , respectively, as well as small perturbations $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.01$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.49$ away from these configurations. Figure 6(a) depicts the synchronization order parameter $Q$ versus $M_{b}$ shown in blue for $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.01$ and in red for $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.5$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.49$ . When starting at $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.5$ , the filaments always synchronize anti-phase. However, this anti-phase synchronization becomes unstable for large $M_{b}$ (dashed red line) because the filaments shift to in-phase synchronization under a small perturbation in the initial conditions ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.49$ ). On the other hand, when starting at $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.01$ , the filaments synchronize in-phase for small $M_{b}$ , shift to anti-phase synchronization as $M_{b}$ increases, and shift back to in-phase synchronization as $M_{b}$ increases further. For $0.6\leqslant M_{b}\leqslant 0.9$ and $2.1\leqslant M_{b}\leqslant 2.4$ , the filaments exhibit both stable in-phase and stable anti-phase synchronization depending on initial conditions. To better understand the sensitivity of these synchronization modes to perturbations in the initial conditions, we perform Monte Carlo simulations with initial conditions randomly chosen from a uniform distribution function $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)\in {\mathcal{U}}(-0.5,0.5)$ . Statistical results of the synchronization modes based on 20 Monte Carlo simulations are shown in figure 6(b). Dotted lines and error bars depict the mean $\langle Q\rangle$ and standard deviation SD $(Q)$ of the synchronization order parameter, respectively. Overlaid blue and red dots are the distributions of the Monte Carlo simulations, coloured in blue and red according to the emergent synchronization modes (blue for in-phase and vice versa). The colour intensity of the dots represents the fraction of simulations corresponding to a particular $Q$ value – lighter colour means fewer simulations out of 20 total number of simulations. The bistable regions where both red and blue dots coexist are consistent with the results in figure 6(a). In the bistable regions, the synchronization mode is sensitive to initial conditions.

To explore the effect of the separation distance $d$ between the filaments on the emergent synchronization modes, we fix the magnitude of the bending moment at $M_{b}=2$ and plot the synchronization order parameter $Q$ versus $d$ in figure 6(c,d). For small $d$ , all initial conditions lead to anti-phase synchronization. As $d$ increases, both in-phase and anti-phase synchronizations coexist, depending on initial conditions, and as $d$ increases further, only in-phase synchronization are observed. The Monte Carlo simulations shown in figure 6(d) are consistent with these findings. In the limit $d\rightarrow \infty$ , the two filaments maintain their initial phase difference. In other words, as $d\rightarrow \infty$ , the two filaments will take an infinitely long time to synchronize.

In figure 7, we report the values of the emergent beating frequencies for the cases considered in figure 6 and compare these values to the case of a single filament from figure 3, which we show in black dashed lines in figure 7. The beating frequencies for the pair of filaments are either faster or slower than the single filament depending on their synchronization modes: anti-phase filaments beat at lower frequencies because the two filaments ‘work against each other’, while in-phase filaments beat at higher frequencies because they ‘work together’. In particular, in anti-phase beating, the two fibres are compressing and extending fluid elements in the region between them; while for in-phase beating, the fluid and the filaments move together. In the limit $d\rightarrow \infty$ , the beating frequency for the pair of filaments converges to the beating frequency of a single filament.

Figure 8 shows the synchronization order parameter $Q$ over the parameter space $(d,M_{b})$ ; figure 8(a) shows $Q$ for nearly in-phase initial conditions $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.01$ and figure 8(b) for nearly anti-phase initial conditions $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.49$ . Figures 8(c) and 8(d) show the mean and standard deviation, respectively, of $Q$ for 20 Monte Carlo simulations with initial phase differences chosen from a uniform distribution ${\mathcal{U}}(-0.5,0.5)$ . Taken together, these results imply that the parameter space can be divided into three distinct regions: a stable anti-phase region where $\langle Q\rangle \in [0.5,1]$ and $\text{SD}(Q)<0.2$ ; a stable in-phase region where $\langle Q\rangle \in [-1,-0.5]$ and $\text{SD}(Q)<0.2$ ; and a bistable region where $\langle Q\rangle \in [-0.5,0.5]$ and $\text{SD}(Q)>0.2$ , in which the synchronization states are sensitive to the initial phase differences. The three regions are illustrated in figure 8(e).

Figure 8.

Synchronization order parameter of two filaments as a function of bending moment $M_{b}$ and separation distance $d$ for (a) nearly in-phase ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.01$ ), (b) nearly anti-phase ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.49$ ), and (c,d) 20 randomly chosen initial conditions. Here $\langle Q\rangle$ and $\text{SD}(Q)$ are the mean and standard deviation of synchronization order parameter. High $\text{SD}(Q)$ indicates sensitivity to initial conditions. (e) In-phase, anti-phase and bistable regions obtained by overlaying (c) and (d).

6 Shape-dependent oscillators

In the geometric-switch model proposed by Kotar et al. (Reference Kotar, Leoni, Bassetti, Lagomarsino and Cicuta2010) and Bruot et al. (Reference Bruot, Kotar, de Lillo, Lagomarsino and Cosentino2012), a rigid spherical particle is free to move along one direction, say the $x$ -axis, under the influence of a driving force $F$ that switches direction when the particle position reaches predefined target positions. Here, we develop a phenomenological model, based on the geometric switch oscillator, that accounts for the filament’s elasticity and shape changes in terms of a ‘lumped’ shape variable, which we denote by $\text{s}$ (not to be confused with the filament’s arclength $s$ ); see figure 1(c). We propose the coupled position–shape system of equations

The elastic ‘particle’ is subject to a configuration-dependent force $\unicode[STIX]{x1D6FC}F$ , where the magnitude $F$ is constant whereas $\unicode[STIX]{x1D6FC}$ switches between $\{-1,1\}$ as the particle position reaches a predefined oscillation amplitude $|x|=X_{s}$ . For a rigid particle, the drive $F$ is balanced by a hydrodynamic drag equal to $\unicode[STIX]{x1D709}{\dot{x}}$ , where $\unicode[STIX]{x1D709}$ is a constant (positive) drag coefficient. Elasticity introduces an internal restoring force that competes with the driving force and couples the shape $\text{s}$ to the orientation dynamics. The elastic force is modelled via a spring with stiffness coefficient $k$ that represents a ‘lumped’ elastic modulus of the filament. The drag coefficient $\unicode[STIX]{x1D709}$ also depends on shape; it should be maximum when $\text{s}=0$ , that is, when the filament is straight and moving transversely to itself, and minimum when the filament reaches its maximum deformation. It should also be symmetric under reflections from $\text{s}$ to $-\text{s}$ . We therefore set $\unicode[STIX]{x1D709}(\text{s})=\max (\unicode[STIX]{x1D709}_{o}-b\text{s}^{2},\unicode[STIX]{x1D716})$ , where the quadratic function $\unicode[STIX]{x1D709}_{o}-b\text{s}^{2}$ is maximum at and symmetric about $\text{s}=0$ and the parameter $b$ characterizes the dependence of $\unicode[STIX]{x1D709}$ on shape. The lower bound $\unicode[STIX]{x1D716}>0$ ensures that the drag coefficient $\unicode[STIX]{x1D709}$ remains positive at all time.

The shape of the filament changes under the influence of the driving force but relaxes to its original shape when it is not actuated. In (6.1), we assume that $\unicode[STIX]{x1D6FC}F$ drives the shape directly and that the shape $\text{s}$ relaxes to the original shape $\text{s}_{o}$ with constant relaxation parameter $\unicode[STIX]{x1D70F}$ . For a fixed value of $\unicode[STIX]{x1D6FC}$ , the solution to the shape equation is of the form $\text{s}=\unicode[STIX]{x1D6FC}F+(\text{s}_{o}-\unicode[STIX]{x1D6FC}F)e^{-t/\unicode[STIX]{x1D70F}}$ , where $\text{s}$ relaxes to $\text{s}_{o}$ when $F=0$ . For non-zero $F$ , the force switches sign at $\pm X_{s}$ , thus coupling position and shape. An alternative form of the shape equation in (6.1) could be written by using ${\dot{x}}$ to drive the shape dynamics instead of directly driving it by $\unicode[STIX]{x1D6FC}F$ . Then, the shape equation becomes nonlinear. We chose the linear form in (6.1) because we are mainly interested in reducing the complexity of the dynamical system, while identifying the main physical mechanisms at play.

Figure 9.

Velocity evolution ${\dot{x}}$ versus time $t$ of the single shape-dependent oscillator for: (a) shape-dependent oscillator, $b=0.5$ , $k=0.5$ ; (b) no restoring force, $b=0.5$ , $k=0$ ; (c) constant drag coefficient, $b=0$ , $k=0.5$ ; and (d) traditional oscillator, $b=0$ , $k=0$ . In all cases, $\unicode[STIX]{x1D70F}=1$ , $\unicode[STIX]{x1D709}_{o}=1$ , $X_{s}=1$ and $F=1$ .

Figure 9(a) shows the typical evolution of ${\dot{x}}$ for the shape-dependent oscillator. The dynamics resembles qualitatively the dynamics of the full filament model shown in the second row of figure 2. Specifically, after each switch, the velocity ${\dot{x}}$ first experiences a sharp decrease and remains small until the next switch. This is a joint effect of the internal restoring moment $-k\text{s}$ due to the filament elasticity and the shape-dependent drag coefficient $\unicode[STIX]{x1D709}(\text{s})$ . If the restoring moment is eliminated ( $k=0$ ), the velocity profile changes such that it first decreases then increases (figure 9 b), and if the drag coefficient is held constant ( $b=0$ ), the decrease in velocity after each switch is more gradual (figure 9 c). In the case of the standard geometric switch $b=k=0$ , the velocity remains constant after each switch because the applied force considered here is constant (figure 9 d).

We now consider two hydrodynamically coupled, shape-dependent oscillators,

Subscripts $i,j=1,2$ are the oscillator indices and $v_{j}(x_{i})={\dot{x}}_{j}/|x_{i}-x_{j}|$ (with $i\neq j$ ) is the far-field approximation of the flow velocity generated by the motion of oscillator $j$ at $x_{i}$ . The state variable $\unicode[STIX]{x1D6FC}_{1}$ switches between $\{-1,1\}$ once $|x_{1}|=X_{s}$ , while $\unicode[STIX]{x1D6FC}_{2}$ switches between $\{-1,1\}$ once $|x_{2}-d|=X_{s}$ , where $d$ is the separation distance between the centres of the two oscillator trajectories.

Figure 10.

Synchronization modes of two shape-dependent oscillators as a function of applied force $F$ and separation distance $d$ for (a) nearly in-phase ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.01$ ), (b) nearly anti-phase ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.49$ ), and (c,d) 20 randomly chosen initial conditions. Here $\langle Q\rangle$ and $\text{SD}(Q)$ are the mean and standard deviation of synchronization order parameter. High $\text{SD}(Q)$ indicates sensitivity to initial conditions. (e) In-phase, anti-phase and bistable regions obtained by overlaying (c) and (d). In all simulations, $\unicode[STIX]{x1D70F}=1$ , $b=2$ , $k=0.1$ , $X_{s}=1$ and $\unicode[STIX]{x1D709}_{o}=1$ .

If $b=k=0$ , the first equation in (6.2) is consistent with the geometric switch oscillators in Kotar et al. (Reference Kotar, Leoni, Bassetti, Lagomarsino and Cicuta2010) and Bruot et al. (Reference Bruot, Kotar, de Lillo, Lagomarsino and Cosentino2012) with one major distinction: here the applied force has constant magnitude. In their model, the applied force depends on the particle position, which determines the type of synchronization: the two oscillators synchronize anti-phase if the force magnitude decreases as the oscillator approaches $X_{s}$ and in-phase if the force magnitude increases. For constant force, the two oscillators do not synchronize. Nowhere are bistable synchronizations observed. Consistent with their findings, when the shape changes are not accounted for ( $b=k=0$ ), the model in (6.2) exhibits no synchronization – the two oscillators maintain their initial phase difference for all time. However, when shape changes are considered ( $b\neq 0$ and $k\neq 0$ ), multiple synchronization modes can arise depending on the parameter values and initial conditions, as shown in figure 10.

Figures 10(a) and 10(b) show the synchronization order parameter $Q$ as a function of the force magnitude $F$ and separation distance $d$ for two sets of initial conditions $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.01$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(0)=0.49$ , respectively. For both initial conditions, in-phase synchronization is favoured as $F$ increases, albeit for different values of $F$ . This tendency to synchronize in-phase at larger $F$ is consistent with the trend observed in the full filament model. Figures 10(c) and 10(d) show the results of 20 Monte Carlo simulations with random initial conditions taken from a uniform distribution function: a stable anti-phase region is observed for $F<0.3$ , a stable in-phase region for $F>0.7$ , and a bistable region for $0.3<F<0.7$ . The three regions are illustrated in figure 10(e).

The simple model in (6.2) captures some of the main features of the full filament model. In particular, it shows the presence of regions where in-phase and anti-phase oscillations are both stable, depending on initial conditions. This bistability is the product of the coupling between hydrodynamic interactions and shape changes. In fact, if the restoring force due to elasticity is eliminated ( $k=0$ but $b>0$ ), the oscillators always synchronize in phase. Meanwhile, if the dependence of drag on shape is eliminated ( $b=0$ but $k>0$ ), the oscillators always synchronize anti-phase. The two types of synchronization modes are observed only when the two effects of $k$ and $b$ are present. These findings imply that the two different synchronization modes observed in the simplified shape-dependent oscillators and in the full filament model are due to the interplay between elasticity, shape-dependent drag and hydrodynamic coupling.