# From flagellar undulations to collective motion:

predicting the dynamics of sperm suspensions

###### Abstract

Swimming cells and microorganisms are as diverse in their collective dynamics as they are in their individual shapes and propulsion mechanisms. Even for sperm cells, which have a stereotyped shape consisting of a cell body connected to a flexible flagellum, a wide range of collective dynamics is observed spanning from the formation of tightly packed groups to the display of larger-scale, turbulence-like motion. Using a detailed mathematical model that resolves flagellum dynamics, we perform simulations of sperm suspensions containing up to 1000 cells and explore the connection between individual and collective dynamics. We find that depending on the level of variation in individual dynamics from one swimmer to another, the sperm exhibit either a strong tendency to aggregate, or the suspension exhibits large-scale swirling. Hydrodynamic interactions govern the formation and evolution of both states. In addition, a quantitative analysis of the states reveals that the flows generated at the time-scale of flagellum undulations contribute significantly to the overall energy in the surrounding fluid, highlighting the importance of resolving these flows.

*Keywords:* sperm locomotion collective dynamics active suspensions fluid-structure interactions

Swimming cells and microorganisms encompass the entire range of cell types and exhibit a great variety of cell geometries and swimming strategies used to move in viscosity dominated environments [1, 2]. Some organisms, such as Volvox colonies, can be nearly spherical and utilize for propulsion thousands of relatively short flexible flagella distributed along their surface. Others, like sperm cells, have a single, long flagellum that propagates bending waves, leading to a large time-periodic shape change of the cell. Though stereotyped, the propagated waveforms, frequencies, and head shapes of sperm cells vary appreciably across different species [3].

There is not only great variation amongst individual cells, but also in the patterns that they form and in how populations are organized. This can range from algal plume formation in bioconvection [4, 5] to turbulence-like swirling in bacterial baths [6, 7, 8, 9]. Great variation is exhibited by sperm suspensions of different species that have been observed to form vortices near planar boundaries [10], aggregate into coherent sperm trains [11], or exhibit turbulence-like swirling as first documented nearly years ago by Lord Rothschild [12, 13, 14, 15]. What is not clear, however, is how the individual differences in sperm cells, such as geometry, flagellum length, waveform and flexibility, give rise to the differences in the way sperm populations organize themselves. Mathematical modeling and simulation provide a route to explore the connection between individual and collective dynamics where experiments might otherwise be very challenging.

In the past to years, there has been much work on the development of mathematical models [8, 16, 17, 18, 19, 20, 21, 22, 23, 24] of microorganism suspensions, especially to understand the turbulence-like state found in bacterial baths. To deal with large numbers of cells, the models rely on a reduced description of the swimmers, often treating them as simple, rigid objects (e.g., rods, dumbbells, spheres, ellipsoids) that interact through steady, dipolar flows and steric repulsion. These models have been effective in reproducing the large-scale motion of the population and relating its formation to the sign of the coefficient of dipolar flow induced by individuals.

While such models now broadly shape our thinking about suspensions of swimming microorganisms, they reduce the diversity of the microscopic world into a handful of parameters and a limited number of degrees of freedom. They ignore the complexity and time-dependence of cell shapes, as well as the time-dependent and beyond-dipolar features of the flows generated by the shape changes. Indeed, many of these details have been included in models of individual or small collections of swimmers. In the case of sperm cells, the addition of flagellar motion leads to hydrodynamically-induced attraction, synchronization, and phase-locking of planar, flagellar waveforms [25, 26, 27]. In simulations of to swimmers [28, 29] these effects were found to induce aggregation and clustering of cells. This tendency to aggregate may aid sperm in forming coherent groups such as sperm trains, however, it remains unclear how a turbulence-like state could then be reached.

In this paper, we perform detailed simulations of up to interacting sperm, resolving the coupled flagellar dynamics along with the flows generated at the undulation time-scale and at sub-flagellum length scales. We report that variation in the individual dynamics across the population, here controlled by the undulation frequency, can suppress the tendency to aggregate and instead lead to a density dependent turbulence-like state with hydrodynamics still playing a key, but now different role. Additionally, an analysis of the clustered and turbulence-like states reveals the strong influence of flagellar undulation on the quantities used to explore large-scale motion in swimmer suspensions. These results suggest that only minor variations in sperm behavior across species are needed to produce very distinct collective dynamics.

### Mathematical model

Our swimmer model is closely related to several models [30, 31, 32, 33, 34] employed to describe the dynamics of sperm cells in a viscous fluid and extends directly from a model previously used to capture undulatory locomotion through structured media [35]. We briefly summarize it here in the context of a single swimmer and provide a more detailed description in the SI.

Our swimmers consist of two elements, the flagellum and the cell head. The flagellum is treated as an inextensible, yet flexible beam of length and bending modulus , while an oblate spheroid of semi-major and minor axes and , respectively, represents the cell head. The overall swimmer length is taken to be , where the distance accounts for a linkage between the flagellum and the head. The flagellum is parametrized by arclength such that the position of a point along the flagellum is and the unit tangent at that point is . The flagellum is driven internally by the moments per unit length, , that arise due to the preferred curvature,

(1) |

where is the undulation frequency, is the wavenumber, is the phase, and is the amplitude. The linear decay in the preferred curvature amplitude near the free-end is introduced to better reproduce observed flagellar waveforms. Allowing the flagellum to also be subject to external applied forces, , and torques, , per unit length due to viscous stresses, the force and moment balances along the flagellum are

(2) | ||||

(3) |

where is the tension and is the bending moment. At one end , the flagellum is attached via a clamped-end condition to the cell head, while the other end remains free.

The continuous beam equations, (2) and (3), are discretized to obtain force and torque balances on segments of the flagellum, as well as the cell head (see SI for details). For the -th segment of the flagellum, the balances are

(4) | ||||

(5) |

