Highly elastic fibers in a shear flow can form double helices

The long-time behavior of highly elastic fibers in a shear flow is investigated experimentally and numerically. Characteristic attractors of the dynamics are found. It is shown that for a small ratio of bending to hydrodynamic forces, most fibers form a spinning elongated double helix, performing an effective Jeffery orbit very close to the vorticity direction. Recognition of these oriented shapes, and how they form in time, may prove useful in the future for understanding the time history of complex microstructures in fluid flows and considering processing steps for their synthesis.

The fiber flexibility increases with an increase of the fiber aspect ratio l/d, the fluid dynamic viscosity η, and the shear rate γ, or a decrease of Young's modulus E. Recently, it has been shown [19] theoretically and experimentally, that for slender, moderately elastic fibers morphological transitions between rigid-like tumbling, springy rotation (C-buckling), and snake orbits (U-turn) take place at certain critical values of a single parameter, the so-called elasto-viscous number μ = 2 9 (l/d) 4 η γ/E.
However, a systematic exploration of the bending patterns of more flexible fibers has not appeared.Examples of their complex dynamics are known, and the formation of compact structures has been shown [13,16,29,31].Hence, in this work, we study the dynamics of highly elastic fibers in a shear flow.Our main result is to show experimentally and numerically that very flexible fibers can form elongated and thin double helices that are oriented close to the vorticity direction.Such shapes and ordering of their orientation have not been reported in the literature so far.We observe that in the simulations, around one-third of sufficiently flexible fibers after a long time bend into such a double-helix shape.In the experiments, there are even more (one-half) double-helix shapes at long times.
The plan of the paper is the following.The experimental system and results are shown in section 2.Here we demonstrate the double helices and describe their shape and dynamics.The theoretical model and numerical method are described in section 3, after which the numerical results are presented in section 4.
Here we quantify the double-helix shape and describe its rotation around the main axis and Jeffery orbits.In section 5 we demonstrate how elastic fiber forms a double helix in the simulations and experiments.In section 6 we discuss the results and present the conclusions in section 7. The appendices contain a description of the supplementary movies, a specification of the initial configurations of the simulations, and a presentation of other shapes of highly elastic fibers.

Experimental setup
We employed a standard microfluidic method for fiber fabrication, as detailed in [34,42,43].This process involves generating a uniform cylindrical jet of oligomer through flow focusing within a microfluidic channel.The oligomer solution primarily comprises water, poly(ethylene glycol) diacrylate (PEG-DA), and a photoinitiator.Exposure of the jet to ultraviolet (UV) light pulses induces gelation, resulting in fiber formation.This technique enables precise control over fiber properties, including length, diameter, and Young's modulus E.
The experimental parameters for fiber fabrication align with those in [34,43], with two notable deviations: the concentration of 2-hydroxy-2-methylpropiophenone (photoinitiator; Sigma-Aldrich) in the oligomer solution is 4%, as opposed to 8%; also, the exposure time for UV exposure is 10 ms, rather than 20 ms.The resulting fibers have a diameter of d = 34.8± 3.3 µm and a length of l = 1.31 ± 0.07 mm; the corresponding aspect ratio is l/d ≈ 40.Due to the reduced photoinitiator concentration, the fibers exhibit a lower Young's modulus E of O (10 4 ) Pa and at rest tend to be slightly curved and spindle-like near the tips.
To study the dynamics of fibers under shear flow, we used a rheo-optical setup similar to that described in [34,43].A dilute suspension of fibers (approximately 0.5% by volume) in PEG-DA (viscosity η = 55 mPas) was placed between two transparent glass parallel plates (radius R = 21.5 mm) of a rheometer (MCR 702; Anton Paar).The gap size of the parallel plates was h = 0.5 mm.The fibers, which are neutrally buoyant in this solution (density ρ = 1.12 g cm −3 ), were illuminated by a white light-emitting diode (LED) positioned above the plates.The fibers are slightly confined between the plates, and this narrow gap size is due to the necessity to keep a small Reynolds number and achieve good-quality imaging.The dynamics of the fibers were captured by a high-speed camera (v7.3;Phantom; sampling frequency 500 Hz) equipped with a 5 × objective (Mitutoyo) and positioned beneath the glass plates.The field of view was 2.57 mm × 3.43 mm (600 pixels × 800 pixels) while its center is 15 mm away from the center of the plates.
A constant shear flow was applied between the counter-rotating plates, with a shear rate of γ = 523 s −1 at the center the field of view (±48 s −1 across the field of view).The streamline is slightly curved, leading to a slow migration of the fibers toward the center of the plate, with a long dimensionless time scale of O(10 5 ).The resulting bending stiffness ratio was A = E/ (64η γ) = O (1).The effective Reynolds number, defined as suggests that inertial effects are negligible.Alternatively, we can estimate the particle-scale Reynolds number by using the diameter of the fiber d.The corresponding Reynolds number, is also small.The flows generally remained stable throughout the experiments, lacking the high-frequency torque fluctuations characteristic of 'elastic instability' in sheared polymer solutions [44][45][46][47][48][49].
Experimental observations below are reported in terms of fiber dynamics projected onto the xy plane, along with dimensionless time t, calculated by dividing real-time by 1/ γ.

