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

Abstract

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 experimentally and theoretically, 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. We analyse a series of models 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.

1. 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, Franks & Worley 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, Hori & Chihara 1980). The photosynthetic activity of the algae in hospite under the worm epidermis provides the nutrients required to sustain the host. The worms propel themselves through the metachronal beating of many cilia (Bailly et al. 2014). They 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 1b). 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).

Figure 1. The plant-animal worm Symsagittifera roscoffensis. (a) Magnified view of adult. (b) S. roscoffensis in situ on the beach.

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, Kruse & Howard 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 leading-order 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 {({\rm \pi} z/2H)}$, where H is the depth of fluid in the dish. 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 {({\rm \pi} z/2H)}$, deriving a Brinkman-like equation for the vertically averaged velocity flow field. In general, further analytic progress cannot 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, consistent with reversibility arguments for Stokes flow, 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 orbit on a circle with centre at 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.

2. 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^\circ 29^\prime 45.3^{\prime \prime }$N $2^\circ 33^\prime 14.3^{\prime \prime }$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. Approximately 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$ f.p.s. using a Canon Eos 5d Mark II camera equipped with a Canon macro lens MP-E $65$ mm f/$2.8$ providing a 1–5× magnification 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.

Figure 2. Field experiments. (a) Set-up used to film milling behaviour in Guernsey. (b) Montage of still images capturing streaklines produced by the flow.

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 anti-clockwise 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 are available at https://doi.org/10.1017/jfm.2020.1112.

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 set-up 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 orbit of a mill centre within a Petri dish and the formation of binary mill systems.

3. 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 in a fluid of kinematic viscosity $\nu$ is given by the swimming speed $U$, length $\ell$ and diameter $\rho$ as $U(\ell \rho ^2)^{1/3}/2\nu \approx 0.36$, so worms swim in the viscous-dominated 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 ({\rm \pi} z/2H)$.

3.1. Locomotion of an individual worm

The individuals of S. roscoffensis studied in the present experiments have a broad distribution of sizes; their length $\ell$ is in the range $\approx 1.5\text {--}6$ mm, with mean $\bar \ell =3$ mm, and diameters $\rho$ falling in the range $\approx 0.2\text {--}0.6$ mm, with mean $\bar \rho \approx 0.35$ mm. Worm locomotion arises from the collective action of carpets of cilia over the entire body surface, each $\approx 20\ \mathrm {\mu }$m long, beating at $\approx 50$ Hz. Muscles within the organism allow it to bend, and thereby alter its swimming direction (Bailly et al. 2014). In unbounded fluid, the average swimming speed of individuals is $U\approx 2$ mm s$^{-1}$.