where are torques arising due to bending, and are constraint forces and torques associated with tension, are the torques due to the preferred curvature, and and are the hydrodynamic forces and torques on the segment. The force and torque balances yield a low Reynolds number mobility problem for the translational and angular motion of the segments that is identical to that for a collection of rigid particles. We solve this mobility problem using a regularized multipole approach for Stokes flows known as the *force-coupling method* (FCM) [36, 37, 38, 35, 39]. In the limit of large , FCM provides an approximation of the hydrodynamics consistent with a regularized version of slender body theory [40, 41, 42, 43, 44, 31, 33]. After solving the mobility problem for the velocities and angular velocities of the flagellum segments and cell head, the differential equations for their positions and orientations are integrated in time using a second-order backwards differentiation scheme. Broyden’s method is then used to solve the resulting system of equations for the updated positions and orientations, and the Lagrange multipliers associated with the forces and torques arising from the inextensibility constraint.

We choose the parameters in our simulations to reproduce a cell geometry, flagellar waveform, and swimming speed close to that reported for ram sperm [45, 14] and in the range of other mammalian sperm [46, 47]. The head size is such that and . The wavenumber is , and the curvature amplitude is . Finally, we set the dimensionless parameter known as the *Sperm number*, which provides a measure of the ratio of the viscous to elastic forces [48, 49, 2, 35, 50], to , where is the dynamic viscosity. With these parameters, 23.7 undulation periods are required for the sperm to swim its flagellum length.

We perform three dimensional simulations in domains of size corresponding to 2D periodic thin films with finite thickness . The corresponding fluid flow is resolved on a grid of size . Free surface conditions at and are established by restricting the motion of the swimmers to the mid-plane and applying periodic boundary conditions in all three directions. This particular choice in boundary conditions and film thickness corresponds to those commonly employed in experiments on swimming microorganisms [6, 16, 51, 52].

When simulating multiple swimmers, a short-ranged repulsive force between nearly touching segments is included to capture steric repulsion between swimmers (see SI). The strength of this force at contact is determined by the parameter . This and all simulation parameters, including their values in simulation units, are summarized in Table 1.

Parameter | Description | Value in simulation units |
---|---|---|

Viscosity | ||

Sp | Sperm number | |

Curvature amplitude | ||

Bending modulus | ||

Segment radius and head height | ||

Head half axis | ||

Swimmer length | ||

Steric reference force | ||

Grid dimensions (in-plane) | ||

Grid dimension (normal) | ||

Linear domain size | ||

Domain height | ||

Number of flagellum segments per swimmer | ||

Number of swimmers | ||

Time-steps per undulation | ||

Relative std. dev. of frequencies |

### Time-dependent and higher-order features of the swimmer flow field

Examining the flow field produced by a single swimmer and the force-moments associated with this flow, we see clear differences between their instantaneous and period-averaged values. The flow at the mid-plane of the thin film and the swimmer at an instant in time are shown in Fig. 1A, while Fig. 1B shows the flow averaged over one undulation period. These flow fields are comparable to those given by other computational studies of individual swimmers propelled by an elastic filament [30, 28, 29, 53, 27, 54]. While the period-averaged flow closely resembles the dipolar flow generated by a so-called pusher, we find that the instantaneous flow field at any point in time is markedly different, containing features of higher-order forces singularities.

To analyze this further, we compute (see SI and [49]) the in-plane dipole and quadrupole force moments for the swimmer as functions of time. The symmetric, trace-less dipole tensor , the symmetric (in the first two indices), trace-less quadrupole tensor , and the anti-symmetric quadrupole tensor , are plotted over two undulation periods in Fig. 2. The time periodic nature of the force moments, especially their change in sign over the course of a period, are consistent with experimental measurements [51, 54] on flagellate microswimmers. The anti-symmetric dipole moment is always zero, as there is no net-torque on the swimmer. From these time series, we can extract the non-vanishing time-averages, the largest of which are the dipole coefficients , and the quadrupole coefficients , where is the in-plane drag on the cell head moving at the average swimming velocity. We note that for both the dipole and quadrupole these average values are much smaller than the maximum values attained during an undulation period. Quadrupole components ( and ) are indeed expected to contribute substantially to the fluid flow around spermatozoa [55, 33, 33] and have also been linked to swimmer velocity correlations in bacterial suspensions [56]. The maximum value of the dipole coefficient in our model is more than a factor of greater than its mean value, while for the quadrupole this increases to a factor of . Additionally, the relatively large values of the quadrupole coefficients produce strong flows that, although they decay faster than the dipolar flow, are non-negligible for significant distances from the swimmer. We estimate that the instantaneous quadrupolar flow dominates the dipolar flow up to separations of approximately .

From this, a picture emerges that not only are the time-dependent aspects of the dipolar flow much stronger than the average flow, but that higher-order contributions are non-negligible when time dependence is included. Both of these features can impact the swimmer-swimmer hydrodynamic interactions, especially in semi-dilute to dense suspensions. As we will see in the proceeding sections, these features, linked to the dynamics at the level of individual cells, do indeed propagate to larger scales and can significantly affect the observed coherent structures and the processes by which they form.

### Suspensions with a monodisperse distribution of undulation frequencies form clusters

Using the methods described in the previous section, we simulate the dynamics of a suspension of swimmers with for undulation periods. In this simulation, all swimmers are driven by a prescribed preferred curvature, Eq. 1, with the same undulation frequency, , and wave-number . Each swimmer has a different phase, , drawn randomly from a uniform distribution. The swimmers are initially distributed uniformly in the domain and are aligned to swim in the -direction. Swimmer-swimmer hydrodynamic interactions are incorporated by solving the full mobility problem that couples the motion of all swimmer segments and heads, while steric interactions are included through short-ranged, pairwise, repulsive forces between all flagellum segments and cell heads. The effective area fraction is . Based on the film thickness, the effective volume fraction, as used in [22], is approximately . For suspensions of pusher dipoles at this effective volume fraction, one would expect rapid decay of the initial polar order followed by the onset of large-scale, fluctuating motion [22, 19].

Fig. 3A shows the suspension after and undulation. periods. The corresponding video showing the complete evolution is provided as ESM. While the suspension quickly evolves from the initial polar state, the dynamics are very different from the long-wavelength bending instabilities predicted by both simulations of rigid, steady pushers [18, 22] and kinetic or continuum models [17, 19, 20, 21, 24]. Instead, we find that the swimmers have a strong tendency to aggregate and form clusters within which flagella are aligned. A quantitative analysis of this cluster growth is included in the SI. Once small clusters have formed, they remain intact and can attract other clusters in their vicinity. As two attracting clusters approach each other, they either rotate to swim in the same direction and merge to form a bigger, synchronized cluster, or instead swim in opposing directions and move apart at an enhanced relative velocity. Earlier simulations of sperm [28, 29] showed similar aggregation dynamics in two dimensional suspensions, but here we see that aggregation also dominates in larger, denser suspensions where long-range hydrodynamics could be expected to produce large-scale swirls and jets.

