Dynamics and Topology of Flexible Chains: Knots in Steady Shear Flows

We use numerical simulations of a bead-spring model chain to investigate the evolution of the conformation of long and flexible elastic fibers in a steady shear flow. In particular, for rather open initial configurations, and by varying a dimensionless elastic parameter, we identify two distinct conformational modes with different final size, shape, and orientation. Through further analysis we identify slipknots in the chain. Finally, we provide examples of initial configurations of an"open"trefoil knot that the flow unknots and then knots again, sometimes repeating several times. These changes in topology should be reflected in changes in bulk rheological and/or transport properties.


Introduction
One area of complex fluids concerns the motion and topology of elastic filaments in fluid flows, which is inspired by both natural and industrial phenomena, such as flagellar motion [1,2,3], fluid-structure interactions, e.g. [4,5], and polymer processing, e.g. [6,7,8,9,10]. Such long flexible filaments can have non-trivial dynamics, e.g. drift across streamlines can occur, e.g. [7,11,12], and topologies. For example, knots are found in Brownian systems such as bacterial DNA, protein structures, and polymer chains, and non-Brownian systems such as ordinary string, elastic fibers and chains of linked beads [13,14]. In this paper we report the time-dependent shapes of flexible non-Brownian filaments in a steady shear flow and identify conditions that allow an "unknotting-knotting" transition, where an unknotted filament is later observed to form an (open) knot as it rearranges continually in the flow.
For rigid non-Brownian particles in shear flow, when the Reynolds number is low, a straight elongated particle has a Jeffery orbit, which is the basis for many studies, and helps organize those initial states that have orientation and rolling about the vorticity axis from those that primarily align with and tumble periodically about the flow direction. For flexible or deformable objects, characterizing the shape of the particle in flow, and the corresponding dynamics, is more challenging. For example, experimental measurements of fiber dynamics commonly use DNA as a model polymer [15], and, on a larger length scale where Brownian effects are less significant, linked colloidal chains are also used as experimental models of a worm-like chain [16,17] and the cellular flagellum [18]. In all of these cases, dynamics in a single plane are measured, but it is difficult experimentally to simultaneously track dynamics in both the shear and vorticity directions [19].
On the other hand, numerical simulations make possible the visualization of multiple viewing planes and the detailed analyses of shape and topology of flexible filaments. Previous work, some of which includes Brownian effects and some of which does not, has identified and tracked the tumbling of flexible filaments (e.g. polymers) in detail, relating the tumbling period and shape to polymer length, flexibility and shear rate, as well as rates of rotation both in and out of plane [12,20,21,22,23]. However, due to computational constraints, these simulations of fibers are generally bead-chains with about N = 50 beads, with the longest ones containing N = 100 beads; an alternative approach to some of these questions is to use slender-body theory, e.g. [24,25].
In recent years, knots in long chains has been studied more widely, including the observation of thermally driven knotting and unknotting of a polymer chain [26], the influence of chain tension on the knot, e.g. [27,28], and the influence of electric fields on knots [29]. Also, the untying of a knot in an extensional flow has been studied recently [30], and has the spirit of the present paper where we ask about the effects of flow. In contrast, for non-Brownian systems other kinds of forcing, such as vibration, can produce knotting of chains of linked spheres [13,14]. We are not aware of any similar studies, experimental or numerical simulations, highlighting conditions that lead to long non-Brownian chains forming knots in flows.
Thus, while significant steps have been taken towards understanding the detailed dynamics of flexible non-Brownian fibers, to understand topology is much more difficult: it requires three-dimensional detail and very flexible, long and thin fibers as, in theory (when there is no flow), the probability of knotting in a random walk increases exponentially with the length of the walk [31]. Moreover, a smart choice of initial conditions is needed and an operational definition is needed for a knot in an open fiber.
In this paper, we focus on the dynamics and topology of long, flexible, non-Brownian chains of beads in a steady shear flow and use numerical simulations to investigate whether knots in the fiber occur as a result of the flow. We will find that unknotted fibers are capable of forming (open) knots, and in some cases we document a sequence of unknotting-knotting transitions. The numerical solution of this kind of problem requires tracking the motion of N beads, where hydrodynamic interactions between the beads mean that the dynamics of the topology of the object involve 3N coupled nonlinear ordinary differential equations. Not surprisingly, such dynamics should be expected to be chaotic though we have not attempted in this paper to link our study of complex shapes and the unknotting-knotting transition to the underlying chaotic dynamics. Our results highlight new aspects of the dynamics of flexible filaments in flow.

