Self-propelled Worm-like Filaments: Spontaneous Spiral Formation, Structure, and Dynamics

Worm-like filaments that are propelled homogeneously along their tangent vector are studied by Brownian dynamics simulations. Systems in two dimensions are investigated, corresponding to filaments adsorbed to interfaces or surfaces. A large parameter space covering weak and strong propulsion, as well as flexible and stiff filaments is explored. For strongly propelled and flexible filaments, the free-swimming filaments spontaneously form stable spirals. The propulsion force has a strong impact on dynamic properties, such as the rotational and translational mean square displacement and the rate of conformational sampling. In particular, when the active self-propulsion dominates thermal diffusion, but is too weak for spiral formation, the rotational diffusion coefficient has an activity-induced contribution given by $v_c/\xi_P$, where $v_c$ is the contour velocity and $\xi_P$ the persistence length. In contrast, structural properties are hardly affected by the activity of the system, as long as no spirals form. The model mimics common features of biological systems, such as microtubules and actin filaments on motility assays or slender bacteria, and artificially designed microswimmers.


I. INTRODUCTION
Its importance in biology and its enormous potential impact in technical applications makes active soft matter a field of rapidly growing interest and progress. [1][2][3] Flexible slender bodies are of particular importance. The majority of natural swimmers propel themselves using flexible, hair-like structures like cilia and flagella. 1 Another important example are actin filaments and microtubules, major constituents of the cytoskeleton, whose capability to buckle decisively controls the mechanical properties of the cell body. 4 Flexibility is the crucial ingredient for the formation of small-scale spirals 5 and possibly also for large-scale swirls 6 of microtubules on motility assays. Even the structure of slender bacteria can be dominated by their flexibility. 7 Elextrohydrodynamic convection can propel colloid chains because they are flexible, 8 just as the swimming mechanism of assembled magnetic beads in an oscillating external magnetic field is possible because of the swimmer's flexibility. 9 Flexibility is of course also the feature that allows for the instabilities leading to cilia-like beating in artificially bundled microtubules. 10,11 Despite its importance, the number of theoretical studies of active agents that incorporate flexibility is still relatively small, and can roughly be subdivided into works that focus on buckling phenomena and on freeswimming agents. Symmetry breaking instabilities leading to rotation and beating motion of active filaments on motility assays can be described with a phenomenological ordinary differential equation for the filaments. 12 The propulsion force of motor proteins has been predicted based on a Langevin model for buckled, ro-tating actin filaments and microtubules. 13 Numerical studies with Lattice-Boltzmann simulations and Brownian or multi-particle collision dynamics have demonstrated that clamped of pinned filaments composed of stresslets or propelled beads can show cilia-like beating or rotation. 14,15 The behaviour of free-swimming actin filaments on motility assays was reproduced in early numerical studies using the Langevin equation. 16 However, it was only recently that theoretical study of flexible, active filaments that can move freely has received significant attention. Lattice-Boltzmann simulations reveal that spontaneous symmetry breaking in chains of stresslets can lead to rotational or translational filament motion. 17 Brownian dynamics simulations of short self-propelled filaments suggest that different types of motion occur for single filaments 18 and that spontaneous rotational motion can arise for pairs of filaments. 19 A combination of Brownian dynamics simulations and analytic theory shows that shot noise in worm-like filaments leads to temporal superdiffusive filament movement and faster-decaying tangent-tangent correlation functions. 20 Finally, chains of active colloids connected by springs have the same Flory exponent but a different prefactor of the scaling law compared to chains of passive colloids, as shown recently both analytically for beads without volume exclusion and numerically with Brownian dynamics simulations for beads with volume exclusion. 21 The free-swimming behaviour of a worm-like filament that is tangentially propelled with a homogeneous force is still unexplored and is the subject of this work. The model is introduced in Section II. Results for the structural and dynamic properties over a wider range or propulsion forces and filament flexibilities are presented in Section III. We find that the filament can spontaneously form spirals, which is the mechanism that dominates the behaviour for large propulsion forces. The rele-vance of our observations for natural and artificial active agents is discussed in Section IV. We present our conclusions in Section V.