The aggregation process is driven by the hydrodynamic attraction and phase-locking behavior that has been found previously with models of undulatory swimmers [26, 28, 29, 57, 34, 27]. We measure this for our system by examining in detail the pairwise interactions between two swimmers moving in the -direction with phase differences and . The vector fields in Fig. 4 show the displacement over four undulation periods of one swimmer centered at relative to another placed at the origin. When the swimmers are far apart, all three vector fields resemble the displacement field associated with pusher dipoles. At smaller separations, however, the fields differ and the influence of the phase difference can be readily seen. We see that there is clear attraction, if the swimmers are in phase. When the phase difference is , the swimmers attract and simultaneously move relative to each other as to align the crests of the flagellar waves. If they are out of phase, the swimmers again move relative to one another attempting to align wave crests, however, the relative shift is much larger. The large displacements observed near the swimmer’s ends are due to the steric forces between the head of one swimmer and the flagellum of the other.

### Suspensions with stochastically varying frequencies reach a turbulence-like state

In the previous section, we saw that flagellum undulations can produce flows that lead to hydrodynamically-induced aggregation and cluster formation. While this provides a mechanism for the formation of coherent groups, such as wood mouse sperm trains [11], the patterns differ markedly from the large-scale swirling found in experiments on ram sperm suspensions [14, 15]. In real samples, sperm cells have distributions of waveforms, undulation frequencies or geometries. Accordingly, the parameters appearing in our model should not necessarily be the same for all individuals and should be chosen to accurately describe variations from individual to individual across the population.

We explore how variations in individual dynamics across the population impact large-scale suspension dynamics by introducing stochastic modulations of the undulation frequency into the model, allowing flagellar waves to differ from individual to individual and to vary over time [46, 58]. As we solve for the flagellum dynamics, changing the frequency also leads to changes in the waveform. The individual frequencies are drawn from a log-normal distribution after time intervals of length as described in the SI. We set the mean frequency to and, in view of matching experimental data [45, 58], set the standard deviation to be .

We repeat our simulations of swimmers, but now allowing for stochastic variations in the frequencies. Snapshots of the suspension evolution are shown in Fig. 3C and corresponding videos can be found as part of the ESM. The stochastic variations in the waveforms suppress the ability of neighboring flagella to synchronize. Rather than forming clusters, the swimmers now form waves, swirls and vortices, as seen in concentrated, quasi-two-dimensional ram sperm suspensions [14]. By running an equivalent simulation starting from an isotropic initial configuration, it was verified that the qualitative features of the fully non-linear state are independent of the initial conditions. We also note that a sufficiently large is needed to successfully suppress aggregation and observe large-scale motion of the suspension. Simulations of smaller ensembles (see SI) show that lower values of still exhibit cluster formation (e.g., or ). In these cases, stochastic frequency changes are rarely large enough to break the synchronization of adjacent flagella.

While hydrodynamic interactions no longer lead to aggregation, they still play a strong role in the evolution of the suspension. We perform simulations where we have removed the hydrodynamic interactions by solving the mobility problem using resistive force-theory (RFT) instead of FCM. As describe in the SI, the drag coefficients appearing in the RFT are chosen as to closely match the swimming speeds and flagellar waves obtained with the full FCM model. Fig. 3B shows snapshots from two different RFT simulations with and where the initial swimmer directions are either all aligned or distributed isotropically. The corresponding videos are provided as part of the ESM. The time-evolution of these suspensions is very different from those with hydrodynamic interactions. Most strikingly, for the initially polar suspension, long bending-waves are entirely absent and at long times we observe “laning” states similar to those found in simulations of self-propelled rods without hydrodynamic interactions [29, 59, 9]. When the initial state is isotropic, laning is is not observed (Fig. 3B), indicating further differences between the RFT and full simulations. Additionally, in thin films, the hydrodynamic interactions between the swimmers are longer ranged than they would be in bulk. As a result, we find that as we increase the film thickness, though the instability of the polar state takes longer to develop, we still observe (see SI) the eventual onset of large-scale motion qualitatively similar to that found in the thinner film cases.

### Order parameters exhibit density dependence

The polar and nematic order of the suspension is dependent on the effective area fraction, . We perform a series of smaller stochastic simulations () ranging from to swimmers run to and compare the order parameters

(6) | ||||

(7) |

which measure the global alignment with the initial swimming direction () and the global nematic order with respect to the initial direction () [22]. Here, is the -th swimmer’s mean orientation and denotes the average over all swimmers. Videos of the simulations are provided as ESM.

We find that the directional order, given by , decreases with time and the decay rate increases with swimmer density (Fig. 5A). The time-scale of the decay is comparable to that found for rod models [22, 23]. The nematic order, given by , initially decreases with time. However, for higher densities, we see large fluctuations emerge after the initial decay (Fig. 5B). The appearance of these fluctuations during the relatively gentle decay of reflects periodic folding and rotation of large patches of aligned swimmers. Fig. 5 also shows and for initially polar RFT simulations with the same effective area fractions. With hydrodynamic interactions removed, we find that higher density suspensions retain their alignment for longer times due to enhanced local caging.

### Flagellar undulations contribute significantly to the energy and velocity spectra

To further quantify the large-scale motion, we compute the fluid energy spectrum,

(8) |

where is the spatial Fourier transform of the fluid velocity and denotes time-averaging while signifies averaging over all wave vectors of magnitude . Fig. 6A shows from simulations with both fixed and stochastically varying frequencies and for different values of . With the exception of the pronounced peak near the undulation wavelength, the spectra are similar to those found in experiments on concentrated sperm suspensions [14]. In particular, we also find that the spectra decay like for large (Fig. 6A). Given that this power-law decay appears over sub-swimmer length scales, we associate it with random forcing of the fluid by the flagella. A calculation presented in the SI demonstrates that an identical decay is obtained when the fluid is driven by spatially uncorrelated forcing. This is different from the decay observed at small scales in suspensions of self-propelled rods [23]. This decay was attributed to the sharp jump in the forcing along the length of each swimmer. Thus, the details regarding how the swimmers propel themselves can lead to changes in the spectrum at length scales comparable to the swimmer size.

