2. Theoretical description and parameters of the simulation

The fibre is modelled as a chain of $N$ identical spherical beads of a diameter $d$, with the centres of the consecutive beads connected by springs. We employ a finitely extensible nonlinear elastic (FENE) model (Warner Reference Warner1972; Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1977) as the stretching potential energy,

(2.1)\begin{equation} U^{FENE} ={-} \frac{1}{2} k (l_{0}-d)^{2} \sum_{i=1}^{N-1}\ln \left[ {1 - \left( \frac{l_{0}-l_{i}}{l_{0}-d} \right)^{2}} \right]. \end{equation}

Here, $k$ is the spring constant, $l_{i}=|\boldsymbol {r}_{i+1}-\boldsymbol {r}_{i}|$ is the time-dependent distance between centres of two consecutive beads at the positions $\boldsymbol {r}_{i+1}$ and $\boldsymbol {r}_{i}$ and

(2.2)\begin{equation} l_0=1.01d \end{equation}

is the length of each spring at elastic equilibrium. In our simulations, the distance $l_i(t)$ between the centres of the consecutive beads stays relatively close to the equilibrium value $l_0$ for all the times $t$: typically, $\max _{t,i}l_i(t) \lesssim 1.017$ and $\min _{t,i}l_i(t) \gtrsim 1.003$.

Each triplet of the consecutive beads $j-1, j, j+1$ is straight at the elastic equilibrium, and it resists bending with the potential energy,

(2.3)\begin{equation} U^{b}_j = \frac{A}{2 l_{0}} \beta_{j}^2, \end{equation}

where $A$ is the bending stiffness and $\beta _{j}$ is the bending angle between $\boldsymbol {r}_{i}-\boldsymbol {r}_{i-1}$ and $\boldsymbol {r}_{i+1}-\boldsymbol {r}_{i}$.

The harmonic potential is commonly used in the literature to represent the bending potential energy in molecular modelling, in particular of proteins and DNA, as described, e.g. by MacKerell et al. (Reference MacKerell1998), Storm & Nelson (Reference Storm and Nelson2003) and Frenkel & Smit (Reference Frenkel and Smit2023). For sedimenting elastic fibres, the Kratky–Porod potential is commonly used by most authors, including Schlagberger & Netz (Reference Schlagberger and Netz2005), Manghi et al. (Reference Manghi, Schlagberger, Kim and Netz2006), Marchetti et al. (Reference Marchetti, Raspa, Lindner, du Roure, Bergougnoux, Guazzelli and Duprat2018), Llopis et al. (Reference Llopis, Pagonabarraga, Cosentino Lagomarsino and Lowe2007) and Saggiorato et al. (Reference Saggiorato, Elgeti, Winkler and Gompper2015). However, Bukowicki & Ekiel-Jeżewska (Reference Bukowicki and Ekiel-Jeżewska2018) showed that sometimes for larger bending angles, the Kratky–Porod potential may lead to spurious dynamics. For highly elastic fibres, large bending angles are expected, and therefore in this work, we use the harmonic bending potential that at large bending angles binds the fibre segments much stronger than the Kratky–Porod potential. It was checked by Bukowicki & Ekiel-Jeżewska (Reference Bukowicki and Ekiel-Jeżewska2018) and Bukowicki (Reference Bukowicki2019) that the harmonic bending potential for highly elastic trumbbells does not lead to spurious dynamics. It was done by comparing with a logarithmic potential that blows up for increasing bending angles.

The bending stiffness

(2.4)\begin{equation} A = \frac{E_Y {\rm \pi}d^4}{64}=\frac{k d^{2} l_{0}}{ 16} \end{equation}

is related to the spring constant $k$ by the model of an elastic cylinder of diameter $d$ with the Young's modulus $E_{Y}$ (Bukowicki & Ekiel-Jeżewska Reference Bukowicki and Ekiel-Jeżewska2018, Reference Bukowicki and Ekiel-Jeżewska2019).

In addition to the elastic forces on a bead $i$,

(2.5)\begin{equation} \boldsymbol{F}_i^e ={-} \frac{\partial \left( U^{FENE} + \displaystyle\sum_{j=2}^{N-1} U^{b}_j \right)}{\partial \boldsymbol{r}_i}, \end{equation}

there is also a constant gravitational force that acts on each bead along the $z$-axis:

(2.6)\begin{equation} \boldsymbol{F}^g ={-} m g\, \hat{\boldsymbol{e}}_z. \end{equation}

Here, $m$ is the bead mass corrected for buoyancy, $g$ is the acceleration due to gravity and $\hat {\boldsymbol {e}}_z$ is the unit vector along the $z$-axis.