A fiber as a string of spherical beads
We consider a long non-Brownian flexible fiber in steady undisturbed shear flow v 0 =γze x of a fluid with viscosity η (figure 1). The vorticity vector is in the ydirection. The fiber deforms owing to elastic and bending forces and is modeled as a chain of N identical spherical beads with diameter d (figure 1a) [12,32]. The parameters of these forces are k, the spring constant normalized by πηdγ, A, the bending stiffness normalized by πηd 4γ , and ℓ 0 , the equilibrium distance between the closest bead centers measured in units of d. All distances are measured in units of d, and time in units of 1/γ. Thus, there are four dimensionless parameters that characterize the chain: k, A, N and ℓ 0 . We work close to the inextensible limit with |ℓ 0 − 1| ≪ 1 and k ≫ 1 and focus on the dynamics for a wide range of values of A.
We assume that the particle-scale Reynolds number is much smaller than unity, i.e. Re = ργd 2 η ≪ 1, where ρ is the fluid density and the fluid flow satisfies the quasi-steady Stokes equations, with no-slip boundary conditions at the bead surfaces. The equations are solved using the multipole expansion method, which accounts for lubrication corrections to speed up the convergence of the friction matrix as the order of the multipole truncation is increased. This solution yields the velocities of the beads, as explained in [12]. The algorithm and its numerical implementation hydromultipole are described in [33,34]. Also, the application of these tools to model dynamics of flexible fibers is outlined in [12]. The time-dependent translational velocities of the bead centers, evaluated by the hydromultipole numerical code, are used by an adaptive fourthorder Runge-Kutta procedure, from which the positions of every bead are updated.
We focus on the dynamics of long fibers (N ≈ 100 −150), which, in flow, tend to fluctuate from elongated shapes to compact shapes with some of the beads coming close together. Spurious overlaps of beads can occur, which can significantly influence the time scales of the dynamics of close particles [35]. Therefore, we analyzed different procedures that allow the dynamics to be continued once an overlap occurs and chose a procedure that caused the least number of overlaps. If the distance between two bead centers r < 1, then the lubrication correction for this pair is evaluated for the rescaled configuration 1 + ǫ lub , where ǫ lub = 10 −4 is the integration accuracy per step, and the dynamics of the original positions are continued. The details of this original simulation shape, which is one of two typical initial configurations studied in this work. The expanded view labels beads i and i − 1. Also, r i is the position vector of each bead center and t i is the vector connecting two adjacent bead centers. (b) A configuration with shading to indicate the chain's ellipsoid of revolution. The shear flow is shown in the background. The major axis of the chain's configuration is shown with an arrow, while the red and blue coloring of beads denote two halves of the fiber for visual clarity. The orientation is measured by the angles θ and φ, such that θ measures the inclination from the vorticity axis (y) and φ measures the orientation in the xz-plane, measured from the z-axis. approach for handling close particle pairs are provided in the next section.
In the simulations reported here, the initial geometry is chosen first from a random rotation of a 'candy-cane' shape with N = 152, i.e. a straight chain that bends into a semi-circle at one end; see figure 1 and Appendix A. The candy-cane shape is convenient for starting with an elongated initial configuration that breaks the symmetries associated with shear flow. Also, we have performed calculations with a trefoil as an initial condition with N = 99 (see Section 4 and Appendix B). For each simulation we fix ℓ 0 = 1.02 and k = 500, and vary A systematically. We evolve the shape in a steady shear flow using the hydromultipole algorithm [33].