Further, we find that the fluid energy spectrum is dominated by contributions at the time-scale of flagellar undulations, even at long length scales. Despite qualitative differences in their dynamics and the quantitative differences in the evolution of the order parameters, is nearly identical for all stochastic simulations, regardless of . The collapse illustrates that scales with the number of swimmers. This is in contrast with large increases with seen at low in simulations of self-propelled rod suspensions [23]. Surprisingly, the curves for the fixed and stochastically varying frequency simulations also do not differ drastically, despite the striking differences in their suspension evolution. We expect that this is due to the energy injected at the short time-scales associated with flagellar undulations being large compared with that due to the large-scale motion of the suspension. To limit the contribution of short time-scales to , we apply a low-pass filter to by computing a running time-average of with a window of before performing the transform. The spectra, , computed using the filtered velocity fields are shown in Fig. 6B. With the short time-scales suppressed, the dependence on becomes visible. We now observe larger values of at low when the sperm density is high, revealing quantitatively the large-scale fluid motion at the time scale of suspension evolution. Additionally, the spectrum for the fixed-frequency simulations is now almost flat at small indicative of the lack of long range correlations.

Following from previous work on bacterial suspensions [9, 59], we also examine the spectrum of the center of mass velocity field, , where and are the center of mass position and velocity, respectively, as defined in the SI. Fig. 6C shows the swimmer velocity spectrum, , for different values of . For each , the lowest spatial modes (lowest ) provide the largest contribution to the total spectrum, with the contribution increasing as increases. From these maximum values, the spectra then quickly decrease with to a flat profile at scales below the flagellum length, corresponding to the level of noise generated in constructing the swimmer velocity field.

Due to flagellar undulations, the swimmers’ center of mass velocities oscillate about their mean values at the time-scale of the undulation frequency. We can again remove contributions of the short time-scale by filtering the swimmers’ center of mass velocities using a running time-average with window size . Fig. 6D shows the spectra associated with the filtered swimmer velocity field. Filtering reduces values of the swimmer velocity spectrum for all , indicating contributions at short time-scales at all length scales. For lower values of , the spectrum is reduced by an approximately constant factor across all , while for higher values of , the reduction is more pronounced at shorter scales.

### Discussion and Conclusions

In this study, we have used numerical simulation to show that variability in sperm flagellum dynamics across a suspension can lead to large changes in collective dynamics. We have incorporated these variations through stochastic changes in the actuation frequency that, in turn, lead to changes in waveform and amplitude. In actual sperm suspensions, variations are also likely to be present in flagellum length and cell geometry, and further, the flagellar waveform is likely to be more complex and fully three-dimensional [47, 60]. We expect these features to have a similar, aggregation-limiting effect as the stochastic variation in frequency. As such, measurements quantifying the distributions of individual sperm properties could not only provide a better classification of the individual cells, but also a better understanding of their collective dynamics.

Our fixed frequency simulations show that flagellar synchronization leads to aggregation and clustering. Certain properties associated with sperm that are found to self-organize into coherent groups, such as wood mouse sperm, might promote synchronization. It is known that these sperm typically have flagella that are longer than those of sperm from other species. By allowing for larger amplitudes, longer flagella may enhance the higher-order, time-dependent flows that we find are responsible for aggregation. Longer, thinner flagella would also bend and flex with greater ease in response to flows generated by neighboring sperm, further reinforcing the tendency to synchronize and aggregate. Additionally, microtubule sliding driving the undulations may depend on the external load experienced by the flagellum. This could allow for a non-trivial coupling between the actuation mechanism and the surrounding flow field, which, in turn, may also promote synchronisation of neighbouring flagella and cluster formation.

The turbulence-like state achieved when there is sufficient variation in the undulation frequency is both quantitatively and qualitatively similar to that observed in suspensions of ram and bull semen. While longer time simulations would be needed to explore this state in more detail, simulations starting from either polar or isotropic states eventually exhibit similar qualitative dynamics and have the same fluid energy spectrum at the final simulation times. Both this and the many previous results obtained using reduced models suggest that the observed jets and swirls are features of a final turbulence-like state. The presence of this state has been empirically correlated with sample fertility [61]. We have shown that its onset depends on sperm density, and its observation could therefore be serving as an indicator of sperm count, a well-known measure of sample fertility.

While in the simulations presented here the fluid domain is fully three-dimensional, the motion of the sperm cells is restricted to a plane. Even though layering in sperm suspensions has been observed in experiments [14] and sperm are attracted to move along planar surfaces [62], the imposed two dimensional motion in our simulations may enhance the role of steric interactions. For example, in a real sample when two sperm cells collide, it is likely that one will pass over the other and they will both experience some change in their swimming directions. In our simulations, the sperm cannot pass by each other as easily and the collisions can dramatically affect their swimming directions, leading to near alignment or anti-alignment. The interactions of self-driven filaments [57] and the trajectories of interacting sperm pairs [34] have recently been explored in simulation. These studies have shown that perturbations from planarity can lead to more complex trajectories, as well as a dependence on the elasticity of the flagella. With further advances in numerical methods for simulating flexible filaments and more accurate models for three-dimensional flagellar waveforms, it will be possible to explore these effects through direct simulation of sperm suspensions. Our theoretical model for the elastic flagellum and the FCM approach to the mobility problem carry over to fully three dimensional simulations. While there will be additional computational costs associated with the increased size of the fluid domain, more limiting is the technical challenge involved in keeping track of the fully three-dimensional deformation of the flagellum and, at the same time, using implicit time integration to handle numerical stiffness. We are currently developing the numerical methods to address this challenge.

Additionally, it would be of interest to ascertain suspension dynamics at even larger length scales and for longer timescales than those that may be accessible through direct simulation. Continuum models incorporating the time dependence of the flows generated by the swimmers have recently been developed [63, 64, 65] and could possibly be adapted to capture the effects of flagellar undulations seen in our simulations and used to explore nonlinear suspension dynamics at these scales.

