Flexible fibers in shear flow approach attracting periodic solutions

The three-dimensional dynamics of a single non-Brownian flexible fiber in shear flow is evaluated numerically, in the absence of inertia. A wide range of ratios A of bending to hydrodynamic forces and hundreds of initial configurations are considered. We demonstrate that flexible fibers in shear flow exhibit much more complicated evolution patterns than in the case of extensional flow, where transitions to higher-order modes of characteristic shapes are observed when A exceeds consecutive threshold values. In shear flow, we identify the existence of an attracting steady configuration and different attracting periodic motions that are approached by long-lasting rolling, tumbling, and meandering dynamical modes, respectively. We demonstrate that the final stages of the rolling and tumbling modes are effective Jeffery orbits, with the constant parameter C replaced by an exponential function that either decays or increases in time, respectively, corresponding to a systematic drift of the trajectories. In the limit of C → 0, the fiber aligns with the vorticity direction and in the limit of C → ∞, the fiber periodically tumbles within the shear plane. For moderate values of A, a three-dimensional meandering periodic motion exists, which corresponds to intermediate values of C. Transient, close to periodic oscillations are also detected in the early stages of the modes.

A basic question in all of these configurations is how the dynamics and shapes of deformable elongated objects in flow depend on their flexibility. This problem has been investigated extensively at macro-and nanoscales for fibers in extensional, cellular, and corner flows [2,8,14,47]. For example, it has been shown that the typical pattern of a fiber's evolution in extensional flow is related to consecutive threshold values of the characteristic ratio of bending to hydrodynamic forces exerted by the fluid flow. When these values are exceeded, higherorder modes of the fiber shape are activated, with shorter characteristic length scales of elastic deformation [2][3][4][5]13,14]. In extensional flow, there exists a family of characteristic, well-defined shapes.
However, in general, the fiber deformation that occurs may depend on the type of flow [5,12,48]. Shear flows typically occur owing to the presence of container walls, and so are of practical and fundamental interest. Therefore, as illustrated in Fig. 1, in this paper we study three-dimensional dynamics of a single flexible fiber immersed in steady shear flow with velocity where e x is the unit vector along x andγ is the shear rate. The Reynolds number of the system is assumed to be much smaller than unity and the fluid flow satisfies the quasisteady Stokes equations. In this article we provide a new perspective on the threedimensional evolution of flexible fiber shapes in shear flow. We show that after a relaxation phase, a flexible fiber is attracted to one of several stationary, periodic, or close to periodic solutions, with different typical sequences of shapes and orientations. Features of these characteristic solutions, their presence or absence, stability or instability, and basins of attraction depend on the ratio A of local bending E π d 2 /64 to hydrodynamic πηγ d 2 forces, where E is the Young's modulus, d is the fiber diameter, and η is the fluid's dynamic viscosity. Therefore, by changing which depends on the fiber's aspect ratio L/d [4]. A rigid straight fiber in a low-Reynolds-number shear flow follows a periodic Jeffery orbit [49]. In this article we demonstrate that for a wide range of parameters and times, the motion and deformation of a flexible fiber can be interpreted as an effective Jeffery orbit that systematically drifts in time, owing to the exponential growth or decrease of the timedependent parameter C. For moderate values of A ≈ 9-12, there exists a range of values of C that correspond to periodic or close to periodic meandering motions.
A significant finding described below is that for a certain range of small values of A ≈ 4-5, fibers typically tend to align with the vorticity of the undisturbed flow. This result indicates a possibility to observe experimentally an ordered, dilute suspension of flexible fibers where all fibers are straight and aligned with the vorticity direction. We suggest that this ordered phase could be produced by adjusting the shear rate to reach the appropriate range of A. We are not aware of any such experimental observations.
As we document in this article, in shear flow (in contrast to extensional flow), the value of the bending-to-hydrodynamic ratio A does not uniquely determine the type of the fiber shape. This conclusion is based on two main features of the dynamics. First, for the same value of A, depending on the initial configuration or orientation, an elastic fiber evolves to a different characteristic sequence of shapes. Moreover, longlasting, chaotic transients are typical: we document close to periodic motions that later spontaneously change into periodic or effective Jeffery motions with different shape sequences.