A simulation procedure for compact structures
As introduced in the previous section, and will be clear when the many simulations in this paper are shown, long flexible fibers can form compact structures with some of the beads coming very close to other beads. As a result, it sometimes happens that surfaces of pairs of beads spuriously overlap. This complication will be faced by all numerical methods and for the long-time dynamics of long flexible fibers, which are the focus of this paper, such artifacts cannot be avoided. The physical problem is that the friction and mobility coefficients are not defined for the distance between the bead centers less than their diameter (equal to one in our units). In the literature, there exist different methods that allow the dynamics to be continued once an overlap takes place (e.g. see the brief review in [35]). In this sub-section we present a new approach for the problem of overlaps; the reader interested only in the results of the fiber dynamics, including the unknotting-knotting transition that we have identified, can skip this sub-section.
Not surprisingly, it is known that numerical procedures for handling the interaction of close pairs can sometimes significantly influence time scales of close-particle dynamics. In particular, examples are known of systems of particles where spurious overlaps increase the period of oscillations, even by a factor of two [35]. Therefore, it is important to choose a numerical treatment of overlaps that minimizes the time-dependent number of overlapping pairs and also prevents the particle surfaces from remaining overlapped for a long time.
We have observed that spurious overlaps of beads sometimes can significantly modify compact shapes of flexible fibers, lead to their orientation along the y (vorticity) axis, and change properties of the dynamics. In particular, a larger number of overlaps leads to a larger fiber curvature and a smaller radius of gyration. Therefore, to address the main questions of this paper we found it essential to use an optimal procedure once an overlap takes place. After analyzing and comparing a few such methods, we proposed a new one. Here we provide a brief description of this procedure (called BWE), which we have applied to continue evolution of the system of equations in case of spurious overlaps of beads. Our new approach gives a small number of overlapping pairs of the beads, which, in addition, usually do not remain overlapped for a long time.
The procedure is based on modification of the lubrication correction δ (2) (r ij ) in the case when the beads i and j overlap. To explain the idea, we first briefly explain the general idea how the lubrication correction is constructed [33,34]. Truncating the multipole expansion at the order L, we evaluate the multipole approximation ζ (N ) L for the friction matrix of the N-bead system. This step can be accomplished even if the distance r ij between the centers of beads i and j is close to but smaller than one (here, i, j = 1, ..., N). We need to add a lubrication correction δ (N ) L to speed up the convergence of the multipole expansion,ζ L . This correction is the sum of all pairwise contributions, k = l, where r kl is the relative position vector of particles k and l. Here, ζ (2) (r kl ) and ζ L (r kl ) are, respectively, the two-body friction matrices -the exact one and its multipole approximation of the order L. Now, consider a pair of overlapping beads i and j, i.e. r ij ≤ 1, for which elements of the matrix ζ (2) are not defined -they diverge at the contact. We overcome this problem by rescaling the lubrication correction for the overlapping pair of beads, and evaluate the dynamics from the unchanged positions of all of the particles. The multipole approximation ζ for the friction matrix of the N-bead system and the lubrication corrections for non-overlapping beads are not modified.
The key idea of the method is to take ǫ lub equal to the integration accuracy per step, which in this work is equal to 10 −4 . Owing to this choice, the time of a spurious overlap is reduced in comparison to smaller values of ǫ lub . Indeed, a typical overlap depth is of the order of the accuracy. Therefore, by assigning to the relative velocity of the beads the new, quite large value, which corresponds to the gap size between their surfaces, i.e. ǫ lub equal to a typical overlap depth, we obtain the required result that the overlapped beads can be disconnected easily, i.e. in relative terms they quickly separate their surfaces from each other.

Results: Elongated initial configurations
Upon starting the shear flow, the fiber, which was initially in the candy-cane configuration, described in Appendix A, aligns approximately with the flow direction and undergoes a 'tumbling' or 'yo-yo' motion (figure 2a,d; the chain halves are colored for clarity), similar to observations in studies of flexible polymers in shear flow [19,36,37,38]. We identify two distinct behaviors for the dynamics, depending on the value of A. To quantify the dynamics, we report the time-dependent approximate shape and orientation of the chain.

