The fluid dynamics of collective vortex structures of plant-animal worms

Circular milling, a stunning manifestation of collective motion, is found across the natural world, from fish shoals to army ants. It has been observed recently that the plant-animal worm $Symsagittifera~roscoffensis$ exhibits circular milling behaviour, both in shallow pools at the beach and in Petri dishes in the laboratory. Here we investigate this phenomenon, through experiment and theory, from a fluid dynamical viewpoint, focusing on the effect that an established circular mill has on the surrounding fluid. Unlike systems such as confined bacterial suspensions and collections of molecular motors and filaments that exhibit spontaneous circulatory behaviour, and which are modelled as force dipoles, the front-back symmetry of individual worms precludes a stresslet contribution. Instead, singularities such as source dipoles and Stokes quadrupoles are expected to dominate. A series of models is analyzed to understand the contributions of these singularities to the azimuthal flow fields generated by a mill, in light of the particular boundary conditions that hold for flow in a Petri dish. A model that treats a circular mill as a rigid rotating disc that generates a Stokes flow is shown to capture basic experimental results well, and gives insights into the emergence and stability of multiple mill systems.


Introduction
From the flocking of birds to the schooling of fish, collective motion, global group dynamics resulting from the interactions of many individuals, occurs all across the natural world. A visually striking example of this is collective vortex behaviour -the spontaneous motion of large numbers of organisms moving in periodic orbits about a common centre. Studied for over a century since Jean-Henri Fabre (1899) first reported the spontaneous formation of continuous loops in columns of pine processionary caterpillars, circular milling has been observed in many species, including army ants (Couzin & Franks 2003), the bacterium Bacillus subtilis (Cisneros et al. 2007;Wioland et al. 2013) and fish (Calovi et al. 2014).
It has been discovered recently that the marine acoel worm Symsagittifera roscoffensis (Bourlat & Hejnol 2009) forms circular mills, both naturally in rivulets on intertidal sand (Sendova-Franks et al. 2018), and in a shallow layer of sea water in a Petri dish (Franks et al. 2016). S. roscoffensis, more commonly known as the 'plant-animal worm' (figure 1(a), Keeble (1910)), engages in a photosymbiotic relationship (Bailly et al. 2014) with the marine alga Tetraselmis convolutae (Norris et al. 1980). The photosynthetic activity of the algae in hospite under the worm epidermis provides the nutrients required to sustain the host. The worms reside on the upper part of the foreshore (regions which are typically underwater for around two hours before and after high tide) of Atlantic coast beaches in colonies of many millions ( figure 1 (b)). It is hypothesised that this circular milling allows worms to self-organise into dense biofilms that, covered by a mucus layer, optimise the absorption of light by the algae for photosynthesis (Franks et al. 2016).
Here, motivated by prior experimental work and in light of new results we consider a range of issues surrounding the fluid dynamical description of mills, with particular attention to the fluid velocity field that is generated by a circular mill and the effect that this flow has on the mill itself. In §2, we describe the parameters of field experiments on mills performed on the Isle of Guernsey and outline the range of questions they pose. Systems such as bacterial suspensions (Woodhouse & Goldstein 2012;Wioland et al. 2013), collections of sperm cells near surfaces (Riedel et al. 2005) and assemblies of molecular motors and biofilaments (Sumino et al. 2012) can spontaneously form vortex-like patterns superficially similar to worm mills. However, their theoretical descriptions (Saintillan & Shelley 2008) treat the constituents as force dipoles (stresslets). The front-back symmetry of ciliated worms would suggest that a stresslet contribution is small if not entirely absent and thus higher-order singularities such as source dipoles ought to appear. Such is the case with the spherical alga Volvox whose flow field has been measured in great detail (Drescher et al. 2010). As there has been little if any work on the collective behaviour of suspensions of singularities beyond stresslets, §3 provides background theoretical considerations on this problem. First, a detailed examination of the relationship between the cilia-generated flow over the surface of an individual worm and the far-field flow behaviour is given within a prolate squirmer model, with which we confirm the absence of a stresslet contribution for a suitably symmetric surface slip velocity and show that the far field is dominated by the source dipole and force quadrupole contributions. Insight into those singularity components that lead to azimuthal flow around a mill composed of such swimmers is obtained by averaging the far-field behaviour over a circular orbit, which is equivalent to considering a ring of swimmers. Within the squirmer model we find that it is the Stokes quadrupole that gives the leadingorder contribution. A model of a complete mill can be constructed from a suitable oblate squirmer, whose far-field behaviour is that of a rotlet dipole. With the particular boundary conditions that hold in a Petri dish the far field from this singularity decays exponentially with z dependence sin (πz/2H). Then in §4, we consider a model of a mill as a rigid disc with varying radius c = c(t) rotating in a Stokes flow. A lubrication analysis of a highly-simplified model of a mill is presented as motivation, and to elucidate the effect of the boundary conditions for this system. Hence, we vertically average the governing equations by setting the z dependence explicitly to sin (πz/2H), deriving a Brinkman-like equation for the vertically averaged velocity flow field. In general, further analytic progress can not be made. However in two particular limits, namely when the mill is close to and when the mill is far away from the centre of the Petri dish, an analytic solution for the fluid velocity field and hence for the force that the flow exerts on the disc can be derived. In §5, we demonstrate the strong agreement between what is predicted by the model and what is observed experimentally. In particular, the viscous force on the disc points in the direction perpendicular to the line between the centres of the mill and of the Petri dish. Hence the centre of the disc will drift on a circle with centre the middle of the Petri dish, precisely as observed experimentally.
Finally, in §6 we extend the analysis to systems with more than one mill, focusing on the simplest binary mill structure. We utilise the knowledge gained from §4 to explain from a purely fluid dynamical viewpoint a large raft of experimental observations, including where a second mill forms and in what direction it rotates, and the conditions for which a second mill will not form. We can also make predictions for the stability of the resulting binary circular mill systems.