The fluid flow is assumed to obey the Stokes equations. The time-dependent positions $\boldsymbol {r}_i$ of the bead centres, $i=1,\ldots,N$, satisfy the set of the first-order ordinary differential equations,

(2.7)\begin{equation} \dot{\boldsymbol{r}}_{i} = \sum_{j=1}^{N} \boldsymbol{\mu}_{ij} \boldsymbol{\cdot} (\boldsymbol{F}_j^e + \boldsymbol{F}^g ). \end{equation}

Here, the 3-× Cartesian mobility matrix $\boldsymbol {\mu }_{ij}$ accounts for hydrodynamic interactions between all the beads and is evaluated from the multipole expansion by the precise numerical codes Hydromultipole (Cichocki, Ekiel-Jeżewska & Wajnryb Reference Cichocki, Ekiel-Jeżewska and Wajnryb1999; Ekiel-Jeżewska & Wajnryb Reference Ekiel-Jeżewska and Wajnryb2009). The multipole truncation order of 3 (Cichocki et al. Reference Cichocki, Felderhof, Hinsen, Wajnryb and Bławzdziewicz1994) was used in these computations, as previous modelling on similar systems showed that the improvement in accuracy was practically negligible when the multipole truncation order of 4 was used (Shashank et al. Reference Shashank, Melikhov and Ekiel-Jeżewska2023).

A dimensionless elasto-gravitation number $B$ is introduced to relate gravitational and bending forces acting on the fibre (Cosentino Lagomarsino et al. Reference Cosentino Lagomarsino, Pagonabarraga and Lowe2005; Llopis et al. Reference Llopis, Pagonabarraga, Cosentino Lagomarsino and Lowe2007; Saggiorato et al. Reference Saggiorato, Elgeti, Winkler and Gompper2015; Bukowicki & Ekiel-Jeżewska Reference Bukowicki and Ekiel-Jeżewska2018, Reference Bukowicki and Ekiel-Jeżewska2019; Marchetti et al. Reference Marchetti, Raspa, Lindner, du Roure, Bergougnoux, Guazzelli and Duprat2018),

(2.8)\begin{equation} B= L_0^2 G / A, \end{equation}

where $G=N m g$ is the gravity corrected for buoyancy acting on the whole fibre and $L_0 = (N-1) l_{0} + d$ is the equilibrium length of the fibre. With (2.2), the fibre aspect ratio is well-approximated as $N$.

In this study, the total number of beads varies in the range $12 \le N \le 40$, and the elasto-gravitation number is in the range $1000 \le B \le 11\,000$. Our choice is related to the results of Llopis et al. (Reference Llopis, Pagonabarraga, Cosentino Lagomarsino and Lowe2007) and Saggiorato et al. (Reference Saggiorato, Elgeti, Winkler and Gompper2015), who assumed a moderate aspect ratio $N=30$. We extended for a range of $N$ around this value. In addition, the choice of $B$ is motivated by the above references. We include the upper limits for the range of values $B\le 5000$ from Llopis et al. (Reference Llopis, Pagonabarraga, Cosentino Lagomarsino and Lowe2007), and $B\le 3000$ from Saggiorato et al. (Reference Saggiorato, Elgeti, Winkler and Gompper2015) and we also go above them to reach highly flexible fibres with much larger values of $B$.

Our choice of the initial conditions is also motivated by previous results. Cosentino Lagomarsino et al. (Reference Cosentino Lagomarsino, Pagonabarraga and Lowe2005) and Llopis et al. (Reference Llopis, Pagonabarraga, Cosentino Lagomarsino and Lowe2007) investigated an elastic fibre initially straight and horizontal, and observed the evolution restricted to a vertical plane. Saggiorato et al. (Reference Saggiorato, Elgeti, Winkler and Gompper2015) pointed out that highly elastic fibres (contrary to moderately elastic ones) are unstable to a perturbation out of this vertical plane, and demonstrated that the dynamics of highly elastic fibres with and without such perturbation are different (in-plane and out-of-plane, respectively).

Therefore, in this work, we follow Saggiorato et al. (Reference Saggiorato, Elgeti, Winkler and Gompper2015) and apply a weak noise on the initially straight and horizontal fibre in the elastic equilibrium, with $l_{i}(0)=l_0$. The amplitude of the random perturbation is equal to $0.0001 d$. The evolution of the fibre is monitored until time $10^6\tau _b/N$, with $\tau _b = {\rm \pi}\eta d^2 / ( mg )$ chosen as the time unit, with $\eta$ being the dynamic viscosity of the fluid. For almost all cases, it takes nearly half of this time before the attracting mode is reached.