A. Elastic fiber
The fiber is modeled [50] as a chain of N = 40 spherical beads of diameter d. The time-dependent position of the center of a bead i is denoted as r i . The centers of the consecutive beads are connected by springs of the equilibrium length 0 = 1.02d and the potential energy Herek is the spring constant and i = |r i − r i−1 | is the distance between centers of beads i and i−1. In this paper we assume that the dimensionless elastic resistance ratio k =k/(πηdγ ) is large, k = 1000, which leads to an almost constant fiber length. At equilibrium, the fiber is straight; its deformation costs energy dependent on the bending resistance, 1 Heret i is the unit vector parallel to the relative positions r i − r i−1 of the centers of beads i and i − 1. The total external force F i acting on bead i is elastic, We assume that the dimensionless bending-tohydrodynamic force ratio (relative bending stiffness) A =Â/(πηd 4γ ), given by Eq. (2), varies in the range of moderate values 4 A 40, where stiffer fibers subsequently deform and straighten while tumbling and very flexible ones remain coiled and do not straighten [51]. The length and time units in the simulations are, respectively, d and 1/γ .

B. Initial fiber configurations
To study the fiber evolution, we analyze the time dependence of the end-to-end vector n(t ) = (δx, δy, δz), shown in Fig. 1(b). We parameterize it by the standard spherical coordinates: the vector length L(t ) = |n(t )|, the angle (t ) between n(t ) and the vorticity direction y, and the angle (t ) between the projection of n(t ) on the xz plane and the x axis. Initially, n(0) = n 0 , (0) = 0 , and (0) = 0 .
To investigate the characteristic features of the flexible fiber dynamics in shear flow, we consider the following family of initial conditions. We assume that at t = 0 the fiber is straight and that all of the springs between the beads are at their equilibrium lengths [48]. The initial orientation of the fiber is given by the orientation of the end-to-end vector n 0 that links the centers of the first and the last beads. The length of this vector is |n 0 | = L 0 ≡ (N − 1) 0 . The direction of n 0 is parameterized by the spherical angles 0 and 0 , as indicated in Fig. 1(a). We consider the whole range of the initial orientations.
In the following, we will systematically investigate how the fiber's dynamics and shape evolution depend on the relative bending stiffness A and the initial orientation ( 0 , 0 ). In this way we will study the three-dimensional dynamics of a 1 For flexible Brownian fibers, the bending resistanceÃ, defined in Eq. (6), is often related to the fiber's persistence length l p =Ã/(k b T ) [2,13,14], where T is the temperature and k B is the Boltzmann constant. For flexible Brownian fibers, the Weissenberg number is defined as Wi =γ τ r [3], where τ r is the Brownian rotational time. Our fibers are non-Brownian, which means that we work in the limit Wi → ∞.

C. Fiber dynamics in flow
The dynamics of an elastic fiber is determined by the external (i.e., elastic) forces exerted on each bead, and the hydrodynamic interactions between them caused by the presence of the shear flow. We assume that the dynamic viscosity η of the fluid that surrounds the fiber is large enough for the Reynolds number to be much less than unity. In this limit, the incompressible fluid velocity v and pressure p satisfy the Stokes equations [52,53], The no-slip boundary condition at the beads' surfaces is assumed. The fluid is unbounded, with the ambient flow velocity v 0 given by Eq. (1); in the presence of the fiber, the fluid velocity v tends to v 0 when the distance from the fiber goes to infinity.
To solve the Stokes equations in the presence of the fiber, we use the advanced theoretical algorithm [54,55], based on the multipole expansion corrected for lubrication [56], implemented in the numerical code HYDROMULTIPOLE ( [56]). The method, similar to Refs. [57,58], is based on the boundary integral representation of the fluid velocity and the boundary integral equations for the surface density of the forces induced at the particle surfaces. These equations are projected on a complete set of elementary solutions of the Stokes equations (spherical multipole functions). The resulting set of linear algebraic equations is truncated at a controlled multipole order. The convergence of the multipole expansion is sped up by applying the lubrication correction [56,59,60].
In this paper, Eq. (9) simplifies, because there are no external torques, T = 0, and we are not interested in the bead rotations . The HYDROMULTIPOLE codes [56] are used to evaluate with high precision μ tt i j , μ tr i j , μ rr i j , F 0 , T 0 for given positions r i of the bead centers, and Eqs. (9) are solved numerically with the adaptive fourth-order Runge-Kutta method. More information about the numerical method and its accuracy is given in Appendix C.