Experimental Methods
Here we describe field experiments done during 12-19 June, 2019 in the Peninsula Hotel Guernsey on the Isle of Guernsey, a channel island near the coast of France. Worms were collected from a nearby beach (49 • 29 45.3 N 2 • 33 14.3 W) just prior to the experiments, minimising the perturbations in the worms' physiology and behaviour resulting from removal from their natural environment. As shown in figure 2(a), Petri dishes of diameter 10 cm were filled with sea water up to a depth of 8 mm. Around ten thousand worms were placed into the dish using 3 ml plastic Pasteur pipettes. The subsequent evolution of the system, including the spontaneous formation of circular mills, was recorded at 25 fps using a Canon Eos 5d Mark II camera equipped with a Canon macro lens MP-E 65 mm f/2.8 and a 1.5× magnifier, mounted above the dish on a copy stand. The system was illuminated uniformly through by a light box below the Petri dishes and LED lights located around them.
In some experiments, small drops of azo dye were injected into the Petri dish using a plastic Pasteur pipette to act as a tracer to track the motion of the fluid. Figure 2(b) is a montage showing the temporal evolution of a red dyed region of fluid, namely streaklines of the flow. As can be seen, the circular mill generates a clockwise flow that is in the opposite direction to the anticlockwise direction of rotation of the worms, that is, a backflow generated by the worms pushing themselves through the fluid.
An instantaneous image of a mill shows that its edge is not well defined. In order to overcome this, every hundred frames (i.e. four seconds of footage) were averaged together to create a coarser time lapse video. This averaging sharpens the mill edges considerably since this process differentiates between worms entering or leaving the mill and worms actually in the mill. Then, the location and radius of the circular mill in each frame were extracted manually, utilising a GUI interface in MATLAB to semi-automate this process.
Appendix A collects relevant information on the many experiments carried out in the field. Selected videos can be found in the Supplementary Material. Circular milling in this system has not previously been studied using the kinds of methods now common in the study of active matter (Marchetti et al. 2013). There are open experimental questions at various levels of organisation in this setup that mirror those that have been successfully answered for bacterial, algal and other microswimmer systems, including measurements of flow fields around individual swimmers, pairwise interactions between them, the temporal dynamics of mill formation from individuals, the flow fields around the mills and the dynamics of the mills themselves within their confining containers. Here our focus experimentally is on the latter; the drift of a mill centre within a Petri dish and the formation of binary mill systems.

From Individual Worms to Mills
We begin with fluid dynamical considerations at the level of individual worms to derive key results that will then be utilised in §4 to motivate a mathematical model for a circular mill. Working in modified prolate spheroidal coordinates, we find that the leading order fluid velocity in the far field produced by the locomotion of a single worm can be expressed in terms of fundamental Stokes flow point singularities as the superposition of a source dipole and a Stokes quadrupole. Among other implications, this result shows that the proper Reynolds number for worm locomotion is given by the swimming speed U , length and diameter ρ as U ( ρ 2 ) 1/3 /2ν ≈ 0.36, so worms swim in the laminar regime. We then consider two possible models for a circular mill. Picturing a mill as the superposition of many rings of worms, we find that the resulting net flow is azimuthal, that is, not in the vertical z direction. Alternately, considering a mill as an oblate squirmer with axisymmetric swirl, we find that away from the mill, the forcing can be expressed as a rotlet dipole and thus the flow has z dependence of the form sin(πz/2H).

