Universal features of the shape of elastic fibres in shear flow

Abstract We present a numerical study of the dynamics of an elastic fibre in a shear flow at low Reynolds number, and seek to understand several aspects of the fibre's motion using the equations for slender-body theory coupled to the elastica. The numerical simulations are performed in the bead-spring framework including hydrodynamic interactions in two theoretical schemes: the generalized Rotne–Prager–Yamakawa model and a multipole expansion corrected for lubrication forces. In general, the two schemes yield similar results, including for the dominant scaling features of the shape that we identify. In particular, we focus on the evolution of an initially straight fibre oriented in the flow direction and show that the time scales of fibre bending, curling and rotation, which depend on its length and stiffness, determine the overall motion and evolution of the shapes. We document several characteristic time scales and curvatures representative of the shape that vary as power laws of the bending stiffness and fibre length. The numerical results are further supported by an interpretation using an elastica model.


Introduction
Physical systems that contain flexible fibres in flow are common in the processing needed to manufacture various textiles, which highlights the properties of fibrous suspensions, in biophysics and cell biology where flagella and cilia are responsible for motion and stirring of fluids and biopolymers constitute the matrix of the structural materials around † Email addresses for correspondence: pjzuk@ippt.pan.pl, mekiel@ippt.pan.pl, hastone@princeton.edu is the attractor for the majority of initial orientations. Therefore, in this paper, we focus on the analysis of fibres that perform their motion entirely in the flow-gradient plane. We use a numerical bead-spring model and the theoretical elastica model to study a single elastic fibre in a low-Reynolds-number shear flow. In particular, we perform extensive bead-spring simulations with the number of beads (fibre aspect ratio) n = 20-300 and two different models of the constitutive relations that determine the resistance of the fibre to bending, i.e. the bead-bead elastic potential energy, and two different models of hydrodynamic interactions. The parameters allow for high aspect ratio, highly flexible fibres. In addition to these bead-spring simulations, we use the elastica and slender-body descriptions of the flexible fibre deformation to rationalize the dynamics.
We characterize the dynamics evaluated numerically from the bead-spring model by identifying typical time scales of the shape deformation and the maximum curvatures that represent the flexible fibre. As one example, we identify a bending time scale intrinsic to the end of a fibre that begins to slowly bend from a straight configuration aligned with the flow direction. The displacement is caused by a transverse force at the end induced by hydrodynamic interactions caused by the rate of strain of the flow. Then, due to this small displacement, in the shear flow a tensile force builds up in the tip region, and eventually a rapid buckling of the tip takes place.
The tumbling motions of a flexible fibre initially aligned with the flow have been analysed in many previous publications, numerically and experimentally, e.g. by Forgacs & Mason (1959a), Yamamoto & Matsuoka (1993) and Skjetne et al. (1997), Lindström & Uesaka (2007), Słowicka et al. (2012Słowicka et al. ( , 2013Słowicka et al. ( , 2015Słowicka et al. ( , 2020, Harasim et al. (2013), Nguyen & Fauci (2014), Farutin et al. (2016), Liu et al. (2018) and LaGrone et al. (2019). This pattern of the dynamics, typical for elongated objects of a non-negligible thickness, is not reproduced by the inextensible infinitely thin Euler-Bernoulli beam (elastica), analysed e.g. by Audoly (2015) and Lindner & Shelley (2015). The elastica does not move out of the stationary configuration perfectly aligned with the flow. Therefore, in this paper, we introduce a modified model that accounts for the dynamics of an elastica initially aligned with the shear flow and allows it to move out of the initial position. The key idea is to extend the Euler-Bernoulli beam model by adding a point force exerted on the end beads of the fibre in the direction perpendicular to the flow. This force is caused by the shear flow, in agreement with the standard theory of the hydrodynamic interactions (Kim & Karrila 1991). Using the elastica model modified in this way, we construct an analytical solution of the early time dynamics, which is in excellent agreement with our numerical simulations.
Moreover, we identify several additional universal scaling laws of the fibre shape and dynamics from the numerical simulations and in some cases are able to rationalize the results using the elastica model. We observe that essential features of the fibre dynamics can be well understood using the parameter space of the fibre's bending stiffness and aspect ratio, which extends the concept of the elasto-viscous number.
This article is organized in the following way. We describe two bead-spring models, M 1 and M 2 , of a flexible fibre in § 2.1. We specify elastic and hydrodynamic interactions in § § 2.1.1 and 2.1.2, respectively. We explain why the fibre made of beads aligned with the flow moves out of this position in § 2.1.3. We evaluate the hydrodynamic force on the tip of the fibre aligned with the flow in § 2.1.4. We present the elastica/slender-body theory in § 2.2. A typical evolution of a flexible fibre, initially aligned with the shear flow, is shown in § 3. Evolution of shape and its maximum curvature are used to introduce typical time scales. The limits of a small and a large ratio A of the bending stiffness to the corresponding viscous stresses associated with the shear flow are discussed. The evolution of the fibre at early times is analysed in § 4. Based on the numerical simulations, in § 4.1 we demonstrate that the bending process originates from the fibre ends, and at early times only the fibre ends are deformed. We define the corresponding bending time τ b and show that it does not depend on the fibre aspect ratio n if n is large enough, and it scales as τ b ∝ A 1/3 . We also provide a scaling of τ b in the whole range of A and n. A similarity solution of small deformations and early times for the elastica is given in § 4.2, and it is used for a comparison with the numerical bead-spring simulations in § 4.3. The essential new input is the addition to the elastica model of a hydrodynamic force exerted on the fibre tip by the rate of strain of the shear flow, in a similar way as follows from the bead-spring models of the hydrodynamic interactions. In § 4.4 we demonstrate that the fibre shapes scale approximately with A 1/3 for times t ≤ τ b , even beyond the range of small deformations, and provide arguments from the elastica model. Highly bent fibres, for times t ≥ τ b , are analysed in § 5. From the numerical simulations based on the bead models M 1 and M 2 we demonstrate in § 5.1 that the maximum fibre curvature κ b2 over time is a local quantity -it does not depend on n if the fibre is long enough, and it satisfies a power-law dependence on A. An explanation for the results in terms of the elastica, and also other numerically found scaling laws, is given in § 5.2. Curling motion of a highly bent fibre is analysed with the M 1 bead model and scaling laws for the curling velocity along the flow are presented in § 5.3.
Characteristic features of the fibre dynamics and curvature in the phase space of the aspect ratio n and the bending stiffness ratio A are analysed in § 6 with the use of the bead models M 1 and M 2 . The universal similarity scaling of the maximum curvature κ b2 is provided in § 6.1. The phase space diagram of the dynamical modes is found in § 6.2. The distinction between the fibres that bend locally (for larger n and smaller A), and the fibres that bend globally (for smaller n and larger A) is demonstrated. The transition between them is shown to take place gradually in a certain range of the phase space. In contrast, within the local bending mode, there exists a sharp transition in the phase space between the fibres that straighten out along the flow while tumbling and the fibres that stay coiled. The transition is described by a simple scaling law. The final conclusions are outlined in § 7. In addition, we give the details of the theoretical and numerical descriptions of the hydrodynamic interactions between the fibre beads in appendix A and we compare the results obtained by the theoretical models M 1 and M 2 in appendix B. Finally, we discuss the universal similarity scaling and the transition between local and global bending in appendices C and D, respectively.

Model of the fibre
We consider the low-Reynolds-number motion of a neutrally buoyant elastic fibre in a fluid of viscosity μ 0 . In particular, the interaction of an elastic fibre with an external shear flow of velocity V ∞ (R) = (γ Z, 0, 0) , (2.1) where R = (X, Y, Z) andγ is the shear rate, is a nonlinear problem and many approaches have been devised to study it theoretically and numerically, e.g. bead-spring models (Warner 1972;Larson et al. 1999;Kuei et al. 2015;Słowicka et al. 2015Słowicka et al. , 2020, cylinder-hinge models (Schmid & Klingenberg 2000;Férec et al. 2009), slender-body and inextensible Euler-Bernoulli beam (elastica) approaches (Becker & Shelley 2001;Tornberg & Shelley 2004;Nazockdast et al. 2017;Liu et al. 2018   bead-spring approach for numerical simulations and the elastica model for rationalization of the numerical results (see figure 1).

Bead model
The bead-spring model illustrated in figure 1(a) describes elastic and hydrodynamic interactions between n numbered i = 1, . . . , n spherical beads of identical radii a (ith bead position is denoted as R i ). In this work, we use three different bead models M i , i = 1, 2, 3 (cf. table 1), which include combinations of two different descriptions of hydrodynamic interactions (HI), described below and in appendix A: the generalized Rotne-Prager-Yamakawa (GRPY) method Zuk et al. 2017) and the multipole method with lubrication correction (HYDROMULTIPOLE) (Cichocki, Ekiel-Jeżewska & Wajnryb 1999;Ekiel-Jeżewska & Wajnryb 2009) with two sets of constitutive laws specifying elastic interactions that are described next.
The results presented in the following sections are based on the numerical simulations from the bead models M 1 and M 2 . Both of them have the same long-distance asymptotics of the hydrodynamic interactions, and for close beads the HYDROMULTIPOLE method is more precise than the GRPY. However, the computations based on the HYDROMULTIPOLE algorithm require more time and memory than the GRPY approach. The GRPY model has been therefore used in simulations of very long fibres. For n ≤ 100 both M 1 and M 2 have been applied.
Qualitative results from the bead-spring models M 1 and M 2 are similar and in regimes of large A and n they are also the same quantitatively. A detailed comparison of the results obtained within both models is given in appendix B. The model M 3 (see table 1) is applied there to explain that some differences between the models M 1 and M 2 are related to different expressions for the bending potential energy, in agreement with the predictions of Bukowicki & Ekiel-Jeżewska (2018).

Elastic interactions
An elastic interaction potential model (constitutive laws) specifies a sum E of all bead-bead interaction energies, which are used to calculate elastic forces F i = −∇ R i E acting on each bead i. We assume that there are no elastic torques acting on the beads. For every pair of beads i, j the connector vector R ij = R j − R i points from the centre of the sphere i towards the centre of sphere j.
For the first set of constitutive laws, set 1, between the centres of every two consecutive beads i and i + 1 we impose the FENE (finitely extensible nonlinear elastic) stretching potential energy (Warner 1972) where j = i + 1, k s is the spring stiffness, R m = 0.2a is the maximum stretching from the equilibrium distance and R ij = |R ij |. Between every triplet of beads i − 1, i, i + 1 we impose a harmonic bending potential energy, where θ i and θ 0 are, respectively, the time-dependent and the equilibrium angles between connector vectors R i,i−1 and R i,i+1 , and A = EI/L 0 is the bending stiffness (per unit length), with the Young modulus E, the moment of inertia I = πa 4 /4 and the distance L 0 between the centres of the consecutive beads. Because a fibre is straight, when in equilibrium, the angle θ 0 = π. In the set 1 of the constitutive laws we assume that L 0 = 2a.
For the second set of constitutive laws, set 2, between centres of every two consecutive beads i and i + 1 we impose a harmonic stretching potential energy with j = i + 1 and the equilibrium distance L 0 between the bead centres usually close to R 0 = 2a but a bit larger. Between every triplet of beads i − 1, i, i + 1 we impose a cosine (Kratky-Porod) bending potential energy This potential energy is a widely used discrete approximation of the elastic bending stiffness, see e.g. Gauger & Stark (2006). Additionally when the GRPY model of hydrodynamic interactions is used we add the repulsive part of the Lennard-Jones potential energy (2.6) between the second nearest or further neighbour beads, where LJ determines the strength of the potential and σ LJ is the characteristic distance. We set σ LJ = 2a and truncate the 914 A31-6

Fibres in shear flow
Lennard-Jones interaction range to 2.5σ LJ . This potential acts to prevent large overlaps of the beads (comparable with 2a − R ij a for R ij < 2a). This is not necessary for the HYDROMULTIPOLE model because the lubrication forces prevent the beads from overlapping.

Hydrodynamic interactions
In this work, we study translational motion of segments of a flexible fibre. In the framework of the bead-spring modelling, the translational motion of the fibre beads is determined by the theory of hydrodynamic interactions between spherical particles. We consider n spherical particles in a fluid of viscosity μ 0 subject to an incompressible external flow V ∞ (R). We investigate the case where the Reynolds number is much smaller than unity and the quasi-steady fluid velocity V (R) and pressure p(R) are described by the Stokes equations (Oseen 1927;Kim & Karrila 1991).
The theory of hydrodynamic interactions is used to calculate the translational velocities U i of the particles, which are in turn necessary to integrate the particle trajectories. In our case the external flows are linear, and there are no torques applied to the particles. Therefore the translational velocities U i satisfy the relations, where F j is the total external force exerted on the particle j and E ∞ = (∇V ∞ + (∇V ∞ ) T )/2 denotes the rate-of-strain tensor of the external fluid flow V ∞ . For the shear flow given by (2.1), There are different methods to evaluate the translational-translational μ tt ij and translational-dipolar μ td ij mobility matrices. The most precise is the multipole expansion, corrected for lubrication, in order to speed up the convergence (Durlofsky, Brady & Bossis 1987;Cichocki et al. 1994;Sangani & Mo 1994;Cichocki et al. 1999;Ekiel-Jeżewska & Wajnryb 2009) through the inverse-power expansion in the inter-particle distance (Kim & Karrila 1991). The analytical Rotne-Prager-Yamakawa approximation is also often used (Rotne & Prager 1969).
In this work we evaluate the mobility matrices as functions of the positions of all the beads using the two methods outlined in appendix A. First, we apply the analytical Rotne-Prager-Yamakawa approximation of the translational-translational mobility μ tt ij (Rotne & Prager 1969), generalized also for the translational-dipolar mobility matrix μ td ij (Kim & Karrila 1991) and implemented in the GRPY numerical program. Second, we use the precise multipole method corrected for lubrication, implemented in the numerical code HYDROMULTIPOLE. The GRPY procedure is less precise, when particle surfaces are closer than the radius of the smaller particle, but computationally much faster than the HYDROMULTIPOLE algorithm. Both methods will be briefly outlined in appendix A.
The equations of motion for the positions R i of the beads arė with U i dependent on the positions R j of all the bead centres j = 1, . . . , n, and given by (2.7).
The equations of motion (2.9) are solved numerically with the use of dimensionless variables. We choose as a characteristic length the bead diameter 2a. The total length of the fibre at equilibrium is L, which in the case of the M 1 model is fixed to L = 2na so that the fibre aspect ratio is n. We choose as a time scale the inverse of the shear rateγ −1 and the forces are normalized with πμ 0γ (2a) 2 . The above introduces the dimensionless stretching stiffness k s = k s /(πμ 0γ (2a)), LJ = LJ /(πμ 0γ (2a) 3 ) and the bending stiffness (2.10) For the GRPY approach, L 0 = 2a. Note that for the HYDROMULTIPOLE treatment of the hydrodynamic interactions, the dimensionless bending stiffness ratio A used here is slightly different from the corresponding parameter EI/(πμ 0γ (2a) 4 ) used by Słowicka et al. (2013Słowicka et al. ( , 2015Słowicka et al. ( , 2020 and denoted there by the same letter. To adjust for this difference, all the numerical values of the bending stiffness based on the HYDROMULTIPOLE codes taken from earlier works were in this paper divided by L 0 /(2a) (typically equal to 1.02 or 1.01, see appendix B).

Why a fibre aligned with the flow moves out of this position
To answer this question, we will use (2.7) to analyse the velocities of the beads for a fibre aligned with the flow and at the elastic equilibrium. We will use the standard pairwise-additive Rotne-Prager-Yamakawa (RPY) approximation for the distinct mobility matrices μ tt ij and μ td ij with i / = j (Kim & Karrila 1991). From the geometric symmetry we can write down the tensorial form of the mobility matrices for a pair of particles i and j (Kim & Karrila 1991), where I is the unit tensor, is a third rank tensor symmetric and traceless in the first and second Cartesian components, i.e.
where the Cartesian components are labelled with the Greek letters. Within the RPY approximation, the translationaltranslational self-mobility matrix μ tt ii = 1 6πμ 0 a I (2.12) and the translational-dipolar self-mobility matrix vanishes, μ td ii = 0. Our goal now is to investigate the initial configuration, when the fibre is parallel to the flow. In this case d ij = ±ê x , with the minus sign for the beads with labels i > j. Since the fibre is at the elastic equilibrium, the external forces vanish, F j = 0, and the only contribution to velocity in the direction perpendicular to the flow comes from the translational-dipolar mobility. From (2.11b) it follows that the contribution to the velocity U i of particle i from the translational-dipolar mobility μ td ij acting on the strain tensor E ∞ 914 A31-8 respectively. Therefore, there exist contributions to the bead velocities perpendicular to the flow, and this is why the fibre moves out of the position aligned with the flow. In the next section, we will show that the largest are perpendicular velocities of the first and last beads, at the initial configuration aligned with the flow, and also later when the fibre is slightly deflected. We will also demonstrate that this effect can be considered as the result of a hydrodynamic force exerted by the shear flow on the fibre.

Hydrodynamic force acting on the tip of the fibre initially aligned with the flow
We now move on to the discussion of the hydrodynamic force exerted by the shear flow on the tip of a fibre aligned with the flow or already slightly deflected from the alignment. In the following, we are going to provide the theoretical explanation for the initial stage of the bending process in terms of the elastica, based on the assumption that a constant hydrodynamic force is exerted on the fibre end by the shear flow. In the standard use of the elastica equations the existence of such a force has not been yet taken into account.
Here, we use the general framework of the theory of hydrodynamic interactions presented in the previous sections to explain the physical origin of this force, and to estimate its value numerically (with the bead model M 1 ).
In the bead models, the tip force can be found rewriting (2.7) where μ tt ii is the translational self-mobility matrix, and defining the hydrodynamic force acting on bead i as The dimensionless form is Our goal is to investigate F Hi at the early stage of the evolution, when the fibre, initially aligned with the flow, slowly moves out of this configuration, but still remains almost parallel to the flow. We will now show that, for the fibre almost aligned with the flow, the hydrodynamic forces F Hi , defined by (2.15), are almost perpendicular to the flow direction e x . We will also provide some theoretical arguments why the value of F Hi is the largest at the ends of the fibre.
The hydrodynamic forces F Hi given by (2.15) are proportional to the shear rateγ . Moreover, the force F Hi is perpendicular to the flow and parallel to the z direction of the flow gradient, F Hi ≈ê z · F Hiêz . Therefore, they displace the fibre beads away from the position aligned with the flow. From the explicit expressions for the functions C and D in (2.11b), given e.g. by Kim & Karrila (1991), it follows that for the first beadê z · F H1 > 0 and for the last beadê z · F Hn < 0. This means that, owing to the hydrodynamic forces (2.15), the fibre follows the rotational component of the shear flow. It is also known that C ∝ R −2 ij and D ∝ R −4 ij , see e.g. Kim & Karrila (1991). Therefore, the major contribution to F Hi comes from relatively close beads j. Additionally, μ td ij is antisymmetric in d ij , which means that the terms in (2.15) corresponding to equally distant left and right neighbours will cancel. Therefore, the total force F Hi is close to zero for i in the middle part of the fibre, and it increases when i is closer to the fibre ends. For longer fibres, the force F Hi is non-negligible only for i close to one of the fibre ends, and it only weakly depends on the total fibre length because it comes from unbalanced local interactions between the bead i and close beads j.
To evaluate f Hi numerically, we use the pairwise-additive GRPY approximation for the mobility matrices. As argued above, in the stage when fibre is only slightly deflected from the straight line, at leading order, f Hi is directed alongê z . In figure 2(a) we plot the dimensionless hydrodynamic forceê z · f Hi as a function of the bead label i for three different fibre lengths n. It is clear that the force is well localized close to the fibre ends.
The orientation of f Hi follows the rotational component of the shear flow. As the fibre gets longer, the force is more localized. Regardless of the fibre length, the end beads support the largest forces, an order of magnitude larger than the forces acting on the other beads. The magnitude of the force acting on the first bead, f H1 ·ê z , initially changes non-monotonically as a function of n (see figure 2b), until it reaches a limiting value f H ≈ 0.16. Indeed, we observe a localized, length-independent tip force perpendicular to the flow. We will use this observation later to construct a modified elastica model, applicable for a fibre initially aligned with the flow. Now it is time to present the standard Euler-Bernoulli beam, based on the local slender-body theory.

The elastica and local slender-body theory
To rationalize the results of numerical simulations from the bead-spring simulations the inextensible elastica model (Duprat & Stone 2015;Lindner & Shelley 2015) is used with the local slender-body theory (SBT) (Cox 1970;Keller & Rubinow 1976;Johnson 1980) to account for the drag forces acting along the fibre. Within the local SBT, in contrast to the bead models, the full long-ranged hydrodynamic interactions are not incorporated, nor is the finite but small thickness of the filament. The last feature is especially important for fibres that are aligned with the flow, as will be discussed in detail later. Similarly as in the bead model, for the elastica we also neglect Brownian motion and buoyancy forces. The fibre has a circular cross-section of radius a and length L = 2na where = a/L = 1/2n 1. We denote as R(S, T) the dimensional position of a fibre segment at the arc length S at time T. The equation of motion of each filament segment as a result of the applied elastic force density F (S, T) per unit length, under the steady undisturbed flow V ∞ can be expressed using the SBT (Cox 1970;Duprat & Stone 2015;Lindner & Shelley 2015 whereṘ = ∂R/∂T, R S = ∂R/∂S and the relative motion of the filament is obtained by applying the mobility tensor, proportional to the anisotropic tensor (I + R S R S ), to the elastic force F (S, T) on the fibre. Here, we consider shear flow For the elastic fibre we use the notation illustrated in figure 1(b), i.e.ê n denotes a unit vector normal to the fibre in the shear plane andê s denotes a unit vector tangent to the fibre. The inextensibility condition |R S | = 1 results inê s = R S and implies the Frenet formulas ∂ê s /∂S = Kê n , ∂ê n /∂S = −Kê s , where K is the local curvature and we have assumed that the fibre shape is planar. In the elastica model the elastic forces acting on the unit segment of the fibre are (Audoly & Pomeau 2000;Audoly 2015) where Σ(S, T) is the tension in the filament (satisfying inextensibility), E is the Young modulus and I is the moment of inertia (I = πa 4 /4), as earlier. Alternatively, the force density per unit length can be expressed as F (S, T) = −EIR SSSS + (Tê s ) S , see e.g. Tornberg & Shelley (2004) and Lindner & Shelley (2015). It is easy to check that With the use of the Frenet formulas it is convenient to write separately the equations of motion in the normal and tangential directions, respectively, We write the dimensionless form (lowercase symbols) of (2.19) by expressing length in the units of 2a and time in the unitsγ −1 , as in § 2.1, to find is a dimensionless compliance. Using the same normalization of EI as for the bead model, EI/(πμ 0γ L 0 (2a) 3 ), we can formally write η = 4(2a)/AL 0 ln( −1 ). A physical comparison between the elastica and bead models will be presented in § 4.3.

A typical bead model simulation
The dimensionless stretching stiffness is fixed to a large value (k s = 2000 in the M 1 model and k s = 1000 in the M 2 model) so that the fibre is close to inextensible. In M 1 , the equilibrium distance between the bead centres corresponds to the touching beads, L 0 = 2a and the dimensionless Lennard-Jones potential coefficient LJ = 5 allows only slight overlaps. In M 2 , lubrication interactions between close particle surfaces prevent overlaps. The equilibrium distance L 0 between the bead centres has to be a bit larger than the bead diameter 2a; here we choose L 0 /(2a) = 1.02. Sensitivity of the M 2 model to the choice of L 0 has been discussed by Słowicka et al. (2015Słowicka et al. ( , 2020. We focus on the fibre dynamics under the influence of the dimensionless bending stiffness A and the number of beads n, indicating the fibre's aspect ratio. The typical shape of a fibre during the evolution is presented in figure 3(a).The simulations (based on the M 1 model) start from a stretched fibre aligned in the flow direction. First, we observe a slow deflection of the fibre tips up to time approximately 30. Later, until the time 35, rapid bending of the tip occurs. Next, a curling motion appears, with the maximum curvature moving to the central part of the fibre, and a typical shape is shown for t = 47. After the kinked parts of the fibre pass over each other (around time 62), the fibre rapidly straightens to a position slightly tilted from the x direction at time 66, after which the fibre slowly stretches and aligns in the x direction until the end beads reach the same z coordinate at time 141.
To characterize the deformation of a fibre, an informative observable is the maximum local curvature κ(t) taken over the fibre length at every time instant, where, similarly to the elastica model, we use lowercase symbols for the dimensionless quantities (see § 2.1). At every time t we inscribe a circle of radius r i−1,i,i+1 (t) on the bead centres (3.1) A typical profile of κ(t) obtained from the simulations is shown in figure 3(b). We identify two characteristic bending curvatures κ b1 , the value of the first plateau, and κ b2 , the maximum value over time. To characterize the shape changes, we introduce characteristic time scales: τ b , the bending time, then τ c , the curling time, and τ s , the stretching time, as indicated in figure 3. Initially, for an almost straight fibre, κ(t) is close to 0. The rapid rise in curvature (starting around t = 30) is connected with rapid bending of the ends until a characteristic curvature κ b1 is reached. We define the time scale τ b as the time needed for a fibre to reach half of its maximum curvature κ b2 starting from a straight fibre.
After the rapid bending, a curling motion occurs. We observe the end beads passing above each other (having the same x coordinate) at a flipping time τ f = 47 (τ f was called the flipping time by Słowicka et al. (2013Słowicka et al. ( , 2015 and is used there to characterize the tumbling dynamics in shear and Poiseuille flows), then the kinks visible in figure 3(a) pass each other. We identify that the last event happens approximately at time τ b2 = 61.4, when the curvature increases to a maximum value κ b2 . Next, there is a rapid decrease of the fibre curvature. In particular, at a turning time τ m = 62.8 the middle beads have the same x coordinate. Later, we observe a rapid relaxation to an almost straight fibre (here at t = 66). We define time scale τ c as the time from the moment τ b when fibre reaches κ b2 /2 for the first time until it reaches κ b2 /2 again after passing the peak of curvature κ b2 . After rapid relaxation, the fibre is close to straight but tilted from the flow direction. The stretching time scale τ s is evaluated from the time of passing κ b2 /2 for the second time until the fibre ends are aligned with the flow direction again (here at time 141). Then, the motion approximately repeats itself periodically although small changes in the times identified in figure 3 are possible. Therefore, the sum τ = τ b + τ c + τ s is the tumbling time scale defined as the half-period of the motion and analysed by Słowicka et al. (2015Słowicka et al. ( , 2020, with typically small variations between the first tumbling and the tumblings observed at long times. With the definitions of τ b , τ c and τ s we seek to capture the time scales of the slow changes between the (much shorter) steep increase and decrease in curvature, which we consider negligible in comparison to τ b , τ c , τ s . Thus, the precise definitions of transitions points between τ b , τ c and τ s can be chosen in a different way and should not have a large influence on the analysis.
We show the changes in the dynamics for different choices of n and A in figure 4. For a small A and n large enough, the end of the fibre bends multiple times and never returns to the straight state again (figure 4a,b), in which case τ c and τ s are not defined. Nevertheless κ remains approximately constant throughout most of the process. A similar qualitative picture was observed with different experimental (Forgacs & Mason 1959b;Harasim et al. 2013;Liu et al. 2018) and numerical (Lindström & Uesaka 2007;Nguyen & Fauci 2014;Liu et al. 2018;LaGrone et al. 2019) methods. The other limit is when, for a small n, A is increased to the point when the fibre bends globally along the whole length.
In the following, we will first analyse the dynamics for 0 ≤ t ≤ τ b when the bending process originates ( § 4), and next we will study the shape evolution in the time range of large deformation, τ b ≤ t ≤ τ b + τ c ( § 5).

Bending time from numerical simulations
In addition to bending, when in a shear flow, a fibre undergoes rotation (it tumbles). It is instructive to compare the bending time τ b , the curling time τ c and the stretching time τ s (see figure 3) with two indicators of a fibre's rotational motion: the flipping time  τ f and the turning time τ m (see figure 5(a) and the insets indicating shapes for τ f and τ m ). We find for A ∈ [1, 10 000) that τ f < τ m , which shows that the ends of the flexible fibre pass above each other earlier than the middle of the fibre rotates. As A increases, both τ m and τ f acquire the interpretation of the half-tumbling time τ/2 or the quarter period T J /4 of the Jeffery orbit (Jeffery 1922), which is understood here as a periodic motion of a certain 'effective' rigid elongated object in shear flow (Słowicka et al. , 2020 with the same period T J = 2τ . The period T J of a Jeffery orbit is approximately proportional to the length of a fibre consisting of n beads (Jeffery 1922;Kim & Karrila 1991;Dhont & Briels 2007;Graham 2018).
For fibres that are very flexible or long enough (e.g. A = 10 and n = 300), τ c and τ s are not defined, because the fibre does not straighten out again , bending multiple times if n is sufficiently large (figure 4a). In the limit of very stiff fibres, all the time scales defined above, τ b , τ f , τ m and τ c + τ s , converge to T J /4, as shown in figure 5(a) for large values of A. More details about the relation between different time scales can be found in appendix B.
The dynamics of bending changes systematically as a function of A. In figures 5(b) and 5(c) we show τ b (see figure 3) as a function of A for different n, using the models M 1 and M 2 , respectively. Three regimes are visible. First, in the small A regime (A 10), the bending dynamics is dominated by large bending angles close to the excluded volume limit. Second, for intermediate A, τ b does not depend on n and it follows a single power law τ b ∝ A 1/3 . Third, in the regime of large A, the bending time systematically deviates from the power law and saturates at a constant value, where larger n have larger limiting Therefore, in figures 5(d) and 5(e) we replot the data from figures 5(b) and 5(c), respectively, using the rescaled bending time, τ b /n. To obtain the universal scaling, we also need to rescale the bending stiffness as A/n 3 . Indeed, after such rescaling, we observe an almost universal curve in the whole range of values of n and A.

4.2.
A similarity solution at early times for the elastica In figure 5(b-e) we show the dependence τ b ∝ A 1/3 , which can be argued with the help of the elastica model, as we will demonstrate in this and the next section. We observe numerically that, in the power-law regime, the bending time does not depend on the fibre length n, which suggests an analysis based on the model of a very long fibre, initially aligned with the flow, with a tip positioned at S = 0. We assume small deflections from the straight line R = Sê x + U(S, T)ê z , which leads to K = U SS . Further, we assume that because of small deflections,ê s =ê x andê n =ê z . Under these assumptions, we rewrite the dimensional linearized equations (2.19) with the Dirac delta δ(S) in the additional term that introduces the hydrodynamic force F E = O(1), perpendicular to the fibre axis, acting on the tip of the fibre at S = 0. This force results from the hydrodynamic interactions of the fibre beads in response to the shear flow (see § § 4.3 and 2.1.4). Alternatively to the delta term, the constant tip force can be formally written as a boundary condition, U SSS (S = 0, T)=F E /(EI). We will use this approach to write (4.1) in the dimensionless form, corresponding to (2.20), In order to solve (4.2a) we apply boundary conditions We seek a similarity solution (Barenblatt 1996;Duprat & Stone 2015;Eggers & Fontelos 2015) and account for the forcing as which leads to the equation with the boundary conditions The solution can be expressed with the help of special (hypergeometric) functions (e.g. use Mathematica) and the gamma Γ (·) function The function U(χ) is shown in figure 6. From U(χ) we calculate (4.10) This result can be expanded around s = 0 and, in particular at s = 0, the end moves according to . (4.12) The complete solution of the similarity ansatz has the magnitude of the deflection u(0, t) ∝ (ηt 3 ) 1/4 . The length scale on which the deflection occurs is s ∝ (t/η) 1/4 ∝ (u(0, t)/η) 1/3 . In time the 'height' of the deflection grows more rapidly than its 'width', which makes the tip more and more steep. The bending stiffness has an opposite effect and makes the deformation less steep with increasing A ∝ η −1 .

Comparing the numerical simulations with the similarity solution
In this section, we present results from the numerical simulations based on the model M 1 and compare them with the predictions from the elastica model. According to the elastica similarity solution, the fibres have features of the shape that follow the scaling laws with t and η presented above. Therefore, we analyse the z coordinate of the relative position of the first bead at time t with respect to its initial position, z 1 (t) =ê z · (R 1 (t) − R 1 (0)), which is calculated from the bead-spring simulations M 1 (figure 7a).We show the same data with the rescaled time t/A 1/3 in figure 7(b). This scaling is suggested by the elastica model, if we identify the height z 1 (t) of the fibre end with the deflection of the elastica tip u(0, t), given by (4.12), and we remember that η ∝ 1/A. We also fit a straight line to the numerical values of log 10 z 1 (t) as a function of log 10 (t/A 1/3 ), in the linear region log 10 (t/A 1/3 ) < −1, where deformations are still small and no deviations from the power law are observed. While fitting, we used data from all the simulations where n ≥ 60 and A ≥ 50. The calculated slope 0.787 is very close to 3/4 theoretically predicted from the dynamics of the elastica. In order to further compare simulations with the theoretical results we will use the best fit of the tip height in the form z 1 (t) = C(t/A 1/3 ) 3/4 , which is suggested by the elastica, that results with C ≈ 10 −0.81 .
In figure 7(c) we present the numerical shapes of the fibres with different A and n taken at different times but still within the range of the similarity solution, with t/A 1/3 10 0.6 ≈ 4. The ends of these shapes can be to a certain extend superimposed onto each other by scaling the coordinates asx = x/(tA) 1/4 andz = z/(t 3 /A) 1/4 , respectively, in accord with predictions from the elastica, and translating by a shift x 0 , which is different for each fibre, as shown in figure 7(d). The rescaled shapes are plotted together with two plots of the similarity solution as a function a k U(xb k ), k = I, II, which correspond to our two different approaches to compare the hydrodynamic forces, f H1 and f E , exerted on the fibre tip in the bead ( § 2.1) and elastica ( § 4.2) models, respectively. (Actually, we will be comparing the limiting value f H for n → ∞ rather than f H1 .) In both approaches, we assume the identical tip heights z 1 (t) and u(0, t) in the bead and elastica models, respectively. In the first approach (I), both forces are assumed to be the same, f H = f E . Therefore, F is related to f H by (4.4). On the other hand, η ≈ 4/(A ln −1 ), as shown below (2.21), and (4.12) links F to the height u(0, t) = z 1 (t) of the fibre. Using the numerical fit of z 1 (t), shown in figure 7(b), we find F, and from this result, using (4.4) we determine the magnitude f E of the dimensionless tip force in the elastica model. The numerically obtained shapes are superimposed onto theoretically calculated shapes a k U(xb k ), k = I, II, with U given by (4.9) and the coefficients a k , b k given in terms of F and η which result from approaches I and II to compare between the bead model and the elastica, given in (4.13)-(4.15) and (4.16), (4.17), respectively.
(4.15) From (4.14) we find −1 ≈ 9, which is rather far from the typical aspect ratios used in the numerical simulations. We use the above values to compare shape of the fibre made of beads with the elastica. Starting from (4.6) we find that a I = 0.1 and b I = 1.16 and plot the corresponding elastica shape (solid line) in figure 7(d). In the second approach (II), we assume that velocities of the fibre segments are the same in the bead and elastica models. In this way, mobilities times forces are equal to each other. The mobility for the elastica comes from the SBT, while in the equations of motion for a fibre made of beads, there appear the single-sphere Stokes mobility. Therefore in approach (II), we match the elastica and bead quantities as follows: which is not very far from the numerical value f H = 0.16 obtained for large n. In this approach the parameter −1 is not used at all. With the use of this set of values and (4.6) we find that a II = 0.1 and b II = 1.32. We plot the corresponding elastica shape (dotted line) in figure 7(d).

At early times beyond small deformations
In this section, we will focus on an early phase of bending for times t τ b . We will analyse the simulations with the bead model M 1 and compare them with the scaling laws following from the elastica model. In the early phase, a flexible fibre aligned with shear flow slowly starts to bend its ends while the middle part of the fibre remains straight. The characteristic length scale of the deformed fibre segments at both ends remains approximately constant in time until a significant, rapid bending is developed at the bending time τ b , associated with a fast increase of both the local curvature κ (figure 3b) and the tip deflection z 1 (t) along z ( figure 7a,b). A corresponding sequence of consecutive fibre shapes, found numerically with the model M 1 , is shown in figure 8. The significant bending from the middle to the last shape occurs on a very short time scale.
In figure 3(b), the bending time τ b was defined as the time when the maximal local fibre curvature reaches one half of its largest value, κ(τ b ) = κ b /2. We now add a physical interpretation: at a time close to τ b , the deflection z 1 (t) of the fibre tip approaches the first local maximum, visible in figure 7(a,b). The corresponding fibre shape at t = τ b is shown in figure 8 as the last one.
The linear regime of small deformations is limited to short times. For example, in figure 8, all shown fibres are already beyond this regime. However, to understand the dynamics in the early phase, it is worthwhile to begin the analysis from the linear regime where the universal scaling of shapes follows directly from the self-similar solution, specified by (4.6) and (4.9). One of characteristic features of the self-similar solution is that the same deflection z 1 (t) = u(0, t) of the fibre tip is reached at the same rescaled time tA −1/3 ∝ tη 1/3 (and with the same rescaled length sA −1/3 along the fibre). The numerical results shown in figure 7(b) illustrate that z 1 (t) is determined by tA −1/3 not only in the range of small deformations, but also beyond it. The tip deflection z 1 (t) remains a universal function of tA −1/3 until its argument reaches value slightly smaller than, but very close to t max A −1/3 ≈ 7.2, with t max defined as the time when z 1 has a local maximum u max . In the log-log scale, as shown in figure 7(b), the universal curve is close to a straight line for tA −1/3 4 (in the linear regime). For tA −1/3 4, a significant deviation from the straight line is observed caused by nonlinear effects.
The deviations from the linear regime in the numerical simulations can be interpreted by the elastica evolution in the range where the nonlinear terms in (2.20) become important, i.e. κ ss ∼ σ κ. (4.19) Next, we observe that the change in the elastica dynamics occurs away from the small deflection state so we can use the relation (4.2b), which implies the scaling, As η ∝ A −1 , these balances suggest that the time scale τ b , when the nonlinear terms become important, follows the power law τ b ∝ A 1/3 , which is in a good agreement with results presented in figure 5(b-e). For times close to τ b and t max , the tip deflection z 1 (t) leaves the universal curve, as shown in figure 7(a,b), and the maximum deflection u max depends on A. We observe that in the range of A, in which the evolution follows a power law, u max ∝ A 0.33 , as determined numerically from the bead model M 1 and shown in figure 9(a). This scaling might reflect a memory of the initial phase of the fibre movement as analysed using the elastica. A linearly deflected fibre changes its shape over a length scale s ∝ A 1/4 t 1/4 and time scale t ∝ A 1/3 (4.6) thus, we find that s ∝ A 1/3 at time t. However, as illustrated in figure 8, for early times t τ b , the typical length scale of bending remains almost constant in time, what allows us to expect that u max ∝ A 1/3 . The numerical results suggest that during the rapid bending the scaling in the x direction becomes comparable to the scaling in the z direction. To demonstrate this feature, we collect several fibre shapes for different values of n and A, at times t max of the maximum deflection (figure 9b), and we replot them in figure 9(c) by rescaling both axes by A 1/3 , which allows for the approximate overlapping of the shapes after translating them in the x direction.
The scalings with A 1/3 , typical for the early phase of the bending process, do not depend on n. For t τ b , bending of the fibre ends is a local process. However, the typical bending length scale increases with A, and therefore for a larger stiffness, the scalings are satisfied by a sufficiently long fibre only.

Bead model simulations
We now move on to discuss the dynamics for times τ b t τ b + τ c , when the fibre is significantly bent, with the maximum local bending curvature κ b2 /2 ≤ κ(t) ≤ κ b2 , where κ b2 is the largest maximum local curvature during this time period (see figure 3). In this range, the main feature of the dynamics is its maximum local curvature κ(t). Therefore, we first discuss if (and how) the characteristic features of the dynamics depend on a specific choice of the time instant when the curvature is determined.
In figure 3(b) we have illustrated that there exists a typical plateau of the curvature, κ b1 , and the largest value κ b2 . Comparison with figures 4(b) and 4(d) indicates that both values vary systematically with A. We analyse this dependence in figure 10(a). On a log-log plot, κ b1 is systematically below κ b2 , and the inset (for n = 140) illustrates that the ratio κ b1 /κ b2 slowly decays with increasing A, but this effect is not large. The numerical data show that, over a few decades of A, the ratio κ b1 /κ b2 changes only by 30 %, and it tends to κ b1 /κ b2 ∼ 0.7 for large A. This observation suggests that κ b2 depends on A in a similar way as κ b1 . The fibre is first in the state of a typical bend and tightens to a maximum curvature for a short time afterwards. However, the plateau in the κ(t) that allows determination of κ b1 , for a given n, occurs only for a finite range of A, while κ b2 is well defined for any value of A. For example, in case of n = 40, in figure 10(a) there are no data points above A = 100 indicating κ b1 while for n = 140, κ b1 can be observed up to A = 2000. Therefore, in the following we will focus on the analysis of κ b2 .
Depending on the values of A, three regimes of fibre bending can be identified, as shown in figure 10(b,c) for the bead models M 1 and M 2 , respectively. The schematics in figure 10(b) show three typical fibre shapes with n = 20 for each of the regimes. First, in the regime of a very flexible fibre (A 10), the maximum curvature is close to the excluded volume (EV) of the beads, with log 10 κ b2 log 10 ( √ 3) ≈ 0.24, which is independent of n. Second, there is a regime A 10, where κ b2 as a function of A continues with a power-law dependence until it deviates from the slope, which happens for different A depending on n. The larger n, the larger range of A that exhibit the power-law dependence. Inside this regime, for a given A, all fibres that are long enough have the same κ b2 , which is independent of n. We interpret this response as local bending. Third, there is a large A regime, which starts after κ b2 departs from the power law. This corresponds to κ −1 b2 comparable to or larger than the fibre length, which we interpret as global bending. This classification of κ b2 is valid even for very long fibres (having multiple loops), and also for very stiff ones (with no pronounced plateau of the fibre curvature κ(t)). For each n, the regimes of A where the power laws are observed for κ b2 agree with the corresponding regimes identified for τ b (compare the ranges of A in figure 10(b,c) with the ranges in figure 5(b,c), respectively).
Comparison of figures 10(b) and 10(c) indicates that, if the bending stiffness ratio A is not very small (A 10) and not too large (with the upper bound dependent on n), the power-law scalings of κ b2 predicted by the M 1 and M 2 models are in a reasonable agreement with each other, and the curves only weakly depend on n. However, for 10 A 100, the maximum curvature κ b2 in the M 2 model decays more rapidly with A than in M 1 , with approximately κ b2 ∝ A −1/3 rather than κ b2 ∝ A −1/4 , respectively (for the M 1 model, κ b2 ∝ A −0.253±0.003 as determined numerically). For A 10, the maximum curvature κ b2 determined from the M 2 model saturates at the EV value while in model M 1 this effect is seen for more flexible fibres with A 1. Although the treatment of hydrodynamic interactions is more precise within the bead model M 2 , it seems that the main reason for some differences between the maximum bending curvature κ b2 obtained by the M 1 and M 2 models is the use of different expressions for the bending potential energy, as discussed in detail in appendix B.

Comparing with the elastica model
We are going to show now that the scaling of the fibre maximum curvature κ b ∝ A −1/4 , independent of n and characteristic of the local bending, can be argued with the elastica model (2.20). We propose that in the local bending process there is only one length scale κ −1 b representative of the deformed fibre. This is consistent with the models of Harasim et al. (2013), LaGrone et al. (2019 and Liu et al. (2018) and with our findings that κ(t) is in the curling motion close to a typical constant value κ b1 .
Next, from the linearity of shear flow, we argue that the magnitude of the flow velocity incident on the fibre and the fibre velocity scale linearly with the length scale (ê Comparing the magnitudes of terms in (2.20b), we find that the dimensionless tension scales as This dependence, together with (2.20a), gives resulting in κ b ∝ η 1/4 ∝ A −1/4 . It is also true that σ ∝ A −1/2 and σ s ∝ A −3/4 . Note that these arguments apply to the maximum curvature in the whole range of the curling motion with a large shape deformation, in particular for κ b = κ b1 and κ b = κ b2 . The scalings obtained from the elastica model can be compared with the results from the bead-spring simulations with the model M 1 . The force F i acting on each bead i as the result of the elastic constitutive laws can be decomposed into the force components normal N i and tangential T i to the fibre centreline. In figure 11(a) we show shapes of locally bent fibres with n = 100 for three different values of A. The colour-coded representations of N i and T i are included in the following way. Each bead is depicted by a hemisphere, which has an orientation that indicates the direction of T i (inset), while the colour coding shows the ratio of the magnitudes of the forces normal and tangential to the fibre, |N i |/|T i |. In order to compare the simulation data quantitatively with the scalings deduced above from the elastica, it is sufficient to choose any time from the curling motion of the fibre. As we compare between different A, we introduce the transformation of the bead numbering i = (i − i 0 )A −1/4 (i is a discrete analogue of the arc length s), where a shift i 0 is chosen (for each fibre separately) to overlap the extrema. In the figure 11(b,c) we show the profile of local curvature κ i over half of the fibre (i = 1 . . . 50) for the shapes presented in figure 11(a). In figure 11(b) the raw data are plotted and in figure 11(c) κ i is multiplied by A 1/4 to show the scaling suggested by the elastica model. From the bead-spring simulations we have direct access to T i ·ê s acting at the centre of bead i, which is the analogue of the derivative of the tension σ s (s) for the elastica. The value of T i ·ê s is shown in figure 11(d),  figure 11(e) we demonstrate that the tangential forces scale as A −3/4 , which has been also suggested by the elastica model.

Curling velocity and curling time
In § 3 we introduced the curling motion and the associated curling time τ c (see figure 3). During the curling motion, the first bead travels from left to right approximately over the distance L = 2na with respect to the fibre's centre. Snapshots illustrating the shape evolution are shown using the schematics at the top of figure 12(a) for the fibre with n = 100 and A = 100. We define the curling velocity v x (t) as the x-component of the velocity of the first bead. At the bottom of figure 12(a), we plot v x versus time, for the fibres with A = 100 and different values of n. Initially, v x is close to zero. Then, it rises significantly and we observe the first peak, the plateau and the second peak. Numerical simulations show that for a given A, the profile of v x (t) in the initial phase of the motion is almost the same for different n. In figure 12(a), the plots of v x (t) for different n are almost superimposed for a long time. The changes between v x (t) with different n occur when the fibre stops to undergo curling motion due to its limited length. The peaks of v x are observed at the times of the steepest changes in κ(t), as illustrated in figure 12(b). The first peak takes place at the time close to the bending time τ b , and the second one approximately after the curling time τ c (compare with figure 3). Therefore, our definitions of τ b and τ c seem to well separate three different phases of the fibre dynamics.
We have found empirically that τ c as a function of A and n can be collapsed on a single universal line when plotted versus A/n 3.5 , as shown in figure 12(c). For small values of A/n 3.5 , the curling time τ c tends to a power law with the exponent −1/3, as determined empirically. Thus we observe that, in the limit of long fibres, we can approximate the curling time as τ c ∝ A −1/3 n 1.17 . That is, τ c is almost linear in n. The deviations might be related to a larger average resistance during the curling motion of longer fibres. Indeed, figure 12(a) illustrates that, for longer fibres, the contribution to the average curling velocity from the initial and final peaks is smaller, and therefore the average curling velocity is smaller, which leads to the curling time increasing with n a little faster than linearly.
The dynamics analogous to the curling motion was investigated in the literature experimentally, numerically and by the elastica model, with and without the Brownian motion (Forgacs & Mason 1959b;du Roure et al. 2019). We have shown here that on the onset of curling motion the characteristic length scale (∝A 1/3 ) is different from the one observed later in the highly bent state (∝A 1/4 ), which leads to the analogous scalings for v x at earlier and later times, respectively. Therefore, we emphasize the importance of the time evolution during the curling motion but at the same time we benefit from the previous studies of Harasim et al. (2013) and Liu et al. (2018) who reported the linear dependence between the local radius of curvature of the bent tip and its track velocity (analogous to our curling velocity), both approximately constant in time.
We have found that, in the regime of the local bending, during the curling motion, the curling velocity and the local curvature are mostly determined by the bending stiffness A and practically do not depend on the fibre aspect ratio n, This finding agrees well with the results of another numerical model of LaGrone et al. (2019) where only minute changes in the local radius of curvature of the fibre tip and its snaking (analogous to our curling) velocity have been observed in a wide range of relatively large fibre aspect ratios.
6. Universal scaling and phase diagram

Shapes of fibres with different n and A
The transition from the locally to globally bent fibres, observed for the increasing values of the bending stiffness A and illustrated in figure 10(b,c), motivated us to search for κ b2 n as a universal function of A/n γ , with a certain value of the exponent γ . We use here κ b2 n because in the global bending mode we expect bending along the whole fibre length. Indeed, in figure 13(a,b), plotted in log-log scale, we find the universal scaling of κ b2 n, based on the numerical simulations M 2 (in a) and M 1 (in b), respectively. We added to figure 13(a) also the results of the M 2 simulations reported by Słowicka et al. (2015), with the parameters n = 10, L 0 /(2a) = 1.01 and k s = 2000. From the numerical data for the model M 2 we obtain the exponent γ = 3.25 and we find the slopes −0.3 and −5 of the two straight lines for the local and global bending regimes, for log 10 (A/n 3.25 ) −2.9 and log 10 (A/n 3.25 ) −2.3, respectively. The fits agree very well with the results of the M 2 simulations, and reasonably well with the results of the M 1 simulations, as shown in figures 13(a) and 13(b), respectively. The deviations from the universal curve are observed only owing to the EV effects seen for very flexible fibres, with the EV value of the maximum local curvature κ b2 = log 10 √ 3 ≈ 0.24. The deviations correspond to the first (small A) regime of the fibre bending described in § 5.1 and shown in figure 10. For the local bending, the relation log 10 (κ b2 n) ∼ −0.3 log 10 (A/n 3.25 ) + 0.48, fitted to the M 2 numerical data, gives the approximate scaling of the maximum curvature κ b2 ∼ A −0.3 independent of n, in agreement with the previous discussion of the local character of the dynamics of very elastic fibres. The exponent −0.3 is close but not identical to the −1/4 fitted to the M 1 numerical data in figure 10(b). In the global bending regime, we find that κ b2 n ∼ (A/n 3.25 ) −5 .
The fitting of the exponent γ in the relation A/n γ is based on the choice of 2an to represent the fibre length L, and it is sensitive to a choice of L. For example, γ ≈ 3 if the fibre length L = (n − 1)L 0 is chosen. Such a shorter fibre length was proposed by Farutin et al. (2016) as the result of comparing shapes of flexible fibres and deformable vesicles in Poiseuille flow and partially accounts for the rigidity of the beads at the fibre ends. On the other hand, in shear flow, a matching of the tumbling period with the half-period of the Jeffery orbit could be used to determine L. In the bead model M 2 for stiffer fibres, the effective aspect ratio L defined in this way is greater than 2an, and it could lead to a γ closer to 4, which means also closer to the scaling κ b2 ∼ A −1/4 of the local bending proposed in figure 10(b).
We expect that also the shape of the whole fibre is a universal function of A/n 3.25 . Indeed, as illustrated in figure 13, for the models M 2 (left) and M 1 (right), the fibre shapes depend on n and A approximately through the ratio A/n 3.25 . We show it separately for the global bending in figure 13(c,d) and for the local bending in figure 13(e, f ). The corresponding values of A/n 3.25 are explicitly indicated below each fibre shape, with approximately the same values for all the similar shapes. . The dynamical modes of the fibres initially aligned with the flow are the following: the fibres that are coiled and do not straighten out (mode 1, triangles); the fibres that straighten out along the flow while tumbling periodically and bend locally (mode 2, rhombus) or globally (mode 3, circles). A sharp transition between the fibres that straighten out while tumbling and the fibres that stay coiled is marked by a dashed line and the stars, taken from Słowicka et al. (2015), for the M 2 model and by a solid line for the M 1 model. In contrast, the transition between fibres bent locally and globally is gradual (grey area). The sizes of symbols for the M 2 model discriminate between data from this work with l 0 = 1.02 and k s = 1000 (large open symbols) and the data of Słowicka et al. (2015) with l 0 = 1.01 and k s = 2000 (small open symbols). The symbols + indicate such points that τ b2 = τ f .
Comparison of the power-law scaling of the fibre shapes with an attempt to find another similarity solution, based on a logarithmic dependence on n and a generalized elasto-viscous number, standard in the SBT and elastica approaches, is presented in appendix C. We show there that, although such a possibility cannot be excluded, it seems to be quite complicated to construct. Using simple arguments, we are able to find such a scaling function only for the local bending mode.

Phase diagram of the dynamical modes
The analysis of the fibre dynamics can be summarized on a phase diagram in the space of the fibre aspect ratio n and the bending stiffness A. In figure 14 we show the numerical results, with essentially the same features for the bead models M 1 and M 2 . The elastic fibre initially aligned with the shear flow has three characteristic modes of motion, depending on values of n and A: (Mode 1) The fibre does not straighten out again. The curvature κ does not return to zero after the first bending event. (Mode 2) The fibre bends locally, curls and then stretches; correspondingly curvature grows, reaches a plateau and then returns to zero in a periodic way. (Mode 3) The fibre periodically bends globally along the whole length. Curvature maxima are observed but the plateau vanishes.
At early times, the fibre is bent only at the ends. During the curling motion, as shown in figures 3 and 12 and earlier by Harasim et al. (2013), Liu et al. (2018) and LaGrone et al. (2019), the range of the most curved segments shifts towards the central part of the fibre. The fibre ends become straight and almost aligned with the flow, and the length of straight ends increases with time. Therefore, in general, we might expect that the end of such a long fibre will behave in a similar way as a fibre of a comparable length aligned with the flow. Therefore, if the curling continues long enough, with τ c τ b , the fibre may bend its end again, even several times, and it will not straighten out. Indeed, such a scenario sometimes happens for very long or very flexible fibres, as shown in figure 4, and earlier by Nguyen & Fauci (2014) and LaGrone et al. (2019). Using our scalings, τ b ∝ A 1/3 and τ c ∝ A −1/3 n 1.17 , a dynamical transition could be expected around A ∝ n 1.75 . However, the physical origin of the transition between the coiled and straightening out modes is more complicated. Shorter fibres cannot bend several times, but still they do not straighten out along the flow when their bending stiffness is small enough.
The transition between the coiled and locally bent modes for shorter fibres with n ≤ 40 and a wide range of values of the bending stiffness A has been analysed by Słowicka et al. (2015). The dynamics of flexible fibres was evaluated over a long time, starting from the initial configuration aligned with the flow. A characteristic value A CS (n) was found for the transition between the fibres that remain coiled and the fibres that straighten out along the flow while tumbling, with A CS (n) ∝ n 3/2 . Moreover, the dynamics was shown to be very sensitive to a small change of A close to A CS . For A slightly below the critical value, fibres often straightened out a smaller or larger number of times before changing to the coiled mode. Słowicka et al. (2015) sorted the data for the modes 1 and 2 based on the long-time behaviour. We present in figure 14 (small open symbols) some of the results obtained by Słowicka et al. (2015). The inset illustrates high precision of the critical values A CS determined there and marked by stars in figure 14. The results of the model M 2 applied in this work (large open symbols) also support the A CS (n) ∝ n 3/2 scaling of the transition between the coiled and straightening out modes. The numerical simulations in the M 1 model also agree well with the above scaling, with a different factor which could be interpreted as the result of different bending potentials in both models.
In contrast to the transition between the modes 1 and 2, the transition between the modes 2 and 3 is not sharp. It takes place in a range of the phase space (n, A) marked in grey in figure 14. This stripe corresponds to −2.3 log 10 (A/n 3.25 ) −2.9, i.e. the range between the local and global bending found numerically with the bead model M 2 and shown in figure 13(a), in agreement with the findings of the model M 1 presented in figure 13(b). The different symbols are the locally and globally bent fibres that indicate just an approximation, based on a comparison of the time instant of the maximum κ b to the flipping time τ f , as described in appendix D. In the regime of the local bending, for a smaller A or larger n, the bending time scales as τ b ∝ A 1/3 and the maximum local curvature scales as κ b ∝ A −1/4 , independently of n. In the regime of the global bending, τ b ∝ n independently of A, and κ b n ∝ (A/n 3.25 ) −5 .
The transition between the local and global bending could be interpreted as a competition between bending and rotation. If the fibre bends before it manages to rotate in shear flow, it belongs to the local bending mode while if it rotates before it bends, it belongs to the global bending mode. Approximating the rotation time as T J /4 ∝ n, and equating τ b ≈ T J /4, we obtain A ∝ n 3 , which is an approximation of the transition between the local and global bending shown in figure 14. Another way of looking at the transition between the local and global bending is to compare the typical length scales. The length of the bent fibre end at τ b scales as A 1/3 . In the local bending mode, it needs to be smaller than half of the fibre length, proportional to n, which again estimates the transition roughly as A ∝ n 3 .

Conclusions
In this paper, we analyse the evolution of an elastic thin fibre that is initially straight and aligned with an ambient shear flow. We consider a wide range of the fibre aspect ratios n and many different values of the bending stiffness ratio A (i.e. the ratio of the bending forces to the hydrodynamic forces caused by the flow rateγ ). We use two theoretical descriptions of the fibre: the bead-spring model with elastic potential energy and hydrodynamic interactions, and also a generalized elastica model. These two approaches complement each other and allow to rationalize analytically many of the observed numerical results for the bead-spring model.
To quantify evolution of the fibre shapes, we introduce and evaluate numerically three main characteristic time-dependent quantities: the deflection of the fibre tip u(0, t) in the direction perpendicular to the flow, with the first maximum at u max , the maximum local curvature κ(t), with the largest value κ 2b = max t κ(t), and the curling velocity v x (t), with the maximum value v xm . Their behaviour allows us to identify three characteristic time scales of the dynamics: the bending time τ b , the curling time τ c and the tumbling time τ equal to the half-period T J /2 of the effective Jeffery rotation.
Accordingly to the time scales, we identify three characteristic stages of the time evolution of flexible fibres initially aligned with the flow: bending of the fibre tips for 0 ≤ t ≤ τ b , curling of the deformation towards the centre of the fibre for τ b ≤ t ≤ τ b + τ c and stretching of the fibre for τ b + τ c ≤ t ≤ T J /2 with an effective Jeffery period T J . In the bending stage, we find the scaling u(0, t) ∝ (t 3 /A) 1/4 , with the maximum u max ∝ A 1/3 at τ b ∝ A 1/3 , all independent of n, in agreement with the local character of the early stage of the fibre dynamics for all the modes. In the curling stage, the maximum curvature κ(t) and the curling velocity v x (t) are approximately independent of n (except for short final time intervals), and for a sufficiently large n change in time only a little (except for short initial and final time intervals), as argued by Harasim et al. (2013) and Liu et al. (2018).
We demonstrate that τ b /n, κ b2 n and τ c depend on n and A approximately through certain universal functions A/n α . Based on the numerical simulations, we determine the exponents α which are equal to 3, 3.25 and 3.5, respectively (close to but different than 4 as in case of the elasto-viscous number). In particular, the shapes of fibres (and the maximum 'global' curvature κ b2 n) are shown to depend on n and A approximately through A/n 3.25 . A referee suggested trying another similarity function, dependent on log n, for the same reason that SBT depends on the logarithm of the aspect ratio. An (unsuccessful) attempt to replace a power law with the exponent 3.25 by a logarithmic dependence is described in appendix C. In figure 16 we present an analogue of figure 13, but with an elasto-viscous number log 10 [A(ln n + ln 2 + 1/2)/n 4 ] on the horizontal axis. The constants in the numerator follow from (8.8) for the SBT transverse motion, derived by Batchelor (1970a). Different constants in the logarithmic expressions are also used e.g. by Becker & Shelley (2001) and Young & Shelley (2007). We find it interesting that the plots of κ b2 n versus the elasto-viscous number in figure 16 seem to indicate that the fibre shape (at the time of its maximum curvature) might be a universal function of the elasto-viscous number in the local bending mode (left part of the plot) but not in the global bending mode (right part of the plot). The difficulty of matching a logarithmic expression might be related to relatively small values of n in our simulations. A scaling which involves ln n might require very large aspect ratios n. It seems logical that comparison of κ b2 n for fibres with different thicknesses and the same length may depend not only on the elasto-viscous number, the parameter adequate for asymptotically large values of n. Moreover, it is known from SBT that the constants added to ln n are sensitive to fibre shape. It is also worth remembering that for moderate values of n a constant added to the logarithm has a significant influence. This constant for a flexible, deformed fibre depends on both n and A, and it is difficult to evaluate it theoretically. However, it is clear that its value is different from Batchelor's result for a straight rigid rod. Therefore, it is clear that an additive constant for a flexible fibre should depend on shape, and therefore on both n and A.
Based on the numerical simulations, we classify the dynamics of flexible fibres in the phase space of n and A, according to the essential features of the motion and shape deformation. We find three different modes of the fibre motion: coiled, locally bent and globally bent, and we identify the characteristic ranges of n and A for each of them. The classification refers to the fibres initially aligned with the flow. In the coiled mode, found for larger n or smaller A, the fibres later do not straighten out along the flow, in contrast to the other two modes. Global bending of fibres takes place at smaller n or larger A and it corresponds to coherent deformation along the whole fibre length. Local bending means that only a part of the fibre is curved, and it is typical for intermediate values of n and A. Essentially, basic features of these three scenarios were identified already in experiments performed by Forgacs & Mason (1959b) who called them a coiled orbit, springy rotation and snake turn, and then analysed e.g. by Lindström & Uesaka (2007), Harasim et al. (2013), Nguyen & Fauci (2014), Liu et al. (2018) and LaGrone et al. (2019), with differences between shapes observed under different physical conditions (e.g. with or without Brownian motion). Here, for the first time, a systematic analysis of these three modes is performed.
In particular, all three stages of the evolution are observed for the local bending dynamical mode. In the global bending mode, the curling stage is absent, and for the coiled mode there is no stretching stage and the curling motion is much more complicated than in the case of the local bending mode. For the local bending mode, we find the approximate scaling κ b2 ∝ A −0.3 independent of n, with the exponent close to −1/4 found by Harasim et al. (2013) and Liu et al. (2018) (in the M 1 model we can also deduce from the numerical results that κ b2 ∝ A −1/4 , see figure 9(b). In this case, our data seem to agree with both scalings). The dependence of the global bending on A has not been analysed. We find that the maximum 'global' curvature κ b2 n ∝ (A/n 3.25 ) −5 decays rapidly with A, which is much faster than in the local bending mode.
Our analysis of the dynamics for different n and A indicates that the transition between the local and global bending modes takes place for −2.3 log 10 (A/n 3.25 ) −2.9. Therefore, it is close but not exactly equal to a certain universal value of the elasto-viscous numberη = 8πμ 0γ (2a) 4 n 4 /EI ln( −1 ) (or effective viscosity (effective flow forcing) μ = 8πμ 0γ (2a) 4 n 4 /EI), which scales as n 4 /A (Becker & Shelley 2001;Tornberg & Shelley 2004;Harasim et al. 2013;Nguyen & Fauci 2014;Liu et al. 2018;LaGrone et al. 2019). Moreover, we have found that a second transition, between the coiled fibres and the fibres that straighten out, is given as log 10 (A/n 3/2 ) = C and it takes place at different values of the elasto-viscous number when n or A are changed. Therefore, we find it beneficial to extend the concept of the elasto-viscous number and analyse the dynamics in the phase space of n and A. Certain features of the dynamics depend on n and A in a more complex way than the elasto-viscous number predicts.
We also analyse the elastica model and rationalize some of the scalings described above. We provide a self-similar exact solution of the linear elastica equations when the fibre is almost aligned with the flow. The main new idea is to assume as the boundary condition the existence of a constant hydrodynamic force exerted on the fibre tip by the rate of strain of the ambient flow. This allows tracing of the early stage of the fibre bending from the initial position aligned with the flow which, is not possible within the standard elastica approach. Moreover, we derive such a hydrodynamic force from the theory of hydrodynamic interactions and evaluate it numerically. These findings indicate that the standard elastica model in some cases may be too simple to predict the dynamics, and cannot always serve as a source of a theoretical explanation.

A.2. HYDROMULTIPOLE -model 2
Consider now a general system of n spherical particles immersed in an incompressible fluid flow with velocity V (R) and pressure p(R) that satisfies the quasi-steady Stokes equations with the boundary condition at infinity, where V ∞ (R) is an arbitrary external fluid flow. Assume the no-slip boundary conditions at the bead surfaces, S i , The integral representation (Pozrikidis 1992) and the method of induced forces (Cox & Brenner 1967;Mazur & Bedeaux 1974;Felderhof 1976) can be used to express the fluid velocity in terms of the Oseen tensor T 0 (R −R), given e.g. by Kim & Karrila (1991), applied to the density f j (R) of the forces exerted by the surface of the particle i on the fluid. Application of the boundary conditions (A2) results in the boundary integral equation for the force density f j (R), which is then projected onto a complete set of elementary (spherical multipole) solutions of the Stokes equations (Cichocki, Felderhof & Schmitz 1988;Felderhof 1988). As a result, an infinite set of algebraic equations is obtained. This set is truncated at a certain multipole order L and solved for the vector of the force multipoles. Converting from the spherical to the Cartesian representation, we obtain a linear relation between (i) the forces F i , torques T i , stresslets S i and higher-order force multipoles exerted on the fluid by the particles i = 1, . . . , N, and (ii) the translational and rotational velocities, U j and Ω j of particle j = 1, . . . , N, and the multipoles of the external velocity field V ∞ (R). This relation is written using the grand friction matrix ζ , ⎛ ⎜ ⎜ ⎝F T S · · · ⎞ ⎟ ⎟ ⎠ = − ⎛ ⎜ ⎜ ⎝ ζ tt ζ tr ζ td · · · ζ rt ζ rr ζ rd · · · ζ dt ζ dr ζ dd · · · · · · · · · · · · · · · ⎞ ⎟ ⎟ ⎠ · ⎛ ⎜ ⎜ ⎝Ṽ In the above the 3N dimensional vectors areF = (F 1 , F 2 , . . . , F N ),T = (T 1 , T 2 , . . . , T N ), U = (U 1 , U 2 , . . . , U N ),Ω = (Ω 1 , Ω 2 , . . . , Ω N ). The velocity multipoles are evaluated at the centres R i of the particle i = 1, . . . , N from the external flow velocity and its derivatives. In particular,Ṽ ∞ = (V ∞ (R 1 ), . . . , V ∞ (R N )). Similarlyω ∞ is the vector of vorticities, with ω ∞ (R i ) = 1 2 (∇ × V ∞ )| R i . Next, we introduce the tensor of strain rates E ∞ = (E ∞ (R 1 ), . . . , E ∞ (R N )) with E ∞ (R i ) = 1 2 (∇V ∞ + (∇V ∞ ) T )| R i . The second rank strain tensors E ∞ (R i ) are symmetric and traceless and thereforeẼ ∞ has 5N independent components. Finally, the symmetric tensorS = (S 1 , . . . , S N ) represents the particle stresslets. To speed up the convergence of the multipole expansion, the lubrication correction is applied to friction matrices, as described by Durlofsky et al. (1987), Sangani & Mo (1994) and Cichocki et al. (1999).  Figure 15. Comparison of three bead models of a flexible fibre. The maximum curvature κ b2 vs. bending stiffness A, evaluated with the use of the models M 1 (filled symbols), M 2 (empty symbols) and M 3 (stars). The EV limit is marked as the horizontal dashed line.
in figure 15. For A 20-30 there is a difference between M 1 and M 2 , and also between M 1 and M 3 . However, the results from M 2 , M 3 are close to each other. Therefore the form of the bending potential at large angles seems to be essential for the dynamics of more flexible fibres, in agreement with the same conclusion for sedimenting flexible fibres given in Bukowicki & Ekiel-Jeżewska (2018). Indeed the cosine bending potential from the models M 2 and M 3 is more flexible than harmonic bending potential from M 1 and for large bends it leads to higher curvatures. For such large values of the bending stiffness that the radius of curvature is three or more times larger than the bead radius, the maximum curvatures obtained with the models M 1 , M 2 , M 3 are the same because the bending potential energies behave alike. This is expected since both are intended to approximate the elastic bending potential energy in the limit of large A (it was estimated by Bukowicki & Ekiel-Jeżewska (2018) that a difference smaller than 5 % is expected for the maximum bending angles π − θ i = κ 0.7).
The agreement between the M 1 , M 2 and M 3 models of more stiff fibres is illustrated in figure 15 where we also present the maximum curvature κ b2 for much longer fibres with n = 100, evaluated with M 1 and M 2 models. For A 100 a good agreement is obtained between all the computations performed with the use of M 1 and M 2 models, regardless of the fibre length. Actually, all the properties of flexible fibres that were discussed in the main text for the model M 1 are analogous for the model M 2 in the range of intermediate and large A, where the constitutive laws are manifesting similar behaviour.

Appendix C. Discussion of the universal scaling
The idea of figures 13(c-f ) and 16(b,c) is to compare properties of fibres of the same length, but different thickness and different bending stiffness. To this goal, on vertical axis we plotted the maximum curvature normalized by the inverse fibre length (rather than by its inverse width), i.e. κ b2 n. The universal scaling of the maximum curvature κ b2 n, provided in § 6.1, is based on the similarity solution as a function of A/n 3.25 . Here we check if a universal scaling can be based on a certain modification of the standard elasto-viscous number, including a logarithmic dependence on n. Therefore in 16(a) we plot in log-log scale κ b2 n versus the elasto-viscous number A(ln n + ln 2 + 1/2)/n 4 , with the constant modified following SBT of Batchelor (1970a). The scaling works reasonably well for the  Figure 16. The elasto-viscous number A(ln n + ln 2 + 1/2)/n 4 , based on the SBT by Batchelor (1970a), is not a universal scaling function of the numerical results obtained with the model M 2 . (a) The maximum curvature κ b2 n, scaled by the inverse of the fibre length, is plotted in log-log scale versus the elasto-viscous number A(ln n + ln 2 + 1/2)/n 4 . The similarity scaling of shapes is observed in (b) for the local bending mode, but it does not work in (c) for the global bending, with values of the elasto-viscous number as indicated. In

Appendix D. Time scales close to the transition between local and global bending
In § 6, the distinction between the locally and globally bent fibres was approximately estimated by comparing the time instant τ b2 of the maximum curvature κ b2 with the flipping time τ f Farutin et al. 2016;Słowicka et al. 2020), i.e. the time when two end beads have the same x coordinate. This procedure is illustrated in figure 17 using the numerical data from the model M 2 . If flipping occurs for τ f < τ b2 , the corresponding point (n, A) is marked by a square in the phase diagram ( figure 14), while if flipping happens for τ f = τ b2 or τ f > τ b2 , by a plus or a circle, respectively. However, it should be kept in mind that the change between the locally and globally bent fibres takes place in a certain range of A. The values marked by the symbols + in figure 14 correspond already to the global mode scaling κ b2 n ∝ (A/n 3.25 ) −5 but they are too large to satisfy κ b2 ∝ A −1/4 typical for the local bending. The shapes at the maximum curvature shown in figure 17(b,c) are bent globally, but the shape presented in figure 17(a) is not locally bent. In the local bending mode, at the moment of the maximum local curvature, the fibre ends are almost parallel to the flow.