III. ATTRACTING MODES OF THE DYNAMICS
One of our most significant findings is that, depending on the initial orientation and bending stiffness, the fiber is attracted to one of three distinct stationary or periodic solutions. The characteristic properties of these evolution patterns are illustrated in Fig. 2 and Supplemental Material Movie 1 [61], for A = 10 and three different initial orientations of the end-toend vector n 0 , ( 0 , 0 ) = (5 • , 10 • ), (45 • , 30 • ), (10 • , 10 • ). The colors blue, red, and green indicate time-dependent modes attracted to different stationary or periodic solutions. In the following these attractors will be called "rolling," "tumbling," and "periodic meandering," respectively.
The rolling mode (also called spin-rotation [28], or log rolling) is shown at early times in Fig. 2(a) and Supplemental Material Movie 1 [61]. In this mode, the fiber end-to-end vector tends at long times to the vorticity (y) direction, i.e., → 0. Moreover, all of the fiber beads rotate around y, and tend to align along y. The tumbling mode is shown at early times in Fig. 2  In this mode, the fiber end-to-end vector tends at long times to the xz plane ( → 90 • ). The convergence is confirmed by much longer simulations with 0 t 15 000.
The attracting character of the rolling and tumbling solutions have been identified previously [26][27][28][29]48,62] as the generic feature of the flexible fiber dynamics, and the corresponding dynamical modes have been called by the same names as their attractors [23,28,29].
To our surprise, we discovered that there exists a third dynamical evolution, illustrated in Fig. 2(c) and Supplemental Material Movie 1 [61]-an approach to an attracting three-dimensional "periodic meandering" solution, with oscillating periodically for 1800 t 10 5 . To the best of our knowledge, this behavior has not been observed before.
The results in Fig. 2 are focused on the relaxation phase after which the attractors are reached. The relaxation time can be very long, of the order of thousands of dimensionless units, as can be seen in Fig. 2. Moreover, long-lasting, close to periodic, transient motions (different from tumbling and meandering), which we will call "squirming," can appear in the relaxation phase, as shown in Fig. 2(c) and Supplemental Material Movie 1 [61]. Squirming motions are typically present for a wide range of initial orientations and values of A, and in all of the modes-those that tend to meandering, rolling, or tumbling solutions.
The approach to three different attracting solutions is also well illustrated by plotting the trajectories drawn by the tip of the end-to-end vector n(t ) = (δx, δy, δz). This concept is explained in Supplemental Material Movie 3 [61]. The results obtained for A = 10, and the same initial orientations as in Fig. 2, are shown in Fig. 3. We use different scales for different modes to illustrate the detailed characteristic features of each evolving mode, which is attracted by a different solution: stationary rolling, periodic tumbling and periodic meandering, respectively. Approaching these solutions takes a long time. We highlight in Fig. 3(a)(ii) the regular, anisotropic spiraling towards the steady state. We demonstrate in Fig. 3(b) a fast convergence of n to the shear plane. We illustrate in Fig. 3(c) the existence of a transient, almost periodic squirming trajectory in the relaxation phase of the meandering mode. The transition between the different modes is illustrated in Fig. 4. The columns (a)-(c) show trajectories of the fiber end-to-end vector for A = 10 and three different initial orientations, which are very close to each other, but represent different modes; we keep the same scale in each row of the figure. It is clear that the dynamics can be sensitive to a small change of the initial conditions. We next analyze the trajectories in Figs. 4(b) and 4(c) during the relaxation phase. By comparing them, we observe that at short and moderate times there appear a range of times with the characteristic squirming trajectory, which repeats almost periodically. This squirming end-to-end trajectory is shown separately in Fig. 5.
Summarizing, the transient periodic squirming motion is visible not only in the meandering mode (which ends up at the periodic meandering solution), but also in the tumbling mode (which ends up at the tumbling solution). For smaller values of the relative bending stiffness A, squirming motion is also often visible before the rolling solution is reached. The attractors-the rolling, tumbling and meandering solutionswill be shown in the next section.