Typical simulations and geometric characterization
First, for a given configuration, we calculate the moment of inertia tensor [39] with the Cartesian coordinates n, m = 1, 2, 3, and the labels i, j of the vertices of each of the N beads. The lengths a, b, and c denote the three semi-axes of the ellipsoid of inertia associated with a configuration (figure 1b). These lengths are related to the three eigenvalues λ n of T nm , i.e. a = √ 3λ 1 , b = √ 3λ 2 , and c = √ 3λ 3 . The largest eigenvalue a gives the primary axis of the ellipsoid, and the orientation angles of the associated eigenvector, θ and φ, are measured from the vorticity and z axes, respectively, as shown in figure 1b. We use a, b, and c to visualize the fiber shape: if a ≈ b ≈ c, the configuration is roughly spherical, whereas if a ≫ b ≈ c, the shape is prolate ellipsoidal, and if a ≈ b ≫ c, the shape is oblate spheroidal.
Also, we calculate the radius of gyration R g of each configuration, i.e.
, by comparing the coordinates of each bead r i to the center of mass r cm . Since R 2 g = trace T nm then R 2 g = λ 1 +λ 2 +λ 3 [40]. With these measures {R g , a, b, c, θ, φ}, we track a given chain's size, shape, and orientation as it evolves in the flow.
For a given initial orientation, specified in Appendix A, chains with lower values of A, e.g. A < 2, behave differently than chains with higher values, A ≥ 3; qualitatively, the former form tight coils aligned with and rolling in the vorticity direction, while the latter forms more open tumbling configurations, by analogy to the results of Ref. [41]. We present detailed results from each category, starting from the same initial 'candycane' geometry.
First, we focus on a simulation of a more flexible chain where A = 1; the evolution is shown in figure 2a and movie 2a and the dynamics are quantified in figure 2b. Initially, the chain extends and orients in the flow direction, resulting in a large radius of gyration and highly skewed axes lengths (a ≫ b, c); the orientation angles of the main axis are θ ≈ π/2 and φ ≈ π/2, i.e. the flow direction. As time continues, the two ends fold towards the center, R g decreases significantly, while a decreases and becomes comparable to b and c. At still later times, the configuration remains relatively constant in the form of a prolate ellipsoid (a ≫ b, c) parallel to the vorticity axis; also, φ systematically increases, which is consistent with a 'rolling' motion.
We next consider a simulation of a stiffer chain, where A = 3; see figure 2c,d and movie 2d. The time dynamics is different than the case of A = 1, e.g. the collapse of the initially extended chain leads to a less compact shape that does not resemble a prolate spheroid. For most times, the orientation angles of the main axis remain close to the flow direction.
Finally, we examine the number of close pairs for both A = 1 and A = 3. For each time-step, we check the distance between each pair of beads; we consider beads a close pair when d ij , the distance between beads i and j, is less than some threshold, and the beads are not the closest neighbors, i.e. i = j ±1. Here we set the threshold d ij /d = 1.05, slightly larger than the dimensionless equilibrium inter-bead spacing ℓ 0 = 1.02. We notice that the number of close pairs in the more flexible fiber (A = 1) is around 4 times greater than the number in the less flexible fiber (figure 2b, 2c). These observations are consistent with A = 1 yielding a tightly packed coil and A = 3 resulting in a less compact shape.
Both visually, as seen in figure 2a,d and movies 2a,d, and quantitatively, through figure 2b,c, a change in A can result in considerably different dynamics. A closer inspection of the configurations shows that the chains may become self-entangled. For example, in panel (V) of figure 2d, the blue end has passed through a loop formed at the opposite red end. Visual inspection shows that this is a slipknot (see section 3.3), but such shapes raise the question of whether or not simple flows can cause chains to tie themselves into knots, and how to identify and classify such configurations.

Knot classification
For a given configuration at a given time, we then study the topology of a chain. A knot is a topological state of a closed curve, where no set of simple line deformations, e.g. 'pushing' a strand, without passing one strand through another, can change its state to that of an open circle, or an 'unknot'. To classify knots, polynomials known as knot invariants, which uniquely identify each curve to a knot type, can be calculated [42]. Here, we calculate the Alexander polynomial of each configuration (see Appendix C).
Any open curve, i.e. a curve with free ends, can always be 'untied' by line deformations and cannot be considered knotted in a rigorous sense. However, these 'open knots' still have characteristics that can be classified, and affect the properties of a chain, such as its shape, size, and relaxation modes [39]; the classification of open knots remains a subject of study. In general, an open curve must have its ends joined together to form a closed curve, from which the topological state can be calculated. There are several well-studied closure schemes [43], but as these methods can result in different topological states for the same open configuration, there is inherently ambiguity.
We use the stochastic closure scheme [44]. For each open configuration, we draw a sphere centered at the center of mass of the chain, with radius ten times the radius of gyration. A point is randomly chosen on this sphere [45], and we connect by a straight line the two ends of the open curve to this point to create a closed curve that can be analyzed for its Alexander polynomial. Repeating this procedure with different random points on the sphere (here, 1000 randomly chosen points), we build a spectrum of knots; i.e. each configuration of the flexible chain can be described as a superposition state of various knot types.
This approach effectively provides a statistical measure of the entanglement of each chain, which we report with a knotted fraction k f , i.e. the fraction of the stochastic closures which result in a knot of any type. For an open flexible fiber, "a knot" corresponds to k f close to 1.
As might be expected from visual analysis of the configurations shown in figure 2, for the initial candy-cane configuration, no knots have been observed for the time of the simulation even though the fiber shapes were often very compact at long times. The calculated values of k f were always close to zero for all time steps in the simulation.
Because the chains we study are highly flexible, it is also necessary to understand the topology of sections of the chain. Thus, for each subchain, the knotted fraction k f is calculated using the stochastic closure scheme and the resulting values are combined into a knot matrix K(i, j), where i, j = 1 · · · N. Each point (i, j) of the matrix corresponds to starting and ending beads of the subchain, i and j, respectively. For example, the origin, at (1, 152), refers to a subchain starting at i = 1 and ending at j = 152, i.e. the entire chain. Another point on the matrix, such as (70, 152), refers to a subchain starting at i = 70 and ending at j = 152, which is the latter half of the full chain. Therefore, every cell of the matrix K(i, j) corresponds to a subchain of the original configuration.
The knot matrix is symmetric, K(i, j) = K(j, i), because the subchain that starts from i and ends at j is just the same subchain that starts at j and ends at i. Therefore, it is sufficient to show values of the knot matrix only for i ≤ j.