Locomotion of an individual worm
The individuals of S. roscoffensis studied in the present experiments have a broad distribution of sizes; their length ranges from ≈ 1.5 − 6 mm, with mean¯ = 3 mm, and diameters ρ falling in the range ≈ 0.2 − 0.6 mm, with meanρ ≈ 0.35 mm. Worm locomotion arises from the collective action of carpets of cilia over the entire body surface, each ≈ 20 µm long, beating at ≈ 50 Hz. Muscles within the organism allow it to bend, and thereby alter its swimming direction (Bailly et al. (2014)). In unbounded fluid, their average swimming speed of individuals is U ≈ 2 mm/s. To model the fluid velocity field produced by a worm, we follow Pöhnl et al. (2020) and consider a spheroidal, rigid and impermeable squirmer (Lighthill 1952) with semi-minor axis b x and semi-major axis b z swimming at speed U = U e z so the z-axis lies along the major axis. The squirmer moves through prescribing a tangential slip velocity u s at its surface Σ. Neglecting inertia, the fluid flow in the swimmer frame u satisfies together with the force-free condition where σ is the stress tensor. We now switch to the modified prolate spheroidal coordinates (τ, ξ, ϕ), utilising the transformations τ = 1 2c x and the squirmer boundary is mapped to the surface τ = τ 0 = b z /c = constant. In this coordinate system, assuming an axisymmetric flow u = u τ e τ + u ξ e ξ and axisymmetric tangential slip velocity u s = u s e ξ , the Stokes streamfunction ψ satisfies Taking from Dassios et al. (1994) the general separable solution for the stream function in prolate spheroidal coordinates and applying the boundary conditions given in (3.2) and (3.3), as in Pöhnl et al. (2020), we obtain where the g n satisfy g 2 (τ ) = C 4 H 4 (τ ) + D 2 H 2 (τ ) − 2c 2 U G 2 (τ ), (3.7a) where G n and H n are Gegenbauer functions of the first and second kind respectively and P 1 n = − √ 1 − x 2 dP n /dx are the associated Legendre polynomials. The integration constants C n and D n are set by the boundary conditions g n (τ 0 ) = 0 and dg n dτ τ =τ0 = τ 0 c 2 n(n − 1)B n−1 for n = 2, 3 · · · (3.8) Here, the {B n } n 1 are the coefficients in the series expansion of u s = τ 0 n 1 B n V n (ξ) using the set of functions {V n = (τ 2 0 − ξ 2 ) −1/2 P 1 n (ξ)} n 1 , which is a basis over the space of L 2 functions satisfying f (ξ = ±1) = 0 together with the inner product <> w = 1 −1 w dξ and the weight function w = τ 2 0 − ξ 2 . Furthermore, the swimming speed U can be expressed in terms of the {B n } n 1 using 9) namely only odd enumerated modes contribute to the squirmer's swimming velocity. Hence, from now on we only consider the case when the forcing is only a linear combination of the odd modes i.e. B 2n+2 = 0 → g 2n+1 = 0 → C 2n+1 = D 2n+1 = 0 ∀n ∈ Z + . When the prescribed forcing arises purely from the first mode, i.e. u s = τ 0 B 1 V 1 , C 2n:n 1 and D 2n:n 1 simplify to become D 2 = −2B 1 c 2 τ 0 (τ 2 0 − 1) and C 2n:n 1 = D 2n:n 2 = 0. (3.10) When the forcing arises from a higher order mode, i.e. u s = τ 0 B 2n+1 V 2n+1 where n 1, D 2 and C 4 simplify to become i.e. the only n dependence arises from the U . Moving into the lab frame, in the far field (τ 1) the dominant term in the expansion for ψ comes from H 2 (τ ) = 1/3τ + · · ·, so ψ = 1 3τ