3. Phase diagram of different dynamical modes

For moderately elastic fibres, a vertical mode, with a stationary, symmetric V-like (or U-like) shape located in a vertical plane, is the only mode observed (Cosentino Lagomarsino et al. Reference Cosentino Lagomarsino, Pagonabarraga and Lowe2005; Schlagberger & Netz Reference Schlagberger and Netz2005; Saggiorato et al. Reference Saggiorato, Elgeti, Winkler and Gompper2015; Marchetti et al. Reference Marchetti, Raspa, Lindner, du Roure, Bergougnoux, Guazzelli and Duprat2018). For highly elastic fibres, i.e. for large values of $B$, the tilted and rotation modes, with stationary out-of-plane shapes, were reported by Saggiorato et al. (Reference Saggiorato, Elgeti, Winkler and Gompper2015) (see figure 1 therein). Here, we confirm this finding and find new families of attracting dynamical modes (we call them crawling and rotation-crawling), using a different theoretical and numerical method. We also show that for very flexible fibres, the dynamics often become irregular.

Physically, there are two different reasons for large fibre flexibility: small ratio of the local elastic to the local gravitational forces, or large aspect ratio. In general, the dynamics are expected to depend on both these parameters. For moderately elastic fibres, the dynamics are often reported in the literature to depend only on a single parameter: the elasto-gravitation number that depends on both physical parameters. Indeed, in the phase diagram in figure 1, we observe that the transition between vertical and tilted mode corresponds to the same value of $B$ in the wide range of the aspect ratios $N$.

Figure 1. Attracting dynamical modes achieved by an elastic fibre with aspect ratio $N$ and elasto-gravitation number $B$.

However, the important finding of our manuscript is that the dynamics of highly elastic fibres depend on both physical parameters. Figure 1 is the map of the attracting dynamical modes of the highly elastic fibre, based on a systematic study of the dynamics in a wide range of the aspect ratio $N$ and the elasto-gravitation number $B$. The borders between the modes in the phase diagram in figure 1 indicate that there is no single parameter to determine the dynamics: two are needed.

It is interesting to study the transition from the tilted to the crawling mode, shown in figure 1, and extract the transition value $B_T$ of the elasto-gravitation number $B$ as a function of the transition value $N_T$ of the fibre aspect ratio $N$. Figure 2 depicts the transition points $B_T(N_T)$, which were identified employing the procedure described in Appendix A.

Figure 2. Transition elasto-gravitation number $B_{T}$ and corresponding aspect ratio $N_T$ for the transition from the tilted to the crawling mode. Data points represent values obtained from numerical simulations. The solid line depicts the fit based on (3.1). Error bars are comparable to the size of the symbol.

We expect a power law (Ekiel-Jeżewska Reference Ekiel-Jeżewska2014) dependence of the form

(3.1)\begin{equation} B_T - B_c = \alpha (N_T - N_c)^p. \end{equation}

We predict that there might exist a critical value of $N_c$ with qualitatively different dynamics for $N< N_c$ and $N>N_c$. This hypothesis is based on the dynamics of elastic loops, investigated numerically (Gruziel-Słomka et al. Reference Gruziel-Słomka, Kondratiuk, Szymczak and Ekiel-Jeżewska2019). It was shown there that for a wide range of the bending stiffness $A$, sedimenting elastic loops made of $N=15$ overlapping beads (with the aspect ratio 9.4) orient horizontally and practically do not deform, in contrast to the loops with $N\ge 16$ and the same $A$, which tend to orient vertically, stay inclined or deform significantly.

The nonlinear least-square fit to the data presented in figure 2 was performed based on (3.1). The identified parameters $B_{c} = 3250 \pm 140$, $N_{c} = 11.8 \pm 0.4$, $\alpha = 44\,000 \pm 14\,000$ and $p = -1.32 \pm 0.14$ yield a function that adequately describes the data, as illustrated in figure 2. To determine the parameters in the power law (3.1) with greater precision, an estimation of $B_{T}$ with significantly higher accuracy is needed, which would require a significantly larger number of simulations near the transition values, and also simulations with the fibres having smaller numbers of beads.

4. Basic properties of the attracting dynamical modes

4.1. Comparison of typical shapes

Typical shapes of highly elastic fibres in different sedimentation modes are shown in figure 3. In the tilted and rotation modes, the shapes do not depend on time. In the crawling and rotation-crawling modes, the shapes change in time with a period $T$.