Slipknots
In figure 2d(V), we observe a configuration that has the blue end of the chain threaded through a loop formed from the red half of the chain. As, intuitively, a knot requires the pulling of an 'end' through a 'loop', the way we might tie a knot in a shoelace, these types of configurations bear further scrutiny. By closing the entire chain with the stochastic closure scheme, these configurations are shown to be unknotted -the fraction of nontrivial closures is low, k f close to 0.
With this visual motivation, next we studied the results of simulations at different values of A, keeping an eye out for possible knots. For A = 1.4, the knot matrix K(i, j) at an intermediate time when an entangled configuration occurs is plotted in figure 3. The white cells have zero knotted fraction, while the darker cells have a higher knotted fraction. We note that the knot matrix is not symmetric with respect to reflection in the line y = 152 − x, because, in general, (i, j) and (152 − j + 1, 152 − i + 1) In plots such as figure 3, dark regions can be identified as knot cores [46]. In these regions, the corresponding subchains have k f approaching unity and are considered knotted, even though k f is close to 0 for the chain as a whole. For example, as evident in figure 3, there is a darker block corresponding to a subchain that starts at the bead 70 and ends at the bead 152. This knotted subchain is shown highlighted in red in the accompanying image in figure 3. In the right panel of figure 3, the entire chain is shown, with the knotted subchain highlighted in red, and the other part highlighted in blue. Thus, we have identified, qualitatively and quantitatively, the formation of a slipknot in a flexible chain subjected to a steady simple shear flow.

Results: Trefoils in steady shear flow
Now that we have explored one extended type of initial condition, we explore a more "knot-like" initial condition. Thus, we ask the question: can shear flow untie a trefoil? We set an initial condition in a trefoil (3 1 knot) shape (with one bead missing), see Appendix B, and vary the initial orientation of this shape relative to the shear flow as shown in Fig. 4. Then, we evolve the chain for different relative bending stiffness A. By using the stochastic closure scheme and counting invariants, we can calculate and plot the fraction knotted k f at each time step, i.e. how many of the stochastic closures result in non-trivial knots. A knot corresponds to k f close to 1.
We used five different orientations of this shape, shown in figure 4a-e. A trefoil with N = 99 is shown in figure 4a and denoted as I 0 = a. Also, the orientations I 0 = b and I 0 = c (see figure 4b and c) were obtained by rotating the configuration I 0 = a along the x-axis by π/2 and π, respectively; the orientations I 0 = d and I 0 = e (see figure 4d and e) were obtained by rotating the configuration along the z-axis by π/2 and 3π/2, respectively.