II. MODEL AND METHODS
We study a single, active, worm-like filament, which is modelled as a sequence of N + 1 beads connected via stiff springs. The overdamped equation of motion is given by where r i are the coordinates of bead i, γ is the friction coefficient, U is the configurational energy, F (i) k B T is the thermal noise force, and F (i) p is the active force that drives the system out of equilibrium. The configurational potential energy is composed of a bond contribution between neighbouring beads a bending energy and an excluded volume term modelled with repulsive Lennard-Jones interactions where r i,j = r i − r j is the vector between the position of the beads i and j, k S is the spring constant for the bond potential, r 0 is the equilibrium bond length, κ is the bending rigidity, and and σ are the characteristic volume-exclusion energy and effective filament diameter (bead size). The drag force γṙ i is the velocity of each bead times the friction coefficient γ. The thermal force F (i) k B T is modelled as white noise with zero mean and variance 2k B T γ/∆t as described in Ref. 22. Note that hydrodynamic interactions (HI) are not included in our model. The model is thus in particular valid for (i) neutral swimmers, for which HI are known to be of minor importance, [23][24][25] (ii) swimmers near a wall, where HI is of less importance, 26,27 and (iii) microorganisms that glide on a surface, such as nematodes like C. elegans. 28,29 FIG. 1. Filament model: Beads are connected via stiff springs. The active force acts tangentially along all bonds. Colour gradient indicates the force direction.
Without propulsion force, F (i) p = 0, the model matches the well-known worm-like chain model for semi-flexible polymers. 30,31 For active filaments, we use a force per unit length f p that acts tangentially along all bonds, i.e., as illustrated in Fig. 1. The force along each bond is distributed equally onto both adjacent beads. We consider systems with parameters chosen such that (i) k S is sufficiently large that the bond length is approximately constant r 0 , that (ii) the local filament curvature is low such that the bead discretization does not violate the worm-like polymer description, and that (iii) the thickness of the chain has negligible impact on the results. When these requirements are met, the system is fully characterized by two dimensionless numbers, where L = N r 0 and ξ P are the length and persistence length of the chain, respectively. ξ P /L is a measure for the bending rigidity of the filament. The Péclet number P e is the ratio of convective to diffusive transport and measures the degree of activity. For its definition, we use that the filament has a contour velocity v c = f p /γ l , and that the translational diffusion coefficient D t = k B T /γ l L, where we have introduced the friction per unit length γ l = γ(N + 1)/L. The ratio of these numbers which we call the flexure number, provides a ratio of activity to bending rigidity. Previous studies showed that this number is decisive for buckling instabilities of active filaments. 12,15 It will be shown below that this is also a determining quantity for spiral stability and rotational diffusion.
Simulations were performed in two dimensions, where volume exclusion interactions have major importance. Equations of motions were integrated using an Euler scheme. Simulation parameters and results are reported In our simulations we used k S = 4000 k B T /r 2 0 , r 0 = σ = L/N , and = k B T if not stated otherwise. A large parameter space for P e and ξ P /L was explored by varying f p , N , and κ. N was varied in the range from 25 to 200 from the highest to the lowest ξ P /L. Almost all simulations were run for more than 5 τ . An initial period of the simulation output is discarded in the analysis. The timestep ∆t was adjusted to the remaining settings to ensure stable simulations. Unless explicitly mentioned, results refer to simulations that were started with a perfectly straight conformation.
All simulations were performed using the LAMMPS molecular simulation package 32 with in-house modifications to describe the angle potential, the propulsion forces, and to solve the overdamped equations of motion.

III. RESULTS
The characteristic filament behaviour depends on its bending rigidity and activity and can be divided into three regimes (see Fig. 2). At low P e or high ξ P /L, the "polymer regime", the active filament structurally resembles the passive filament with P e = 0. The main difference compared to the passive filament is that the active force drives the filament along its contour, leading to a directed translational motion -we name this characteristic movement "railway motion". At high P e and low ξ P /L, the filament spontaneously winds up to a spiral. The "spiral state" is characterized by ballistic rotation but only diffusive translation. Spirals can spontaneously break up. Their lifetime determines whether spiral formation has a minor impact on the overall filament behaviour, the "weak spiral regime" at intermediate P e, or whether spirals are dominating, the "strong spiral regime" at large P e. Because spiral formation has a major impact on both the structure and the dynamics, features related to spiral formation are addressed first. Black lines are a guide to the eye. Green area: threshold for spiral stability against break-up by widening. Spirals above this threshold will unfold by widening, spirals below will not. Purple: parameter space that can be obtained by actin filaments on a myosin carpet at T = 300 K using parameters for fp and κ from Ref. 12. Structural and dynamic properties of the elongated and spiral state are presented afterwards.