IV. THE ATTRACTORS: STATIONARY AND PERIODIC SOLUTIONS
The three modes presented in Figs. 2-4 tend to three different attracting solutions that are compared with each other in Figs. 6 and 7. As illustrated in Fig. 6, the ordering of the end-to-end vector is different for each of these attractors, with the corresponding mean values δy/L 0 well-separated from each other.
In the rolling solution, the fiber is straight and it spins along the vorticity direction y. In the tumbling periodic solution (with the period T t = 181 for A = 10), the fiber end-to-end vector is perpendicular to y. It tumbles in the shear plane around the vorticity direction while the fiber straightens along the flow and then becomes coiled, with large-amplitude oscillations of length of its end-to-end vector.
The periodic meandering orbit (with the period T m = 735 for A = 10) has been observed even in very long simulations with more than 130 periods. Details about the periodic meandering motion are given in Appendix A. The evolution of shapes in the meandering and tumbling periodic solutions  [61]. The time range is 600 t 1224. Trajectories of the fiber end-to-end vector are close to periodic; shapes are a little different, but almost repeating after t ≈ 305.
FIG. 6. Three attractors of the fiber dynamics for A = 10: rolling stationary configuration (blue diamonds), together with tumbling (red), and meandering (green) periodic orbits. The xy, xz, and yz projections of the end-to-end periodic trajectories are plotted for t end − 800 < t < t end , with t end = 10 5 for meandering, and t end = 15 000 for tumbling trajectories. are compared with each other in Figs. 7(a) and 7(b) and Supplemental Material Movie 2 [61]. In both periodic orbits, the xz projections of the fibers have small amplitudes δz along the flow gradient direction, but the corresponding shapes are essentially different from each other. Moreover, in the tumbling motion, all of the fiber beads stay in the shear plane (i.e., δy = 0), while the meandering motion is three-dimensional, with large values of the end-to-end projection δy along the vorticity direction. Therefore, the meandering and tumbling periodic motions differ significantly from each other.

V. DEPENDENCE OF THE DYNAMICS ON THE BENDING STIFFNESS RATIO A
In the proceeding sections, we have systematically investigated how the evolution pattern depends on the fiber's initial orientation ( 0 , 0 ) for different values of the bending stiffness ratio A. A summary of the most important findings is illustrated in Figs. 8(a)-8(c) for 4 A 40. Flexible fibers with A = 4, shown in Fig. 8(a), belong to the rolling mode (blue dots) for all the initial orientations other than in the shear plane xz ( 0 = 90 • ). For rather stiff fibers with A = 40, shown in Fig. 8(c), the tumbling mode (red dots) dominates for most of the initial orientations, except some of those that are very close to the vorticity direction and lead to the rolling mode (blue dots).
For fibers of a moderate stiffness A = 10, shown in Fig. 8(b), there exists a range of the initial orientations that correspond to the periodic mode, with the same attracting periodic meandering orbit. This range, marked in Fig. 8 by green dots, is approximately contained between two black solid lines. Their meaning will be explained in the next section by comparing with the Jeffery orbits [49]. Close to the boundaries that separate different dynamical modes in Fig. 8, some irregularities appear and a relatively small change of initial orientations can result in a different dynamical mode. An example of such behavior for three close initial orientations was shown in Fig. 4. This sensitivity is related to chaotic properties of the transition between the modes.
We have analyzed the fiber dynamics for a wide range of values of the bending stiffness A. The periodic meandering solutions are observed for a certain range of A. In addition, for values of A not far from this range, close-to-periodic meandering motions are observed. Examples of such motions are shown and discussed in Appendix B.
In Fig. 8(d), the dependence of the attracting solutions on the bending stiffness A is shown for the initial orientations with 0 = 10 • and different values of 0 . The attracting periodic or close to periodic meandering motions (green dots) are visible for 9 A 12. Examples of close to periodic meandering solutions are shown in Appendix B. We observe that for A = 7.3, some of the initial orientations are attracted to another three-dimensional periodic orbit, corresponding to the squirming solution (marked by green diamonds). This finding illustrates that the squirming periodic motion can be transient or attractive, depending on the value of the bending stiffness ratio A.
The results shown in Fig. 8 indicate that there exist welldefined ranges of the fiber bending stiffness in which different solutions dominate. This result is significant for dilute suspensions of flexible fibers where hydrodynamic interactions between fibers are weak. In particular, our results [ Fig. 8(a)] predict the existence of an ordered phase with all (or almost all) the fibers straightened out and parallel to each other (and to the vorticity direction), for a narrow range of values of the ratio A of bending to hydrodynamic forces (4 A < 6 for N = 40). This range can be reached (or avoided) by adjusting the shear rate. Therefore, it should be possible to control the ordering of flexible fibers in dilute systems.