To model the fluid velocity field produced by a worm, we follow (Pöhnl, Popescu & Uspal 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 $\boldsymbol {U} = U\boldsymbol {e_z}$ so the $z$-axis lies along the major axis. The squirmer moves through prescribing a tangential slip velocity $\boldsymbol {u_s}$ at its surface $\boldsymbol{\varSigma}$. Neglecting inertia, the fluid flow in the swimmer frame $\boldsymbol {u}$ satisfies

(3.1)\begin{equation} \mu \nabla^2 \boldsymbol{u} = \boldsymbol{\nabla} p, \end{equation}

with boundary conditions

(3.2a,b)\begin{equation} \boldsymbol{u} |_{\boldsymbol{\varSigma}} = \boldsymbol{u_s}, \quad \boldsymbol{u}(|r| \rightarrow \infty ) \rightarrow - \boldsymbol{U}, \end{equation}

together with the force-free condition

(3.3)\begin{equation} \int_{\varSigma} \boldsymbol{\sigma} \boldsymbol{\cdot} \boldsymbol{n} \, \textrm{d} \boldsymbol{\varSigma} = 0, \end{equation}

where $\boldsymbol {\sigma }$ is the stress tensor, $p$ is the pressure, and $\mu$ is the dynamic viscosity of water. We now switch to the modified prolate spheroidal coordinates $(\tau , \xi , \varphi )$, utilising the transformations

(3.4a)\begin{gather} \tau = \frac{1}{2c}( \sqrt{x^2 + y^2 + (z + c)^2} + \sqrt{x^2 + y^2 + (z - c)^2} ), \end{gather}

(3.4b)\begin{gather}\xi = \frac{1}{2c}( \sqrt{x^2 + y^2 + (z+c)^2} - \sqrt{x^2 + y^2 + (z - c)^2} ), \end{gather}

(3.4c)\begin{gather}\varphi = \arctan\left(\frac{y}{x}\right), \end{gather}

with $c = \sqrt {b_z^2 - b_x^2}$ and the squirmer boundary is mapped to the surface $\tau = \tau _0 = b_z/c =$ constant. In this coordinate system, assuming an axisymmetric flow $\boldsymbol {u} = u_{\tau }\boldsymbol {e_{\tau }} + u_{\xi }\boldsymbol {e_{\xi }}$ and axisymmetric tangential slip velocity $\boldsymbol {u_s} = u_s\boldsymbol {e_{\xi }}$, the Stokes stream function $\psi$ satisfies

(3.5a)\begin{gather} u_{\tau} = \frac{1}{h_{\xi}h_{\varphi}}\frac{\partial \psi}{\partial \xi} = \frac{1}{c^2 \sqrt{(\tau^2 - \xi^2)(\tau^2 - 1)}}\frac{\partial \psi}{\partial \xi}, \end{gather}

(3.5b)\begin{gather}u_{\xi} = -\frac{1}{h_{\tau}h_{\varphi}}\frac{\partial \psi}{\partial \tau} = - \frac{1}{c^2 \sqrt{(1 - \xi^2)(\tau^2 - \xi^2)}}\frac{\partial \psi}{\partial \tau}. \end{gather}

Taking from Dassios, Hadjinicolaou & Payatakes (1994) the general separable solution for the stream function in prolate spheroidal coordinates and applying the boundary conditions given in (3.2a,b) and (3.3), as in Pöhnl et al. (2020), we obtain

(3.6)\begin{equation} \psi(\tau,\xi) = \sum_{n = 2}^{\infty} g_n(\tau)G_n(\xi), \end{equation}

where the $g_n$ satisfy

(3.7a)\begin{gather} g_2(\tau) = C_4H_4(\tau) + D_2H_2(\tau) - 2c^2U G_2(\tau), \end{gather}

(3.7b)\begin{gather}g_3(\tau) = -\frac{C_3}{90}G_0(\tau) + C_5H_5(\tau) + D_3H_3(\tau), \end{gather}

(3.7c)\begin{gather}g_{n \geq 4}(\tau) = C_{n + 2}H_{n + 2}(\tau) + C_nH_{n - 2}(\tau) + D_nH_n(\tau), \end{gather}

where $G_n$ and $H_n$ are Gegenbauer functions of the first and second kind respectively and $P^1_n = -\sqrt {1 - x^2}\, \textrm {d}P_n/{\textrm {d} x}$ are the associated Legendre polynomials. The integration constants $C_n$ and $D_n$ are set by the boundary conditions

(3.8)\begin{equation} g_n(\tau_0) = 0 \mbox{ and } \left. \frac{\textrm{d} g_n}{\textrm{d}\tau}\right|_{\tau = \tau_0} = \tau_0 c^2 n(n-1)B_{n-1} \quad \mbox{for } n = 2,3 \ldots . \end{equation}

Here, the $\{B_n\}_{n \geq 1}$ are the coefficients in the series expansion of $u_s = \tau _0 \sum _{n \geq 1} B_n V_n(\xi )$ using the set of functions $\{ V_n = (\tau _0^2 - \xi ^2)^{-1/2}P_n^1(\xi ) \}_{n \geq 1}$, which is a basis over the space of $L^2$ functions satisfying $f(\xi = \pm 1) = 0$ together with the inner product $\langle \rangle _w = \int ^1_{-1} w \, \textrm {d} \xi$ and the weight function $w = \tau _0^2 - \xi ^2$. Furthermore, the swimming speed $U$ can be expressed in terms of the $\{B_n\}_{n \geq 1}$ using

(3.9)\begin{equation} U = -\frac{\tau_0}{2}\int^1_{-1} \frac{\sqrt{1 - \xi^2}}{\sqrt{\tau_0^2 - \xi^2}} v_s(\xi) \, \textrm{d} \xi = \frac{\tau_0^2}{2} \sum_{n \mbox{ odd}} B_nU_n \quad \mbox{where } U_n = \int^1_{-1} \frac{P_1^1 P_n^1}{\tau_0^2 - x^2}\, {\textrm{d} x}, \end{equation}

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 solely a linear combination of the odd modes i.e. $B_{2n+2} = 0 \rightarrow g_{2n+1} = 0 \rightarrow C_{2n+1} = D_{2n+1} = 0 \, \forall n \in \mathbb {Z}^{+}$. When the prescribed forcing arises purely from the first mode, i.e. $u_s = \tau _0 B_1V_1$, $C_{2n \colon n \geq 1}$ and $D_{2n \colon n \geq 1}$ simplify to become

(3.10a,b)\begin{equation} D_2 = -2B_1c^2\tau_0(\tau_0^2 - 1) \quad \mbox{and} \quad C_{2n \colon n \geq 1} = D_{2n \colon n \geq 2} = 0. \end{equation}

When the forcing arises from a higher-order mode, i.e. $u_s = \tau _0 B_{2n+1}V_{2n+1}$ where $n \geq 1$, $D_2$ and $C_4$ simplify to become

(3.11)\begin{align} \left( \begin{array}{ll} C_4 \\ D_2 \end{array} \right) &= \frac{B_{2n+1}c^2U}{H_4(\tau_0)H_2'(\tau_0)-H_2(\tau_0)H_4'(\tau_0)} \nonumber\\ &\quad \times \left( \begin{array}{ll} 1 \\ \dfrac{2}{3} + \dfrac{5\tau_0^4}{4} - \dfrac{25\tau_0^2}{12} - \dfrac{5\tau_0}{8}(1 - \tau_0^2)^2\log{\left(\dfrac{\tau_0+1}{\tau_0-1}\right)} \end{array} \right), \end{align}

i.e. the only $n$ dependence arises from the $U$. Moving into the laboratory frame, in the far field ($\tau \gg 1$) the dominant term in the expansion for $\psi$ comes from $H_{2}(\tau ) = 1/3\tau + \cdots$, so

(3.12a)\begin{align} \psi &= \frac{1}{3\tau} \left( \frac{D_2}{2}(1-\xi^2) - \frac{C_4}{8}(1 - 6\xi^2 + 5\xi^4) \right) + \cdots, \end{align}

(3.12b)\begin{align} u_{\xi} &= -\frac{1}{\tau c^2 \sqrt{1 - \xi^2}}\frac{\partial \psi}{\partial r} + \cdots \nonumber\\ &= \frac{1}{3\tau^3c^2\sqrt{1 - \xi^2}}\left( \frac{D_2}{2}(1-\xi^2) - \frac{C_4}{8}(1 - 6\xi^2 + 5\xi^4) \right) + \cdots, \end{align}

(3.12c)\begin{align} u_{\tau} &= \frac{1}{\tau^2c^2}\frac{\partial \psi}{\partial \xi} + \cdots = -\frac{\xi}{3\tau^3c^2}\left( D_2 + \frac{C_4}{2}(5\xi^2 - 3) \right) + \cdots. \end{align}

Converting this back to vector notation yields

(3.13)\begin{equation} \boldsymbol{u} = -\frac{1}{2c^2}\left[ \left(D_2 + \frac{C_4}{6}\right)\boldsymbol{u_D} + \frac{5C_4}{36}\boldsymbol{u_Q} \right] + \cdots ,\end{equation}

where $\boldsymbol {u_D}$ and $\boldsymbol {u_Q}$, the flows generated by a source dipole and a Stokes quadrupole respectively, satisfy

(3.14a)\begin{gather} {u_{Di}} = \left( \frac{c}{r} \right)^3\left( \frac{x_ix_3}{r^2} - \frac{\delta_{i3}}{3} \right), \end{gather}

(3.14b)\begin{gather} {u_{Qi}} = \frac{\partial}{\partial x_3^2}\left( \frac{x_ix_3}{r^3} + \frac{\delta_{i3}}{r} \right) = 3\left(\frac{c}{r}\right)^3\left(\frac{5x_ix_3^2}{r^2} - (2 + \delta_{i3})\frac{x_ix_3}{r^2} - \left( \frac{x_ix_3}{r^2} - \frac{\delta_{i3}}{3} \right) \right). \end{gather}

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 $\tau _0 \rightarrow 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 $\tilde {a} = c^3\tau _0(\tau _0^2 -1)$. When $\tau _0 \rightarrow \infty$, namely the spherical limit, as expected $\tilde {a} \rightarrow b_x = b_z = a$, the radius of the sphere. When $\tau _0 \rightarrow 1$, the elongated limit, $\tilde {a} \rightarrow ( b_z b_x^2 )^{1/3}$. Thus, at the scale of an individual worm, the Reynolds number $U \tilde {a}/\nu$ in water ($\nu =1\ \textrm {mm}^2\ \textrm {s}^{-1}$) is $\approx 0.36$ where $\tilde {a} = ( b_z b_x^2 )^{1/3} = (\bar \ell \bar \rho ^2)^{1/3}/2$ is the correct length scale for locomotion of an individual worm. Inertial effects are modest and individual worms swim in the viscous-dominated regime.

3.2. Ring of spheroidal squirmers

Given the results 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 {\theta }, c \sin {\theta },0)$ and orientated in the direction $\boldsymbol {p} = (\sin {\theta }, -\cos {\theta },0)$, utilising a Cartesian coordinate system $(x,y,z)$ with origin at the centre of the circle. If each squirmer generates a source dipole $\boldsymbol {u_{D}}$, the fluid velocity $\boldsymbol {u_{DX}}(c,\theta )$ at $\boldsymbol{X} = (R,0,0)$ is

(3.15)\begin{equation} \boldsymbol{u_{DX}} = \frac{c^3}{3r^5}\left( \begin{array}{l} (2R^2-c^2)\sin{\theta} - cR\sin{\theta}\cos{\theta} \\ (R^2 + c^2)\cos{\theta} - cR(3 - \cos^2{\theta}) \end{array} \right), \end{equation}

where $r = (R^2 + c^2 - 2cR\cos {\theta })^{1/2}$. The total velocity $\boldsymbol {u_{D\kern0.05em ring}}\,(R)$ at $\boldsymbol{X}$ due to a ring of clockwise swimming worms of radius $c$ with line density $\lambda _{ring}$ is then

(3.16)\begin{equation} \boldsymbol{u_{D\kern0.05em ring}}\,(R) = c \lambda_{ring} \int^{\pi}_{-\pi} \boldsymbol{u_{DX}}(c,\theta) \, \textrm{d} \theta = 0. \end{equation}

Hence, a ring of uniformly distributed source dipole swimmers generates no net flow field outside of the ring. By contrast, the flow field $\boldsymbol {u_{Q\kern0.05em ring}}\,(R)$ due to a ring of clockwise swimming worms, each generating a Stokes quadrupole, is finite:

(3.17)\begin{equation} \boldsymbol{u_{Q\kern0.05em ring}}\, (R) = c\lambda_{ring}\left( \frac{c}{R} \right)^3 \varUpsilon\left( \frac{c}{R} \right)\boldsymbol{e_y},\end{equation}

where

(3.18)\begin{align} \varUpsilon(x) &= \int^{2\pi}_0 \frac{\cos{\theta}( 4-23x^2-x^4 +\cos^2{\theta}(17x^2-5))}{(1 + x^2 - 2x\cos{\theta})^{7/2}} \, \textrm{d} \theta \nonumber\\ & \quad - \int^{2\pi}_0 \frac{x( 12 - 13x^2 + \cos^2{\theta}( 9x^2-31 ) + 15\cos^4{\theta} )}{(1 + x^2 - 2x\cos{\theta})^{7/2}} \, \textrm{d} \theta. \end{align}

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 $\varphi$ direction with imposed surface flow $u_{\varphi 0}(\xi ) \, \boldsymbol {e_{\varphi }}$ in free space. Assuming that the generated fluid flow is purely in the $\varphi$ direction with no $\varphi$ dependence, the $\varphi$ component of the Stokes equations, $\mu (\nabla ^2 \boldsymbol {u})_{\varphi } = 0$, becomes

(3.19)\begin{equation} (\tau^2 - 1)\frac{\partial^2 u_{\varphi}}{\partial \tau^2} + 2\tau\frac{\partial u_{\varphi}}{\partial \tau} - \frac{u_{\varphi}}{\tau^2 - 1} + (1 - \xi^2)\frac{\partial^2 u_{\varphi}}{\partial \xi^2} - 2\xi\frac{\partial u_{\varphi}}{\partial \xi} - \frac{u_{\varphi}}{1 - \xi^2} = 0. \end{equation}

This admits the general separable solution that tends to zero at infinity

(3.20)\begin{equation} u_{\varphi} = \sum_{n = 1}^{\infty} C_{pn} P^1_{n}(\xi) Q^1_{n}(\tau), \end{equation}

where $C_{pn}$ are constants and as before $P^1_n(\xi )$ and $Q^1_n(\tau )$ are associated Legendre polynomials. Furthermore, since the squirmer is force and torque free, $C_{p1} = 0$. Decomposing $u_{\varphi 0}$ using the basis $\{ V_n(\xi ) \}$, i.e. $u_{\varphi 0} = \sum _{n = 2}^{\infty } C_{n0} V_n(\xi )$, we find that $C_{pn} = C_{n0}/W_{pn}(\tau _0)$ where $\tau _0 = 1/\sqrt {1 - (r_e)^2}$ and $r_e = b_x/b_z \leq 1$ is the aspect ratio of the spheroid. Note that, in the spherical limit ($r_e \rightarrow 1$, $b_x = b_z = a$) (3.20) simplifies to become

(3.21)\begin{equation} u_{\varphi} = \sum_{n = 1} a \bar{C}_n \frac{a}{r^{n+1}} V_n(\xi), \end{equation}

where $\bar {C}_n$ are constants, and we recover the general form for a spherical squirmer with swirl (Pak & Lauga 2014; Pedley, Brumley & Goldstein 2016).

Returning to the general case, the dominant term in the far field arises from mode $2$,

(3.22)\begin{equation} u_{\varphi} = -\frac{2C_{p2}}{5\tau^3}\xi\sqrt{1-\xi^2} + \cdots \rightarrow \boldsymbol{u} = \frac{2C_{p2}}{5}\left( \frac{c}{r} \right)^3\frac{\boldsymbol{x} \times \boldsymbol{x_k}}{r^2}, \end{equation}

where $\boldsymbol{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 $(\tau ,\xi ,\varphi )$ to the oblate spheroidal coordinate system $(\lambda ,\xi ,\varphi )$ using the substitutions

(3.23a,b)\begin{equation} \tau = i\lambda, \quad c = -i\bar{c}, \end{equation}

where $\bar {c} = \sqrt {b_x^2 - b_z^2}$ and $(\lambda ,\xi ,\varphi )$ can be expressed in terms of the Cartesian coordinates $(x,y,z)$ using

(3.24a)\begin{gather} \lambda = \frac{1}{2\bar{c}}( \sqrt{x^2 + y^2 + (z - i\bar{c})^2} + \sqrt{x^2 + y^2 + (z + i\bar{c})^2} ), \end{gather}

(3.24b)\begin{gather}\xi = \frac{1}{-2i\bar{c}}( \sqrt{x^2 + y^2 + (z - i\bar{c})^2} - \sqrt{x^2 + y^2 + (z + i\bar{c})^2} ), \end{gather}

(3.24c)\begin{gather}\varphi = \arctan{\left( \frac{y}{x} \right)}. \end{gather}

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,

(3.25)\begin{equation} \boldsymbol{u} = \frac{2C_{ob2}}{5}\left( \frac{\bar{c}}{r} \right)^3\frac{\boldsymbol{x} \times \boldsymbol{x_k}}{r^2}. \end{equation}

This result is intuitive; in the absence of a net torque on the object there cannot 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 non-vanishing 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, Lauga & Goldstein 2021) 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

(3.26)\begin{equation} \boldsymbol{u_p} = \frac{2{\rm \pi}^2 C_{ob2}}{15}\left( \frac{\bar{c}}{H} \right)^3\frac{\boldsymbol{x} \times \hat{\boldsymbol{k} }}{\rho}\cos{\left( \frac{{\rm \pi} h}{2H} \right)}\sin{\left( \frac{{\rm \pi} z}{2H} \right)}K_1{\left( \frac{\rho {\rm \pi}}{2H} \right)}, \end{equation}

where $\rho ^2 = x^2 + y^2$ and $K_1$ is a modified Bessel function of the second kind. The flow decays exponentially away from the squirmer with decay length $2H/{\rm \pi}$.

4. Mathematical model for a circular mill

4.1. Background

We now proceed to develop a model for 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\text {--}20$ mm and rotates roughly as a rigid body with period $T = 2{\rm \pi} c/U$ in the range $15\text {--}60$ s and angular frequency $\omega = U/c$ in the range $0.4\text {--}0.1$ s$^{-1}$ whereas in § 3.1 the average worm swimming speed $U \approx 2$ mm s$^{-1}$. 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 a 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 {({\rm \pi} z/2H)}$. Hence, we then vertically average the governing equations by setting the $z$ dependence of $u$ to be precisely $\sin {({\rm \pi} z/2H)}$ i.e. we employ the factorisation

(4.1)\begin{equation} \boldsymbol{u} = \frac{\rm \pi}{2}\sin{\left(\frac{{\rm \pi} z}{2H}\right)} \boldsymbol{U}(r), \end{equation}

where $\boldsymbol {U} = \boldsymbol {U(r)}$ is independent of $z$ and the factor ${\rm \pi} /2$ is for convenience.

Thus, a suitable rotational Reynolds number on the scale of a mill is $Re \sim UX/\nu \approx 10$ where $X = 2H/{\rm \pi}$. Moreover, the dominant velocities are azimuthal, with gradients in the radial direction. This suggests that the fluid dynamics of milling is certainly in the viscous-dominated regime, and the neglect of inertial terms is justified.

4.2. Defining notation

As shown in figure 3 and in supplementary movie 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 $-\varOmega$, as a rigid disc of radius $c$ and height $d$ with imposed angular velocity $\varOmega$, generating a flow in a cylindrical domain with cross-sectional radius 1 where $d_0 \ll d \ll c , b , H$. Let the centre of the disc M have instantaneous position $(-b,0,d_0 + d/2)$.