III.1. Spiral Formation
The processes that lead to the formation and break-up of spirals are depicted in Fig. 3. Spontaneous spiral formation (cf. Fig. 3a) results from the leading tip colliding with a subsequent part of the chain. Volume exclusion then forces the tip to bent. By further forward movement, the chain winds to a spiral. Two spiral break-up mechanisms occurred in our simulations. The first is the thermally activated mechanism in Fig. 3b. The leading tip of the wound-up chain spontaneously changes direction and the spiral deforms. This break-up mechanism requires strong local bending and is therefore predominant for small ξ P . The second mechanism is spiral breakup by widening and is depicted in Fig. 3c. The bending potential widens the spiral until the leading tip looses contact to the filament end. This break-up mechanism is predominant when ξ P is too large for spontaneous spiral break-up. Because high stiffness is also unfavourable for spiral formation, spiral break-up by widening was almost exclusively observed in simulations that started with a spiral configuration.
To understand spiral formation more quantitatively, we introduce the spiral number where φ(s) is the bond orientation at position s along the contour of the filament, as measure for the instantaneous chain configuration. The definition is illustrated for three sample structures in Fig. 4a. It effectively measures how often the filament wraps around itself. The time evolution of s is depicted in Fig. 4b for the same P e and ξ P /L as in Fig. 2. At P e = 200, s is always close to zero. At P e = 1000, s behaves similarly, except that peaks with larger values for |s| occur occasionally, i.e., when spirals with a short lifetime form. At P e = 5000, extended plateaus develop at large |s|. The spirals are wound up stronger and have a much longer lifetime. Probability distributions p(|s|) are depicted in Fig. 4c. For the simulation without spirals (P e = 200), the his-togram resembles the right half of a Gaussian distribution. For the simulation in the weak spiral regime, p(s) is similar for low |s|, but also has a small peak at |s| ≈ 2−3. For strong spiral formation at P e = 5000, p(|s|) has only a small peak at low |s|, which corresponds to the elongated state, and a large peak at large |s|, the predominating spiral state. It turns out that the different regimes can be well distinguished by the kurtosis where . . . denotes the ensemble average and σ s is the standard deviation of s. Results for the kurtosis are shown in Fig. 4d for selected ξ P /L. β 2 ≈ 3 in the polymer regime, as expected for Gaussian distributions. In the weak spiral regime, the small peak at larger values increases the numerator in Eq. (13) and has only a weak impact on σ s , leading to an increase of β 2 . When the spiral state is dominating, σ s grows drastically, resulting in a much smaller kurtosis β 2 . Note that to reduce statistical uncertainties we symmetrized the s-distribution in the computation of β 2 by counting each measured |s| as +s and −s and only used data from the spiral state for simulations that do not show spiral break-up.
With the kurtosis as measure to characterize spiral formation, a phase diagram can be constructed as depicted in Fig. 4e. Low filament rigidity ξ P and high propulsion P e is beneficial for spiral formation. In particular, for a fixed propulsion strength per unit length, any chain will form spirals if it is sufficiently long, because increasing the chain length without modifying any other parameter corresponds to moving to the lower right in the phase diagram.
The dimensionless numbers ξ P /L and P e completely characterize the system if the filament diameter -or the filament aspect ratio -is of minor importance. This is the true in the entire polymer regime, where volumeexclusions interactions hardly come into play because of the elongated chain structure. For the spiral regimes, the aspect ratio has an impact on the structure of the spiral and does in this way influence the results. Which features of the spiral regime can be approximated well by the dimensionless numbers can be understood from the spiral formation and break-up mechanisms. The aspect ratio is hardly relevant for spiral formation and spiral break-up by widening, where the decisive moments are when the filament tip collides with subsequent parts of the chain, or when it looses contact to the chain end, respectively. That break-up by widening is characterized well by the dimensionless numbers is also confirmed by a series of simulations that we start from a spiral configuration in which we vary N , f p , κ, and k B T . It turns out that spirals will break up by widening if In contrast, spontaneous spiral break-up by a change of orientation of the leading tip is dependent on a strong local curvature close to the tip and the structure of the spiral, which in turn is dependent on the filament diameter. The dimensionless description does therefore not provide a full characterization of the strong spiral regime, where spontaneous spiral break-up is the only mechanism to escape the spiral state. This is also confirmed by results for spirals that never broke up (cf. dark red squares in Fig. 4e), results for ξ P /L = 0.2 and ξ P /L = 0.14 show non-monotonic behaviour in the direction of ξ P /L. This is a result of a combination of that the dimensionless description is only partially valid in this regime and that we chose N = 200 for ξ P /L < 0.2 but N = 100 for 0.2 ≥ ξ P /L ≥ 2.0 in our simulations, i.e., the aspect ratio L/σ is halved in our simulations for filaments with ξ P /L < 0.2.
Finally, the bead discretization with the chosen parameters favours a staggered arrangement of beads of contacting parts of the filament, 33,34 which implies an effective sliding friction between these parts. To study the importance of this effect, we increase the diameter of the beads at fixed bond length so that neighboring beads are heavily overlapping, which leads to a strongly smoothened interaction potential. Results for an increased diameter σ = 2L/N are shown in Fig. 4b and d. We find that the spiral formation frequency is hardly affected by smoothening the filament surface. In contrast, spontaneous spiral break-up is largely alleviated for the smoother filament, leading to decreased spiral life-time, as can be seen from the evolution of s for P e = 5000 in Fig. 4b. Smoother filaments thus show a qualitatively similar phase behaviour with slightly moved phase boundaries.