12c)
Converting this back to vector notation yields where u D and u Q , the flows generated by a source dipole and a Stokes quadrupole respectively, satisfy Thus, the far field fluid velocity field decays like 1/r 3 , consisting of a combination of a source dipole and a Stokes quadrupole. The far field fluid velocity generated by a higher order than one odd mode squirmer contains both source dipole and quadrupole components with the quadrupole component dominating as τ 0 → 1. Similarly, the far field fluid velocity generated by a mode one squirmer is purely a source dipole. Using Lauga (2020), this is the same as an efficient spherical squirmer (forcing only arising from the first mode) with effective radiusã = c 3 τ 0 (τ 2 0 −1). When τ 0 → ∞, namely the spherical limit, as expectedã → b x = b z = a, the radius of the sphere. When τ 0 → 1, the elongated limit,ã → (b z b 2 x ) 1/3 . Thus, at the scale of an individual worm, the Reynolds number Uã/ν in water (ν = 1 is the correct length scale for locomotion of an individual worm. Inertial effects are modest and individual worms swim in the laminar regime.

Ring of Spheriodal Squirmers
Given the results in above, it is natural to ask which singularities associated with individual worms contribute to the azimuthal flow around a mill. This can be investigated by averaging over the contributions from a swimmer in a circular orbit, as has been done in the stresslet case (Michelin & Lauga 2010). Hence, consider a spheroidal squirmer swimming clockwise horizontally in a circle of radius c, instantaneously located at the point P = (c cos θ, c sin θ, 0) and orientated in the direction p = (sin θ, − cos θ, 0), utilising a Cartesian coordinate system (x, y, z) with origin at the centre of the circle. If each squirmer generates a source dipole u D , the fluid velocity u DX (c, θ) at X = (R, 0, 0) is where r = (R 2 + c 2 − 2cR cos θ) 1/2 . The total velocity u Dring (R) at X due to a ring of clockwise swimming worms of radius c with line density λ ring is then Hence, a ring of uniformly distributed source dipole swimmers generates no net flow field outside of the ring. By contrast, the flow field u Qring (R) due to a ring of clockwise swimming worms, each generating a Stokes quadrupole, is finite: This decays in the far field as 1/R 4 . We conclude that, viewing the mill in terms of its individual constituents, it is the Stokes quadrupole from individual swimmers that drives the dominant azimuthal flow.
3.3. Oblate Squirmer with Swirl Further insight into the flow field generated by a mill can be obtained by viewing it as a single, self-propelled object with some distribution of velocity on its surface arising from the many cilia of the constituent worms. With a shape like a pancake, it can be modelled as an oblate squirmer with axisymmetric swirl. First, consider a prolate squirmer with aspect ratio r e = b x /b z rotating in the ϕ direction with imposed surface flow u ϕ0 (ξ) e ϕ in free space. Assuming that the generated fluid flow is purely in the ϕ direction with no ϕ dependence, the ϕ component of the Stokes equations, µ(∇ 2 u) ϕ = 0, becomes This admits the general separable solution that tends to zero at infinity where C pn are constants and as before P 1 n (ξ) and Q 1 n (τ ) are associated Legendre polynomials. Furthermore, since the squirmer is force and torque free, C p1 = 0. Decomposing u ϕ0 using the basis {V n (ξ)}, i.e. u ϕ0 = ∞ n=2 C n0 V n (ξ), we find that C pn = C n0 /W pn (τ 0 ) where τ 0 = 1/ 1 − (r e ) 2 and r e = b x /b z 1 is the aspect ratio of the spheroid. Note that in the spherical limit (r e → 1, b where C n are constants, and we recover the general form for a spherical squirmer with swirl (Pak & Lauga 2014;Pedley et al. 2016).
Returning to the general case, the dominant term in the far field arises from mode 2, where x k points in the z direction. This is a rotlet dipole. Now, using Dassios et al. (1994), to compute the velocity fluid for an oblate squirmer with swirl of aspect ratio r e > 1, we translate the results from the prolate spheroidal coordinate system (τ, ξ, ϕ) to the oblate spheroidal coordinate system (λ, ξ, ϕ) using the can be expressed in terms of the Cartesian coordinates (x, y, z) using Similarly to the prolate case, in the far field the second mode dominates, giving a fluid velocity field also in the form of a rotlet dipole, This result is intuitive; in the absence of a net torque on the object there can not be a rotlet contribution, so analogously to the case of a single bacterium whose body rotates opposite to that of its rotating helical flagellum, the rotlet dipole is the first nonvanishing rotational singularity.
We close this section by asking how a free-space singularity of the type in (3.25) is modified when placed in a fluid layer with the boundary conditions of a Petri dish. Here we quote from a lengthy discussion (Fortune et al. 2020) of a number of cases that complements earlier work (Liron & Mochon 1975) on singularities bounded by two no-slip walls; the leading order contribution in the far field to the fluid velocity from an oblate squirmer with swirl at z = h between a no-slip lower surface at z = 0 and an upper free surface at z = H is where ρ 2 = x 2 + y 2 and K 1 is the Bessel function. The flow decays exponentially away from the squirmer with decay length 2H/π.

Mathematical Model for a Circular Mill
4.1. Background We now proceed to develop a model the collective vortex structures observed experimentally in §2. A laboratory mill of the kind studied here typically has a radius c in the range 5 − 20 mm and rotates roughly as a rigid body with period T = 2πc/U in the range 15 − 60 s and angular frequency ω = U/c in the range 0.4 − 0.1 s −1 whereas in §3.1 the average worm swimming speed U ≈ 2 mm/s. Almost all the worms swim just above the bottom of the Petri dish in a layer typically only one worm thick, with even in the densest regions of the mill at most two or three worms on top of each other.
We observe minimal variation in the height H of the water in the Petri dish. Furthermore, by tracking dye streaklines we also observe minimal fluid flow in the vertical (z) direction. These observations are consistent with the considerations in §3; we can picture a circular mill as a superposition of many rings of worms, each of which lies in a horizontal plane. Combining (3.13), (3.16) and (3.17), the far field flow field for each ring is azimuthal, not in the vertical direction, and thus the net flow for the circular mill is horizontal as well.
As a final reduced model, we make the further simplification of considering a mill as a rotating disc with a defined centre b(t), radius c(t) and height d(t), quantities that are allowed to vary as a function of time. To maximise swimming efficiency, isolated worms away from the mill will propel themselves mostly from the first order mode V 1 given in §3.1. Since this fluid velocity field decays rapidly away from a worm and using §3.2 is zero across a full orbit, we neglect the fluid flow generated by worms away from the mill. Since a mill contains a high density of worms together with the interstitial viscoelastic mucus, we assume that the disc is rigid. Finally, since the locomotion of the worms generates a fluid backflow in the opposite direction to their motion, the disc is assumed to rotate in the opposite angular direction to that of the worms.
A common mathematical tool for solving problems in a Stokes flow is to express the forcing as the sum of a finite set of fundamental Stokes flow point singularities (Jeong & Moffatt 1992;Crowdy & Or 2010). Given in §3.3, by considering the circular mill as an rigid oblate squirmer with swirl, the dominant contribution from the forcing can be approximated as a rotlet dipole. However from (3.26), the leading order contribution in the far field for a rotlet dipole trapped between a lower rigid no-slip boundary and an upper free surface which deforms minimally has z dependence of the form sin (πz/2H). Hence, we then vertically average the governing equations by setting the z dependence of u to be precisely sin (πz/2H) i.e. we employ the factorisation where U = U (r) is independent of z and the factor π/2 is for convenience. Thus, a suitable rotational Reynolds number on the scale of a mill is Re ∼ U X/ν ≈ 10 where X = 2H/π. Moreover, the dominant velocities are azimuthal, with gradients in the radial direction. This suggests that the fluid dynamics of milling is certainly in the laminar regime, and the neglect of inertial terms is justified.

Defining Notation
As shown in Figure 3 and in Supplementary Video 1, we define a coordinate basis (x, y, z) with origin at the centre of the Petri dish P, where the bottom of the dish is at z = 0 and the free surface is at z = H, a constant. We model an established circular mill, rotating a distance d 0 above the bottom of a circular Petri dish with angular velocity −Ω, as a rigid disc of radius c and height d with

Lubrication Picture
Insight into the mill dynamics comes first from an extremely simplified calculation within lubrication theory in which the disc (mill) has a prescribed azimuthal slip velocity on its bottom surface and a simple no-slip condition on its top surface, as if only the bottom of worms have beating cilia. This artiface allows the boundary conditions at z = 0 and z = H to be satisfied easily. For the thin film of fluid between the bottom of the mill and the bottom of the dish, namely {(x, y, z): 0 z d 0 ; , x 2 + y 2 c 2 }, let the imposed slip velocity at z = d 0 be u slip = u slip (r)ê θ . In the absence of any pressure gradients, the general results of lubrication theory dictate a linear velocity profile for the flow u b in the film, where Ω is the as yet unknown angular velocity of the disc. The flow in the region above the mill is simply u(r, z) = Ωr, independent of z; it rotates with the disc as a rigid body. The torque T b on the underside of the mill is while there are no torques from the flow above because of its z-independence.
Since the mill as a whole is torque free, T b = 0 and we deduce If u slip has solid body like character, u slip = u 0 r/c, then Ω = −u 0 /c and u b = 0. This "stealthy" mill generates at leading order no net flow in the gap between the mill and the bottom of the dish and it is analogous to the stealthy spherical squirmer with swirl that generates no external flow (Lauga 2020). Any slip velocity other than the solid body form will generate flow in the layer, and one notes generically that it is in the opposite direction to the slip velocity. This is consistent with the phenomenology shown in figure 2 involving the backwards advection of dye injected near a mill. Hence, from now on, we will assume that the mill effectively imposes a constant velocity boundary condition at the edge of the mill, (x + b) 2 + y 2 = c 2 .

Full Governing Equations
Assuming Stokes flow with fluid velocity u = (u x , u y , w = 0), pressure p = p(x, y, z) and viscosity µ, the governing equations are µ∇ 2 u = ∇p , ∇ · u = 0. (4.5) Employing no slip boundary conditions (Batchelor (1967)) at both the outer edge and the bottom of the Petri dish yield u = 0 at z = 0 , u = 0 at r = x 2 + y 2 = 1. (4.6) On the surface of the fluid, z = H, the dynamic boundary condition is where σ f and σ a = −p 0 I are stress tensors for the fluid and the air respectively. Finally, the boundary conditions on the surface of the mill become 8c) where the tangent and normal vectors e t and e n satisfy e n = 1 (4.9b) 4.5. Vertically-averaged Governing Equations Defining r as the in-plane coordinates (x, y), as set out in §4.1, we employ the factorisation where κ = π/2H plays the role of the inverse Debye screening length in screened electrostatics, and the z-dependent prefactor guarantees both the lower no slip and the upper stress free vertical boundary conditions. Vertically averaging, i.e considering H −1 H 0 · · · dz, gives the Brinkman-like equation µ ∇ 2 − κ 2 U = ∇p , ∇ · U = 0, (4.11) U = U x e x + U y e y = U n e n + U t e t has a corresponding Stokes streamfunction ϕ satisfying together with boundary conditions u n = 0 at r = x 2 + y 2 = 1, (4.13a) u t = 0 at r = x 2 + y 2 = 1, (4.13b) The fluid dynamics of collective vortex structures of plant-animal worms 13 u t = Ωc on Γ = {(x, y) : (x + b) 2 + y 2 = c 2 }, (4.13d) i.e. no-slip is imposed at the edge of the Petri dish while a constant azimuthal velocity boundary condition is imposed at {(x, y) : (x + b) 2 + y 2 = c 2 }. Using separation of variables, the general series solution in polar coordinates to 4.12 is cos nθ A n r n + B n r −n + C n I n (κr) + D n K n (κr) + ∞ n=1 sin nθ Ã n r n +B n r −n +C n I n (κr) +D n K n (κr) , (4.14) where {I n (r) , K n (r)} are solutions of the first and second kind respectively for the modified Bessel equation f + f /r − f (1 + n 2 /r 2 ) = 0. In general, this system does not admit an analytic solution. However, significant analytic progress can be made in two particular limits, namely when the mill is close to and when the mill is far away from the centre of the Petri dish.

The Drift of the Circular Mill Centre
The flow exerts a force F on the mill where F = F xx + F yŷ satisfies Here the bipolar basis vectorsη andξ satisfyη = fx + gŷ andξ = gx − fŷ where f = 1 − cosh η cos ξ cosh η − cos ξ and g = − sin ξ sinh η cosh η − cos ξ , (5.2) while σ ηη and σ ηξ are components of the stress tensor σ ij = −pδ ij + µ (∂u i /∂x j + ∂u j /∂x i ). Since this is not a standard result given in the literature (Wakiya (1975) is the closest reference which can be found), for completeness appendix B.1 gives the full form of ∂u i /∂x j when expressed in bipolar coordinates for general u.
This system, in a domain symmetric about the line θ = 0, is forced by a fluid flow even in θ. Hence, since it admits a general separable form where each term is either even or odd in θ (4.14), ϕ is even in θ and hence p and σ rr are also even in θ. Similarly, σ rθ is odd in θ and hence from rewriting (5.1) in terms of cylindrical polar coordinates, F x = 0. This force causes the mill centre to slowly drift on a larger timescale than the period of rotation of a mill, maintaining a constant distance from the centre of the Petri dish.
In general, F y does not admit an analytic form. However, as in §4.6 and §4.7, further progress can be made analytically for circular mills both close to and far away from the centre of the Petri dish.

Far-field Circular Mill
Since (4.11) reduces in this case to the Stokes equations, F x = 0 follows immediately by utilising the properties of a Stokes flow. Reversing time and then reflecting in the x axis returns back to the original geometry but with the sign of F x flipped i.e. F x = −F x → F x = 0. Substituting (4.36) into (4.11) and then integrating gives the pressure p = 2µ a 2 ( E sinh η sin ξ + F sinh 2η sin 2ξ − 2F sinh η sin ξ + G cosh 2η sin 2ξ − 2G cosh η sin ξ). (5.10) Shifting the basis vectors back to Cartesian coordinates, the force can be expressed in the form where f x and f y are explicit functions of η, ξ and {A, B, C, E, F, G}. However, f x is odd with respect to ξ at η = η 2 since +2G cosh 2 η sinh η) = 0 at η = η 2 .
(5.12) Therefore as expected F x = 0. f y can be similarly simplified, removing the terms odd in ξ, to give where g i = g i (η) : i ∈ {0, 1, 2, 3, 4} are given for completeness in Appendix B.2 while I n satisfies I n = π π cos (n ξ) (cosh η − cos ξ) 3 dξ. (5.14) To investigate this force more quantitatively, we numerically calculate F y from (5.13) as a function of b and c by utilising MATLAB's symbolic variable toolbox. Figure 6a plots F y as a function of c for a range of values for b (chosen to demonstrate the full phase space of behaviour of F y (b, c)). Large mills have positive F y (in the grey region), i.e. the mill centre drifts in the same angular direction as the worms. Small mills have negative F y (in the white region), i.e. the mill centre drifts in the opposite direction to the worms. Hence, we define the critical radius c = c (b) as the mill radius at which F y = 0 (plotted as a function of b in Figure 6b). Note that when b is large the critical geometry is a lubrication flow since b + c → 1. Furthermore, c has a maximum of 0.222 at b = 0.70.