VI. COMPARISON WITH JEFFERY ORBITS
In this section, we will demonstrate that the rolling and tumbling modes can be interpreted as effective Jeffery orbits of the flexible fiber end-to-end vector n [see Fig. 1(b) for the notation], with an effective amplitude that depends exponentially on time. We will also show that the meandering mode is essentially different from the Jeffery solution. For the sake of clarity, we will focus on a fixed value of the bending stiffness A = 10. The results can be easily generalized for a wider range of values.

A. Jeffery orbits
To compare, we will first remind a reader of the classical Jeffery equations of motion of a rigid prolate spheroid with an aspect ratio r, oscillating periodically in shear flow [49,63]. We will keep the same parametrization of the spheroid orientation [ J (t ), (t )] as for the end-to-end vector of the flexible fiber, and assume for simplicity that tan J (0) = C 0 , (0) = 90 • . The period of the motion, in units ofγ , is and the evolution of the orientation angles satisfies We will use the relation, which follows from Eq. (12), cos 2 (t ) = r 2 sin 2 rt 1+r 2

B. Rolling and tumbling modes
In this section, we will briefly recall basic features of Figs. 2(a) and 2(b) and then analyze them both qualitatively and quantitatively, in a more general context. The characteristic variables and parameters in the rolling and tumbling modes will be labeled by (r) and (t), respectively. For both modes, the time dependence of the fiber end-to-end orientation angle (t ), shown in Figs. 2(a) and 2(b), illustrates two important features. First, the characteristic oscillation time is almost constant in time. Second, the amplitude of the oscillations changes monotonically with time.
We highlight important features of the dynamics in Fig. 9. In the inset we confirm that the oscillation time is practically time independent, with T ≡ T r = 459 for the rolling mode, T ≡ T t = 181 for the tumbling mode, respectively. In the main panel of Fig. 9, we plot ln | tan | versus time and obtain linear dependence of the maxima and minima on time. In this way we demonstrate that for the rolling mode, the amplitude of tan decays exponentially while for the tumbling mode, the amplitude of tan grows as illustrated in Fig. 9 by straight lines of the corresponding inclinations. We have checked that the values of the characteristic oscillation times T and the relaxation times τ are practically the same for other initial conditions leading to the same mode. Moreover, we will now demonstrate that with the periodic function J (t ), which will be matched to the Jeffery solution given by Eq. (12), and negative or positive values of τ given by Eq. (15a).
Matching the period T of J (t ), given by Eqs. (14b), with the Jeffery period of a spheroid with the aspect ratio r, given by Eq. (10), we obtain (for A = 10) r = r r = 73.0 as "the effective aspect ratio" of the flexible fiber in the rolling mode. Note that this value is significantly larger than the geometrical aspect ratio of the deformed fiber (the spring constant k was chosen to be so large that the fiber is practically inextensible). Moreover, the shape of the fiber significantly changes in time, while the period of the damped oscillations and 'the effective aspect ratio' remain almost constant in time, decreasing only by a few percent. For the tumbling mode, the effective aspect ratio following from Eq. (10) is much shorter, r = r t = 28.8. 2 2 The tumbling motion in the shear plane of the end-to-end vector of a flexible fiber was matched with an effective Jeffery orbit in Ref. [51]. In Fig. 10 we compare J (t ) = exp(−t/τ ) tan (t ) and cos 2 (t ) evaluated numerically (solid lines) with the effective Jeffery Eqs. (11) and (13) (dashed lines), in which time is shifted back (by 5 and 13 units) to match the initial conditions and to compensate for transients, including the small decrease of T with time in the first stage of the evolution. Also, for the same reason, we use the fitted values of C 0 given in the caption of Fig. 10. They slightly differ from the corresponding initial values: C(0) = 0.0152 for the rolling mode and C(0) = 0.501 and for the tumbling mode.
For the rolling mode, the numerical and analytical curves are superimposed. For the tumbling mode, there appear some differences, but still J (t ) can by approximated remarkably well by the Jeffery solution.
In this way, we arrive to one of the main conclusions of this paper. The rolling and tumbling modes can be interpreted as "effective Jeffery orbits." The word "effective" means that the standard Eqs. (11)-(13) for the Jeffery trajectories still hold, but now with a constant value of the parameter C 0 replaced by the time-dependent amplitude C(t ), with the negative and positive values of τ corresponding to two different attractors of the dynamics. The decay of C(t ) to zero for the rolling mode and the growth to infinity for the tumbling mode describe, respectively, the rates of convergence to the rolling stationary state and to the periodic two-dimensional tumbling restricted to the shear plane.
To analyze the two different patterns of the evolution, and inspect the transition between them, we investigate other initial conditions. In particular, in Fig. 11 we choose the initial values of the amplitude C much closer to each other, i.e., C 0 = 0.0163 and C 0 = 0.268. The different time dependence of C in both modes is clearly visible in Fig. 11 where we use the simulation results to plot versus . In this plot, the change of C can be traced along the trajectories as the consecutive values of arctan at = 90 • . For different initial conditions with C 0 > 0.25 and C 0 < 0.0195, the evolution in time is qualitatively the same as shown in Fig. 11. In the next subsection, we will discuss the dynamics for 0.0195< C 0 < 0.25, which corresponds to the empty space that separates the rolling and tumbling modes in Fig. 11.