Figure 3. System containing a single mill. (a) Experimental view. (b,c) Plan and front views for the corresponding schematic showing a disc, rotating with angular velocity $\varOmega$, which has radius $c$, thickness $d$ and centre $M$ a distance $b$ away from the centre $P$ of a circular Petri dish of unit radius.

4.3. 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 ventral surfaces of worms have beating cilia. This artifice 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) \colon 0 \leq z \leq d_0;, x^2 + y^2 \leq c^2 \}$, let the imposed slip velocity at $z = d_0$ be $\boldsymbol {u_{slip}} = u_{slip}(r)\hat {\boldsymbol {e}_{\theta }}$. In the absence of any pressure gradients, the general results of lubrication theory dictate a linear velocity profile for the flow $\boldsymbol {u_b}$ in the film,

(4.2)\begin{equation} \boldsymbol{u_b}(r,z) = \frac{z}{d_0}( u_{slip} + \varOmega r )\hat{\boldsymbol{e}_{\theta}}, \end{equation}

where $\varOmega$ is the as yet unknown angular velocity of the disc. The flow in the region above the mill is simply $u(r,z) = \varOmega r$, independent of $z$; it rotates with the disc as a rigid body. The torque $T_b$ on the underside of the mill is

(4.3)\begin{equation} T_b = \frac{2{\rm \pi}\mu}{d_0} \int^c_0 r^2 \, \textrm{d} r (u_{slip} + \varOmega r), \end{equation}

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

(4.4)\begin{equation} \varOmega = - \frac{4}{c^4} \int^c_0 r^2 \, \textrm{d} r u_{slip}(r). \end{equation}

If $u_{slip}$ has solid-body-like character, $u_{slip} = u_0 r/c$, then $\varOmega = -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$.

4.4. Full governing equations

Assuming Stokes flow with fluid velocity $\boldsymbol {u} = (u_x, \, u_y, \, w = 0)$, pressure $p = p(x,y,z)$ and viscosity $\mu$, the governing equations are

(4.5a,b)\begin{equation} \mu \nabla^2 \boldsymbol{u} = \boldsymbol{\nabla}p ,\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0. \end{equation}

Employing no-slip boundary conditions (Batchelor 1967) at both the outer edge and the bottom of the Petri dish yields

(4.6a,b)\begin{equation} \boldsymbol{u} = 0 \mbox{ at } z = 0 , \quad \boldsymbol{u} = 0 \mbox{ at } r = \sqrt{x^2 + y^2} =1. \end{equation}

On the surface of the fluid, $z = H$, the dynamic boundary condition is

(4.7)\begin{equation} (\boldsymbol{\sigma}_{\boldsymbol{f}} - \boldsymbol{\sigma}_{\boldsymbol{a}}) \boldsymbol{\cdot} \hat{\boldsymbol{z}} = 0,\end{equation}

where $\boldsymbol {\sigma _f}$ and $\boldsymbol {\sigma _a} = -p_0\boldsymbol{\mathsf{I}}$ are stress tensors for the fluid and the air respectively. Finally, the boundary conditions on the surface of the mill become

(4.8a)\begin{gather} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{e}_{\boldsymbol{t}} = \varOmega \tilde{c} \mbox{ on } \boldsymbol{\varGamma} = \{ (x,y,z) \colon \tilde{c}^2 = (x + b)^2 + y^2 \leq c^2, z = d_0 \}, \end{gather}

(4.8b)\begin{gather}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{e}_{\boldsymbol{t}} = \varOmega c \mbox{ on } \boldsymbol{\varGamma} = \{ (x,y,z) \colon (x + b)^2 + y^2 = c^2, d_0 \leq z \leq d_0 + d \}, \end{gather}

(4.8c)\begin{gather}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{e}_{\boldsymbol{t}} = \varOmega \tilde{c} \mbox{ on } \boldsymbol{\varGamma} = \{ (x,y,z) \colon \tilde{c}^2 = (x + b)^2 + y^2 \leq c^2, z = d_0 + d \}, \end{gather}

where the tangent and normal vectors $\boldsymbol {e_t}$ and $\boldsymbol {e_n}$ satisfy

(4.9a)\begin{gather} \boldsymbol{e_n} = \frac{1}{( y^2 + (x+b)^2 )^{1/2}}((x+b)\boldsymbol{e_x} + y\boldsymbol{e_y}), \end{gather}

(4.9b)\begin{gather}\boldsymbol{e_t} = \frac{1}{( y^2 + (x+b)^2 )^{1/2}}(-y\boldsymbol{e_x} + (x+b)\boldsymbol{e_y}). \end{gather}

4.5. Vertically averaged governing equations

Defining $\boldsymbol {r}$ as the in-plane coordinates $(x,y)$, as set out in § 4.1, we employ the factorisation

(4.10)\begin{equation} \boldsymbol{u} = \frac{\rm \pi}{2}\sin{\left( \frac{{\rm \pi} z}{2H} \right)}\boldsymbol{U}(\boldsymbol{r}) \longrightarrow \nabla^2 \boldsymbol{u} = \frac{\rm \pi}{2}\sin{\left( \frac{{\rm \pi} z}{2H} \right)}(\nabla^2 - \kappa^2 ) \boldsymbol{U}, \end{equation}

where $\kappa ={\rm \pi} /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}\int ^{H}_{0}\cdots \, \textrm {d} z$, gives the Brinkman-like equation

(4.11a,b)\begin{equation} \mu ( \nabla^2 \boldsymbol - \kappa^2 )\boldsymbol{U} = \boldsymbol{\nabla}p , \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0, \end{equation}

where $\boldsymbol {U} = U_x \boldsymbol {e_x} + U_y \boldsymbol {e_y} = U_n \boldsymbol {e_n} + U_t \boldsymbol {e_t}$ has a corresponding Stokes streamfunction $\varphi$ satisfying

(4.12)\begin{equation} \left \{ U_x = \frac{\partial \varphi}{\partial y}, U_y = - \frac{\partial \varphi}{\partial x} \right \} \longrightarrow \nabla^4 \varphi = \kappa^{2} \nabla^2 \varphi, \end{equation}

together with boundary conditions

(4.13a)\begin{gather} u_n = 0 \mbox{ at } r = \sqrt{x^2 + y^2} =1, \end{gather}

(4.13b)\begin{gather}u_t = 0 \mbox{ at } r = \sqrt{x^2 + y^2} =1, \end{gather}

(4.13c)\begin{gather}u_n = 0 \mbox{ on } \boldsymbol{\varGamma} = \{ (x,y) \colon (x+b)^2 + y^2 = c^2 \}, \end{gather}

(4.13d)\begin{gather}u_t = \varOmega c \mbox{ on } \boldsymbol{\varGamma} = \{ (x,y) \colon (x+b)^2 + y^2 = c^2 \}, \end{gather}

i.e. no-slip is imposed at the edge of the Petri dish while a constant azimuthal velocity boundary condition is imposed at $\boldsymbol{\varGamma}=\{ (x,y) \colon (x+b)^2 + y^2 = c^2 \}$. Using separation of variables, the general series solution in polar coordinates to 4.12 is

(4.14)\begin{align} \varphi &= A_0 + B_0 \log{r} + C_0 I_0(\kappa r) + D_0 K_0(\kappa r) \nonumber\\ & \quad + \sum_{n = 1}^{\infty} \cos{n \theta}( A_n r^n + B_n r^{-n} + C_n I_n(\kappa r) + D_n K_n(\kappa r) ) \nonumber\\ & \quad + \sum_{n = 1}^{\infty} \sin{n \theta}( \tilde{A}_n r^n + \tilde{B}_n r^{-n} + \tilde{C}_n I_n(\kappa r) + \tilde{D}_n K_n(\kappa r) ), \end{align}

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.