Observation of double-helix fibers
Under shear, flexible fibers displayed complex shapes and dynamics, as shown in figure 1.The imaging shows the projection of the fibers on the xy plane, where x is along the flow and y is along the vorticity direction.
Interestingly, a subset of fibers exhibit double-helix shapes, as highlighted in figure 1.These double-helix configurations were relatively stable over time under shear.We tracked three individual double-helix fibers to analyze their dynamics under shear, as illustrated in figure 2 and corresponding movies 2, 3, and 4.
A time sequence for one such fiber is shown in figure 2(a) and highlights changes in both shape (from t = 1065 to 1074) and orientation (at t = 1065 and 1087).To quantify these dynamics, we measured the distance from a helix tip to its first cross-point, l r , and the angle of the fiber's tip-to-end direction relative to the y-axis, θ.Figures 2(b)-(d) show the detailed measurements for the three individual double-helix fibers.These data reveal similar dynamic patterns for the different fibers.Typically l r increased over time until the formation of a new cross-point near the tip.Both l r and θ exhibited periodic changes, with a faster frequency for l r and a slower frequency for θ.
Our observations reveal two motions of double-helix fibers: rapid 'waving and rolling' [period of O( 10)] and slower Jeffery-orbital movements [period of O(100)].Notably, the deformation of the double helix propagates as a wave from tip to end.

Theoretical description
We consider an elastic fiber in a shear flow, of a fluid with dynamic viscosity η.We assume that the Reynolds number is much smaller than unity and Brownian motion is negligible.An elastic fiber is modeled as a chain of N identical, spherical, solid beads with diameter d.The elastic interactions of the beads are determined by harmonic stretching and bending forces.At elastic equilibrium, the distance between centers of two consecutive beads is equal to ℓ 0 = 1.02d, and in most cases, the fiber is straight.To model the experimental conditions, we also evaluated the evolution of fibers that form a half-circle (C-shape) at the elastic equilibrium (see appendix B).
Consecutive beads are connected by springs with the spring constant k.The potential stretching energy of the fiber has the form where ℓ i = |r i − ri−1 | is the distance between the centers of consecutive beads i and i−1.Here, the position of the center of bead i is denoted as ri , as shown in figure 3.
In addition, each triplet of the consecutive beads resists bending, with the bending stiffness Â equal to the product of Young's modulus E and the second moment of area, ( The bending potential energy of the fiber is given as with the bending angle β i determined by the relation cos ), and shown in figure 3.
To avoid long-lasting spurious overlaps, we introduce a small repulsion interaction between non-consecutive overlapping beads, with the potential energy where rij = |r i − rj | and B is the repulsion strength.The external force F i acting on bead i has the form, In the following, we will use dimensionless variables.We choose d and 1/ γ as the length and time units, respectively.The dimensionless positions of the bead centers are denoted as r i = ri /d.The normalized spring constant k, bending stiffness ratio A, and repulsion coefficient B are respectively.To model an almost inextensible fiber, we choose a large spring constant k = 1000.We assume a weak repulsion, with B = 0.1.We consider the number of beads N = 40 or 60 and the bending stiffness 0.
The positions r i of the bead centers satisfy the system of first-order ordinary differential equations [32][33][34]39], where is the rate-of-strain tensor of the simple shear flow v 0 given by equation ( 3).
The translational-translational and translational-dipolar mobility matrices, µ tt ij and µ td ij , respectively, are long-ranged functions of positions of all the bead centers.In general, they depend on the boundary conditions.Here the no-slip boundary conditions on the surfaces of the beads are assumed.The mobility matrices are evaluated by the multipole expansion of the Stokes equations, [50][51][52][53], corrected for lubrication to speed up the convergence [52][53][54], and implemented in the HYDROMULTIPOLE numerical programs of a controlled, high precision [52,53].
In [34], simulations of moderately elastic filaments in a shear flow were performed for initially straight shapes, for the whole 3D range of the orientations [34].As a result, in [34] only odd shapes were obtained.(Here 'odd' means that the fiber deformation is an anti-symmetric function of position along the main fiber axis, relative to its center.)To observe even shapes, an even initial perturbation is needed, as in the numerical   simulation of [19], and in the spectral analysis and numerical results of [55].Also, in [34] (see figure 9 there) it was tested that small deviations from straight initial configurations can lead to different, non-planar buckled shapes, in agreement with the results for buckled shapes of slender fibers in the compressional flow [13].Therefore, in this work, the fibers are initially straight at different orientations (Θ 0 , Φ 0 ) with a small random perturbation, as illustrated in figure 4(a), or form planar bent shapes shown in figure 4(b).The details of the initial conditions are given in appendix B.
The equations of motion (10) are solved numerically by the adaptive fifth-order Runge-Kutta procedure for dimensionless times up to 8000.The same numerical codes were applied in our previous studies of flexible filaments in a shear flow [32-34, 39, 56].