C. Meandering modes
The essential property of the periodic meandering solution is that it differs significantly from the Jeffery orbits [49]. In Fig. 12, we plot (in green) (t ) versus (t ) for the meandering mode, evaluated numerically for 0 = 10 • , 0 = 10 • . The corresponding initial value of C is C 0 = 0.0307, with r = r m = 117, as evaluated from Eq. (10) for T m = 735. It is clear that the meandering trajectory is much more complicated than the periodic Jeffery orbit. In this case there is no systematic drift of the trajectory, and therefore the concept of "an effective Jeffery orbit" cannot be applied.
In Fig. 12, we also plot two black solid lines that follow from the Jeffery relation Eq. (11) as in Fig. 11. The meandering trajectory is almost everywhere located inside the region determined by these lines.
We investigated numerically the evolution of flexible fibers with different initial orientations. For most of the initial conditions with 0.0195 < C 0 < 0.25, we have found the meandering mode, which fills the empty space that separates the rolling and tumbling modes in Fig. 11.

D. Physical interpretation
The analysis of the effective Jeffery orbits and the meandering modes, performed in the previous subsections, provides a simple physical explanation for the dependence of the modes on the initial orientation, determined in Fig. 8(b) for a moderate bending stiffness ratio A = 10. In this diagram, the phase space of the initial orientation angles separates into three distinct regions, each of them leading to a different dynamical mode, and finally to a different dynamical attractor. The borders between these regions are separated by the two special Jeffery solutions with C 0 = 0.0195 and C 0 = 0.25, plotted as the black lines in Figs. 8, 11, and 12. The effective Jeffery solutions, convergent to the tumbling and rolling modes exist for a sufficiently large and a sufficiently small value of C 0 , respectively. However, they do not exist for a range of intermediate values of C 0 . In this range, periodic or close to periodic motions dominate.
This generic classification of the modes is perturbed by the existence of transient, close to periodic motions for some of the orientations and bending stiffness ratios A, as illustrated in Figs. 3 and 4.
In future work we hope to provide analytical justification for these new observations and quantitative results.