Figure 3. Snapshots of the fibre made of $N=36$ beads in different dynamical modes, reached after a long time. Positions of the fibre beads are normalised by their diameter. (a) The notation used. (b,c) Fixed shapes in the tilted ($B=3500$) and rotation ($B=6500$) modes, respectively. (d,e) Periodic evolution of shapes in the crawling ($B=7500$) and rotation-crawling ($B=8500$) modes, respectively. The reference frames are given in (4.2a–c) for (b), (4.3) for (c), (4.4a–c) for (d) and (4.5) for (e). See also the supplementary movies available at https://doi.org/10.1017/jfm.2024.729.

To quantify properties of the sedimentation modes, it is useful to compute the moment of inertia tensor (alternatively, the gyration tensor $({1}/{N})\sum _{j=1}^{N} r'_{j\alpha } r'_{j\beta }$ is often evaluated),

(4.1)\begin{equation} {I}_{\alpha \beta}=\sum_{j=1}^{N} (\delta_{\alpha \beta} {{r}'_j}^{2} - r'_{j\alpha} r'_{j\beta}), \end{equation}

where $\alpha = x, y, z,\ \beta = x, y, z$ and $\boldsymbol {r}'_j=(r'_{jx},r'_{jy},r'_{jz})$ is the position of $j$th bead centre in the centre-of-mass reference frame, $\boldsymbol {r}'_j=\boldsymbol {r}_j-\boldsymbol {r}_{CM}$, with the centre-of-mass position $\boldsymbol {r}_{CM}=(x_{CM},y_{CM},z_{CM})$ and $r_j'=|\boldsymbol {r}'_j|$. Knowledge of the moment of inertia tensor allows evaluation of its eigenvectors and eigenvalues. A unit vector $\boldsymbol {n}$, which is parallel to the eigenvector associated with the largest eigenvalue, can be introduced. For a planar shape, $\boldsymbol {n}$ would be perpendicular to this plane. For shapes slightly out of plane, $\boldsymbol {n}$ is perpendicular to a certain ‘average plane’ and determines its orientation. The unit vector $\boldsymbol {n}$ can be expressed in spherical coordinates with the polar angle $\theta$ that it makes with the vertical axis $z$ (parallel to gravity) and the azimuthal angle $\phi$ that its projection onto $xy$-plane makes with the $x$-axis (see figure 3a). The $x$ and $y$ axes are chosen arbitrarily and are not related to the initial fibre orientation that usually changes during the transient time of the evolution. In the following, the angles $\theta$ and $\phi$ are determined as functions of time and are used to characterise the dynamical modes.

The time-dependent end-to-end vector $\boldsymbol {r}_{1N} = \boldsymbol {r}_{N} - \boldsymbol {r}_{1}$ provides more information about the fibre shape and motion in different modes. Here we are interested in the radius $\rho$ and the altitude $\Delta z$, i.e. the cylindrical coordinates of $\boldsymbol {r}_{1N}$, with the cylindrical axis parallel to gravity and the origin at the point $\boldsymbol {r}_1$. The projection of the trajectory of the tip of the end-to-end vector $\boldsymbol {r}_{1N}$ on the plane of the cylindrical coordinates $\rho$ and $\Delta z$ will be analysed in § 5 to show the difference between the modes and bifurcations they undergo.

4.2. Tilted mode

The fibre shape in the tilted mode does not change with time, and it is not planar, as shown in figure 3(b) (compare with the similar numerical results in figure 1 of Saggiorato et al. Reference Saggiorato, Elgeti, Winkler and Gompper2015). The mode name, tilted, indicates that the ‘averaged plane’ of the fibre shape is inclined with respect to gravity. In figure 3(b), the horizontal axes of the side and front planes and vertical axis correspond respectively to

(4.2a–c)\begin{equation} (x-x_{CM})/d,\quad (y-y_{CM})/d,\quad (z-z_{CM})/d. \end{equation}

More information about the fibre shape, orientation and motion is provided in figure 4 (blue dashed lines). Its orientation is fixed: the polar angle $\theta$ and the azimuthal angle $\phi$ are constant in time, as shown in figure 4(a,d). The shape is symmetric with respect to reflection in the vertical plane parallel to the side plane of the snapshot shown in figure 3(b). In particular, the fibre ends are at the same horizontal level, with $\Delta z =0$, as visible in figure 4(b). Due to the inclined shape, i.e. the non-zero inclination angle $\theta$, the fibre moves horizontally along $x$. Therefore, the top view of the centre-of-mass trajectory is a straight line, shown in figure 4(c).

Figure 4. Basic features of the regular (a–d) and irregular (e,f) modes of fibre with $N=36$ (one example of each mode). Time dependence of (a) polar angle $\theta$, (b) vertical coordinate $\Delta z$ of the fibre end-to-end vector $\boldsymbol {r}_{1N}$, and (d) azimuthal angle $\phi$. Top view (c) of the centre-of-mass trajectory (the ranges of times are $1150 \tau _{b}$ for tilted, $1545 \tau _{b}$ for rotation, $1800 \tau _{b}$ for crawling and $3700 \tau _{b}$ for rotation-crawling modes). Evolution of an irregular mode: (e) azimuthal angle $\phi$ (with the results for crawling and rotation-crawling modes shown for comparison) and ( f) top view of the centre-of-mass trajectory.

4.3. Rotation mode

The rotation mode is another mode in which the shape of the fibre does not change with time. The shape is chiral and, therefore, it rotates around the gravity direction $z$ with a certain period $T_R$. In the reference frame with the origin on the rotation axis of the fibre centre of mass, $x_{CM}=R \cos ({2 {\rm \pi}t}/{T_{R}})$ and $y_{CM}=R \sin ({2 {\rm \pi}t}/{T_{R}})$, where $R$ is the radius of rotation. Projections of the fibre shape in the reference frame rotating and translating with the fibre and the origin on the rotation axis are shown in figure 3(c). Here, the horizontal axes of the side and front planes and the vertical axis correspond to

(4.3) \begin{equation} \left(\begin{matrix} \cos\left(\tfrac{2 {\rm \pi}t}{T_{R}}\right) & \sin\left(\tfrac{2 {\rm \pi}t}{T_{R}}\right) & 0 \\ -\sin\left(\tfrac{2 {\rm \pi} t}{T_{R}}\right) & \cos\left(\tfrac{2 {\rm \pi}t}{T_{R}}\right) & 0\\ 0 & 0 & 1 \end{matrix}\right)\boldsymbol{\cdot} \left(\begin{matrix} \left(x-R\cos\left(\dfrac{2 {\rm \pi} t}{T_{R}}\right)\right)/d \\ \left(y-R \sin\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right)\right)/d \\ (z-z_{CM})/d \end{matrix}\right), \end{equation}