Comparison with Experiments
We now compare these predictions with experimental data for b(t) and c(t), generated using the methodology given in §2. Unlike the simple circle considered in the model, the shape of a real circular mill is complicated. Not only does a mill at any one time consist of thousands of individual worms but also, as the mill evolves, this population changes as worms enter and leave. Hence, mills typically have constantly varying effective radii and are not simply connected. Furthermore, the edges of a circular mill are not well-defined, leading to a greater uncertainty in measuring the mill radius. However, despite these complications, the experimental results agree well with the predictions made above in §5. Within experimental uncertainty, b is constant i.e. the centre of the mills do indeed drift on circles centred at the middle of the arena. Furthermore, the direction of drift also matches with the theory given in §5.2 for the force on a mill in the far field.
To illustrate this, consider Figure 7 which presents graphically the experimental data for two representative experiments which sit at either end of the phase space of mill centre trajectories. In the first experiment (Figures 7(a) and 7(b)), the circular mill radius c decays slowly over time, always being less than the critical radius c (b) (in Figure 7(b) the green points lie below the red dashed line). Hence we are in the white region of the phase space in Figure 6. Since the worms are moving clockwise, the model predicts that the mill centre should drift anticlockwise, increasing in angular speed as time progresses. This is indeed what we see in Figure 7(a) with the darker later time blue points less clustered together than the lighter earlier time points.
In contrast, we see much more variation in c in the second experiment (Figures 7(c) and 7(d)) with points both above and below c (in Figure 7(d) green points lie either side of the red dashed line). The predicted sign of F y thus oscillates i.e. the model predicts that the net drift of the mill centre should be minimal. This is indeed what we see in Figure 7(c) with light and dark blue points equally scattered.