Along with variations amongst the cells, collective dynamics are likely to be influenced by more complex features of the environment, such as nearby boundaries that could curve sperm trajectories [58], or particles or other micro-scale structures in the surrounding fluid that are known to enhance swimming speeds of undulatory swimmers [35]. Non-Newtonian features of the surrounding fluid, such as viscoelasticity, typically associated with biological fluids have also been been found to promote clustering [66]. These, and other important, perhaps more complex effects, such as chemotaxis, could further contribute to the richness and diversity of sperm collective dynamics.

### Competing interests

The authors declare no competing interests.

### Authors’ Contributions

S.F.S. and E.E.K. designed research, performed research, analyzed data, and wrote the paper.

### Acknowledgements

The authors would like to thank Pierre Degond and Demetrios Papageorgiou for helpful discussions.

### Funding

S.F.S. gratefully acknowledges funding by an Imperial College President’s PhD scholarship. E.E.K. acknowledges support from EPSRC grant EP/P013651/1.

### Supporting Information

see below

## References

- [1] E. M. Purcell. Life at low reynolds number. American Journal of Physics, 45(1):3–11, 1977.
- [2] Eric Lauga and Thomas R Powers. The hydrodynamics of swimming microorganisms. Rep. Prog. Phys., 72(9):096601, Aug 2009.
- [3] Scott Pitnick, David J. Hosken, and Tim R. Birkhead. Sperm morphological diversity. In Tim R. Birkhead, David J. Hosken, and Scott Pitnick, editors, Sperm Biology, pages 69 – 149. Academic Press, London, 2009.
- [4] T J Pedley and J O Kessler. Hydrodynamic phenomena in suspensions of swimming microorganisms. Annu. Rev. Fluid Mech., 24:313–358, 1992.
- [5] N.A. Hill and T.J. Pedley. Bioconvection. Fluid Dyn. Res., 37(1–2):1 – 20, 2005. Biofluiddynamics.
- [6] Xiao-Lun Wu and Albert Libchaber. Particle diffusion in a quasi-two-dimensional bacterial bath. Phys. Rev. Lett., 84:3017–3020, Mar 2000.
- [7] Christopher Dombrowski, Luis Cisneros, Sunita Chatkaew, Raymond E. Goldstein, and John O. Kessler. Self-concentration and large-scale coherence in bacterial dynamics. Phys. Rev. Lett., 93:098103, Aug 2004.
- [8] D. L. Koch and G. Subramanian. Collective hydrodynamics and of and swimming microorganisms and living fluids and. Annu. Rev. Fluid Mech., 43:637–59, 2011.
- [9] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans. Meso-scale turbulence in living fluids. Proc. Natl. Acad. Sci. USA, 109(36):14308–14313, 2012.
- [10] Ingmar H. Riedel, Karsten Kruse, and Jonathon Howard. A self-organized vortex array of hydrodynamically entrained sperm cells. Science, 309(5732):300–303, 2005.
- [11] H. Moore, K. Dvoráková, N. Jenkins, and W. Breed. Exceptional sperm cooperation in the wood mouse. Nature, 418(6894):174–177, Jul 2002.
- [12] Lord Rothschild. The activity of ram spermatozoa. J. Exp. Biol., 25(3):219–226, 1948.
- [13] Lord Rothschild. Measurement of sperm activity before artificial insemination. Nature, 163:358–359, 1949.
- [14] Adama Creppy, Olivier Praud, Xavier Druart, Philippa L. Kohnke, and Franck Plouraboué. Turbulence of swarming sperm. Phys. Rev. E, 92(032722), Sep 2015.
- [15] Adama Creppy, Franck Plouraboué, Olivier Praud, Xavier Druart, Sébastien Cazin, Hui Yu, and Pierre Degond. Symmetry-breaking phase transitions in highly concentrated semen. Journal of The Royal Society Interface, 13(123), 2016.
- [16] Igor S. Aranson, Andrey Sokolov, John O. Kessler, and Raymond E. Goldstein. Model for dynamical coherence in thin films of self-propelled microorganisms. Phys. Rev. E, 75(040901), Apr 2007.
- [17] R. A. Simha and Sriram Ramaswamy. Hydrodynamic fluctuations and instabilities in ordered suspensions of self-propelled particles. Phys. Rev. Lett., 89(5), Jul 2002.
- [18] Juan P. Hernandez-Ortiz, Christopher G. Stoltz, and Michael D. Graham. Transport and collective dynamics in suspensions of confined swimming particles. Phys. Rev. Lett., 95:204501, Nov 2005.
- [19] David Saintillan and Michael J. Shelley. Instabilities and pattern formation in active particle suspensions: Kinetic theory and continuum simulations. Phys. Rev. Lett., 100(178103), Apr 2008.
- [20] Christel Hohenegger and Michael J. Shelley. Stability of active suspensions. Phys. Rev. E, 81(046311), Apr 2010.
- [21] Aparna Baskaran and M. Cristina Marchetti. Statistical mechanics and hydrodynamics of bacterial suspensions. Proc. Natl. Acad. Sci. USA, 106(37):15567–15572, 2009.
- [22] David Saintillan and Michael J. Shelley. Orientational order and instabilities in suspensions of self-locomoting rods. Phys. Rev. Lett., 99(058102), Jul 2007.
- [23] D. Saintillan and M. J. Shelley. Emergence of coherent structures and large-scale flows in motile suspensions. J. R. Soc. Interface, 9:571–585, Aug 2011.
- [24] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, Madan Rao, and R. Aditi Simha. Hydrodynamics of soft active matter. Rev. Mod. Phys., 85:1143–1189, Jul 2013.
- [25] Lisa J Fauci. Interaction of oscillating filaments: A computational study. J. Comput. Phys., 86(2):294 – 313, 1990.
- [26] Gwynn J. Elfring and Eric Lauga. Hydrodynamic phase locking of swimming microorganisms. Phys. Rev. Lett., 103(088101), Aug 2009.
- [27] Sarah D. Olson and Lisa J. Fauci. Hydrodynamic interactions of sheets vs filaments: Synchronization, attraction, and alignment. Phys. Fluids, 27(121901), Dec 2015.
- [28] Yingzi Yang, Jens Elgeti, and Gerhard Gompper. Cooperation of sperm in two dimensions: Synchronization, attraction, and aggregation through hydrodynamic interactions. Phys. Rev. E, 78(061903), Dec 2008.
- [29] Yingzi Yang, Vincent Marceau, and Gerhard Gompper. Swarm behavior of self-propelled rods and swimming flagella. Phys. Rev. E, 82(031904), Sep 2010.
- [30] Lisa J Fauci and Charles S Peskin. A computational model of aquatic animal locomotion. Journal of Computational Physics, 77(1):85 – 108, 1988.
- [31] D. J. Smith. A boundary element regularized stokeslet method applied to cilia- and flagella-driven flow. Proc. R. Soc. A, 465(2112):3605–3626, 2009.
- [32] Sarah D. Olson, Susan S. Suarez, and Lisa J. Fauci. Coupling biochemistry and hydrodynamics captures hyperactivated sperm motility in a simple flagellar model. Journal of Theoretical Biology, 283(1):203 – 216, 2011.
- [33] Montenegro-Johnson, Thomas D., Smith, Andrew A., Smith, David J., Loghin, Daniel, and Blake, John R. Modelling the fluid mechanics of cilia and flagella in reproduction and development. Eur. Phys. J. E, 35(10):111, 2012.
- [34] Julie Simons, Lisa Fauci, and Ricardo Cortez. A fully three-dimensional model of the interaction of driven elastic filaments in a stokes flow with applications to sperm motility. J. Biomech., 48(9):1639–1651, Jun 2015.
- [35] T. Majmudar, E. E. Keaveny, J. Zhang, and M. J. Shelley. Experiments and theory of undulatory locomotion in a simple structured medium. J. Roy. Soc. Interface, 9(73):1809–1823, Feb 2012.
- [36] M.R. Maxey and B.K. Patel. Localized force representations for particles sedimenting in stokes flow. International Journal of Multiphase Flow, 27(9):1603 – 1626, 2001.
- [37] Sune Lomholt and Martin R. Maxey. Force-coupling method for particulate two-phase flow: Stokes flow. Journal of Computational Physics, 184(2):381 – 405, 2003.
- [38] D. Liu, E.E. Keaveny, M.R. Maxey, and G.E. Karniadakis. Force-coupling method for flows with ellipsoidal particles. J. Comput. Phys., 228(10):3559–3581, Jun 2009.
- [39] Martin Maxey. Simulation methods for particulate flows and concentrated suspensions. Annual Review of Fluid Mechanics, 49(1):171–193, 2017.
- [40] G. K. Batchelor. Slender-body theory for particles of arbitrary cross-section in stokes flow. Journal of Fluid Mechanics, 44(3):419–440, 1970.
- [41] R. E. Johnson and C. J. Brokaw. Flagellar hydrodynamics – a comparison between resisitive-force theory and slender-body theory. Biophys. J., 25:113–127, January 1979.
- [42] James Lighthill. Reinterpreting the basic theorem of flagellar hydrodynamics. Journal of Engineering Mathematics, 30(1):25–34, Mar 1996.
- [43] Anna-Karin Tornberg and Michael J. Shelley. Simulating the dynamics and interactions of flexible fibers in stokes flows. J. Comput. Phys., 196(1):8 – 40, 2004.
- [44] R. Cortez and M. Nicholas. Slender body theory for stokes flows with regularized forces. Comm. App. Math. and Comp. Sci., 7(1):2012, 2012.
- [45] Margaret Anne Denehy. The propulsion of nonrotating ram and oyster spermatozoa. Biology of reproduction, 13(1):17–29, 1975.
- [46] Robert Rikmenspoel. The tail movement of bull spermatozoa: Observations and model calculations. Biophys. J., 5(4):365 – 392, 1965.
- [47] D. M. Woolley. Evidence for “twisted plane” undulations in golden hamster sperm tails. J. Cell. Biol., 75(3):851–865, 1977.
- [48] C. P. Lowe. Dynamics of filaments: modelling the dynamics of driven microfilaments. Philos. T. Roy. Soc. B, 358(1437):1543–1550, Sep 2003.
- [49] E. E. Keaveny and M. R. Maxey. Spiral swimming of an artificial micro-swimmer. J. Fluid Mech., 598:293–319, Feb 2008.
- [50] Blaise Delmotte, Eric Climent, and Franck Plouraboué. A general formulation of bead models applied to flexible fibers and active filaments at low reynolds number. J. Comput. Phys., 286:14–37, Apr 2015.
- [51] Jeffrey S. Guasto, Karl A. Johnson, and J. P. Gollub. Oscillatory flows induced by microorganisms swimming in two dimensions. Phys. Rev. Lett., 105(168102), Oct 2010.
- [52] Hüseyin Kurtuldu, Jeffrey S. Guasto, Karl A. Johnson, and J. P. Gollub. Enhancement of biomixing by swimming algal cells in two-dimensional films. Proc. Natl. Acad. Sci. USA, 108(26):10391–10395, 2011.
- [53] E.A. Gaffney, H. Gadêlha, D.J. Smith, J.R. Blake, and J.C. Kirkman-Brown. Mammalian sperm motility: Observation and theory. Annual Review of Fluid Mechanics, 43(1):501–528, 2011.
- [54] Kenta Ishimoto, Hermes Gadêlha, Eamonn A. Gaffney, David J. Smith, and Jackson Kirkman-Brown. Coarse-graining the fluid flow around a human sperm. Phys. Rev. Lett., 118:124501, Mar 2017.
- [55] David J. Smith and John R. Blake. Surface accumulation of spermatozoa: a fluid dynamic phenomenon. Mathematical Scientist, 34:74–87, Dec 2009.
- [56] Qian Liao, Ganesh Subramanian, Matthew P. DeLisa, Donald L. Koch, and Mingming Wu. Pair velocity correlations among swimming escherichia coli bacteria are determined by force-quadrupole hydrodynamic interactions. Physics of Fluids, 19(6):061701, 2007.
- [57] I. Llopis, I. Pagonabarraga, M. Cosentino Lagomarsino, and C. P. Lowe. Cooperative motion of intrinsic and actuated semiflexible swimmers. Phys. Rev. E, 87(3), Mar 2013.
- [58] B. M. Friedrich, I. H. Riedel-Kruse, J. Howard, and F. Julicher. High-precision tracking of sperm swimming fine structure provides strong test of resistive force theory. J. Exp. Biol., 213(8):1226–1234, Mar 2010.
- [59] H H Wensink and H Löwen. Emergent states in dense systems of active rods: from swarming to turbulence. J. Phys. Condens. Matter, 24(46):464130, 2012.
- [60] Anton Bukatin, Igor Kukhtevich, Norbert Stoop, Jörn Dunkel, and Vasily Kantsler. Bimodal rheotactic behavior reflects flagellar beat asymmetry in human sperm cells. Proc. Natl. Acad. Sci. USA, 112(52):15904–15909, 2015.
- [61] Ingrid David, Philippa Kohnke, Gilles Lagriffoul, Olivier Praud, Franck Plouarboué, Pierre Degond, and Xavier Druart. Mass sperm motility is associated with fertility in sheep. Anim. Reprod. Sci., 161:75 – 81, 2015.
- [62] Rothschild. Non-random distribution of bull spermatozoa in a drop of sperm suspension. Nature, 198(488):1221, 1963.
- [63] Sebastian Fürthauer and Sriram Ramaswamy. Phase-synchronized state of oriented active fluids. Phys. Rev. Lett., 111(238102), Dec 2013.
- [64] M. Leoni and T. B. Liverpool. Synchronization and liquid crystalline order in soft active fluids. Phys. Rev. Lett., 112(148104), April 2014.
- [65] Tommaso Brotto, Denis Bartolo, and David Saintillan. Spontaneous flows in suspensions of active cyclic swimmers. Journal of Nonlinear Science, 25(5):1125–1139, Oct 2015.
- [66] Chih-kuan Tung, Chungwei Lin, Benedict Harvey, Alyssa G. Fiore, Florencia Ardon, Mingming Wu, and Susan S. Suarez. Fluid viscoelasticity promotes collective swimming of sperm. Sci. Rep., 7(1):3152, 2017.
- [67] Uri M Ascher and Linda R Petzold. Computer methods for ordinary differential equations and differential-algebraic equations, volume 61. SIAM, 1998.
- [68] C. G. Broyden. A class of methods for solving nonlinear simultaneous equations. Math. Comp., 19(92):577–593, 1965.