respectively.

More information about the fibre shape, orientation and motion is provided in figure 4 (orange dotted lines). The polar angle $\theta$ remains constant in time, as visible in figure 4(a). The azimuthal angle $\phi$ increases with a constant slope that defines the characteristic period of rotation, with $T_{R}=1545 \tau _{b}$ for the orange dotted line in figure 4(d). A non-zero, time-independent value of the difference $\Delta z$ between vertical positions of the fibre ends can be seen in figure 4(b). The fibre centre of mass follows a helical trajectory with a circular cross-section, shown in figure 4(c). Another example of the fibre motion in the rotation mode was evaluated numerically by a different method (Saggiorato et al. Reference Saggiorato, Elgeti, Winkler and Gompper2015).

4.4. Crawling mode

The crawling mode has not been reported in the literature so far. The shape of the fibre in the crawling mode does not remain fixed but changes periodically with a period $T$. The snapshots of the fibre shape, taken at time intervals equal to $T/8$, are shown in figure 3(d), where the horizontal axes of the side and front planes and the vertical axis correspond respectively to

(4.4a–c)\begin{equation} (x-x_{CM})/d,\quad y/d,\quad (z-z_{CM})/d. \end{equation}

The fibre appears to be making crawling motions with its arms, hence the name of the mode. First, for half of the period $T$, the left arm is positioned above the right one and, then, for the next half of the period $T$, the right arm becomes above the left one. Moreover, the fibre shape at time $t +T/2$ is the mirror reflection of the shape at time $t$, with the reflection plane parallel to the side plane, as seen in figure 3(d).

More information about the fibre shape, orientation, and motion is provided in figure 4 (green solid lines). It is clear already in figure 3(d) that the shape is inclined, with the inclination angle $\theta > 0$, displayed in figure 4(a). This inclination causes a lateral drift of the fibre centre of mass along $x$, with small periodic side movements along $y$, as seen in figure 4(c). The azimuthal angle $\phi$ and the difference $\Delta z$ between the vertical positions of the fibre ends are periodic functions with the same period $T= 440 \tau _{b}$ (see the solid green lines in figure 4b,d). The polar angle $\theta$ is a periodic function with the period $T/2$ (see figure 4a).

4.5. Rotation-crawling mode