Binary Circular Mill Systems
During the evolution of the system, multiple mills can emerge at the same time ( Figure 8 and Supplementary Video 2). This is to be expected since the worms can only interact locally with each other and hence can not coordinate globally to produce a single mill. Here for simplicity we will only consider the most common example of this phenomenon, namely a pair of circular mills. Since the radii of the mills is of the same order of magnitude as the distance between them, a perturbation expansion in terms of c is not possible. Hence the method utilised in §4 can not yield an analytic solution here.
However, using the insight revealed from §4 regarding the flow field produced by an individual mill, we can explain the experimentally observed behaviour from a fluid dynamical viewpoint. In particular, we can explain both the location where the second mill forms and the direction in which it rotates and predict the stability of the binary system. Experimentally, we observed a total of nine binary circular mill systems (summarised in appendix A). Figure 9(a-c) gives a snapshot from three of these experiments. Using the theoretical model, streamlines for the flow produced by each of the two mills if they existed in isolation were generated and superimposed on the same plot (Figures 9(d-f)). The first important observation is that secondary mills only appear when the flow produced by the first mill has a stagnation point, forming in the corresponding stagnation point region. All nine observed binary circular mill systems obey this hypothesis while all observed circular mills which do not generate a stagnation point are stable to the emergence of secondary mills. Note that this relation is not a one-to-one correspondence between having stagnation points and secondary mills emerging. Many other factors can prevent secondary mills forming e.g. a low density of worms swimming in the stagnation point region.
Furthermore, the worms forming the secondary mill tend to swim in the direction of the flow around the stagnation point i.e. the two mills tend to rotate in the same angular direction. In seven of the nine binary circular mill systems examined, the mills rotate in the same direction while the second mill in one of the other systems is seeded by a single flotilla of worms who were tracking around the edge of the Petri dish.
Finally, we can gain a qualitative understanding of the stability of the binary system from looking at the streamlines produced by the second mill. If these streamlines do not have a stagnation point in the vicinity of the first mill (Figures 9(a) and 9(d)), the system is unstable as the first mill breaks up. Alternatively if a stagnation point exists and aligns with the first mill (Figures 9(b) and 9(e)), the system will be stable. Figures 9(c) and 9(f) show the intermediate regime where the first mill is partly (but not fully) inside the second mill's stagnation region. The system is unstable over a much longer time scale. In this particular case, since the first mill is much larger than the second mill, it dominates and the second mill breaks up.