III.2. Structural Properties
The structural properties of the filament conformations can be best understood from the end-to-end vector r e , as depicted in Fig. 5. As long as no spirals form, simulation results are in good agreement with the Kratky-Porod model (valid for worm-like, non-active polymers without volume-exclusion interactions) that predicts 30,31 in two dimensions. At low ξ P /L, volume-exclusion interactions lead to slight deviations between the Kratky-Porod model and the simulation results. Strong deviations between the Kratky-Porod model and simulation results only occur in the strong spiral region in the phase diagram. The same trend was observed for the tangenttangent correlation function, the radius of gyration, and the static structure factor, but is not reported here to avoid unnecessary repetition.

III.3. Dynamic Properties
The characteristic filament motion can be understood from the mean square displacement (MSD) of a bead j +i relative to bead j, as shown in Fig. 6. For comparison, the MSD of the reference bead j of a passive filament is also shown. Note that the displacement of this bead is subdiffusive at the short lag times shown here. 35 Displacement functions of the propelled beads are always larger than in the passive case. The curves show three distinct regimes. At small lag times, the MSDs of active filaments display plateaus due to the average distance of the two beads j + i and j along the filament. At large lag times, the increased motion caused by activity effects the MSDs of the propelled beads to grow more rapidly than that of the passive bead. The relevant part of the MSD that shows that the characteristic filament motion is movement along its contour is at intermediate lag times, where the MSDs of the propelled beads pass through minima that touch the reference MSD for thermal motion. At that lag time, the bead j + i has moved approximately to the position of bead j at zero lag time. The beads have moved along the chain contour, similar to the movement of a train on a railway. The deviation from the exact starting position of bead j exactly matches the thermal motion, which results in the MSDs of the propelled beads touching the MSD of the passive filament. Thus, the characteristic movement of the filament is motion along its contour superimposed with thermal noise, as depicted in Fig. 7a. Note that ξ P /L = 0.3 was selected in Fig. 6 because this corresponds to a rather flexible filament, where stronger deviations from the characteristic railway motion might be expected. This type of motion was observed in all simulations in the polymer regime and weak spiral regime. In the strong spiral regime, the MSDs of the propelled beads even fall below the reference line for purely thermal motion.   The rotational diffusion can be accessed from the orientation of the end-to-end vector θ. Its mean square rotation (MSR) is given in Fig. 8a. Note that complete rotations around the axis are accounted for in our computations. θ(t) can therefore be much larger than 2π. In both the spiral and the elongated state, there is a regime at short lag times in which the MSR is dominated by the internal filament flexibility. For the spiral state, this regime is followed by a ballistic regime with MSR ∝ t 2 . For simulations in which the spirals break up spontaneously, a subsequent regime at high times with MSR ∝ t is expected but could not be detected in our simulations because of finite simulation time and strong noise at large lag times in the MSR. For the elongated state, the regime dominated by internal flexibility is followed by a diffusive regime with MSR ∝ t. The rotational diffusion coefficient D r can be extracted by fitting MSR = 2D r t to the regime of the MSR with gradient unity on a double-log scale. Measured D r are given in Fig. 8b as a function of the flexure number F. The diffusion coefficients collapse to a single curve, which has a plateau at low F and then grows linearly. Strong deviations from this trend are only observed at high F when the filament is in the weak spiral regime, and at F = P e = 0 for flexible filaments, where strong deviations from a rod-like shape increase D r .
The rotational diffusion coefficient D r can be predicted from the characteristic railway motion in Fig. 7 and the relation of the rotational diffusion coefficient to the au-tocorrelation function of the end-to-end tangent vector t e t e (t) · t e (0) = e −Drt , which is valid for lag times t that are sufficiently large such that variation of t e is not dominated by nondiffusive behaviour at early lag times caused by the filament flexibility (cf. Fig. 8a). With where t(s, t) is the tangent vector at position s of the filament at time t, the left hand side of Eq. (16) becomes where the order of summations has been changed to arrive at the right-hand side of Eq. (18). As a representation of the characteristic railway motion (cf. Fig. 7b), we write t(s, t) = t(s + v c t, 0). (19) Note that this equation disregards the passive equilibrium rotation D r,p . With Eq. Integrating Eq. (18) provides t e (t) · t e (0) = −ξ 2 P /L 2 e −vct/ξ P 2 − e L/ξ P − e −L/ξ P .
(21) A second order Taylor expansion in (small) L/ξ P then gives so that a comparison with Eq. (16) finally yields the activity-induced rotational diffusion Note that v c /ξ P = F/4τ . Assuming that uncorrelated activity-induced and thermal rotation D r,a and D r,p contribute to the overall rotation, we write where D r,p depends on ξ P /L and has the lower bound D r,p = (9/4) τ −1 for rod-like filaments. 36 As can be seen from Fig. 8b, Eq. (24) matches the simulated rotational diffusion coefficient accurately. The characteristics of the center-of-mass MSD are shown in Fig. 9. For the polymer regime, the typical S-shape of subsequent short-time diffusive, intermediatetime ballistic, and long-time effective diffusive behaviour develops; 37 stronger propulsion increases the MSD. An important difference compared to rigid bodies is that the transition time τ r = 1/D r to long-time diffusive behaviour is dependent on the propulsion strength.
When spiral formation becomes important, the general trend of the MSD changes, as shown in Fig. 9b for a flexible filament. In the polymer regime or weak spiral regime, increasing P e leads to a larger displacement. In the strong spiral regime, however, the MSD decreases. For very stable spirals, the MSD is only weakly affected by the propulsion and almost matches the case of purely diffusive motion.
The MSD for active point particles, spheres, or stiff rods is given by 1 where D t is the translational diffusion coefficient and v 0 is a ballistic velocity. It turns out that Eq. (25) can be used to describe the MSD for active filaments, when the three coefficients D t , v 0 , and D r are chosen properly. The translational diffusion coefficient is D t = L 2 /4τ = k B T /γ l L. We predict the rotational diffusion coefficient D r with Eq. (24). Finally the effective velocity can be expressed via as a balance of the net external force |F p | with the total friction force γ l Lv 0 . |F p | can conveniently be expressed as the propulsive force per bond f p times the end-to-end vector, thus leading to v 0 = f p r 2 e γ l L .
As shown in Fig. 9, using these correlations for the coefficients provides an accurate prediction of the MSD. The last item we address is the effect of propulsion on conformational sampling. Figure 10a shows results for the dynamic structure factor (28) averaged over different directions of q. In the phase without spirals, S(q, t)/S(q, 0) decays more rapidly with increasing P e, indicating a faster change of conformations with increasing propulsion. When spirals form, S(q, t)/S(q, 0) is indepenent of P e and larger than S(q, t)/S(q, 0) at P e = 0, indicating a slow change of conformations, which agrees with the observation of hardly any internal motion of the chain in this regime in our simulation output. Note that for the strong spiral regime, the data is from simulations where spirals formed spontaneously and did not break up. The depicted data is a result of averaging over the spiral states only.
To better quantify the behaviour of S(q, t), we compute the characteristic decay time of the dynamic structure factor Results for τ S (q) at q ≈ 5π/L, a q-vector large enough to capture the behaviour of mainly the internal degrees of freedom, are given in Fig 10b. τ S (q) decays slowly at low P e. At high P e, τ S decays inversely proportinal to P e when no spirals form This is consistent with the picture that instantaneous conformations are essentially identical to those of passive filaments, but they are traversed with velocity v c , corresponding to τ ∝ P e −1 . In the strong spiral regime, τ S is large and independent of P e, which is a sign for that conformational changes are irrelevant and that τ S is determined by the quasi-diffusive center of mass movement. Note that the measured τ S at different ξ P /L collapse to a single line for both the polymer and the strong spiral regime.