The rotation-crawling mode has not been reported in the literature so far. It is a more complex mode than those discussed previously. The rotation-crawling mode seems to be a quasi-periodic motion with two incommensurable characteristic time scales, $T_R$ and $T$. The shape is chiral, and the fibre rotates around $z$ with a period $T_R$, which is the time it takes for the time-averaged azimuthal angle $\phi$ to change by $2 {\rm \pi}$. In addition, similar to the crawling mode, the shape of the fibre changes with time, allowing to introduce a characteristic time $T$ required for the fibre to return to its initial conformation after completing a full cycle of deformation. To identify the period $T$, the time dependence of either $\theta$ or $\Delta z$ is used.

The period $T_R$ of the rotation-crawling mode is typically larger than the rotation period of the rotation mode with similar values of $B$ and $N$. Therefore, the shape oscillations influence the rotational dynamics.

The sequence of snapshots of the fibre shape, taken at the time intervals $T/8$, are shown in figure 3(e). The moving frame of reference is chosen in such a way that the horizontal axes of the side and front planes and the vertical axis are specified using the following equation:

(4.5)\begin{equation} \left(\begin{matrix} \cos\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right) & \sin\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right) & 0 \\ -\sin\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right) & \cos\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right) & 0 \\ 0 & 0 & 1 \end{matrix}\right)\boldsymbol{\cdot} \left(\begin{matrix} (x-x_{CM})/d \\ \left(y-R \sin\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right)\right)/d \\ (z-z_{CM})/d \end{matrix}\right), \end{equation}

where $R$ is the radius of rotation of the time-averaged centre-of-mass trajectory $y_{CM}(x_{CM})$. Therefore, $y_{CM}$ oscillates around $R \sin ({2 {\rm \pi}t}/{T_{R}})$.

The snapshots in figure 3(e) are strikingly reminiscent of the sequence of snapshots from the crawling mode, shown in figure 3(d). In contrast to the crawling mode, for the rotation-crawling mode, the shapes at $t+T/2$ differ from the mirror reflection of shapes at time $t$.

The centre-of-mass trajectory in this mode shows the superposition of three motions: settling down, rotation around gravity and periodic sideways motion (see figure 4c). In the reference frame settling down with velocity $\dot {z}_{CM}$, two independent characteristic time scales describe this motion: a longer period $T_{R}$ that is associated with the rotation and a shorter period $T$ that characterises the smaller variations. For the example shown by red solid lines in figure 4, $T_R\approx 3690 \tau _{b}$ and $T\approx 390 \tau _{b}$. Both time scales are also seen in figure 4(d) as the azimuthal angle $\phi$ increases with the superimposed small periodic oscillations. The difference $\Delta z$ between the vertical positions of the fibre ends is also a periodic function with the period $T$, although its positive and negative amplitudes have different magnitudes, as opposed to the crawling mode (see figure 4b).

5. Time scales, asymmetry and bifurcations of the modes

In this section, we study in detail the properties of the modes. We first focus on the regular modes and determine the dependence on $B$ and $N$ of two characteristic periods: the rotation period, $T_R$, and the oscillation period, $T$.

The rotation period $T_{R}$ is the longest characteristic time that is present for rotation and rotation-crawling modes. $T_{R}$ as a function of $B$ and $N$ is shown in colours in figure 5(a). For both modes, $T_{R}$ varies in the comparable range of values. There is a clear dependence on $N$ and $B$ observed in the case of the rotation mode: $T_{R}$ decreases with an increase of $B$ for a fixed $N$, and it decreases also with growing $N$ for a fixed $B$. More flexible fibres form stationary shapes that rotate faster.

Figure 5. Characteristic features of the modes. (a) The rotation period $T_{R}$ for the rotation and rotation-crawling modes. (b) The oscillation period $T$ for the crawling and rotation-crawling modes. (c) Maximum difference between vertical positions of the fibre ends ${| \Delta z |}_{max}$ during the mode. (d) Average time $T_{100d}$ for a fibre to sediment by $z_{CM}=100d$. In (a) and (b), the white colour corresponds to the absence of rotation or oscillation, and the hatched area to the irregular mode; in (c), the white colour corresponds to $\Delta z = 0$ (exactly) for the tilted mode. In (a) $T_{R}=93\,280 \tau _{b}$ for $N=34$ and $B=5500$, and $T_{R}=20\,800 \tau _{b}$ for $N=22$ and $B=9000$.