Double-helix shape, orientation, and rotation
A double helix shape of a highly elastic fiber, observed in our experiments, is also detected in our numerical simulations.It forms spontaneously after a sufficiently long dimensionless time (typically, around t = 1000), and keeps this shape until the end of the simulation (at t = 8000).
An example with N = 40, A = 0.6, Θ 0 = 70 • , Φ 0 = 90 • , and t = 3500 is shown in figure 5. Here, we took approximately the same aspect ratio N as in the experiments.
Figure 6 illustrates that the elongated double helix shape, almost aligned along the vorticity direction, is maintained for a long time.Here, N = 40, A = 0.6, Θ 0 = 80 • and Φ 0 = 175 • .The colored positions of 4 beads at the snapshots of different times in figure 6(c) indicate that the helical shape rotates.The rotation is also visible as almost periodic curves in figures 6(a) and (b) where we investigate time dependence of the end-to-end vector, r e = r N − r 1 , which links the centers of the first and the last bead, as illustrated in figure 5.
In figure 6(a) we plot the length of the end-to-end vector, L e = |r e |, normalized by its value L 0 = (N − 1)ℓ 0 at the elastic equilibrium, and in figure 6(b), the polar angle Θ e , defined in figure 5.Both L e /L 0 and Θ e are almost periodic functions of time.Figure 6(a) illustrates that the double helix shape is not rigid-it fluctuates almost periodically at the time scale equal to half of the rotation period τ r = 22.The small shape oscillations found in the numerical simulations are in agreement with the experimental results from section 2.
Rotation and shape oscillations during half of the period for the same simulation are studied further in figure 7. The xz projection of the fiber shape is shown in figure 7(a).The double helix projection resembles an ellipse, with the smaller axis parallel to the flow gradient, i.e. to z.The aspect ratio of the ellipse fluctuates in time.