## Supporting Information

### Swimmer model

As each swimmer is identical, for clarity, we describe the model for a single swimmer and subsequently describe how it generalizes when considering a suspension. Recall that a swimmer is composed of two elements, its cell head and its flagellum, where the cell head is treated as a rigid oblate spheroid with semi-major and minor axes and , respectively. Using the arc-length parametrization where , the force and moment balance equations are

where is the tension and is the bending moment. The external force per unit length, , and torque per unit length, , arise due to the viscous stresses along the flagellum and, in the case of multiple swimmers, steric interactions between neighboring swimmers. The torques per unit length arise due to the time-dependent preferred curvature. At the free end , the boundary conditions are such that the tension and moment are zero, while at the attachment point to the cell head, we require a clamped-end condition be satisfied.

To solve these equations numerically, we first discretize the flagellum into segments of length . The segments have positions and orientations for . The orientations are the discrete representation of the centerline tangents at the segment positions. Considering the tension and bending moment at the midpoints between adjacent segments and applying central differencing, we obtain the discretized force and torque balances

(9) | |||

(10) |

where . Multiplying through by and introducing and as the total applied force and torque, respectively, on segment , we may write Eqs. [9] and [10] as

(11) | ||||

(12) |

where , , and . The driving torques, , arising from the preferred curvature are given by