For the rotation-crawling mode, the dependence of $T_{R}$ on $N$ and $B$ is not so simple. In most cases, a general tendency is observed, but there are exceptions. When $N=26$ is considered, $T_{R}$ for the rotation-crawling mode decreases with an increase of $B$, similar to the behaviour of $T_{R}$ for the rotation mode. This is also observed for $N=38$, though on two data points only. However, exceptionally, for $N=34$, the value of $T_{R}$ for $B=10\,000$ is much larger than for $B=8500$. In general, for the rotation-crawling mode, the dependence of $T_{R}$ on $N$ for a fixed $B$ is as for $B=8500$: a gradual increase of $T_{R}$ with increasing $N$ is observed, in contrast to the rotation mode. This difference is related to the time-dependent shape and time-dependent angular velocity of the fibre in the rotation-crawling mode. However, for $B=9000$, there is an exception from this general tendency: the value $T_{R} = 20\,800 \tau _{b}$ for $N=22$ is much larger than $T_{R}$ for a slightly larger $N=24$.

The characteristic oscillation period $T$ for the crawling behaviour as a function of $N$ and $B$ is shown in colours in figure 5(b). It is seen that, in general, $T$ increases with increasing $N$ for a fixed $B$, and $T$ decreases with increasing $B$ for a fixed $N$. But again, there are exceptions. Three regions stand out from the general tendency described previously. Two of them correspond to the exceptional behaviour of $T_R$. Namely, the value of $T$ for $N=22$ and $B=9000$ is much larger than the values of its neighbours. The values of $T$ for $N=32$, $34$ and $B=9500$, $10\,000$ are also larger than $T$ for smaller values of $B$. The third region is spanned by $30 \leq N \leq 34$ and $4500 \leq B \leq 6000$. In this range, extremely large values of $T$ are observed for the largest value $N=34$, close to the border with the rotation mode. It seems that the exceptions from the general tendency occur for values of $B$ and $N$ close to transitions between different modes.

It is useful to compare $T_R$ and $T$ to a characteristic time of the gravitational settling, chosen as the time $T_{100d}$ needed for a fibre to move along gravity the distance equal to $100d$. The time $T_{100d}$ is computed separately for each pair $B$ and $N$, based on the mean vertical centre-of-mass velocity in the respective attracting mode. In figure 5(d), we show $T_{100d}$ as a function of $N$ and $B$ with the use of colours.

It is clear that $T_{100d}$ is much smaller than the rotation period $T_{R}$ and also much smaller than the oscillation period $T$. The settling time $T_{100d}$ only weakly depends on $B$, but it systematically decreases with the increase of $N$. The increase of $N$ results in the increase of the fibre mass, but also in the increase of its drag. The competing nature of these two effects leads to a small increase of the fibre settling velocity and, therefore, a small decrease of $T_{100d}$ (for larger values of $N$, the settling velocity scales as $\ln N +C$ Adamczyk et al. Reference Adamczyk, Sadlej, Wajnryb, Ekiel-Jeżewska and Warszyński2010). Therefore, even though the fibre's mass, as well as the fibre's flexibility, change by up to a factor of three or more, figure 5(d) shows that the sedimentation time $T_{100d}$ changes by no more than 1.4 times. Recall that for a single bead, $T_{100d}=300 \tau _{b}$.

For the tilted mode, there is a clear dependence of $T_{100d}$ on $N$ and $B$: the sedimentation time $T_{100d}$ increases with increasing $B$ for a fixed $N$, but for a fixed $B$ it decreases with increasing $N$ as the fibre becomes heavier and heavier. In the case of the rotation mode, a similar decreasing tendency is observed with increasing $N$ for a fixed $B$. For this mode, on the other hand, $T_{100d}$ practically does not change with increasing $B$ for a fixed $N$. For the other sedimentation modes, the fibre shape changes over time, which significantly affects the mean sedimentation velocity and also $T_{100d}$.

More information about the time-dependent sedimentation velocity is provided in Appendix B. In figure 9, we plot the maximum sedimentation velocity and the difference between the maximum and the minimum sedimentation velocities for each mode with given values of $B$ and $N$.

The regular and irregular modes can be compared with each other by ‘an asymmetry parameter’: the maximum difference ${| \Delta z |}_{max}$ between the vertical positions of the fibre ends achieved during the observation of the mode with given values of $B$ and $N$. In figure 5(c), ${| \Delta z |}_{max}$ is normalised by $L_0$: the fibre length in the elastic equilibrium. In general, the increase in the fibre flexibility, i.e. the increase of $B$ or $N$, leads to an increase in ${| \Delta z |}_{max}/L_0$. However, due to the symmetry of the tilted mode, $\Delta z = 0$ for any values of $N$ and $B$. The rotation mode exhibits an increase in ${| \Delta z |}_{max}/L_0$ with increasing $N$ (or $B$) for a fixed $B$ (or $N$). For the crawling and rotation-crawling modes ${| \Delta z |}_{max}/L_0$ behaves similarly to the rotation mode, except the three regions of $N$ and $B$ which exhibit exceptional behaviour of $T_R$ or $T$, as described earlier.