VII. DISCUSSION AND CONCLUSIONS
One of the main results of this paper is that there exist periodic and close to periodic three-dimensional motions of a flexible fiber in shear flow. The meandering periodic solution is an attractor for a certain range of the relative bending stiffness A. The squirming motion is a transient in this range, and an attractor in a narrow range of smaller values close to A = 7.3. Our results indicate that a change of the relative bending stiffness A may trigger a transformation between different periodic (or close to periodic) solutions, and the corresponding change of the characteristic sequence of three-dimensional fiber shapes. We have demonstrated that the complexity of a flexible-fiber shape (and in particular the number of local maxima of the curvature) may change significantly with time. A more detailed study of periodic and close to periodic motions of flexible fibers in shear flow, for a wide range of the bending stiffness A, will be presented elsewhere.
We have shown also that the time-dependent rolling and tumbling modes of a flexible fiber in shear flow can be interpreted as effective Jeffery solutions, with the constant C 0 replaced by an exponential function of time C(t ) given by Eq. (17). Therefore, the effective Jeffery orbits drift towards one of two attractors: the fiber aligned with the vorticity direction and the fiber performing periodic motions entirely in the shear plane. There exist two thresholds for the initial value of C(t ): small C(0) lead to the rolling mode, intermediate C(0) to the meandering mode, and large C(0) to the tumbling one. These thresholds are sensitive to the bending stiffness ratio A.
Unlike the slender body theory [1,4], the HYDROMUL-TIPOLE method used in this study takes into account the fiber thickness. This feature is important to study the threedimensional dynamics in shear flow, and in particular, the effective Jeffery orbits and the tumbling motion, in contrast to the elastica models [13,14] that predict an infinite tumbling time of infinitely thin fibers. Finally, we comment on possible comparisons with experiments. In Refs. [3,4], quasi-2D trajectories of actin were investigated, with the focal plane of the motion perpendicular to the vorticity direction. The elastoviscous numberμ = 8(L/d ) 4 /A ≈ 2.2×10 6 , i.e., the largest value ofμ analyzed in Ref. [4], corresponds to A = 10 for our system. Indeed, for A = 8 we recover a shape evolution similar to Movies 5 and 6 from Ref. [4], with U turns and S turns that appear irregularly.
However, we observe that more stiff fibers, after a long time, approach the periodic tumbling motion, with the S-turns only, as illustrated in Fig. 13. These shapes are different from the C-turns observed in Ref. [4]. A more detailed study is needed to understand the reason for such a difference. The obvious differences are that the fibers we study do not perform Brownian motion and are much shorter than actin. Moreover, in Ref. [4], the simulations were performed in 2D and the experiments were carried out in a bounded geometry. Finally, we note that even a very small curvature of the external flow can have a drastic influence on flexible fiber shapes [20].
As a last remark, we suggest searching for threedimensional, periodic motions in experiments. Our estimates indicate that they appear in the range of the bending stiffness accessible in experiments. With this goal in mind, we highlight the need to perform video recordings in the xy (flowvorticity) plane where the characteristic repeatable evolution of shapes should be visible. We hope to report the results of such experiments in a future communication.

APPENDIX A: PERIODIC MEANDERING ATTRACTORS
In this section, we provide more details about the periodic meandering solution. In Fig. 2(c) the maximal values of the relative length of the end-to-end vector are close to one. In the figure, by looking at the orientation (t ) we observe that during the periodic meandering motion, the fiber straightens while tumbling at an orientation that is not along the flow. This property is also visible in Fig. 7 where fiber shapes are shown. In particular, for the periodic meandering solution with A = 10, the fiber is almost straight for = 0 • and ≈ 78 • . In contrast, for the tumbling solution, the fiber straighten along the flow, i.e., at = 90 • .
In the meandering motion, every bead performs a periodic orbit superposed with translation along the shear flow with a velocity v a , equal to the mean velocity of the fiber centerof-mass, averaged over the period. The periodic meandering motion of the center-of-mass for A = 10 is shown in Fig. 14 and Supplemental Material Movie 4 [61], in the frame of reference translating along the shear flow with velocity v a . Here, A = 10, and the initial conditions are the same as in Figs. 2, 3, 6, 7, and Supplemental Material Movie 2 [61]. The time range is 99 188 t 99 988, with the period T m = 735. Owing to symmetry, trajectories of the end-to-end vector, displayed in Fig. 6, close after T /2. Our results show that the meandering period is longer than the transient squirming period and much longer than the tumbling period.