4.6. Near-field perturbation analysis

Motivated by perturbative studies of screened electrostatics near wavy boundaries (Goldstein, Pesci & Romero-Rochin 1990), we consider a small perturbation of the mill centre away from the middle of the Petri dish, i.e. $b = c\epsilon$ where $\epsilon \ll 1$. Expanding in powers of $\epsilon$, i.e. for a given function $f$ considering $f = f^0 + \epsilon f^1 + \epsilon ^2 f^2 + \ldots$, (4.12) becomes

(4.15)\begin{equation} \nabla^4 \varphi^{i} = \kappa^2 \nabla^2 \varphi^{i} \quad\mbox{for } i = 1,2,3 \ldots \end{equation}

with corresponding boundary conditions at the edge of the Petri dish

(4.16)\begin{equation} \frac{\partial \varphi^i}{\partial r} = \frac{\partial \varphi^i}{\partial \theta} = 0 \quad \mbox{for } i = 1,2,3 \ldots. \end{equation}

Furthermore, since $\boldsymbol{\varGamma}$ can be expressed in polar coordinates as

(4.17)\begin{align} \boldsymbol{\varGamma} &= \left\{ (r,\theta) \colon r = R(\theta) = -b \cos{\theta} + \sqrt{c^2 - b^2 \sin^2{\theta}} \vphantom{\frac{c \sin^2{\theta}}{2}} \right.\nonumber\\ &= \left.c - \epsilon c \cos{\theta} - \epsilon^2 \frac{c \sin^2{\theta}}{2} + {O}(\epsilon^3) \right\}, \end{align}

while (4.9a) and (4.9b) expand to become

(4.18)\begin{gather} u_n = \frac{1}{r}\frac{\partial \varphi}{\partial \theta} + \epsilon \frac{c \sin{\theta}}{r}\frac{\partial \varphi}{\partial r} - \frac{\epsilon^2 c^2}{2r^2}\left( \sin{2\theta}\frac{\partial \varphi}{\partial r} + \frac{\sin^2{\theta}}{r}\frac{\partial \varphi}{\partial \theta} \right) + {O}(\epsilon^3), \end{gather}

(4.19)\begin{gather}u_t = -\frac{\partial \varphi}{\partial r} + \epsilon \frac{c \sin{\theta}}{r^2}\frac{\partial \varphi}{\partial \theta} + \frac{\epsilon^2 c^2}{2 r^2}\left( \sin^2{\theta}\frac{\partial \varphi}{\partial r} - \frac{\sin{2 \theta}}{r}\frac{\partial \varphi}{\partial \theta} \right) + {O}(\epsilon^3), \end{gather}

(4.13c) and (4.13d) reduce to

(4.20)\begin{align} &\frac{1}{c}\frac{\partial \varphi^0}{\partial \theta} + \epsilon \left( \frac{1}{c}\frac{\partial \varphi^1}{\partial \theta} + \sin{\theta}\frac{\partial \varphi^0}{\partial r} + \frac{\cos{\theta}}{c}\frac{\partial \varphi^0}{\partial \theta} - \cos{\theta}\frac{\partial^2 \varphi^0}{\partial r \partial \theta} \right) \nonumber\\ & \quad + \epsilon^2 \left( \frac{1}{c}\frac{\partial \varphi^2}{\partial \theta} + \sin{\theta} \frac{\partial \varphi^1}{\partial r} + \frac{\cos{\theta}}{c}\frac{\partial \varphi^1}{\partial \theta} - \cos{\theta}\frac{\partial^2 \varphi^1}{\partial r \partial \theta} + \frac{ \cos^2{\theta}}{c}\frac{\partial \varphi^0}{\partial \theta} \right. \nonumber\\ & \quad - \left.\frac{1 + \cos^2{\theta}}{2}\frac{\partial^2 \varphi^0}{\partial r \partial \theta} - \frac{c \sin{2 \theta}}{2}\frac{\partial^2 \varphi^0}{\partial r^2} + \frac{c \cos^2{\theta}}{2}\frac{\partial^3 \varphi^0}{\partial r^2 \partial \theta} \right) \nonumber\\ & \quad + {O}(\epsilon^3) |_{r = c} = 0. \end{align}

(4.21)\begin{align} &-\frac{\partial \varphi^0}{\partial r}+ \epsilon \left( -\frac{\partial \varphi^1}{\partial r} + \frac{\sin{\theta}}{c}\frac{\partial \varphi^0}{\partial \theta} + c \cos{\theta}\frac{\partial^2 \varphi^0}{\partial r^2} \right) \nonumber\\ & \quad + \epsilon^2 \left( -\frac{\partial \varphi^2}{\partial r} + \frac{\sin{\theta}}{c}\frac{\partial \varphi^1}{\partial \theta} + c \cos{\theta}\frac{\partial^2 \varphi^1}{\partial r^2} + \frac{\sin^2 \theta}{2}\frac{\partial \varphi^0}{\partial r} - \frac{\sin{2 \theta}}{2}\frac{\partial^2 \varphi^0}{\partial r \partial \theta} \right. \nonumber\\ & \quad + \left.\frac{\sin{2\theta}}{2c}\frac{\partial \varphi^0}{\partial \theta} + \frac{c \sin^2{\theta}}{2}\frac{\partial^2 \varphi^0}{\partial r^2} - \frac{c^2 \cos^2{\theta}}{2}\frac{\partial^3 \varphi^0}{\partial r^3} \right) \nonumber\\ & \quad + {O}(\epsilon^3) |_{r = c} = c \varOmega. \end{align}

Hence, at ${O}(1)$, namely when the circular mill is concentric with the Petri dish, we find

(4.22)\begin{gather} \left.\frac{\partial \varphi^0}{\partial r}\right|_{r = 1} = \left.\frac{\partial \varphi^0}{\partial \theta}\right|_{r = 1} = \left.\frac{\partial \varphi^0}{\partial \theta}\right|_{r = c} = 0 , \left. \frac{\partial \varphi^0}{\partial r}\right|_{r = c} = -\varOmega c \Longrightarrow \end{gather}

(4.23)\begin{gather}\varphi^0 = -(\alpha_0 K_0(\kappa r) + \beta_0 I_0(\kappa r)) \Longrightarrow \end{gather}

(4.24)\begin{gather}u_{\theta} = c\varOmega \frac{I_1(\kappa r)K_1(\kappa) - K_1(\kappa r)I_1(\kappa)}{I_1(\kappa c)K_1(\kappa) - K_1(\kappa c)I_1(\kappa)}, \end{gather}

where

(4.25a)\begin{gather} \alpha_0 = \frac{c}{\kappa}\varOmega \frac{I_1(\kappa)}{I_1(\kappa c)K_1(\kappa) - K_1(\kappa c)I_1(\kappa)}, \end{gather}

(4.25b)\begin{gather}\beta_0 = \frac{c}{\kappa}\varOmega \frac{K_1(\kappa)}{I_1(\kappa c)K_1(\kappa) - K_1(\kappa c)I_1(\kappa)}. \end{gather}

For comparison, the corresponding Couette solution is

(4.25c)\begin{equation} u_{\theta} = \frac{c^2 \varOmega}{1 - c^2}\left( \frac{1}{r} - r \right). \end{equation}

Figure 4(a), plots $u_{\theta }(r)$ for both the Brinkman and Couette solutions when $c = 10/45$ and $H = 8/45$ i.e. the Brinkman fluid velocity field decays much faster away from the mill than the Couette fluid velocity field.