where for , and

as well as, .

The tension enforces flagellum inextensibility at the level of each segment. In the discrete setting, they are the Lagrange multipliers associated with the constraints,

(13) |

In the discretized system, the boundary conditions at the free end are satisfied by taking . The clamped-end condition on the cell head is recovered by treating the head as another segment of size with the bending moment and inextensibility constraint computed with respect to the attachment point rather than the head center.

Until now, we have discussed the model in the context of a single swimmer. To allow for swimmers, we must have force and moment balances, one pair for each of the flagella. As the driving torques, tension, and bending moments are internal to each flagellum, the force and moment balances are coupled only through the external forces and torques. In our simulations, the coupling arises due to steric repulsion between the segments and heads and hydrodynamic interactions. Accordingly, the total external force on segment is given by where is the steric force and is the hydrodynamic force on the segment. The only external torque on the segment is the viscous torque,

The steric force on segment is given by the generic repulsive barrier force

where the sum runs over , the set of all segments/heads, , with and . The parameter specifies the center-to-center distance for particles and at contact and sets the range over which the barrier force acts. If and are both segments, the contact distance is , where as if and are both heads, we have . Lastly, if and are a head-segment pair then . The parameter sets the strength of the repulsion at contact, which in our simulations is .

### Force-coupling method

The force and torque balances, Eqs. (11) and (12) respectively, establish a low Reynolds number mobility problem whose solution is the coupled motion of the particles through the fluid. To solve this mobility problem, we utilize the force-coupling method (FCM) [36, 37, 38], which, for the sake of clarity, we describe here in the case of spherical particles of equal radius, however, it readily extends to spheroidal particles [38], such as our cell heads, and allows for polydispersity in particle sizes.

In FCM forces and torques on the particles are projected onto the fluid using a truncated and regularized force multipole expansion. This leads to the incompressible Stokes flow

for pressure field, , and fluid velocity field, . The Stokes flow is generated by the forces and torques each particle exerts on the fluid, which, in the case of our swimmers, will be , and . The force and torques are projected to the fluid via the Gaussian functions, which for particle are given by

The motion of the particles is obtained by volume averaging the resulting fluid velocity using the same Gaussian functions. Specifically, the translational and angular velocities are given by

If the particles have radius , FCM reproduces the correct Stokes drag and viscous torque with and . In our simulations, we preserve these ratios between the particle and the Gaussian envelope sizes. The segments are taken to have radius , while the semi major-axes of the spheroidal heads are and the semi minor-axes are . All simulation parameters are summarized in table 1 and the datasets provided have the same units.

### Updating positions and orientations

As our simulations are quasi-2D in that the motion of the swimmers is restricted to a plane, we have for all and can introduce angle such that . Therefore, after computing the and for each segment and head, the positions and orientations are obtained by integrating in time

(14) | ||||

(15) |