Dynamics: An unknotting-knotting transition
We performed 29 simulations starting from a trefoil shape at different orientations I 0 = a, b, c, d, e shown in figure 4, and several values of the bending stiffness A. For certain I 0 and A, we observed an unknotting-knotting transition, by which we mean Table 1: The existence of an unknotting-knotting transition as characterized by a trefoil unknotting then knotting (1) or not (0), as a function of the stiffness A and initial orientation I 0 . The table summarizes the 29 different simulations that were performed.  . Here k f denotes the knotted fraction based on the stochastic closure scheme described in the paper.
that the trefoil unties and then ties again. In Table 1, we specify values of the bending stiffness A and the initial orientation I 0 for which such a transition was (or was not) observed. In most cases, after a transient period, the trefoil unties but, in some cases, at later times the chain knots again, as characterized by the knotted fraction, k f , decreasing to very small values and then rising up to almost one. We have found 13 examples of the unknotting-knotting transition: for fibers with A = 1.4, initially oriented as I 0 = b, and for fibers with A = 2.1, 2.2, 2.3, initially . Here k f denotes the knotted fraction based on the stochastic closure scheme described in the paper. oriented as I 0 = c. In most cases, the unknotting-knotting transition is observed more than once during evolution of the same fiber. In figures 5-8 we plot the time-dependent knotted fraction k f based on the stochastic closure scheme, and display some unknotted and some knotted configurations. Also, in the figures we use different colors to help highlight the complex topology of the complicated chain.
Two examples of the unknotting-knotting transition with A = 2.2, visible in the time evolution of k f in figure 5, are shown in movie 5. In particular, in figure 5 it is also illustrated how the chain first unknots as k f decreases below 0.1, and then knots again, as evidenced by k f then increasing to approach unity. Eventually, however, the chain is visibly unknotted and k f remains near zero. The other examples of the unknottingknotting transitions that we have uncovered are shown in movies 6-8 and the plots of the time-dependent k f in figures 6-8. Specifically, for A = 1.4, as displayed in figure 6 and movie 6, a knotted fiber is unknotted after time t ≈ 330 and knots again at time t ≈ 380. At the end of the simulation, the chain remains in the knotted configuration, with k f close to one. The evolution of a fiber with A = 2.1 is shown in figure 7 and movie 7. Again, the initially knotted chain unknots and knots again several times. Finally, the chain is knotted and stays in this form to the end of the simulation.
In yet another example, the dynamics of a fiber with A = 2.3 is presented in the figure 8

Characterizing configurations that display the unknotting-knotting transition
Comparing the different simulations with each other, e.g. figures 5-8, we conclude that the existence of the unknotting-knotting transition is not correlated with a final knotted or unknotted state. The question then can be asked if the existence of the unknottingknotting transition for a trefoil is correlated with some specific mode of the flexible fiber dynamics, similar to those observed in the simulations that started from the initial candy-cane configuration, which were shown in figure 2. For example, the characteristic features of the two main conformational modes where the fiber primary axis tends to be aligned along the vorticity axes (mode I) or the flow direction (mode II), discussed earlier in connection to figure 2, are summarized in table 2.
close pairs many a few Therefore, in figure 9, we compare the dynamics of trefoils with different bending stiffnesses A and initial orientations I 0 in terms of their characteristic modes of conformation. We present evolution of the same parameters of the fiber shapes as in figure 2, but this time for fibers that have an initial trefoil shape. For two sets of the parameters, corresponding to the dynamics shown in columns (ii) and (iii), the unknotting-knotting transition is observed. For two other sets of the parameters, related to columns (i) and (iv), the unknotting-knotting transition is not observed.
We first compare two cases for which the unknotting-knotting transition is observed, i.e. columns (ii) and (iii) in figure 9. For example, the evolution displayed in column (iii) for larger times approximately corresponds to mode I, which is characterized by an elongated shape with relatively large aspect ratios, and the tendency to align with the vorticity direction and roll around it. For comparison, at long times, the dynamics visible in column (ii) approximately corresponds to mode II, which is characterized by a smaller elongation, and no special features of the dynamics related to the vorticity direction.
On the other hand, comparing columns (i) and (ii) in figure 9, we conclude that the same mode II configuration is reached by fibers for which the unknotting-knotting transition, respectively, did not and did take place. These examples of two topologically different behaviors correspond to different initial orientations of the trefoil characterized by a single value of A. The complementary comparison of columns (iii) and (iv) in figure  9 demonstrates that fibers with and without the unknotting-knotting transition can both later exhibit the same mode I configuration. In this case, examples of two topologically different behaviors correspond to the same initial orientation of the trefoil but different values of A. Therefore, the examples displayed in figure 9 illustrate that the existence of the unknotting-knotting transition is not correlated with the final configuration mode I or mode II of the dynamics.
Finally, we investigate if the process of creating a knot is correlated with any simultaneous specific features of the fiber shape. Two examples shown in figure 10 contribute to a negative answer. In both cases, the fiber starts from very small knotted fraction k f , which increases up to almost one. However, in the first case (shown in the left part of figure 10), the fiber is in mode II, and in the second case (shown in the right part of figure 10), the fiber is in mode I. To conclude this discussion, we note that the process of knot creation takes place both for fiber shapes that are compact (figures 6, 7) and loose (figures 5, 8).