More detailed description
In general, the orientation of the fiber in the numerical simulations is close to, but not identical to the vorticity direction, and it changes with time, as illustrated in figure 8 and movie 7.This effect is also visible in the experiments (see figures 1-2 and movies 1-4).
Therefore, to describe the fiber dynamics quantitatively, we evaluate the time-dependent inertia tensor Î in itcenter of mass frame [57], where α = x, y, z, β = x, y, z, and r ′ j = r j − r cm is the position of j-th bead in the center-of-mass reference frame, with the center-of-mass position r cm = N i =1 r i /N, and r ′ j = (r ′ jx , r ′ jy , r ′ jz ).We calculate the eigenvalues I 1 , I 2 , I 3 of Î, and the corresponding normalized eigenvectors n 1 , n 2 , n 3 , which are perpendicular to each other.We assume, without loss of generality, that I 1 < I 2 < I 3 .Then, n 1 represents the normalized eigenvector with the smallest eigenvalue, and it determines the principal axis of the fiber.The double helix is elongated along n 1 , which depends on time; n 1 (t) is described by the spherical angles Θ(t) and Φ(t), as shown in figure 8 6 .
The double helix is a long, thin, three-dimensional, non-planar structure.To determine how much it deviates from a planar structure, we project the positions r ′ i of all the beads i = 1, . .., N on n 3 , which is the unit eigenvector corresponding to the largest eigenvalue of the inertia tensor, and then evaluate the so-called non-flatness parameter as The time-dependent non-flatness δ of the double helix is plotted in figure 7(b).It is small, just a fraction of the bead diameter d, and changes periodically with time, according to the change of the aspect ratio of the   13) ; (c) bead-to-bead projection b, given by equation ( 14); (d) azimuthal angle ϕ i , given by equation (15).See movie 7 for the time evolution of the xy projection of the double helix.double helix cross-section, shown in figure 7(a), both triggered by the rotation.The oscillations of δ are very small, with the amplitude of the oscillations 0.05.Now we move on to a more detailed analysis of the double helix dynamics, focusing on a highly elastic fiber with a larger aspect ratio N = 60, and a similar bending stiffness ratio A = 0.5.
The xy projections of shapes are shown in figure 9(a) at approximately equal time intervals τ r /7.As for N = 40, the shape is thin and elongated, oriented close to the vorticity direction y.The shape rotates with small, almost periodic changes.But the 16th and 17th of the beads at the top (i.e. with the largest value of coordinate y), and the 41st and 42nd of the beads at the bottom (i.e. with the smallest value of coordinate y) do not change with time.For the initial orientation Θ 0 = 20 • and Φ 0 = 90 • , shown in figure 9, the fiber rotation period τ r = 14.
To investigate the motion and shape in detail, for each bead i = 2, . .., N − 1 we evaluate the time-dependent local curvature, In figure 9(b), colors are used to represent the local curvature vs the bead label i and time.Clearly, κ i is the largest at the beads i = 16, 17 and i = 41, 42 with the smallest and the largest y coordinates where the chain bends significantly to allow for the bundling of its other parts.
The regular pattern of the colors in figure 9(b) indicates that the local curvature at each bead changes with time periodically.This can be understood by taking into account that the double helix is systematically squeezed along the direction of the flow gradient while rotating, as illustrated in figure 7(a) at the yz projection.Therefore, the period of the oscillations of the local curvature of each bead is synchronized with the half-period of the rotation of the whole fiber.Actually, the rotation period τ r can be determined as the time separation between every second stripe of a given color in figure 9(b).The inclination of these stripes carries the information of the twist, and therefore, the number of coils in the double helix, in this case around 6, in agreement with the projection of shape shown in figure 9(a).
A similar regularity of colors is observed in figure 9(c) for a different local characteristic parameter of the fiber shape, i.e. the bead-to-bead projection b i , i = 2, . .., N, of vector r i − r i−1 on the unit vector n 1 , The bead-to-bead projection rapidly changes sign at the turning points of the fiber, where it bends significantly to bundle its arms, as shown in figure 9(c).The labels of the beads at the turning points do not change with time.
To trace the rotation, for each bead i = 1, . .., N one can introduce an azimuthal angle ϕ i between the projection of r ′ i on the xz plane and the x axis, where e x and e y are the unit vectors along the x and y axes, respectively.The color map of ϕ i versus time t and the bead label i is shown in figure 9(d).The same period of rotation and the same inclination of the color stripes can be detected from this plot.However, the color pattern is significantly less regular than for the local curvature and bead-to-bead projection in figures 9(b) and (c).The irregularity of colors in figure 9(d) is related to a small difference between the vorticity direction y and the orientation n 1 of the fiber principal axis, as illustrated in figure 9(a).
The fiber rotates around n 1 rather than y.To show it, in figure 10, we introduce the Cartesian coordinates x and z and the azimuthal angle φ in the plane perpendicular to n 1 .We select one bead, with the label 14, and trace its evolution in time.In figure 10(a), we plot the azimuthal coordinate of this bead, φ14 , vs time, and obtain a periodic function.The perpendicular to the n 1 projection of the trajectory of the bead 14 is shown in figure 10(b).The rotation is visible remarkably well, with only small trajectory deviations from a circle.Such deviations are in agreement with the small deformations of the fiber shape discussed before.
Their elongated shape suggests that double helices should perform an effective Jeffery motion [58].Therefore, we will now focus on the time dependence of the double-helix orientation, i.e.Θ and Φ.In figure 11 we compare it with Jeffery orbits of a spheroid with the aspect ratio r [51,58,59], Table 1.Characteristic parameters of double helices formed in our simulations for fibers with N = 60, A = 0.5 and different initial orientations Θ0, Φ0: average length Lav, average thickness Dav, rotation period τ r , and parameters of the Jeffery orbits: the period τ of the Jeffery orbit, effective aspect ratio R eff , the constants C b , Ce at the beginning of the periodic motion and at the end of simulations, respectively.By matching the periods, T = 2π(r + 1/r), we determine the effective aspect ratio of the spheroid, r = R eff = 13.5.Figures 11(a) and (b) illustrate that for the double helix, Θ(t) changes almost periodically as for a rigid spheroid, but with a very slow drift: the Jeffery constant C slowly changes with time, from C = 0.007 in figure 11(a) to C = 0.01 in figure 11(b).The slow change of C results in a relatively thick stripe of the points of the Θ(Φ) trajectory of the double helix in figure 11(c).Most of the points are approximately located between the solid lines for a rigid spheroid with C = 0.007 and C = 0.01.A slow drift of effective Jeffery orbits of flexible fibers with different shapes was reported earlier in [32].
So far, we discussed the basic features of the dynamics and shape of a double helix using two examples with different aspect ratios.In table 1, we compare the characteristic parameters of 6 regular double helices formed in our simulations with N = 60, A = 0.5, and different initial orientations.
First, we characterize the double helix shape.We evaluate the time-dependent length L and thickness D of the double helix.The length is defined as the maximum distance between the bead centers along the main axis n 1 : and the thickness is calculated as the average distance of a bead center from the main axis, Then, we evaluate L av = ⟨L⟩ t and D av = ⟨D⟩ t as the averages of L and D over time 7000 ⩽ t ⩽ 8000 (i.e. the last 1000 time units of the simulation), and list the resulting values in table 1.We also list values of the rotation period τ r -the period of spinning around the double helix main axis, parallel to n 1 , and characteristic values of the Jeffery orbits, i.e. the period τ of the Jeffery orbit, effective aspect ratio R eff of a rigid spheroid with the same period, and Jeffery parameters C b and C e , corresponding to different time ranges-at the beginning of the periodic motion and at the end of simulations, respectively.The results in table 1 indicate that the drift of Jeffery orbits is sometimes towards the vorticity direction (C e < C b ), and otherwise away from it (C e > C b ).
For regular double helices, we recover approximately the same ratio of the Jeffery period to the rotation period, τ /τ r ≈ 6, as in the experiments.Table 1 illustrates that some of the Jeffery periods are very long.They correspond to irregular double helices, with the double helix shape restricted to only one or two parts of the fiber.Examples of irregular double helices performing Jeffery's motion are shown in figure 12.The irregular double helices found in the numerical simulations correspond to the irregular double helices observed in the experiments (compare with figure 1 and movie 1).