Figure 4. (a) Azimuthal fluid velocity profile, plotted as a function of distance from the Petri dish centre $r$, for both the Brinkman solution (black) and the corresponding Couette solution (blue) when $c = 10/45$ and $H = 8/45$. (b) Perturbation fluid velocity profile, plotted as a function of distance from the Petri dish centre $r$, showing both the radial ($u^1_r/\sin {\theta }$, black) and tangential flow ($u^1_{\theta }/\cos {\theta }$, blue) when $c = 10/45$ and $H = 8/45$.

Similarly at ${O}(\epsilon )$, we obtain

(4.26a)\begin{gather} \left. \frac{\partial \varphi^1}{\partial r}\right|_{r = 1} = 0, \quad \left.\frac{\partial \varphi^1}{\partial r}\right|_{r = c} = c \cos{\theta} \left.\frac{\partial^2 \varphi^0}{\partial r^2}\right|_{r = c}, \end{gather}

(4.26b)\begin{gather}\left. \frac{\partial \varphi^1}{\partial \theta}\right|_{r = 1} = 0, \quad \left. \frac{\partial \varphi^1}{\partial \theta}\right|_{r = c} = - c \sin{\theta} \left. \frac{\partial \varphi^0}{\partial r}\right|_{r = c} \Longrightarrow \end{gather}

(4.27)\begin{gather} \varphi^1 = - \cos{\theta}\left( \alpha_1 K_1(\kappa r) + \beta_1 I_1(\kappa r) + \gamma_1 r + \frac{\delta_1}{r} \right), \end{gather}

where $\{ \alpha _1 , \beta _1 ,\gamma _1 , \delta _1 \}$ are known functions of $c$ and $X$ which satisfy the following set of simultaneous equations

(4.28a)\begin{gather} \alpha_1 K_1(\kappa) + \beta_1 I_1(\kappa) + \gamma_1 + \delta_1 = 0, \end{gather}

(4.28b)\begin{gather}\alpha_1 K_1(\kappa c) + \beta_1 I_1(\kappa c) + \gamma_1 c + \frac{\delta_1}{c} = \varOmega c^2, \end{gather}

(4.28c)\begin{gather}-\alpha_1( \kappa K_0(\kappa) + K_1(\kappa) ) + \beta_1(\kappa I_0(\kappa) - I_1(\kappa) ) + \gamma_1 - \delta_1 = 0, \end{gather}

(4.28d)\begin{gather} -\alpha_1\left( \frac{K_0(c/X)}{X} + \frac{K_1(c/X)}{c} \right) + \beta_1\left( \frac{I_0(c/X)}{X} - \frac{I_1(c/X)}{c} \right) + \gamma_1 - \frac{\delta_1}{c^2} \nonumber\\ \quad = \kappa\varOmega c^2\left( -\frac{1}{\kappa c} + \frac{I_0(\kappa c)K_1(\kappa) + K_0(\kappa c)I_1(\kappa)}{I_1(\kappa c)K_1(\kappa) - K_1(\kappa c)I_1(\kappa)} \right). \end{gather}

Figure 4(b) plots $u^1_r/\sin {\theta }$ and $u^1_{\theta }/\cos {\theta }$ as functions of $r$ when $c = 10/45$ and $H = 8/45$. This perturbation flows also decays exponentially away from the mill i.e. the Brinkman term still plays a key role. As will be shown in § 5.1, this perturbation flow leads to the centre of the mill orbiting clockwise on a circle centred on the middle of the Petri dish. That is, the stationary point where the mill and the Petri dish are concentric is unstable.

4.7. Far-field solution

When the circular mill is away from the centre of the Petri dish ($b = {O}(1)$), the boundary conditions at the edge of the mill can no longer be expressed straightforwardly in terms of the polar coordinates $(r,\theta )$. Instead we switch to the bipolar coordinates $(\eta , \xi )$ utilising the transformations

(4.29a,b)\begin{equation} x = a \frac{\sinh{\eta}}{\cosh{\eta} - \cos{\xi}} + d \quad \mbox{and}\quad y = a \frac{\sin{\xi}}{\cosh{\eta} - \cos{\xi}}. \end{equation}

In particular, we can map the outer boundary, $r = 1$, to $\eta = \eta _1$ and the disc boundary to $\eta = \eta _2$ by defining the constants $a,d,\eta _1$ and $\eta _2$ as satisfying

(4.30a–d)\begin{align} d = -a \coth{\eta_1} , \quad b = a(\coth{\eta_1} - \coth{\eta_2}) , \quad 1 = \frac{a}{\sinh{\eta_1}} \quad \mbox{and}\quad c = \frac{a}{\sinh{\eta_2}}, \end{align}

i.e.

(4.31a,b)\begin{gather} \eta_1 = \ln{( a + \sqrt{a^2 + 1} )} , \quad \eta_2 = \ln{\left( \frac{a + \sqrt{a^2 + c^2}}{c} \right)}, \end{gather}

(4.32)\begin{gather}a = \frac{1}{2b}(\sqrt{(1 + c^2 - b^2)^2 - 4c^2}). \end{gather}

In this basis, the system becomes

(4.33a,b)\begin{align} \nabla^2 ( \nabla^2 - \kappa^2 )\varphi = 0 \quad \mbox{where } \nabla^2 = \frac{1}{h^2}\left( \frac{\partial^2}{\partial \xi^2} + \frac{\partial}{\partial \eta^2} \right) \quad \mbox{and}\quad h = \frac{a}{\cosh{\eta} - \cos{\xi}}, \end{align}

with boundary conditions

(4.34a,b)\begin{equation} \frac{\partial \varphi}{\partial \xi} = \left\{ \begin{array}{@{}ll} 0, & \eta = \eta_1 \\ 0, & \eta = \eta_2. \end{array} \right. \quad \mbox{and}\quad \frac{1}{h}\frac{\partial \varphi}{\partial \eta} = \left\{ \begin{array}{@{}ll} 0, & \eta = \eta_1 \\ -\varOmega c, & \eta = \eta_2. \end{array} \right. \end{equation}

Now, in general, this does not admit an analytic solution. However, for large mills away from the Petri dish centre, namely $1/a > \kappa$, the biharmonic term dominates and thus using Melesko & Gomilko (1999), (4.33a,b) reduces to

(4.35)\begin{equation} \nabla^4 \varphi = 0 \longrightarrow \frac{\partial^4 \varPhi}{\partial \xi ^4} + 2 \frac{\partial^4 \varPhi}{\partial \xi^2 \partial \eta^2} + \frac{\partial^4 \varPhi}{\partial \eta^4} + 2 \frac{\partial^2 \varPhi}{\partial \xi^2} - 2 \frac{\partial^2 \varPhi}{\partial \eta^2} + \varPhi = 0 \quad \mbox{where } \varPhi = \frac{\varphi}{h}. \end{equation}

From Kazakova & Petrov (2016), this yields the analytic solution

(4.36)\begin{equation} \varphi = N(\eta) + \frac{M(\eta)}{\cosh(\eta) - \cos(\xi)}, \end{equation}

where

(4.37)\begin{align} \left.\begin{array}{c@{}} N(\eta) = A\eta - F\cosh{2\eta} - G\sinh{2\eta},\\ M(\eta) = B\sinh{\eta} + C\cosh{\eta} + E\eta\sinh{\eta} + F\cosh{\eta}\cosh{2\eta} + G\cosh{\eta}\sinh{2\eta}. \end{array}\right\} \end{align}

Here, $A,\, B,\, C,\, E,\, F$ and $G$ are constants which, letting $\alpha = \eta _1 + \eta _2$ and $\beta = \eta _1 - \eta _2$, satisfy

(4.38)\begin{align} \left.\begin{array}{c@{}} A = \dfrac{\varOmega c a \cosh{\beta}}{\sinh{\beta}}\dfrac{\beta \sinh{\eta_2} - \sinh{\beta} \sinh{\eta_1}}{\beta (\cosh{\alpha}\cosh{\beta} - 1) - \sinh{\beta} (\cosh{\alpha} - \cosh{\beta})},\\ E = \varOmega c a \dfrac{\cosh{\beta} \cosh{\eta_1} - \cosh{\eta_2}}{\beta (\cosh{\alpha}\cosh{\beta} - 1) - \sinh{\beta} (\cosh{\alpha} - \cosh{\beta})},\\ C = \dfrac{\sinh{\eta_2}}{2}(E \sinh{\eta_2} + A \cosh{\eta_2} + \varOmega c a) + \dfrac{\sinh{\eta_1}}{2}(E \sinh{\eta_1} + A \cosh{\eta_1}),\\ B = -E \eta_2 - \cosh{\eta_2}(E \sinh{\eta_2} + A \cosh{\eta_2} + \varOmega c a),\\ F = -\dfrac{A}{2}\dfrac{\sinh{\alpha}}{\cosh{\beta}} , \quad G = \dfrac{A}{2}\dfrac{\cosh{\alpha}}{\cosh{\beta}}. \end{array}\right\} \end{align}