Summary and conclusions
We have focused on the question of whether a simple shear flow can tie a knot in a flexible fiber. If yes, then can we predict the future existence of a knot based on the fiber flexibility and its initial shape and orientation? Also, is the knot a transient effect or a stable configuration? To address these questions, the fiber was modeled as a chain of N beads connected by very stiff springs with additional elastic forces acting against the fiber bending. The fiber evolution was determined using the numerical code hydromultipole, based on the accurate multipole algorithm for solving the Stokes equations with no-slip boundary conditions on the bead surfaces. In addition, to characterize the topology a knot on a flexible open fiber was defined as in Refs. [33][34] and associated with the knotted fraction k f close to unity.
In the first part of our study, we investigated the evolution of fibers of different flexibilities, but with the same initial orientation and the same untied shape (with a vanishing knotted fraction). We concluded that for such a simple initial configuration, a change of flexibility is not sufficient to cause formation of a knot, although evidence of slipknots was observed.
Then, in the second part of our study, we constructed a more complicated initial shape: a trefoil with the knotted fraction equal to one. We used the same initial shape at different initial orientations and for fibers of different flexibilities, and we followed their evolution in a shear flow. A few of the more flexible fibers remained knotted, but most of the fibers untied after some time. In particular, we then observed that some of these untied fibers again became knotted, as illustrated in figure 11 for A = 2.3 and N = 99. In this figure we display two shapes with the knotted fraction close to zero, which later on lead to evident knots (with the knotted fraction close to one).
Therefore, indeed, a simple shear flow can tie a knot on a flexible fiber. Sometimes, such an unknotting-knotting transition repeated several times during the evolution of the fiber dynamics. We observed 13 such unknotting-knotting transitions in 4 (out of 29) simulations, for 4 different (moderate) flexibilities and 2 different initial orientations. Based on all of our observations, we conjecture that the formation of knots is a transient effect and may possibly be related to the chaotic properties of the dynamics of this N bead system.
To conclude, we have analyzed basic features of the dynamics of long flexible fibers in unbounded steady shear flow. In particular, we give examples of flexible fibers, initially shaped as open trefoils, which are later untied and then knotted again, sometimes several times, by the flow. This finding provides a new perspective for future studies of topological structures of long flexible non-Brownian chains evolving as they are entrained by ambient flows, and possibly changing the macroscopic transport properties. Finally, as we noted in the introduction, the numerical solution of this problem requires the solution of 3N coupled nonlinear ordinary differential equations for the positions of the N beads. We should expect such a dynamical system to have all of the attributes of chaos. How the observed dynamics of knot formation and unknotting may (or may

Appendix C. Knot Invariants
After our open curve is closed by the stochastic closure scheme, it can be analyzed for its polynomial knot invariant. There are by now many different knot invariants, which have varying ability to distinguish knots; the first one, which was discovered in 1928 by Alexander, is one of the simplest to compute. It has some drawbacks, such as the inability to distinguish handedness (i.e. it cannot distinguish between two mirror images), but as we only need to distinguish between an 'unknot' and a 'knotted' curve to calculate the knot fraction, the Alexander polynomial is sufficient. The invariant is calculated by projecting a curve into 2D, then identifying 'crossing points' and 'regions', while keeping track of the 'overlying' and 'underlying' branches at each crossing. If there are n crossings, there will be n + 2 regions in the knot diagram. Depending on the orientation of the overlying/underlying branches at each crossing, the regions adjacent to the crossing are assigned different indices. A n × (n + 2) incidence matrix is constructed from the indices at each crossing and region, and two columns corresponding to adjacent regions are removed, to form a n×n matrix. The determinant of the matrix is then the Alexander polynomial. More details, and proof that removing two arbitrary columns (provided they refer to adjacent regions) results in the same normalized Alexander polynomial, can be found in his 1928 paper [47].
We have written a custom code to compute the Alexander polynomial of any configuration, given the (x, y, z) coordinates of each bead in the chain, following the algorithm described by Alexander, with the help of functions from other Alexander polynomial calculators [48].