How elastic fiber forms a double helix?
After describing in the previous section the numerical results for the attracting double helix mode of the dynamics of highly flexible fibers in a shear flow, we now move on to show how a double-helix shape is created.Typically, the fiber rather quickly forms a compact structure, resembling a 3D paperclip clump, or irregular twisted hairpin, which then slowly becomes more and more elongated, and oriented approximately along the vorticity direction.This evolution pattern is detected in our numerical simulations of a single fiber, and in our experiments with a dilute suspension of fibers.For the numerical simulations, the formation of a double helix for the initially almost straight fiber with Θ 0 = 20 • and Φ 0 = 90 • is shown in the supplemental movies 8 and 9 (the final stage of this simulation is visible in figure 9(a).
Figure 13 and the corresponding supplemental movie 10 illustrate how the initial U-shaped configuration evolves towards a double-helix shape.In particular, figure 13 shows that for such an initial shape, the long arms parallel to the vorticity direction do not bundle immediately as the fiber rotates.Instead, the fiber quickly nearly straightens close to the shear plane (i.e. the plane spanned by the flow and flow-gradient directions) and then coils into a compact symmetric structure that slowly becomes asymmetric and then even slower forms a double helix.This sequence of shapes is the same for fibers that are initially almost straight.The same stages of the bending process are also seen in our experiments.Figure 14 demonstrates statistically the shape changes of the fibers of the dilute suspension.Under shear flow, most of the fibers quickly elongate close to the shear plane and then undergo complex coiling, form a  time-dependent compact structure that slowly elongates and orients close to the vorticity direction.Eventually, after a long time, it often adopt a regular or irregular double helix configuration.
In the experiment, fibers after some time leave the camera field of view, entrained by the shear flow.However, it is possible to trace the evolution of certain fibers for quite a long time.For example, figure 15(b) shows a fiber of compact, irregular, time-dependent shape slowly evolving to a double helix with a larger and larger aspect ratio.A very similar time sequence of shapes is detected in our numerical simulations, as shown in figure 15(a

Discussion
Double helices are formed for different initial orientations and different initial shapes.However, they are not observed in the numerical simulations for initial orientations that correspond to Jeffrey orbits with a very small constant C; such orientations seem to prevent the fiber ends from bundling.Other shapes formed by highly elastic fibers ('zigzags') are presented and discussed in appendix C. Their analogs are also found in our experiments and the early experiments from [29].
We have checked that only highly elastic fibers form double helices.By 'highly elastic' we mean with a low value of the bending stiffness ratio A for a given fiber aspect ratio N. A low value of A can be achieved by a low value of Young's modulus or a high value of the shear rate, see equation (9).In the numerical simulations with a value of A increased twice (e.g. from A = 0.5 to A = 1 for N = 60), double helices have not been detected.A similar conclusion follows from our experiments: we increased A by decreasing the shear rate while keeping the same low value of Young's modulus.The experiments were performed for a range of shear rates 1 s −1 ⩽ γ ⩽ 523 s −1 .When decreasing the shear rate from γ = 523 s −1 , the double helices were observed less frequently and were typically less regular.For γ < 20 s −1 , they were absent.
This paper focuses on highly elastic fibers.For highly elastic fibers, the straight configuration in the xz plane is not stable.We have chosen values of the bending stiffness 0.1 ⩽ A ⩽ 1, which is much smaller than a critical value A c (N) = CN 3/2 , with C = 0.029, that is determined in [33,39] by the condition that for A < A c fibers initially aligned with the flow do not align with it again, as they are always bent [31].For the chosen aspect ratios N = 40 or N = 60, A c (40) = 7.34 and A c (60) = 13.5.Therefore, the highly flexible fibers considered in this work are bent all the time.
The chosen ranges of the aspect ratio and bending stiffness correspond to the elasto-viscous numbers μ = 8N 4 /A = 2.3 • 10 7 − 2.1 • 10 8 .Our fibers are significantly more flexible than the fibers studied in [13,19].The essential difference is that the moderately elastic fibers in [13,19] straighten along the shear flow.Therefore, there is no reason to expect that the dynamics of highly elastic fibers depends only on a single parameter, the elasto-viscous number μ.
The presence or absence of double helices seems to depend on both N and A rather than only μ.For example, in our numerical simulations, double helices have been detected for N = 60 and A = 0.5 (μ = 2.1 • 10 8 ) and also for N = 40 and A = 0.6 (μ = 3.4 • 10 7 ), but they are absent for an intermediate value μ = 1.0 • 10 8 for N = 60 and A = 1.The dependence of the dynamics of highly elastic fibers on N will be studied elsewhere.
At long times, we have not observed compact tangled structures reported in [16,29,31,33] for very flexible fibers with larger aspect ratios at short times.In principle, the shear flow can untangle flexible fibers [56].Therefore, it would be interesting to study the long-time dynamics of highly elastic fibers with large aspect ratios and check if the compact shapes untangle and form long-lasting double-helices.
In our experiments, the gap size between the plates, 0.5 mm, is smaller than the fiber length, 1.31 mm.Below we argue that this does not qualitatively alter the double-helix orientation and motion found numerically for the unbounded fluid.First, the gap size is larger than the width of the double helix.Indeed, based on figure 2, the double-helix thickness can be estimated as 150 µm.Moreover, the double helix is oriented close to the vorticity direction y.It follows that the Jeffery orbits are thin in the z direction and wider in the x direction.For example, the orbit shown in figure 11, according to equations ( 16)-( 17), has the width along z only around 0.5d, and the width along x around 6d.For the experimental fiber width, d = 35 µm, the thickness of the orbit along z is below 20 µm, much smaller than the gap size.However, the thickness of the orbit along x is around 210 µm, and therefore, the motion along this orbit is easily observed in the experimental plane xy.The above estimates illustrate that the double helix fits in between the plates.Therefore, in our experiments, it is possible to observe the characteristic features of the double-helix dynamics, even though there appear hydrodynamic interactions with the plates.
In this work, we have not observed an irregular flow, typical of elastic turbulence [44][45][46][47][48][49].This observation might be related to the relatively small fiber aspect ratio, lack of Brownian motion, very small volume fraction, or relatively small curvature of the flow.

Conclusions
In this work, we have analyzed experimentally and numerically the long-time behavior of highly elastic fibers in a shear flow.We have found that such fibers often spontaneously form elongated helical structures, folded and bundled; we call them 'double helices' .We have found the double-helix shapes in many simulations; they are even more common in our experiments.(In particular, for N = 60 and A = 0.5, double helices were detected in 36 out of 103 simulations, as listed in table 2.) Typically, the main axis of a double helix is oriented close to the vorticity direction.There are two characteristic time scales of the double helix dynamics: the period τ r of rotation around the main axis, with small oscillations of shape, and the period τ of Jeffery orbits.
The formation of such shapes, their ordered orientations, and dynamics of highly elastic fibers have not been reported so far.Until now, based on a short-time analysis, very flexible fibers in a shear flow have been expected to form irregular, three-dimensional compact shapes [29,31].Here we show that after a sufficiently long time, a compact structure can escape the shear stresses of the fluid by elongating along the vorticity and becoming thinner in the perpendicular directions.For moderate values of the aspect ratio, often a double-helix shape is formed.Our results indicate that dilute suspensions of very flexible fibers in shear flow after a long time become highly ordered.Double helices of highly elastic fibers in a shear flow, and their ordered suspensions, may have various technological, biological, and medical applications.On the other hand, our study demonstrates that highly elastic fibers can be produced, enabling such applications.Moreover, the highly elastic fibers may also be used in other fluid flows, allowing, e.g. for applications in sedimenting systems [60,61].
Once the fiber is perturbed, the distances between the centers of adjacent beads are changed.Therefore, the initial relative positions of the consecutive beads are rescaled to their elastic equilibrium value ℓ 0 /d = 1.02.
For fixed values of N and A, the amplitude ∆ of the perturbations is estimated in the following way: for each initial configuration, the maximum distance of a bead center from the straight unperturbed fiber is calculated, and then averaged over all initial orientations.For N = 60 and A = 0.5, ∆ = 0.28, and for N = 40 and A = 0.6, ∆ = 0.15, so they are much smaller than the bead diameter d.The positions of the bead centers for all the initial configurations used in our simulations are listed in Repository [62].

B.2. Fibers initially significantly curved
For N = 60 and A = 0.5, we also perform simulations of elastic fibers that initially are significantly curved.The planar C-shaped, U-shaped, and S-shaped fibers, shown in figure 4(b), are used as the initial configurations.The C-shape is a half-circle, the S-shape consists of two half-circles connected by a straight line, and the U-shape is a half-circle with two parallel straight arms.
The shapes shown in figure 4(b) are used as the initial configurations for 23 simulations of the fibers straight at the elastic equilibrium.The initial C-shapes have 9 different initial orientations of the end-to-end vector r e (t=0), given by the following angles : Θ e (t=0) = 70 The positions of all the bead centers at each initial configuration of the fiber described above are listed in Repository [62] in the files C-shape2.zip,U-shape.zip and S-shape.zip for the initial C-shapes, U-shapes, and S-shapes, respectively.
The fibers used in our experiments form sort of a C-like shape in the elastic equilibrium.Therefore, we also perform simulations of fibers that form the C-shape (a half-circle) in the elastic equilibrium, using the same C-shape as the initial configuration.We perform 40 simulations for such a fiber, with its end-to-end vector, defined as r e = r N − r 1 , initially oriented at the following angles : Θ e (t=0) = 10 • , 20 • , 30 • , 45 • , 60 • , 70 • , 80 • , 90 • and Φ e (t = 0) = 90 • , 120 • , 160 • , 170 • , 175 • , with Θ e and Φ e defined in figure 5.For each initial configuration of the fiber, positions of all the bead centers are given in Repository [62] in the file C-shape1.zip.

Appendix C. Other shapes of highly elastic fibers
Our simulations indicate that double helices are not the only attracting dynamical modes for N = 60 and A = 0.5.The other typical evolution pattern we call the zigzag mode.Its characteristic features are illustrated in figure 16.This dynamical mode is also almost periodic, but the shape is different.The fiber forms an elongated helical shape, oriented along y, rather thick in the direction of x and thin in the direction of z.Similar shapes and their orientations have been found in the experiments reported in this work in section 2 and in movie 1, and also in [29] in figure 4.
Differences between the shapes of a double helix and a zigzag are illustrated in figures 17(a) and (c) by comparing the xy projection of their shapes, and in figures 17(b) and (d), by comparing their time-dependent ratio L/D of length L and thickness D, defined in equations ( 18) and (19), respectively.In particular, figure 17 shows that double helices are significantly thinner and longer than zigzags.The characteristic time scales of the shape oscillations, seen in the time dependence of L/D, in both cases are associated with the rotation around the main axis.The rotation period for the double helix, τ r = 14, is much smaller than the rotation period for the zigzag, τ r = 102.
In table 2, we list the initial and final shapes of highly elastic fibers, together with the corresponding numbers of the simulations.Almost one-third of all the simulations lead to the formation of double helices, for all the initial conditions considered here.Some of the double helices from table 2 are less regular than those described quantitatively in table 1, with the double helix restricted to only part of the fiber.Examples of such shapes are shown in figure 12.

Figure 1 .
Figure 1.A cropped snapshot of flexible fibers under shear at normalized time t = 862.Fibers exhibiting regular double-helix structures are highlighted with white dashed circles.Many irregular double helices are also visible.See movie 1.

Figure 2 .
Figure 2. The dynamics of double-helix fibers.(a) Time sequence for fiber 1, with zoomed-in views near its tip.We measure the distance from the tip of a helix to the first cross-point, lr, and the angle from the fiber's tip-to-end direction to the y-axis, θ. (b)-(d) Measurements of lr and θ for three double-helix fibers (fibers 1, 2, and 3).Blue circles represent lr and red squares denote θ.

Figure 3 .
Figure 3. Part of an elastic fiber modeled as a chain of beads, with the notation used.

Figure 4 .
Figure 4. Initial configurations of a flexible fiber in a shear flow used in the numerical simulations.(a) A straight fiber at different orientations (Θ0,Φ0) is slightly perturbed.(b) C-shaped, U-shaped, and S-shaped fibers are planar and not perturbed.

Figure 5 .
Figure 5.A double helix from the simulations.Its rotation can be seen, e.g. by tracing the time-dependent vector re, which links the centers of the first (red) and last (green) beads.The orientation of re is given by the angles Θe, Φe.

Figure 6 .
Figure 6.A double helix, spontaneously formed in our simulation of elastic fiber with N = 40, A = 0.6, Θ0 = 80 • , Φ0 = 175 • , oscillates periodically.a,b) Normalized length Le/L0 and polar angle θe of the end-to-end vector versus time.(c) Shape of the flexible fiber projected on the xy plane at the times indicated by * in (b);.yellow: 1st bead, red: 16th bead, green: 35th bead; magenta: 40th bead.See movie 5 for the time evolution of the xy projection of the double helix.

Figure 7 .
Figure 7. Rotation and oscillations of the same double helix as in figure 6 during half of the rotation period.(a) Shape projected on the xz plane at the times indicated by dots in panel b.(b) Non-flatness δ of the fiber, defined in equation (12), versus time.See movie 6 for the time evolution of the xz projection of the double helix.

Figure 9 .
Figure 9. Shape and periodic rotation of a double helix with N = 60, A = 0.5, Θ0 = 20 • and Φ0 = 90 • : (a) shape projection onto the xy plane, with the red dot marking the first ball in the chain; the beads i = 16-17 are at the top, and the beads i = 41-42 are at the bottom; (b)-(d) the dependence (shown in colors) on time t and the bead label i of: (b) local curvature κ i , given by equation (13) ; (c) bead-to-bead projection b, given by equation (14); (d) azimuthal angle ϕ i , given by equation(15).See movie 7 for the time evolution of the xy projection of the double helix.

Figure 10 .
Figure 10.Rotation around n1 of 14th bead of the fiber from figure 9: (a) the time dependence of the azimuthal angle φ14 of 14th bead in the plane xz perpendicular to n1; (b) the xz projection of the 14th bead trajectory in the same time range.

Figure 11 .
Figure 11.Jeffery's orbits of the double-helix from figure 9 (red dots), compared to Jeffery orbits of a rigid spheroid with the aspect ratio R eff = 13.5 and C = 0.007 (black line) and C = 0.01 (blue line).The angle Θ: (a) versus time, earlier; (b) versus time, later; (c) versus angle Φ.In (c), red points are plotted for times 2400 ⩽ t ⩽ 8000.

Figure 13 .
Figure 13.Numerical simulation: a typical pattern of the double helix formation, observed in our simulations: a fiber quickly straightens in the shear plane, and then forms a compact structure, which slowly becomes elongated and oriented close to the vorticity direction.As an example, we present here the snapshots from the evolution of a fiber with N = 60 and A = 0.5, initially U-shaped in a plane that is inclined at an angle of 20 • to the xy-plane.Times of snapshots from the left side: t = 0, 62, 123, 204, 565, 1006, 1645.See movies 8-10.

Figure 14 .
Figure 14.Experiment: a typical evolution of fiber shapes in a dilute suspension under shear flow, leading to the creation of double helices.Cropped snapshots of several flexible fibers at normalized time t = 0, 105, 314, and 942.In a relatively short time, fibers straighten in the shear plane.Later, fibers tend to form compact structures.At large times, the structures become elongated, with their orientation correlated positively with the vorticity direction.See movie 1.

Figure 15 .
Figure 15.Typical process of a double-helix formation from a compact irregular shape that slowly becomes elongated close to the vorticity direction.Normalized times of the snapshots are indicated.(a) Numerical simulation of a fiber with N = 40, A = 0.6, Θ0 = 70 • and Φ0 = 90 • .(b) Experiment shown in movie 1.

Table 2 .
Number of simulations for N = 60 and A = 0.5: total and leading to the double helix or zigzag shapes, starting from different initial shapes of the fibers (C-shaped, S-shaped, U-shaped).C-shaped elastic equilibrium is denoted by ce.Otherwise, the fibers are straight at the elastic equilibrium.