1. Overview

Biological processes are one of many inspirations of elastohydromechanics (Shelley & Ueda Reference Shelley and Ueda2000; Schoeller et al. Reference Schoeller, Townsend, Westwood and Keaveny2021). Slender biological objects emerge in multiple contexts, motivating detailed investigation. Starting from the sub-cellular level, examples include DNA and protein folding dynamics (Goldstein & Langer Reference Goldstein and Langer1995), lipids usually forming cell walls assembling into long filaments (Rudolph, Ratna & Kahn Reference Rudolph, Ratna and Kahn1991) or microtubules helping healing by contracting wounds (Ehrlich, Grislis & Hunt Reference Ehrlich, Grislis and Hunt1977). Another large area of interest is motility – cells moving inside a fluid environment or cells inducing motion of a fluid. An iconic example of such motion is sperm cells. On closer investigation it turns out that flexibility plays an important role in their motion, and the interplay between elasticity and viscous forces causes changes to the beating pattern in response to the changing environment (Cosentino Lagomarsino, Capuani & Lowe Reference Cosentino Lagomarsino, Capuani and Lowe2003; Fauci & Dillon Reference Fauci and Dillon2006; Gaffney et al. Reference Gaffney, Gadêlha, Smith, Blake and Kirkman-Brown2011). Biological ‘optimisation’ for viscosity gradients can also be found in mucus transport inside the lungs where the correct length, stiffness and active deformation of cilia provide the necessary movement of multiple layers of fluid with varying viscosity, essential for healthy respiration (Fulford & Blake Reference Fulford and Blake1986). On larger length scales, bacterial complexes were observed joining into elongated structures exhibiting a complex dynamics because of elastohydrodynamic effects (Mendelson et al. Reference Mendelson, Thwaites, Kessler and Li1995; Goldstein, Powers & Wiggins Reference Goldstein, Powers and Wiggins1998).

Most of these examples are set in a microscale context, and thus the observed dynamics is dominated by the viscous interactions with the surrounding fluid (Lauga & Powers Reference Lauga and Powers2009).

Flexible fibres with free ends have been studied in multiple settings, including their sedimentation, both experimentally (Herzhaft & Guazzelli Reference Herzhaft and Guazzelli1999) and numerically (Li et al. Reference Li, Manikantan, Saintillan and Spagnolie2013). The free-end configuration was investigated first, because methods of producing slender filaments were already developed, and because it is of interest for both industrial applications and biological settings. Further, the one-dimensional structure provides a particularly elegant, treatable and successful way of modelling (Wiggins et al. Reference Wiggins, Riveline, Ott and Goldstein1998).

In the case of low Reynolds number flows in such settings, elastic elongated filaments can be modelled using various simulation methods, for example: the immersed boundary (IB) method (Peskin Reference Peskin2002), regularised singularity methods (Cortez, Fauci & Medovikov Reference Cortez, Fauci and Medovikov2005), bead-spring models (Kuei et al. Reference Kuei, Słowicka, Ekiel-Jeżewska, Wajnryb and Stone2015; Słowicka, Wajnryb & Ekiel-Jeżewska Reference Słowicka, Wajnryb and Ekiel-Jeżewska2015; Schoeller et al. Reference Schoeller, Townsend, Westwood and Keaveny2021) and discrete and continuous variants of slender-body theory (SBT) (Tornberg & Shelley Reference Tornberg and Shelley2004; Saintillan, Darve & Shaqfeh Reference Saintillan, Darve and Shaqfeh2005). The reduced dimensionality of the filament offers a computational advantage, which has been used in variants of the IB technique to study the whirling instability of spinning filaments (Lim & Peskin Reference Lim and Peskin2004) and hydrodynamic bundling of bacterial flagella (Lim & Peskin Reference Lim and Peskin2012). A combination of SBT and the regularised Stokeslet method has also been formulated by Cortez & Nicholas (Reference Cortez and Nicholas2012) and profitably applied e.g. to explain the motion of flagella in dinoflagellates (Hguyen et al. Reference Hguyen, Ortiz, Cortez and Fauci2011). See Nguyen, Cortez & Fauci (Reference Nguyen, Cortez and Fauci2014) for a review of this approach. On the other hand, in methods which treat the filament as fundamentally one-dimensional, such as SBT, one faces problems when the mesh along the filament is too fine (small in comparison with the reduced length scale), even when smoothing the integral kernels, as discussed by Tornberg & Shelley (Reference Tornberg and Shelley2004). Finally, any numerical scheme for the time evolution of elastic filaments must address the stiffness of the equations caused by presence of high-order spatial derivatives in the equations of motion responsible for bending rigidity. Due to the very high rate at which disturbances of small wavelength are damped, the issue of stiffness becomes even more pronounced with finer mesh sizes. Practically all numerical works to date use an implicit integration scheme, while here we present a different approach.

We focus on a different configuration – looped filaments with no free ends. The dynamics of microscale loops in viscous flows was previously analysed in the context of growing smectic-A liquid crystals, which were modelled by Shelley & Ueda (Reference Shelley and Ueda2000) using SBT, and for circular filaments with a non-zero inherent twist and bend, explored using a variant of the IB method by Lim et al. (Reference Lim, Ferent, Wang and Peskin2008) in the context of over- and underwinding of DNA leading to dynamic transitions of shapes. Our work is motivated by the experimental work of Alizadehheidari et al. (Reference Alizadehheidari2015) on circular DNA confined to nanofluidic channels (and in particular its breaking), and that of Koche et al. (Reference Koche2020) linking extrachromosomal circular DNA properties with neuroblastoma, and partly inspired by previous numerical work using bead-spring hydrodynamic models (Gruziel-Słomka et al. Reference Gruziel-Słomka, Kondratiuk, Szymczak and Ekiel-Jeżewska2019). Electrophoretic and ultracentrifugation measurements of mobility pose questions about what constitutes a flexible regime, correct values of drag coefficients and the stress distribution along the filament.