Milling Conclusions
Vortex motions in animal groups have been studied for over a century in many animal species. In this paper, we have demonstrated for the very first time that in order to understand these behaviours in aquatic environments, of which the circular milling of S. roscoffensis is a prime example, one has to understand the underlying fluid dynamics of the system. From the drift of the vortex centre to the formation of secondary vortices and their subsequent stability, it is fundamentally the fluid flow processes that drive these mesmerising and constantly evolving structures. This fluid velocity field may allow nutrient circulation as well as providing an efficient method of dispersal of waste products away from the main body of worms. Furthermore, it exerts a force on the circular mill which causes the mill to slowly drift. In particular, for a single mill in a circular arena, the centre of the mill drifts on a circle whose centre is the middle of the arena.
We present a simple model for the system, (a rigid disc rotating in a Stokes flow), parametrised by only two key variables; namely the radius of the mill c and the distance to the centre of the arena b. This fits the experimental results well, both in terms of the mill centre drift direction but also the predicted streamlines. Utilising this understanding, we are able to shed light on the fluid dynamical stability of circular mills. Secondary circular mills form around stagnation points of the flow. The resulting system evolves to one of two kinds of stable states; namely a single mill with no nearby stagnation points or a set of linked mills where each mill centre is located in the stagnation region of another mills. Although in real life the geometry of the arena is more complicated than our circular model, the same principle remains, namely that stagnation points of the flow occur near a mill when that mill is close to a boundary. This allows the worm population to passively organise towards the arena centre without needing to know the exact extent of the domain. Typically the arena centre will be less shaded and more resource rich. A next step is to estimate the speed of drift of the mill centre. As each worm secretes a layer of mucus around itself, creating a non-Newtonian boundary layer between the mill and the bottom of the Petri dish, both the thickness of this boundary layer and the mechanical properties of the mucus need to be quantified before one is able to calculate this drift. A second line of enquiry results from the fact that a dense core of stationary worms is often experimentally observed to form in the centre of a mill. This core can be unstable and break up or it can take over the whole mill, forming a biofilm. At a more microscopic level, it is of interest to examine the extent to which the formation and breakup of mills can be captured by the kinds of continuum models that have been used successfully to study collective behaviour in bacterial systems (Saintillan & Shelley 2008), where vortex formation is now well-established (Wioland et al. 2013).
Funding. This work was supported in part by the Engineering and Physical Sciences Research Council, through a doctoral training fellowship (GTF) and an Established Career Fellowship EP/M017982/1 (REG), an ERC Consolidator grant 682754 (EL), and by the Schlumberger Chair Fund (REG).
Declaration of interests. The authors report no conflict of interest.
Author ORCID. G.T. Fortune, https://orcid.org/0000-0003-0817-9271; A. Worley, https://orcid.org/0000-0002-7734-7841; A.B. Sendova-Franks, https://orcid.org/0000-0001-9300-6986; N.R. Franks, https://orcid.org/0000-0001-8139-9604; K.C. Leptos, https://orcid.org/0000-0001-9438-0099; E. Lauga, https://orcid.org/0000-0002-8916-2545; R.E. Goldstein, https://orcid.org/0000-0003-2645-0598 Appendix A. Experimental Data on Binary Circular Mill Systems Table 1 gives collocated experimental data of the evolution of circular mills across eighteen distinct experiments. The experimental net drift was obtained by plotting the angle between the line through the centres of the mill and Petri dish and a fixed reference line as a function of time. Linearly interpolating this data, if the magnitude of the gradient of the plotted line is greater than 1.05 rad/h (i.e. changes by more than 10 • during a experiment of typical duration 10 minutes), then we can definitively say that there is a net drift i.e. clockwise if the gradient of the line is negative or anticlockwise if the gradient of the line is positive. Otherwise we write none, since there is no observable net drift within the bounds of experimental error. Similarly, the net drift predicted by net drift was obtained by plotting c − c (b) as a function of time. If the mean of these data points is greater than one standard deviation, then we predict that the mill should drift clockwise. If the mean is less than zero but has magnitude greater than on standard deviation, we predict that the mill should drift anticlockwise. Otherwise, we predict that there is should be no observable net drift within the bounds of experimental error.
While the radius of a circular mill c varies considerably throughout its evolution, its centre remains within experimental error at a constant distance b from the centre of the Petri dish. For fourteen of the eighteen mills, the predicted net direction of drift of the mill centre from the model (assuming that the mill centre drifts in the direction of the force that the flow imposes onto it) matches with the actual net direction. The discrepancy in the other four experiments arises from inertial effects, which particularly come into play for circular mills close to the centre of the Petri dish (experiments 8 and 9) where F from (5.13) is small. Table  2 gives the corresponding data for nine distinct binary circular mill systems.  Table 2: Radius bi, distance from the arena centre ci, and orientation data from the evolution of each of the two mills (i ∈ {1, 2}) in nine experimentally observed binary circular mill systems. CW is clockwise while ACW is anticlockwise.