subject to the inextensibility constraints given by [13]. To do this numerically, we rely on a second-order implicit BDF scheme [67] and Broyden’s method [68] to find the updated positions and orientations, as well as the Lagrange multipliers .

### Swimmer center of mass and director

We define the center of mass of the -th swimmer as

Here, is the position of the -th swimmer’s head, and is that of its -the segment. The center of mass velocity is then given by

The director of swimmer is given by

### Force multipole expansion for a single swimmer

The singular force density for a single swimmer is given by

where the asymmetric dipole is related to the torque through . We have suppressed the superscript index, since we are considering only one swimmer. Expanding this force density about the swimmer’s center of mass, , and defining the position of the center of mass relative to the -th segment/head as , we have

The first term is zero because there is no external force on the swimmer. The remaining two terms are the force dipole and quadrupole,

respectively, which drive the flow observed in the limit . We compute the symmetric, trace-less and anti-symmetric parts of these tensors, given by

as functions of time. These are shown in Fig. 2 in the main text, with the exception of which is zero for all time as there is no net-torque acting on the swimmers.

### Stochastic variation of the undulation frequency

We introduce stochastic changes of the undulation frequencies to explore how variability within a population, as well as the fluctuations of individual frequencies over time can impact collective dynamics. After every fixed average period , a new undulation frequency, , is randomly assigned to a swimmer to replace its previous value, . New frequencies are drawn from a log-normal distribution with parameters () given by

where the mean frequency is and the standard deviation of the distribution is . To enforce continuity in time for each swimmer’s preferred curvature , the phase must also be updated from to . Specifically, at the time when the frequency changes, we require

This can be fulfilled by taking

for each individual swimmer.

We note that an alternative approach to introducing stochasticity is to assign a random frequency to each swimmer at the beginning of the simulation and keep it fixed for all time. This was adopted previously in [28]. We found that on time-scales that are relatively short compared to the evolution of the suspension this leads to qualitatively similar results as those obtained by allowing the frequency to vary over time. Over longer times, however, we expect, based on our monodisperse frequency simulations, that keeping the frequencies fixed could potentially lead to preferential clustering of swimmers with similar frequencies and hence structurally different dynamics.

### Identifying and counting swimmer clusters

To effectively identify swimmer clusters, we must account for relative swimmer orientations in addition to relative distances. Based on this, swimmer is taken to belong to a cluster if the modified distance

with at least one other swimmer in the cluster is . The parameter controls the influence of particle alignment and is chosen to be . With this choice of , swimmers with orientations differing by an angle of must be half as far apart as aligned swimmers to still be considered clustered. Using this notion of distance, we employ the hierarchical cluster algorithm in Matlab for cluster identification. The number of clusters over time for the two simulations with individual swimmers are shown in Fig. 7. Individual swimmers are counted as clusters of size one. We see that the number of clusters for the fixed frequency simulation is much lower than that for the simulation where the undulation frequency is allowed to vary stochastically. Additionally, we see that when stochasticity is introduced, cluster growth is halted very early in the simulation.

### Effect of film thickness on suspension dynamics

As the flows induced by singular force distributions decay more slowly in 2D than in 3D, we expect the hydrodynamic interactions to play a stronger role as and, as a result, the instability to exhibit a faster growth rate for thinner films. Figs. 8 and 9 show results from simulations containing 85 swimmers in domains of linear size and for thicknesses , and . The effective area fraction is and the swimmers are initially in a polar state.

In accordance with the slower formation of the bending waves as shown in Fig. 8, the slower decay (see Fig. 9) of the order parameters and (as defined in the main text) for thicker films indicate that the shorter range of the hydrodynamic interactions lead to slower growth rates of the instability. The smallest value of used in Fig. 9 matches that for the larger simulations presented in the main text. The order parameter decay rate in those simulations, however, is faster still due to the longer modes afforded by the larger system size.

### Effect of undulation frequency variability on aggregation

As mentioned in the main text, low variability in the swimmers’ undulation frequencies does not entirely suppress cluster formation and still allows for neighboring swimmers to synchronize and align their waveforms.

Snapshots from simulations of a system of size with initially aligned swimmers and different values of reveal that for small , cluster formation persists (see Fig. 10). As increases, the number of visually identifiable clusters is reduced and we see instead the emergence of a bending wave.

### Tuning the resistive force theory

To remove the hydrodynamic interactions between the swimmers yet still allow for propulsion via undulation, we solve the mobility problem using a drag-based resistive force theory (RFT) [2, 41] rather than FCM. With the drag-based model, the velocity and angular velocity of segment of a flagellum are related to the force and torque on the segment through

where , , and are mobility coefficients for the segments. These mobility coefficients are initially estimated based on FCM simulations of straight filaments and subsequently adjusted to ensure that the swimmer waveform and free swimming velocity are comparable to those given by the simulations with FCM, see Fig. 11. For the cell heads, we have for all swimmers

where and are mobility coefficients taken from FCM simulations of a single spheroid. The numerical values used in the simulations are

### Energy density spectrum of a randomly forced Stokesian fluid

In this section, we show that a behavior of the fluid energy spectrum can be related to spatially uncorrelated forcing of a 2D fluid. In the following, the Fourier transform and its inverse are defined as

We begin by considering a force density, , restricted to a square domain,

where and is the indicator function defined on the set . The Fourier transform of the restricted force density is then

The Fourier transform of the fluid velocity resulting from the restricted force density is given by

where is the Fourier representation of the inverse Stokes operator, , and .

The fluid energy density spectrum associated with is

where is related to the two-dimensional wave-vector through and denotes the ensemble average. The energy density spectrum of the infinite system will then be given by .

As

the fluid energy density spectrum is directly related to the forcing through

For to hold, we must have independent of . If, in addition, we require statistically independent and identical force components that do not depend on , we have that

where is a constant independent of . The dependence on the system size guarantees that is well-defined and non-zero as and that the total energy will scale with the system size.

Assuming that the forcing is a statistically stationary spatial process, we may express the correlations of the force density as

As the force density is strictly restricted to , we may replace the integral over with one over all of . Thus,

which, from the convolution theorem, may be expressed as

Finally, using the expression for given above, we see that

for any . Allowing then yields

and we see that a decay in the fluid energy density spectrum is given by uncorrelated spatial forcing of the fluid.