APPENDIX B: CLOSE-TO-PERIODIC MEANDERING ATTRACTORS
In Fig. 15 we illustrate that by changing values of the bending stiffness A away from A = 10, we observe close to periodic (rather than periodic) meandering motions, which are long-lasting and vary similarly to the periodic meandering solution for A = 10. For A = 9, departure from periodicity has a different pattern: the system oscillates back and forth between the twin meandering orbits, symmetric with respect to δx → −δx and δz → −δz.
The end-to-end periodic meandering trajectory, shown in Fig. 6, is nonsymmetric. Naturally, after rotation by π around the y axis, one obtains another periodic orbit. Both appear in our simulations (compare Figs. 3, 4, and 6). Moreover, we observed long time effects of flipping from one to the other trajectory, as shown by the example in Fig. 16 (again, for 0 = 0 = 10 • ). The periodic meandering mode for A = 10 seems to be stabilized already for times around 2000. However, as shown in Fig. 16, there appear some long-time effects: the trajectory later flips to another one, symmetric with respect to rotations along the vorticity direction (i.e., with respect to δx → −δx and δz → −δz), and later remains in this shape.

APPENDIX C: ACCURACY OF THE SIMULATIONS
The Stokes flow dynamics of an elastic fiber modeled as a chain of beads, evaluted using the HYDROMULTIPOLE numerical algorithm, is more accurate for a larger multipole truncation order L and a smaller distance l 0 between the centers of consecutive beads. Numerical accuracy of the simulations of flexible fibers in shear flow, performed with the HYDROMULTIPOLE numerical codes [56], was analyzed in detail in Ref. [51] for the periodic tumbling motion of fibers located in the xz plane. The accuracy of the truncation at L = 2 was estimated there as 3%, and for the l 0 = 1.02 approximation as 10-15%. It is essential to keep the gap between consecutive bead surfaces small enough to prevent the beads from rotating under shear independently from each other. For small gaps, the lubrication effects serve this goal quite effectively.
Similar accuracy estimates are expected to hold for the three-dimensional motions considered in this article. However, in this case, the analysis of the accuracy is more subtle because we have already observed that the existence of FIG. 17. The long-time trajectories of the fiber end-to-end vector. Comparison between fibers with different values of l 0 and A (as indicated), for the same initial orientation ( 0 , 0 ) = (10 • , 10 • ). The plots are made for 5 000 t 10 000 (upper row) and for t end − 5000 t t end (lower row), where t end = 5 · 10 4 for A = 9 and 11, and t end = 10 5 for A = 10. See also Table I. periodic meandering solutions is limited to a narrow range of the bending stiffness. A small departure from this range gives rise to motions that are close to periodic (for example, see Fig. 15). Therefore, taking into account the chaotic nature of flexible fibers dynamics in shear flow [51], the basic question is if three-dimensional periodic meandering motions are also observed for more precise (but computationally more demanding) bead models of a flexible fiber, with a smaller l 0 . Therefore, we choose A = 9, 11 as the values of the bending stiffness close to A = 10, which for l 0 = 1.02 corresponds to the meandering attractor with the period T = 735. In Fig. 17 and Table I we compare the long-time solutions obtained for l 0 = 1.005 and l 0 = 1.02. While changing l 0 , we fix values of N and k, and consider the identical initial orientations 0 = 0 = 10 • .
As shown in Fig. 17 and Table I, for l 0 = 1.005, the periodic meandering solution is not recovered at A = 10; instead, for these values of the parameters, we find the close-to-periodic meandering motion, with a smaller maximal time τ = 636 between the consecutive tumblings. The decrease of the tumbling time for smaller values of l 0 was also observed for the tumbling periodic orbit [51], and explained as the influence of the bead rotation. To minimize this effect, the gaps between the bead centers should be kept reasonably small. In practice, l 0 1.02 is sufficient.
The long-time periodic meandering solution, detected for l 0 = 1.02 and A = 10 with the period T = 735, is observed also for the smaller distance l 0 = 1.005, but at a larger value of the bending stiffness A = 11, and with a smaller period T = 653. Supplemental Material Movie 5 [61] illustrates that the evolution of shapes in these two periodic solutions is almost identical (here time is normalized by the two different periods, respectively). We also find that the long-time meandering quasiperiodic solutions for l 0 = 1.02 are similar to those obtained for l 0 = 1.005, but with larger values of the bending stiffness A, as shown in Fig. 17 and Table I.