To emphasise the complexity of the motion and shape deformation of different modes, and to highlight their differences, figure 6 shows the projection of the trajectories of the tip of the end-to-end vector $\boldsymbol {r}_{1N}$ on the plane of the cylindrical coordinates $\rho$ and $\Delta z$, for $N=36$ and $3000 \le B \le 8500$, changed by the step of 500. For the tilted mode, $| \Delta z | = 0$ and, therefore, the radial coordinate $\rho$ is constant in time, and the same is true for the length of the end-to-end vector $|\boldsymbol {r}_{1N} |$. As the shape of the fibre does not change in the rotation mode, the radial coordinate $\rho$ of the end-to-end vector does not change in time. Symmetry with respect to the reflection $\Delta z \rightarrow -\Delta z$ is present in the crawling mode when the $\rho -\Delta z$ projection of the tip of the vector follows a figure of eight: $\Delta z$ is a symmetric periodic function with a period $T$ (see also figure 4b). Finally, a distorted figure of eight is drawn on the $\rho -\Delta z$ plane by the tip of the end-to-end vector in the rotation-crawling mode, as $\Delta z$ is now a non-symmetric periodic function with a period $T$ (see also figure 4b).

Figure 6. Projection of the trajectories of the tip of the end-to-end vector $\boldsymbol {r}_{1N}$ on the plane of the cylindrical coordinates $\rho$ and $\Delta z$ for $N=36$ and $3000 \le B \le 8500$, changed by the step of $B_{step} = 500$. The colours of the curves in the online version correspond to the same colours of the modes as in figures 1 and 4. The values of $B$ are stated for the tilted ($B=3000,\ 3500$), crawling ($B=7500,\ 8000$) and rotation-crawling ($B=8500$) modes. The values $4000\le B \le 7000$ for the rotation mode correspond to the dots (orange online) under an arrow, which marks the increase of $B$.

Overall, the data analysis suggests that sedimentation modes with a non-stationary fibre shape cannot be easily described. There are several indications that this is the case. First, there exists the region $30 \leq N \leq 34$ and $4500 \leq B \leq 6000$, for which the behaviour of $T$ and ${| \Delta z |}_{max}/L_0$ differs from the behaviour for the other values of $N$ and $B$. Second, there are two other regions (with $N=22$ and $B=8500$, and $N=32$, $34$ and $B=9500$, $10\,000$), which break the monotonic dependence on $N$ and/or $B$ for $T_{R}$, $T$ and ${|\Delta z |}_{max}/L_0$. Finally, there exists the irregular mode for different, separated from each other regions of $N$ and $B$.

Therefore, in figure 7 we analyse the differences between the dynamics for $B=5500$ and different values of $N$. We show the projection of the trajectory of the tip of the end-to-end vector on the plane of $\rho$ and $\Delta z$ for $N=28,30$ (crawling mode), $N=32$ (irregular mode) and $N=34$ (rotation-crawling mode), close to the transition to the rotation mode, observed at $N=36$. Between $N=28$ and $N=30$, a period-doubling bifurcation takes place. A sequence of bifurcations between $N=30$ and $N=32$ might lead to a strange attractor and chaotic behaviour of the irregular mode. A sequence of period-halving bifurcations between $N=32$ and $N=36$ may lead to the formation of the rotation mode.

Figure 7. Bifurcations of the dynamics with the increasing aspect ratio $N$ (top-down: $N=28, 30, 32, 34$) for a fixed value $B=5500$. In each frame, points or curves for $N=36$ (rotation mode) are drawn in orange for comparison. (Left column) Cylindrical coordinate $\rho$ as a function of $\Delta z$, (centre) vertical coordinate $\Delta z$ as a function of time $t$, and (right) polar angle $\theta$ as a function of time $t$. For $N=32$, the parts of the curves corresponding to $20\,000 \le t/\tau _{b}\le 23\,000$ are highlighted in red on all three plots.

Figure 7 illustrates that the existence of the irregular mode is an inherent feature of the dynamics of highly elastic fibres. Transitions between the crawling, rotation-crawling and irregular modes, seem to be sensitive to even small changes in values of the parameters $B$ and $N$, as illustrated in Appendix C. To determine the nature of the bifurcations in the phase diagram, the Floquet stability analysis would be useful, but it goes beyond the scope of this paper.