IV. DISCUSSION
The spontaneous formation of spirals is the feature dominating the overall behaviour of self-propelled filaments, both for dynamic and structural properties. Formation of spirals was previously observed for long, slender bacteria surrounded by short bacteria. 7 It was concluded that interaction with other active particles is a prerequisite for spiral formation. In contrast, the study at hand shows that spirals can form even for isolated filaments, as long as (i) the filament is sufficiently flexible, (ii) the propulsion is sufficiently strong, and (iii) excluded volume interactions force the tip of the filament to wind up.
The first two conditions will be met automatically for any real system by choosing L sufficiently large and leaving all other parameters constant (leading to increased P e and decreased ξ P /L, i.e., favouring spirals). Meeting the third condition can in general not be achieved so easily. A free-swimming filament in three dimensions or a filament in two dimensions with low resistance of crossing its own body will not form spirals. This is also one reason why spiral formation has not yet been observed in more experimental studies. Agents that are similar to our model are actin filaments or microtubules on a motility assay. The former have a high crossing probability, 38 formation of spirals is therefore not expected. The area enclosing the actin-filament parameter space in Fig. 4 must thus be understood as that the regime where the flexibilities and propulsion strengths permit spiral formation can in principle be reached in real systems, and not so much as that sufficiently long actin filaments will form spirals. Microtubules on dynein carpets, which have a much lower crossing probability, 6,34 will possibly form spirals if they are grown to sufficient size.
Overall, except for slender bacteria, 7 we are unaware of a microscopic example in which spiral formation was observed. Yet, the formation of spirals is a feature that deserves more attention in the future. First, formation of spirals is an extremely simple non-equilibrium phenomenon that, in contrast to many other phenomena of active matter, arises for a single self-propelled particle and cannot easily be mapped qualitatively to passive systems in which activity is replaced by attractive forces. It can thus be used as a model phenomenon for the study of non-equilibrium thermodynamics. Second, our model is very simple; a realization in experiment seems possible within the near future. Finally, the formation of spirals leads to a sudden, strong change in structural and dynamic properties. The effect can thus potentially be used as a switch on the microscopic scale.

V. SUMMARY AND OUTLOOK
We report an extensive study for the behaviour of dilute, self-propelled, worm-like filaments in two dimensions. The spontaneous formation and break-up of spirals is the feature that dominates the filament behaviour. Spiral formation is favoured by strong propulsion and low bending rigidity. Propulsion has a noticeable impact on structural properties only when spirals are dominating. The Kratky-Porod model 30 is therefore valid for filaments that are weakly propelled or have high bending rigidity. When spiral formation becomes significant, structural properties change drastically.
The characteristic filament motion is what we call the railway behaviour. The chain moves along its own contour superimposed with noise. With the understanding of the structural properties and the characteristic motion, rotational diffusion and the center-of-mass mean square displacement can be predicted to high accuracy when no spirals form. In contrast to rigid bodies, propulsion has an impact on the rotational diffusion coefficient. Finally, propulsion enhances conformational sampling in the regime without spirals.
An obvious next step is understanding the collective motion of such active filaments. We expect that our single filament results will help to understand the collective behaviour, which is nonetheless strongly influenced by the additional interactions. In particular collision with other constituents might enhance spiral formation and lead to swirl-like patterns.