From Kazakova & Petrov (2016), this flow field takes one of two forms. When the mill is relatively close to the centre of the Petri dish, the flow has no stagnation points and the streamlines are circular (figure 5(a) gives a typical example). When the mill is close to the boundary of the Petri dish, the flow has a stagnation point (figure 5(b) gives a typical example). Mathematically, a stagnation point exists when $u_{\eta } = u_{\xi } = 0$ i.e. $\xi = 0$ and $\eta = \eta ^{\star }$ where $\eta _1 < \eta ^{\star } < \eta _2$ satisfies

(4.39)\begin{equation} V(\eta^{\star}) = 0 : V(\eta) = \frac{\textrm{d}}{\textrm{d} \eta}\left( N(\eta) + \frac{M(\eta)}{\cosh{\eta} - 1} \right). \end{equation}

Note that when $\xi = {\rm \pi}$, although $u_{\eta } = 0$, $u_{\xi } \neq 0\ \forall \, \eta \in (\eta _1, \eta _2)$. Without loss of generality, let $\varOmega > 0$ i.e. $V(\eta _1) = 0$ while $V(\eta _2) <0$. If $V'(\eta _1) < 0$, $V$ decreases monotonically and no such $\eta \in (\eta _1, \eta _2)$ exists. Conversely if $V'(\eta _1) > 0$, $V$ achieves positive values in $[\eta _1, \, \eta _2 ]$ and so by the intermediate value theorem, such a $\eta \in (\eta _1, \eta _2)$ exists. Hence, in $\{ b, \, c\}$ phase space, the critical curve separating the two regions satisfies $V'(\eta _1) = 0$. Furthermore from Kazakova & Petrov (2016), a good approximation to the boundary is the interpolation curve

(4.40)\begin{equation} b(\textit{c}) = 0.36(1 - c) + 0.08 c(c-1). \end{equation}

Figure 5. Streamlines of the flow highlighting the two distinct possibilities; namely, no stagnation points in (a), where $b = 0.25$ and $c = 0.2$, and stagnation points in (b), where $b = 0.373$ and $c = 0.298$.

5. The orbit of the circular mill centre

The flow exerts a force $\boldsymbol {F}$ on the mill where $\boldsymbol {F} = F_x \hat {\boldsymbol {x}} + F_y \hat {\boldsymbol {y}}$ satisfies

(5.1)\begin{equation} \boldsymbol{F} = \int^{\pi}_{-\pi} (\sigma_{\eta\eta} \hat{\boldsymbol{\eta}} + \sigma_{\eta\xi} \hat{\boldsymbol{\xi}})_{\eta = \eta_2} h \, \textrm{d} \xi. \end{equation}

Here, the bipolar basis vectors $\hat {\boldsymbol {\eta }}$ and $\hat {\boldsymbol {\xi }}$ satisfy $\hat {\boldsymbol {\eta }} = f \hat {\boldsymbol {x}} + g \hat {\boldsymbol {y}}$ and $\hat {\boldsymbol {\xi }} = g \hat {\boldsymbol {x}} - f \hat {\boldsymbol {y}}$ where

(5.2a,b)\begin{equation} f = \frac{1 - \cosh{\eta}\,\cos{\xi}}{\cosh{\eta} - \cos{\xi}}\quad \mbox{and} \quad g = - \frac{\sin{\xi}\sinh{\eta}}{\cosh{\eta} - \cos{\xi}}, \end{equation}

