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. Reference Bailly2014). 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 Reference Pöhnl, Popescu and Uspal2020) and consider a spheroidal, rigid and impermeable squirmer (Lighthill Reference Lighthill1952) 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 (Reference Dassios, Hadjinicolaou and Payatakes1994) 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. (Reference Pöhnl, Popescu and Uspal2020), 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 (Reference Lauga2020), 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 Reference Michelin and Lauga2010). 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 Reference Pak and Lauga2014; Pedley, Brumley & Goldstein Reference Pedley, Brumley and Goldstein2016).

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. (Reference Dassios, Hadjinicolaou and Payatakes1994), 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 Reference Fortune, Lauga and Goldstein2021) of a number of cases that complements earlier work (Liron & Mochon Reference Liron and Mochon1975) 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.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 Reference Jeong and Moffatt1992; Crowdy & Or Reference Crowdy and Or2010). 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 Reference Lauga2020). 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 Reference Batchelor1967) 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 Reference Goldstein, Pesci and Romero-Rochin1990), 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 (Reference Melesko and Gomilko1999), (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 (Reference Kazakova and Petrov2016), 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 (Reference Kazakova and Petrov2016), 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 (Reference Kazakova and Petrov2016), 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.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.