The inclusion of elasticity is necessary to analyse the aforementioned systems. It was observed experimentally for red blood cells (Jay & Canham Reference Jay and Canham1972) and numerically for flexible chains of beads (Gruziel et al. Reference Gruziel, Thyagarajan, Dietler, Stasiak, Ekiel-Jeżewska and Szymczak2018; Gruziel-Słomka et al. Reference Gruziel-Słomka, Kondratiuk, Szymczak and Ekiel-Jeżewska2019) that high flexibility leads to a change in orientation (and sometimes shape) of sedimenting objects, affecting their sedimentation speed. Independently, in the case of linear filaments, Reichert & Stark (Reference Reichert and Stark2005) observed that including elasticity can change the behaviour qualitatively when looking at bundling vs non-bundling flagella. Even for a fixed shape, a change in orientation can alter the sedimentation velocity by $25\,\%$, as shown by Tchen (Reference Tchen1954) and in further analytical solutions for sedimenting tori (Cox Reference Cox1970; Majumdar & O'Neill Reference Majumdar and O'Neill1977; Johnson & Wu Reference Johnson and Wu1979). Furthermore, as noted by Box et al. (Reference Box, Kodio, O'Kiely, Cantelli, Goriely and Vella2020) and Kodio, Goriely & Vella (Reference Kodio, Goriely and Vella2020), dynamical buckling can occur in similar settings, resulting in significant shape changes of the filaments. In bead-spring models, Gruziel-Słomka et al. (Reference Gruziel-Słomka, Kondratiuk, Szymczak and Ekiel-Jeżewska2019) observed the existence of an elasticity threshold beyond which flexible loops undergo significant changes in sedimentation dynamics. While stiff loops were seen to attain almost planar, oval shapes and sediment vertically or at an acute angle to gravity, depending on their stiffness, more flexible fibres exhibited a complex shape evolution. Our work aims to explore this stability threshold in slender-body dynamics, both analytically and in terms of numerical simulations.

In this work, we analyse the dynamics of slender elastic loops by linear stability analysis in a coupled elastohydrodynamic model, and by numerical simulations introducing a new method based on Fourier expansions. The mathematical elegance of the periodic boundary conditions allows us to simplify the theoretical analysis and gain an analytical insight into the stability question. Our results contribute to the explanation of horizontal sedimentation preference. We also derive explicit expressions for tension distribution along the filament, which comply with the work of Alizadehheidari et al. (Reference Alizadehheidari2015) on DNA loops in microfluidic channels, where typical locations of ruptures correspond to the highest tension in our model.

2. Qualitative description

We focus on a thin, inextensible, looped elastic filament, settling in a viscous fluid under gravity. The filament has a length $L$ and bending stiffness $A$. We consider the dynamics in the Stokesian regime of low Reynolds numbers, where the fluid drag forces are proportional to the local velocity of the filament. Solutions for the terminal velocity for a rigid loop were known even before the development of the SBT (Tchen Reference Tchen1954; Majumdar & O'Neill Reference Majumdar and O'Neill1977) and were tested experimentally in some cases (Amarakoon et al. Reference Amarakoon, Hussey, Good and Grimsal1982). The distribution of shear forces from the fluid acting on toroidal particles in these solutions is not uniform and thus has to be balanced by an equal and opposite force from the particle. For loops that are not perfectly stiff, these forces may partially arise from elastic deformation.

The presence of elastic forces can give rise to a complex dynamics of the sedimenting loop. To understand the stability of sedimenting circular loops, we first consider a simpler case. A classical example of beam instability under external compression has a solution known since the mid-eighteenth century (Euler Reference Euler1759). Buckling under internal forcing (a heavy column buckling under its own weight) was revisited later by Greenhill (Reference Greenhill1881). For a beam of length $L$ under the action of an external force $F$, stability results from a competition of this forcing with the stiffness of the beam. Because the bending rigidity is quantified by the product $A=EI$, with $E$ being the Young's modulus of the beam material, and $I$ being its cross-sectional moment of inertia, the relevant dimensionless quantity is $EI/FL^{2}$, capturing the stiffness to external force ratio, and there exists a critical threshold value of this quantity above which the initial shape becomes unstable.

We thus expect that a large enough compression is sufficient to destabilise an elastic filament. On the other hand, when we compute forces on a straight beam bent into the shape of a circle, we find that bending forces are balanced by the tension (negative compression, or stretching) in the beam of value $T = EI / R^{2}$, where $R$ is the radius of the circle. Hence, in the absence of the fluid, we have two mechanisms which stabilise the shape against perturbations: the negative compression rate in the beam and the forces that resist bending. The presence of a fluidic environment introduces additional hydrodynamic forces that may be responsible for a local compression of the sedimenting loop.

Uniform drag would not lead to any compression, so even for a qualitative explanation of hydrodynamic buckling it is necessary to include the dependence of the drag force density on the location on the loop. The intuition that the drag is larger in the areas where the filament is perpendicular to the flow is exemplified by a known result of local SBT: that the ratio of the drag coefficients of a slender body in the directions parallel and perpendicular to the centreline of the body is $1/2$ (Guazzelli & Morris Reference Guazzelli and Morris2012, p. 79) up to $O(\epsilon )$, where $\epsilon \ll 1$ is the slenderness parameter of the body, or the filament aspect ratio. Consider now a sedimenting loop of circular shape, as in figure 1. The local gravitational force density is uniform on the circle, but we expect the drag forces to be higher at the top and the bottom in comparison with the sides, which results in the front of the circle being hydrodynamically compressed and its back being stretched. The compression can be destabilising, leading to a spontaneous shape change at the fore side of the loop – in sedimentation this is the lowest point of the loop – which can then lead to reorientation and further deformation. For stiff loops, the global balance of hydrodynamic, elastic and gravitational forces suggests that the hydrodynamic forces are proportional to the total weight of the loop only via a geometric, dimensionless factor. Because the flow field around stable configurations of flexible loops has to be identical to that around stiff loops of the same shape, we expect this scaling to hold for flexible loops as well. Recalling the stability threshold form to be $EI/FL^{2}$, and with the force scaling as mass times gravity, this reasoning finally hints at $EI / ((\rho _L g L ) L^{2})$ being the dimensionless number governing this set-up, similarly to simpler buckling examples (here, $\rho _L$ is the fibre linear density corrected for buoyancy).

Figure 1. A three-dimensional sketch of the studied system. The direction of gravity $\boldsymbol {g_0}$, angle of inclination $\theta$ and the arc length variable $s$ are highlighted. Additionally, the colour indicates the tension in the beam, either positive (aft or trailing side) or negative valued (fore or leading side). Such a tension distribution is possible when the dimensionless gravity-to-stiffness ratio and inclination angle are large enough to cause compression due to drag anisotropy.

3. Local slender-body equations

We model the fibre as a slender elastic beam of length $L$ in a viscous fluid. To account for hydrodynamic interactions, we use a local SBT which is a far-field approximation of the Stokes flow due to an obstacle with a very small aspect ratio (Batchelor Reference Batchelor1970; Cox Reference Cox1970; Johnson Reference Johnson1980), allowing for shape parametrisation using the centreline position $\boldsymbol {x}(s)$, where $s \in [0,L]$ is the arc length.

The Stokes approximation is valid if two dimensionless constants are very small: the Reynolds and Stokes numbers, measuring the relative importance of viscous to advective and viscous to inertial terms in the Navier–Stokes equation, respectively. In this case, we neglect the inertial and time-dependent terms in the Navier–Stokes equations, and arrive at the Stokes equations describing the flow field $\boldsymbol {u}$ of an incompressible fluid with viscosity $\mu$ under external body force density $\boldsymbol {g}$

(3.1)\begin{gather} \mu \nabla^{2} \boldsymbol{u} = \boldsymbol{\nabla} p - \boldsymbol{g}, \end{gather}

(3.2)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

where $p$ denotes pressure. These equations are linear and thus admit the Green's fundamental solution, also called the Stokeslet, which reads

(3.3)\begin{equation} \boldsymbol{u}_{S}(\boldsymbol{r}) = \frac{1}{8{\rm \pi}\mu r}\left(\mathbb{I} + \frac{\boldsymbol{r}\boldsymbol{r}}{r^{2}} \right)\boldsymbol{\cdot}\boldsymbol{f}, \end{equation}

where $\mathbb {I}$ denotes a unit tensor and $r=|\boldsymbol {r}|$. The Stokeslet is associated with a point force $\boldsymbol {f}$ acting on the fluid at the origin. Notably, its derivatives are also solutions to the Stokes equations. One of particular importance is the Stokes doublet, which has a dipolar flow character (Blake & Chwang Reference Blake and Chwang1974) and decays faster ($\sim 1/r^{2}$) than the Stokeslet solution.

The SBT solves for the flow around a slender object of radius $r$ and length $L$, with a typical aspect ratio (slenderness parameter) $\epsilon = r/L\ll 1$, by approximating the force density on its surface by a distribution of Stokeslets and Stokes doublets along the centreline.

This is motivated by the idea that distributing the singularities should be sufficient to model the flow at distances large in comparison with the typical radius of the filament. Matching the ‘inner’ expansion of the flow field with the ‘outer’ flow produced by the body as a whole, and taking into account the boundary conditions on the surface of the rod, allows for expression of the centreline velocity of the filament, $\boldsymbol {u}(s)$, in terms of two linear operators $\boldsymbol {\varLambda }$ and $\boldsymbol {{K}}$ acting on the force density applied to the filament as

(3.4)\begin{equation} \boldsymbol{u}(s) ={-} \boldsymbol{\varLambda}[\boldsymbol{f}](s) - \boldsymbol{{K}}[\boldsymbol{f}](s). \end{equation}

The operators take the form

(3.5)\begin{gather} \boldsymbol{\varLambda}[\boldsymbol{f}](s) = \left[(c+1)\mathbb{I} + (c-3)\partial_s\boldsymbol{x}\,\partial_s\boldsymbol{x}\right]\boldsymbol{\cdot} \boldsymbol{f}(s), \end{gather}

(3.6)\begin{gather}\boldsymbol{{K}}[\boldsymbol{f}](s) = \int_0^{L} \left( \frac{\mathbb{I} + \hat{\boldsymbol{R}}(s,s')\hat{\boldsymbol{R}}(s,s')}{|\boldsymbol{R}(s,s')|} \boldsymbol{\cdot}\boldsymbol{f}(s') - \frac{\mathbb{I}+\partial_s\boldsymbol{x}\,\partial_s\boldsymbol{x}}{|s-s'|} \boldsymbol{\cdot}\boldsymbol{f}(s) \right) \,\textrm{d}s', \end{gather}

where $\boldsymbol {R}(s,s') = \boldsymbol {x}(s)-\boldsymbol {x}(s')$, $\hat {\boldsymbol {R}}(s,s')=\boldsymbol {R}(s,s')/|\boldsymbol {R}(s,s')|$ and $c = -2\log (\epsilon )$ is a function of the slenderness parameter. This method was initially developed by Batchelor (Reference Batchelor1970) and later improved by Cox (Reference Cox1970), Keller & Rubinow (Reference Keller and Rubinow1976) and Johnson (Reference Johnson1980). The non-local contribution ($\boldsymbol {K}[\boldsymbol {f}]$ together with the $c$-independent part of $\boldsymbol {\varLambda }[\boldsymbol {f}]$) vanishes on comparison with the local one at a rate $o(1/\log (\epsilon ))$. In this contribution, we shall take advantage of this asymptotic behaviour by neglecting $\boldsymbol {K}[\boldsymbol {f}]$ entirely. For slender fibres, neglecting the non-local term is a great simplification towards an analytical treatment of the resulting equations and leads to the local SBT, also known as resistive-force theory (RFT) introduced by Gray & Hancock (Reference Gray and Hancock1955), in which the local velocity (in a quiescent fluid) is related to the local hydrodynamic drag force on the filament, $\boldsymbol {f}_{h}$, by

(3.7)\begin{equation} \boldsymbol{u}(s) ={-}\frac{c}{8 {\rm \pi}\mu}\left(\mathbb{I} + {\partial_s\boldsymbol{x}}\,{\partial_s\boldsymbol{x}}\right)\boldsymbol{\cdot}\boldsymbol{f}_{h}(s). \end{equation}

More recently, RFT has been increasingly popular as a modelling technique in biological fluid dynamics and the analysis of the motion of slender filaments, leading, e.g. to a general qualitative agreement with experimental observations of deforming flagella (Lauga & Eloy Reference Lauga and Eloy2013), or giving insights into the buckling (De Canio, Lauga & Goldstein Reference De Canio, Lauga and Goldstein2017) and swirling instabilities (Stein et al. Reference Stein, De Canio, Lauga, Shelley and Goldstein2021) in the microtubule cytoskeleton. In fact, high-precision tracking of swimming sperm revealed that RFT can quantitatively predict the complex trajectory of a sperm cell (Friedrich et al. Reference Friedrich, Riedel-Kruse, Howard and Jülicher2010). Even in the case when the slender filaments come close together, RFT has proven to be useful in predicting their bundling behaviour (Man, Koens & Lauga Reference Man, Koens and Lauga2016). In Stokes flow, the forces acting on a suspended body balance out to zero. In our case, this involves elastic forces, gravity and hydrodynamic drag on the filament, so the no-net-force condition can be written as

(3.8)\begin{equation} \boldsymbol{f}_{el} + \boldsymbol{f}_{g} +\boldsymbol{f}_{h} = 0. \end{equation}

We model the elastic forces $\boldsymbol {f}_{el}$ according to the Euler–Bernoulli beam theory, which takes into account only the local curvature of the filament and the longitudinal tension (Tornberg & Shelley Reference Tornberg and Shelley2004; Euler, Fellmann & Mikhai Reference Euler, Fellmann and Mikhai2016). The elastic force density is then given by

(3.9)\begin{equation} \boldsymbol{f}_{el} = EI\,\partial_{ssss}\boldsymbol{x} - \partial_s (T(s) \partial_s\boldsymbol{x}), \end{equation}

where $EI$ is the flexural (bending) rigidity and $T(s)$ is the tension of the filament. The second term imposes a constant length of the filament with $T$ acting as a Lagrange multiplier.

The tension is determined by the inextensibility equation $|\partial _s\boldsymbol {x}| = 1$, which can be rewritten as a condition on the filament velocity by taking the time derivative

(3.10)\begin{equation} 0 = \tfrac{1}{2}\partial_t | \partial_s\boldsymbol{x} |^{2} = \partial_{st}\boldsymbol{x} \boldsymbol{\cdot} \partial_s\boldsymbol{x} = \partial_s\boldsymbol{u} \boldsymbol{\cdot} \partial_s\boldsymbol{x}, \end{equation}

and noting that it is satisfied initially. One problem that arises due to this treatment is the lack of a correcting mechanism in cases when the length changes slightly due to numerical errors. We implement the solution proposed by Tornberg & Shelley (Reference Tornberg and Shelley2004) by introducing a numerical stabilisation term, recasting (3.10) as

(3.11)\begin{equation} 0 = \partial_s\boldsymbol{x} \boldsymbol{\cdot} \partial_s\boldsymbol{u} - w (1 - \partial_s\boldsymbol{x} \boldsymbol{\cdot} \partial_s\boldsymbol{x}), \end{equation}

with $w$ controlling the absolute extension penalty.

The force balance condition, (3.8), governs the dynamics. The hydrodynamic force density on the filament is determined by the sum of the gravitational and elastic forces (as in (3.9)), which is then used to compute the velocity of the filament centreline via (3.7). Equation (3.10) closes the system by imposing the filament inextensibility. We now rescale these equations to arrive at a dimensionless system. Firstly, we choose the dimensionless arc length to have a domain $s \in [0, 2{\rm \pi} ]$ for additional convenience when expanding in Fourier series. This results in the characteristic length $L / 2{\rm \pi}$ scaling for $\boldsymbol {x}$. Secondly, we introduce $L^{3} / (8 {\rm \pi}^{3} EI)$ as the force scale. Finally, we rescale the time by choosing $2 \mu L^{3} / ({\rm \pi} EI c)$ as our velocity scale. This leads to equations in the form

(3.12)\begin{gather} \boldsymbol{f} ={-}\partial_{ssss}\boldsymbol{x} + (\partial_s\boldsymbol{x} \partial_s + \partial_{ss}\boldsymbol{x}) T(s) - g_0 \boldsymbol{\hat{e}_z}, \end{gather}

(3.13)\begin{gather}\boldsymbol{u} = \boldsymbol{f} + (\partial_s\boldsymbol{x} \boldsymbol{\cdot} \boldsymbol{f}) \partial_s\boldsymbol{x} , \end{gather}

(3.14)\begin{gather}0 = \partial_s\boldsymbol{x} \boldsymbol{\cdot} \partial_s\boldsymbol{u} - w (1 - \partial_s\boldsymbol{x} \boldsymbol{\cdot} \partial_s\boldsymbol{x}), \end{gather}

with a single dimensionless parameter $g_0 = L^{3} g \rho _L / (8 {\rm \pi}^{3} EI)$, analogous to (the inverse of) that used by Gruziel-Słomka et al. (Reference Gruziel-Słomka, Kondratiuk, Szymczak and Ekiel-Jeżewska2019) for the bead-spring model. Here, ${\rho _L = {\rm \pi}r^{2} (\rho _{{beam}} - \rho _{{fluid}})}$ is the mass per unit length of the fibre corrected for buoyancy, with $\rho _{beam}$ and $\rho _{fluid}$ being the densities of the beam material and the fluid, respectively. We note that $w$ is merely a numerical stabilisation constant having no influence on the solutions under exact evolution.

4. Linear stability analysis of the planar circle solution

We now use the equations of motion to study the sedimentation dynamics of looped filaments with the initial condition that they are perfectly circular and inclined at an angle $\theta _0$ to the horizontal plane, as in figure 1. A rigid circle solution (including full hydrodynamics with non-local terms) was already known previously (Tchen Reference Tchen1954; Majumdar & O'Neill Reference Majumdar and O'Neill1977; Johnson & Wu Reference Johnson and Wu1979) and serves as the starting point for our stability analysis.

Equations (3.12)–(3.14) admit a translating ($\partial _s \boldsymbol {u} = \boldsymbol {0}$) solution with a single parameter $\theta _0$ as

(4.1)\begin{gather} \boldsymbol{x}_0(s,t) = [\sin s, \cos \theta_0 \cos s, \sin \theta_0 \cos s] + \boldsymbol{u}_0 t, \end{gather}

(4.2)\begin{gather}T_0(s) = 1 + \frac{g_0 \sin \theta_0}{3} \cos s, \end{gather}

(4.3)\begin{gather}\boldsymbol{u}_0 = \left[0, \frac{g_0 \sin 2 \theta_0}{6} , \frac{ g_0 (7 - \cos 2 \theta_0)}{6}\right]. \end{gather}

Note that, for $g_0 (\sin \theta _0) / 3 > 1$, there appears an area of negative tension (compression) in the beam, also the tension is largest on the aft side of the loop, explaining the observations of Alizadehheidari et al. (Reference Alizadehheidari2015) that such flexible loops of DNA tend to break near fore or aft more frequently than in between – we propose that this effect is compounded by the higher curvature, as noted in the mentioned work.

Similarly to a sedimenting slender rod, this solution exhibits a lateral drift due to the friction anisotropy. The maximal settling angle (the angle between the sedimentation velocity and gravity) is $\gamma _{max} = \tan ^{-1}(1/(4 \sqrt {3})) \approx 8.2^{\circ }$ (as compared with $\gamma _{max} \approx 19.47^{\circ }$ for a thin rod (Guazzelli & Morris Reference Guazzelli and Morris2012, p. 83)). It is relatively small, because some parts of the circle contribute to the downwards force, while not contributing to the sideways force at any chosen angle.

We perturb the solution (4.3) by taking

(4.4)\begin{equation} \boldsymbol{x}(s) = \boldsymbol{x}_0 + \boldsymbol{\tilde{x}} + O(\boldsymbol{\tilde{x}}^{2}), \end{equation}

where a tilde over a symbol denotes the perturbation function. A general form of the perturbation turns out to be analytically intractable due to the complexity of the coupled equations for tension perturbation, so further simplifying assumptions are necessary. Here, we consider specific perturbations in the direction perpendicular to the circle's plane, as in figure 2, so that $\boldsymbol {\tilde {x}}$ is of the form

(4.5)\begin{equation} \boldsymbol{\tilde{x}} = \zeta(s) [0, -\sin \theta_0, \cos \theta_0]. \end{equation}

The presented method gives rise to two problems when trying to expand to in-plane perturbations. Firstly, taking a dot product with an in-plane vector instead of a normal vector leads to significantly more complex equations. The second complication is that in-plane perturbations are inherently two-dimensional and cannot be described by just a single scalar function. Choosing only specific normal perturbation is justified by the intuitive insight that comes from the tractable form of the resulting linear stability analysis problem.

Figure 2. Diagram of the considered initial perturbation $\zeta (s)$. Presented linear stability analysis focuses only on the perturbation in the direction perpendicular to the plane in which the unperturbed solution lies.

Assuming that the associated perturbation of tension $\tilde {T}$ is $O(\zeta )$, we neglect quadratic- and higher-order terms in $\zeta$. Then the force $\boldsymbol {f}$ in the perturbed system is of the form $\boldsymbol {f} = \boldsymbol {f}_0 + \boldsymbol {\tilde {f}}$ where

(4.6) \begin{gather} \boldsymbol{f}_0 = \frac{ g_0 \sin \theta_0}{3} \, \left[\sin 2s,\, \cos 2s \cos \theta_0,\, \frac{3}{\sin \theta_0} + \cos 2s \cos \theta_0\right], \end{gather}

(4.7) \begin{gather}\boldsymbol{\tilde{f}} ={-} \partial_{ssss} \boldsymbol{\tilde{x}} + \partial_s T_0 \, \partial_s \boldsymbol{\tilde{x}} + \partial_s \tilde{T} \, \partial_s \boldsymbol{x}_0 + \tilde{T} \, \partial_{ss} \boldsymbol{x}_0 + T_0 \, \partial_{ss} \boldsymbol{\tilde{x}}. \end{gather}

Finally, we get a linear resulting perturbation to the velocity $\boldsymbol {u} = \boldsymbol {u}_0 + \boldsymbol {\tilde {u}}$ with

(4.8)\begin{equation} \boldsymbol{\tilde{u}} = \boldsymbol{\tilde{f}} + (\,\boldsymbol{\tilde{f}} \boldsymbol{\cdot} \partial_s \boldsymbol{x}_0) \partial_s \boldsymbol{x}_0 + (\boldsymbol{f}_0 \boldsymbol{\cdot} \partial_s \boldsymbol{\tilde{x}}) \partial_s \boldsymbol{x}_0 + (\boldsymbol{f}_0 \boldsymbol{\cdot} \partial_s \boldsymbol{x}_0) \partial_s \boldsymbol{\tilde{x}}. \end{equation}

This can be put into the inextensibility condition $\partial _s\boldsymbol {u} \boldsymbol {\cdot } \partial _s\boldsymbol {x} = 0$. With $\partial _s \boldsymbol {u}_0=\boldsymbol {0}$ and ${\partial _s \boldsymbol {\tilde {u}} \boldsymbol {\cdot } \partial \tilde {\boldsymbol {x}}}$ being a second-order term, the inextensibility equation for the perturbed shape is simply

(4.9)\begin{equation} 0 = \partial_s \boldsymbol{\tilde{u}} \boldsymbol{\cdot} \partial_s \boldsymbol{x}_0. \end{equation}

It could be in principle solved for the tension perturbation. Instead, a more convenient way of proceeding is to note that

(4.10)\begin{equation} \partial_t \zeta \!=\! \boldsymbol{\tilde{u}} \boldsymbol{\cdot} [0, -\sin \theta_0, \cos \theta_0] \!=\! \partial_{ss}\zeta \!+\! \frac{g_0 \sin \theta_0}{3} \left[(\partial_s\zeta) \sin s \!+\!(\partial_{ss}\zeta) \cos s \right] \!-\! \partial_{ssss}\zeta = \mathcal{L}[\zeta], \end{equation}

which has all the information needed to analyse the evolution of perturbation in the direction of the initial perturbation. The perturbation dynamics is now governed by a single parameter only

(4.11)\begin{equation} q = \tfrac{1}{3} g_0 \sin \theta_0. \end{equation}

Equation (4.10) can be rewritten as a diffusion-like equation of the form

(4.12)\begin{equation} \partial_t \zeta ={-} \partial_{ssss}\zeta + \partial_s(T_0 \partial_s\zeta). \end{equation}

This highlights the essential role of negative tension in the development of shape instability, which takes the role of the diffusion coefficient in (4.12) and the only other term is $\partial _{ssss}\zeta$, which has an additional stabilising effect. To determine the stability of the linear partial differential equation (4.12), we examine the eigenvalues of the linear operator $\mathcal {L}$ on the right-hand side. The value of the initial tension $T_0(s) = 1 + q \cos s$ gives rise to a simple analytical form of this operator

(4.13)\begin{equation} \mathcal{L} ={-} \partial_{ssss} + \partial_{ss} + q ( (\sin s) \,\partial_s + (\cos s) \,\partial_{ss}). \end{equation}

The periodicity of $\zeta$ can be enforced by analysing $\mathcal {L}$ action on $\{\sin ks,\cos ks\}$ basis on $L^{2}(S^{1})$. We note that

(4.14)\begin{equation} \mathcal{L}[\sin ks] ={-}(k^{2}+k^{4})\sin ks + \tfrac{1}{2}kq(k+1)\sin(k-1)s + \tfrac{1}{2}kq(k-1)\sin(k+1)s, \end{equation}

and similarly for $\cos ks$. Therefore, to find the eigenvalues of $\mathcal {L}$, it is sufficient to consider linear combinations of the sine and cosine parts of the Fourier expansion separately, as $\mathcal {L}$ maps the span of either one to itself. Moreover, on each of the subspaces, the restricted maps are the same and thus have identical eigenvalues.

The operator $\mathcal {L}$ on $\operatorname {span}[\sin (kx)]$ has the following matrix representation:

(4.15)\begin{equation} \mathcal{L} = \begin{bmatrix} -2 & -3q & & & & & \phantom{\ddots}\\ 0 & -20 & -6q & & & & \phantom{\ddots}\\ & -q & -90 & -10q & & & \phantom{\ddots}\\ & & -3q & -272 & & & \phantom{\ddots}\\ & & & -6q & \ddots & -qk(k+1) /2 & \phantom{\ddots}\\ & & & & & -k^{2}-k^{4} & \phantom{\ddots}\\ \phantom{\ddots} & \phantom{\ddots} & \phantom{\ddots} & \phantom{\ddots} & \phantom{\ddots} & -qk(k-1)/2 & \ddots\\ \end{bmatrix}. \end{equation}

We can obtain approximations to $\mathcal {L}$ by truncating at a desired $n$. For a given $n$ the condition that $q$ is critical translates to $\mathcal {L}$ having one eigenvalue equal to zero, which can be expressed as $\det \mathcal {L} = 0$, which is a polynomial equation in $q$. Such equations have fast numerical solvers, allowing for computation of the critical value with high accuracy. We examined this for $n \in (1, 2, \ldots 60)$ to verify that the highest eigenvalue of $\mathcal {L}$ was determined with satisfactory precision – the convergence is extremely fast (at least exponential), as illustrated in figure 3. We find that the critical value of $q = g_0 (\sin \theta _0)/3$ is $14.56105439107$. Above this critical value the largest eigenvalue is positive, as illustrated.

Figure 3. (a) Largest eigenvalues of the linear operator $\mathcal {L}$ depending on the stability parameter $q$ (dark blue lines) describing the linearised stability problem. At $q$ approximately $14.5$ the largest eigenvalue crosses zero, which corresponds to the appearance of an unstable solution of the time-dependent equation. (b) The critical value of $q$ computed for various values of the truncation frequency $n$. For $n<4$ the behaviour of $\mathcal {L}$ is completely different, but for $n\geq 4$ the critical value of the stability parameter changes by a very small fraction. This is possible because the most unstable mode is dominated by low-frequency oscillations.

5. Numerical method

In order to verify the predictions of the theoretical model and the simplified linear stability analysis, we solve the equations of motion numerically. Because all the functions characterising the elastic loop are periodic, we represent them in the form of (truncated) Fourier series. More precisely, an approximation to a function $f(s)$ is numerically represented by a complex valued $2n$-dimensional vector $f_\alpha$, such that

(5.1)\begin{equation} f(s) \approx \sum_{\alpha ={-}n}^{n} f_\alpha \exp(\textrm{i} \alpha s). \end{equation}

For a smooth $f$, such series converge exponentially. We simulate the equations of motion by computing truncated series approximations to the position $\boldsymbol {x}$, velocity $\boldsymbol {u}$, and tension $T$ up to a fixed order $n$.

The governing equations can be coded as three affine maps (of type $\boldsymbol {x} \mapsto \boldsymbol{\mathsf{A}} \boldsymbol {x} + \boldsymbol {b}$) as follows:

(5.2)\begin{equation} T \xrightarrow[(1)]{\text{Euler{--}Bernoulli equation}} \ \boldsymbol{f} \xrightarrow[(2)]{\text{local SBT}} \ \boldsymbol{u} \xrightarrow[(3)]{\text{inextensibility equation}} \ \varepsilon. \end{equation}

These correspond to the Euler–Bernoulli theory (map $(1)$ – (3.9)) allowing for computation of the force density given the tension $T$, local SBT mapping force density $\boldsymbol {f}$ to the local velocity $u$ (map $(2)$ – (3.7)) and the inextensibility equation mapping the local velocity to local length creation, the error term $\varepsilon$, which we try to minimise in the simulation (map $(3)$ – (3.11)). Maps $(1),(2)$ are approximated by matrix equations of dimension corresponding to the truncation order $n$. This is chosen in order to keep $\boldsymbol {x}$ and $\boldsymbol {u}$ expanded to the same order. Nevertheless, it is essential to compute the inextensibility equation map $(3)$ including higher-order ($2n$) terms. These three maps are combined to obtain a relationship between the coefficients of the Fourier expansion of the tension distribution $T$ and the local filament length creation $\varepsilon$, which should be as close to zero as possible. Because there are more terms in the expansion of the local length creation than in the tension distribution, we can only attempt to make them as close to zero as possible. The combined affine map from the tension $T$ to the error term $\varepsilon$ therefore induces an overdetermined system of linear equations. These are solved for $T$ by $L^{2}$ error minimisation, which is the same as solving the ordinary least squares (OLS) problem (Fourier basis is orthonormal) with an additional restriction that $T(s)$ is real valued. It might be tempting to simplify this procedure by using the same truncation order on the error terms as on the tension expansion, leading to an exact solution for $T$ instead of an optimisation problem, but if we leave out the highly oscillatory terms in the error map $(3)$, then the solver is oblivious to the filament length increase due to oscillations with frequencies higher than $n/2$ in the tension expansion, resulting in an exponential explosion of high-frequency vibrations. Such behaviour comes from the terms in the equations where two functions are multiplied (such as $\partial _s T\, \partial _s \boldsymbol {x}$); there two terms of a given wavenumber can combine to one term with double the wavenumber. Additional care needs to be taken to ensure that the trajectories remain real valued (the error can accumulate in the complex-valued $\boldsymbol {x}$ Fourier expansions). This was achieved by projecting the solution onto the allowed subspace at each time step.

Equations (5.2) for the tension $T$ are solved at each evaluation of $\boldsymbol {u}$, and this value of tension is used to compute the velocity in an explicit integration scheme with a variable time step of the Runge–Kutta–Feldberg (fifth-order) method. This algorithm, however, is at best $O(n^{3})$ in the truncation frequency (because the OLS minimiser is $O(n^{3})$) and in reality even slower, as more degrees of freedom necessitate a decrease of the time step. For $n=6$, our implementation was running at a speed of $40$ dimensionless time units per hour for typical values of parameters on one thread of a typical 2.5 GHz processor; for $n=8$, the speed decreased to $5$ dimensionless time units per hour (giving a very crude estimate for complexity of approximately $O(n^{7})$). This makes investigations of large values of $n$ impractical. For our calculations, we choose $n=6$. We discuss this choice further in § 6.1. Most of the simulations were run with the help of GNU Parallel software (Tange Reference Tange2011).

6. Sedimentation modes

Motivated by the linear stability analysis performed above for circular loops, we now explore numerically the evolution of elastic rings starting at arbitrary inclination angles for a range of the elastohydrodynamic parameter $g_0$. We thus choose similar initial conditions to those used by Gruziel-Słomka et al. (Reference Gruziel-Słomka, Kondratiuk, Szymczak and Ekiel-Jeżewska2019). In the simulations, we observe three distinct sedimentation regimes (corresponding to terminal shapes) depending on the stiffness of the loop. They are demonstrated in figure 4 with top and in-plane views of the terminal sedimentation shapes. For very stiff loops, characterised by low values of $g_0$, marked in purple in the figure and referred to as stable, we observe stable sedimentation of a circular shape. The dynamics is then given by the translating solution of (4.3) – we see no change in the shape or the sedimentation angle. When the elasticity of the loop is increased, we observe a different terminal regime referred to as in-plane dynamics, where the shape of the loop evolves in a two-dimensional plane defined by the initial angle (cases marked in green–blue) – we only see changes of the shape, while the inclination angle of the loop remains unchanged. The circular state is then unstable and evolves into a prolate loop. When the stiffness is reduced further, for high values of $g_0$, we observe three-dimensional evolution (marked in yellow and referred to as three-dimensional dynamics) with a change in the angle of the terminal sedimentation plane with respect to the initial inclination. In the three-dimensional dynamics regime the loop leaves the initial plane essentially immediately – it does not go through dynamics similar to the in-plane regime. Deviations from a circular shape are necessary but not sufficient for a change in the inclination angle, as exemplified by the existence of the in-plane dynamics regime. For the three-dimensional dynamics regime in each simulation we observe in-plane and out-of-plane perturbations appearing spontaneously every time, and the inclination angle always changes.

Figure 4. Selection of terminal shapes for varying values of forcing (inverse stiffness) $g_0$ with an initial angle $\theta _0 = 30 ^{\circ }$ and $n = 6$ coloured by trajectory type: stable (dark purple), in-plane dynamics (green-blue), three-dimensional dynamics (yellow). The left panels $(a{,}c{,}e{,}g{,}i{,}k{,}m{,}o{,}q{,}s{,}u{,}w)$ show side views (sedimentation downwards) with a unit circle and initial tilt plane marked by dashed lines; black arrows indicate the direction of gravity. When the terminal shape aspect ratio is 1, regardless of the terminal inclination angle, the projection fits inside the unit circle. The right panels $(b{,}d{,}f{,}h{,}j{,}l{,}n{,}p{,}r{,}t{,}v{,}x)$ show the shape within the final sedimentation plane aligned with the principal axes of the loop. Lowest points on the respective loops are located on the left side of panels (b,d,f,h,j,l,n,p,r,t,v,x). Dashed unit circle is plotted for reference. Note that for sufficiently large values of $g_0$ (highly flexible filaments) the terminal tilt angle changes erratically with initial conditions.

In order to test the stability criterion derived in § 4, we plot the observed sedimentation modes depending on the dimensionless gravity $g_0$ and the initial sedimentation angle $\theta _0$ in figure 5. Clearly, the stability criterion with $q\approx 14.5$ (solid line) divides the regions of absolute stability (with purple markers) and full three-dimensional dynamics (yellow markers), with the planar evolution states in between. All instances of stable behaviour are inside the predicted stability region. The initially unstable behaviour – yellow points to the right of the stability curve – involves a complex transient evolution that finally settles on a stable configuration at a smaller inclination angle. Even though the final shape in figure 4 is planar (prolate or circular), an example of the full shape transition in figure 6 shows significant bending with a complete deviation from the initial plane and the establishment of a new ellipsoidal shape in a different plane, essentially always at an angle smaller than the initial $\theta _0$, and followed by a relaxation to the final shape. In most cases, the trajectory consists of three phases. Firstly, the loop folds in half starting with the fore side of the loop falling faster than the centre of mass and immediately after being ‘blown’ backwards by the drag force. This is related to the loop deforming towards a more prolate shape. With the increasing fore–aft distance, and the centre of mass position remaining symmetrically in the middle of the loop, it is necessary for the fore side to move faster than the centre of mass during this stage. We regard this as the primary reason for the different behaviours of flexible rods and loops. In the initial stages of motion a loop extends in plane, elongating its long axis, while a rod retains its constant length. We note, however, that this effect lasts for a very short period of time, and is present only at the early stages of the evolution. Secondly, the two lobes formed by the fold relax towards the terminal plane. The dynamics then becomes very slow and the loop attains the terminal shape within the terminal plane. Notably, this behaviour is different from the case of a free-end filament, in which the lowest part of the filament initially sediments slower than the centre of mass (Li et al. Reference Li, Manikantan, Saintillan and Spagnolie2013), leading to a different shape evolution path.

Figure 5. Sedimentation regimes showing stable translation (dark purple), in-plane dynamics (green-blue) and three-dimensional dynamics (yellow), observed in simulations for truncation frequency $n=6$. Solid line represents the linear stability analysis prediction ($q \approx 14.5$); dashed lines are empirical in-plane dynamics boundaries ($q \approx 10.5$ and $q \approx 18.5$) presented as eye guides.

Figure 6. Sample of three-dimensional time evolution dynamics. Three-dimensional rendering (a), side view (b), top view (c). Images show shape evolution when starting in an unstable equilibrium. Snapshots are taken at regular time intervals with changing colour: initial light yellow to final dark purple. The loop starts at initial angle $\theta _0 = 60^{\circ }$ subject to a gravitational force of dimensionless value $g_0 = 150$. The fore side of the loop folds upwards leading to the formation of two lobes and an ellipsoidal shape (best visible in top view). At long times, the loop converges to a near-horizontal plane to finally (after a very slow dynamics) relax to a perfect circle (not shown in the figure).

Lastly, the green–blue points span across the stability boundary in figure 5. This type of behaviour is not taken into account by our simplistic linear stability analysis, because of the assumption of a normal direction of perturbation, while here the loop stays in the initial plane. Nevertheless, this type of dynamics is observed in the vicinity of the predicted stability boundary and can be regarded as an intermediate stage between complete instability and the complete lack of shape change. A vast majority of such trajectories are bounded by $10.5< (g_0 \sin \theta _0) /3 < 18.5$. This region closely follows the stability curve in a wide range of the control parameter values and has boundaries that appear to have a similar functional form to the analytical predictions.

To further analyse the relationship between the initial and final sedimentation planes, in figure 7 we plot the initial stability parameter $q=g_0\sin \theta _0/3$ vs the final parameter (calculated as $g_0\sin \theta /3$) for all the cases investigated. The stability threshold is again marked with a solid line, together with the empirical strip of two-dimensional evolving shapes between the dashed lines. For stiff loops, at low values of $q$, the evolution does not affect the sedimentation angle, and thus we see the expected linear correlation, which persists for the semi-stable states which still remain in the initial plane. For unstable loops of high flexibility, we see no apparent correlation between the initial and the final sedimentation angles, as seen clearly in figure 7(b). Although the transition from stable to unstable dynamics is quite accurately grasped by our estimates, we note that the detailed evolution in this mode of motion is sensitive to the initial conditions, as discussed further.

Figure 7. (a) A comparison between the initial and the final values of the stability parameter in all simulations. Linear stability analysis prediction (solid line at $q \approx 14.5$) is marked together with estimated in-plane regime boundaries of $q \approx 10.5$ and $q \approx 18.5$ (dashed lines). The points are coloured by the type of trajectory: stable (dark purple), in-plane dynamics (green-blue), three-dimensional dynamics (yellow). (b) A comparison between the initial and the final values of tilt angle. The stable (darkest) and the in-plane dynamics (medium) points were shifted $1^{\circ }$ in (a,b) respectively, for clarity. The unordered scatter of the three-dimensional dynamics points (yellow) shows that the final tilt is difficult to predict, but in the vast majority of cases it is smaller than the initial tilt angle.

In figure 8, we explore the deviation from a circular shape vs the final stability parameter. To this end, we plot the aspect ratio of the final shape for all our data, defined as the ratio of the highest to the middle eigenvalue of the spatial covariance matrix. The absolutely stable loops remain circular, and once we enter the unstable region, in the vicinity of the stability criterion we see a continuous increase of the aspect ratio. For unstable loops undergoing the full three-dimensional evolution, we do not see a systematic trend, but the final aspect ratios remain within the trend seen for stiffer loops, confirming that these configurations are globally stable solutions of the loop evolution equations. In addition, we see that the planar evolution is governed solely by the stability parameter $q = g_0 \sin \theta /3$, also in the cases where the final value of $\theta$ is very small and the initial three-dimensional dynamics converges to a perfectly circular equilibrium shape, the same shape as in the case of stable sedimentation modes.

Figure 8. Aspect ratio of the terminal centreline shape plotted against the terminal value of stability parameter $g_0 (\sin \theta ) /3$. Large aspect ratios refer to highly elongated loops. A functional relationship confirms that the in-plane dynamics depends only on the stability parameter. The final shape is determined by the final sedimentation angle (and vice versa). Regime boundary between stable and the in-plane dynamics is clearly visible on this graph as the point where medium shade points reach the aspect ratio $=1$ line. When the initial stability parameter is greater than $q^{\ast }\approx 9$, loops move from a stable to an in-plane dynamics regime. Note that for the stable and in-plane dynamics trajectories (darkest and medium points) the initial and the final tilt are the same, but for the three-dimensional dynamics the terminal tilt angle is smaller than the initial one (explaining lightest points with aspect ratios of unity). Regardless of the initial configuration and parameters, it is sufficient for the terminal stability parameter to be smaller than $q^{\ast }$ to ensure the terminal configuration takes the shape of a circle. Dashed line shows the unstable branch of the solutions.

6.1. Influence of truncation order

All the presented summaries of simulations are results of a numerical scheme terminated at $t=10$ of the dimensionless time, which corresponds to a sedimentation distance of the order of $10^{3}$ in terms of loop radius for moderate values of stiffness, upon confirming no further shape evolution. This is long enough for all the simulated shapes to attain the final configuration. In most cases, all rapid changes in shape (each at time scales of approximately $0.5$) are finished when $t$ reaches approximately $3$, then the terminal angle is selected. A much slower relaxation of shape within the terminal plane follows with the characteristic time of $t\approx 1$.

In the section above, we studied the stability of initial configurations and the attraction to a stable shape for different values of the initial angle $\theta _0$ and the truncation orders $n=6$. Numerical investigation of the eigenvalues of the truncated $\mathcal {L}$ operator shows that the stability boundary should be largely independent of the truncation order $n$. Even though the initial rate at which the instability develops will be independent of $n$ (in an unstable equilibrium it is the numerical noise that initiates movement), we should expect quantitatively similar behaviour for the terminal motion of the loops.

In figure 9, we redraw the stability diagram, figure 5, for different truncation order of $n=4$ (a) and $n=8$ (b). Regardless of the initial angle and the truncation order, the stability region is correctly predicted by the simplified linear stability analysis. Moreover, when the loops are initially in an unstable position, they eventually reach the terminal angles, which are constrained by the stability region. However, the truncation order changes the exact behaviour predicted for a particular stiffness and initial sedimentation angle. Notably, for very small values of $n$ the intermediate regime of in-plane dynamics cannot be observed. It is possible that this effect appears because the loop is unable to attain configuration close enough to the terminal shape of the in-plane dynamics, due to the small number of degrees of freedom. For a larger number of truncation modes, we see that the extent of the two-dimensional evolution region decreases, with more unstable states close to the stability criterion, but nevertheless the clear division between the modes remains in place, rendering the simplistic criterion a useful tool for assessing the loop stability.

Figure 9. Comparison of regimes of sedimentation observed in simulation for truncation frequency $n=4$ (a) and $n=8$ (b). Solid line represents the boundary of the stable region as obtained from the linear stability analysis. Dashed lines are eye guides to facilitate the comparison of the two panels. Very small values of $n$ do not exhibit the intermediate regime of in-plane dynamics. However, the proposed stability criteria seem to hold sway.

Far from the stability boundary, the dynamics of the loop follows the qualitative conclusions derived from the analytical treatment, and the lower boundary (stable to in-plane dynamics) is weakly sensitive to the truncation frequency – additional simulations for $n=10$ on a restricted set of initial angles confirm that a lower stability boundary shows no further change with inclusion of higher-order terms in expansions. We note, however, that the upper boundary (in-plane to three-dimensional dynamics) is a fragile one, and the selection of one attractor over the other is sensitive to the details of numerical implementation. The discussed discrepancies can be of three origins: (i) change in stability due to the truncation of high frequencies, (ii) change in perturbation power spectrum due to the change of dimensionality or (iii) change in perturbation due to the change in numerical stability of the OLS minimiser procedure. The presented analysis of truncation order in figure 3 gives us confidence that the high-frequency modes make negligible contributions to the stability problem due to the extremely high damping by elastic $\partial _{ssss}$ terms. Therefore, our expectation is that effects (ii) and (iii) are the primary reasons for the observed differences between the smallest and largest values of the truncation order.