while $\sigma _{\eta \eta }$ and $\sigma _{\eta \xi }$ are components of the stress tensor $\sigma _{ij} = -p \delta _{ij} + \mu (\partial u_i/\partial x_j + \partial u_j/\partial 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 $\partial u_i/\partial x_j$ when expressed in bipolar coordinates for general $\boldsymbol {u}$.

This system, in a domain symmetric about the line $\theta = 0$, is forced by a fluid flow even in $\theta$. Hence, since it admits a general separable form where each term is either even or odd in $\theta$ (4.14), $\varphi$ is even in $\theta$ and hence $p$ and $\sigma _{rr}$ are also even in $\theta$. Similarly, $\sigma _{r \theta }$ is odd in $\theta$ and hence from rewriting (5.1) in terms of cylindrical polar coordinates, we find $F_x = 0$, a result to be expected from reversibility. This force causes the mill centre to slowly orbit on a larger time scale 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.

5.1. Near-field circular mill

Building from § 4.6, substituting (4.23) and (4.27) into (4.11a,b) using standard properties of modified Bessel functions and then integrating yields

(5.3a,b)\begin{equation} p^0 = p_0 , \quad p^1 = \mu \kappa^2\sin{\theta}\left( \frac{\delta_1}{r} - \gamma_1 r \right), \end{equation}

where $p_0$ is a constant. Furthermore, since $\sigma _{rr} = -p + 2\mu \, {\partial u_r}/{\partial r}$ while $\sigma _{r\theta } = \mu ({\partial u_{\theta }}/{\partial r} + ({\partial u_r}/{\partial \theta })/r - u_{\theta }/r)$, we obtain

(5.4a,b)\begin{gather} \sigma^0_{rr} = 0 , \quad \sigma^0_{r\theta} = \alpha_0 \mu \left( \kappa^2K_0(\kappa r) + \frac{2\kappa K_1(\kappa r)}{r} \right) + \beta_0 \mu \left( \kappa^2 I_0(\kappa r) - \frac{2\kappa I_1(r/X)}{r} \right), \end{gather}

(5.5)\begin{gather} \sigma^1_{rr} = \mu \sin{\theta} \left( -2\alpha_1 \left( \frac{\kappa K_0(\kappa r)}{r} + \frac{2K_1(r/X)}{r^2} \right) + 2\beta_1 \left( \frac{\kappa I_0(\kappa r)}{r} -\frac{2 I_1(\kappa r)}{r^2} \right) \right. \nonumber\\ \hspace{-10.7pc} + \left. \gamma_1 \kappa^2 r - \delta_1 \left( \frac{\kappa^2}{r} + \frac{4}{r^3} \right) \right), \end{gather}

(5.6)\begin{align} \sigma^1_{r\theta} &= \mu \cos{\theta} \left( \alpha_1 \left( \frac{2\kappa K_0(\kappa r)}{r} + K_1(\kappa r) \left( \kappa^2 + \frac{4}{r^2} \right) \right)\right.\nonumber\\ &\quad - \left.\beta_1 \left( \frac{2\kappa I_0(\kappa r)}{r} - I_1(\kappa r) \left( \kappa^2 + \frac{4}{r^2} \right) \right)+ \frac{4\delta_1}{r^3} \right). \end{align}

Hence, the flow exerts a force $\boldsymbol {F}$ on the mill, where $\boldsymbol {F} = F_x \hat {\boldsymbol {x}} + F_y \hat {\boldsymbol {y}}$ satisfies

(5.7)\begin{gather} \boldsymbol{F} =\left. \int^{\pi}_{\pi} ( \sigma_{rr} \hat{\boldsymbol{r}} + \sigma_{r\theta} \hat{\boldsymbol{\theta}} ) \right|_{r = c} c \, \textrm{d} \theta \Longrightarrow \end{gather}

(5.8)\begin{gather}F^0_x = F^0_y = F^1_x = 0, \end{gather}

(5.9)\begin{gather}F^1_y = {\rm \pi}\mu\kappa^2 c \left( \alpha_1 K_1(\kappa c) + \beta_1 I_1(\kappa c) + \gamma_1 c - \frac{\delta_1}{c}\right) = {\rm \pi}\mu\kappa^2 c^2\left( \varOmega c - \frac{2\delta_1}{c^2} \right). \end{gather}

Note that, in the front of this expression, we have $(\kappa c)^2$ rather than $c^2$ i.e. the effective radius of the mill is modulated by the screening length $\kappa$. For the values taken in figure 4, $F^1_y$ is positive, i.e. the mill centre orbits clockwise in a circle centred on the middle of the Petri dish.

5.2. Far-field circular mill

Since (4.11a,b) 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 \rightarrow F_x = 0$. Substituting (4.36) into (4.11a,b) and then integrating gives the pressure

(5.10)\begin{align} p &= \frac{2\mu}{a^2}(E\sinh{\eta} \sin{\xi} + F\sinh{2\eta} \sin{2\xi} - 2F\sinh{\eta} \sin{\xi} \nonumber\\ & \quad + G\cosh{2\eta} \sin{2\xi} - 2G\cosh{\eta} \sin{\xi}). \end{align}

Shifting the basis vectors back to Cartesian coordinates, the force can be expressed in the form

(5.11)\begin{equation} \boldsymbol{F} = \frac{\mu}{a^2}\int^{\pi}_{-\pi} (\,f_x \hat{\boldsymbol{x}} + \,f_y \hat{\boldsymbol{y}})_{\eta = \eta_2} h \, \textrm{d} \xi, \end{equation}

where $f_x$ and $f_y$ are explicit functions of $\eta$, $\xi$ and $\{ A,B,C,E,F,G\}$. However, $f_x$ is odd with respect to $\xi$ at $\eta = \eta _2$ since

(5.12)\begin{align} &\frac{1}{2}(\,f_x(\eta,\xi) + \,f_x(\eta,-\xi))\nonumber\\ &\quad =\frac{2\sin^2{\xi}(1-\cosh{\eta} \cos{\xi})}{(\cosh{\eta} - \cos{\xi})^2}(C\cosh{\eta} \!+\! B\sinh{\eta} + E\eta\sinh{\eta} \!+\! F\cosh{\eta}(2\cosh^2{\eta} - 1) \nonumber\\ &\qquad + 2G\cosh^2{\eta}\sinh{\eta}) = 0 \mbox{ at } \eta = \eta_2. \end{align}

Therefore, as expected, $F_x = 0$. Also, $f_y$ can be similarly simplified, removing the terms odd in $\xi$, to give

(5.13)\begin{equation} \boldsymbol{F} = \frac{\mu \hat{\boldsymbol{y}}}{8 a} \sum_{i=0}^{4} g_i(\eta_2) I_{i}, \end{equation}

where $g_i = g_i(\eta ) : i \in \{0,1,2,3,4\}$ are given for completeness in appendix B.2 while $I_n$ satisfies

(5.14)\begin{equation} I_n = \int^{\pi}_{\pi} \frac{\cos{(n \xi)}}{(\cosh{\eta} - \cos{\xi})^3} \, \textrm{d} \xi. \end{equation}

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 6(a) 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 orbits in the same angular direction as the worms. Small mills have negative $F_y$ (in the white region), i.e. the mill centre orbits in the opposite direction to the worms.

Figure 6. (a) The force that the flow exerts on the disc expressed as a function of $c$ for a range of values of $b (0.85, 0.8, 0.65, 0.52, 0.4)$. (b) The critical radius $c^{\star }$ expressed as a function of $b$. As $b \rightarrow 1$, $c^{\star } \rightarrow 1 - b$ (the red dashed line).

Hence, we define the critical radius $c^{\star } = c^{\star }(\textit {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^{\star } \rightarrow 1$. Furthermore, $c^{\star }$ has a maximum of 0.222 at $b = 0.70$.

5.3. 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 centres of the mills do indeed move on circles centred at the middle of the arena. Furthermore, the direction of the movement also concurs 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 7a and 7b), the circular mill radius $c$ decays slowly over time, always being less than the critical radius $c^{\star }(\textit {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 orbit anti-clockwise, increasing in angular speed as time progresses. This is indeed what we see in figure 7(a) with the darker, later-time points less clustered together than the lighter, earlier-time points.

Figure 7. Circular mill data for two different experiments, (a and b) and (c and d). In (a and c) the location of the mill centre is plotted on the ($x, y$) plane. Points shaded a darker blue denote a later time, as quantified by the colour bar. (b,d) Show $b$ (orange filled circles), $c$ (green filled circles) and the critical radius $c^{\star }(\textit {b})$ (red dashed line), plotted as functions of time. Here, $\{ x, \, y, \, b, \, c, \, c^{\star } \}$ have all been normalised by the Petri dish radius (denoted $R$).

In contrast, we see much more variation in $c$ in the second experiment (figures 7c and 7d) with points both above and below $c^{\star }$ (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 orbit 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.

6. Binary circular mill systems

During the evolution of the system, multiple mills can emerge at the same time (figure 8 and supplementary movie 2). This is to be expected since the worms can only interact locally with each other and hence cannot 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 are 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 cannot yield an analytic solution here.

Figure 8. Systems containing multiple mills. (a) Nine circular mills of different sizes observed ex situ (in a tub). (b) Binary circular mill system observed in situ (on a beach). Both images are reproduced from Sendova-Franks et al. (2018).

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 (figure 9d–f).

Figure 9. Snapshots of a binary circular mill system for three distinct experiments together with corresponding streamline plots (first mill is light green with red streamlines, second mill is dark green with blue streamlines). (a) and (d) Unstable with the second mill dominating, (b) and (e) stable, (c) and (f) unstable on a longer time scale with the first mill dominating.

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 9a and 9d), the system is unstable as the first mill breaks up. Alternatively, if a stagnation point exists and aligns with the first mill (figures 9b and 9e), 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.

7. 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 $\textit {S. roscoffensis}$ is a prime example, one has to understand the underlying fluid dynamics of the system. From the orbit 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. Such a 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 orbit slowly. In particular, for a single mill in a circular arena, the centre of the mill orbits on a circle whose centre is the middle of the arena, a result consistent with considerations of reversiblity in Stokes flow.

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, not only in terms of the mill centre orbit 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 mill. 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 passively to 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 orbit 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 movement. A second line of enquiry should focus on the observation 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).