Stability of dancing Volvox

Abstract

Biflagellate algal cells of the genus Volvox form spherical colonies that propel themselves, vertically upwards in still fluid, by the coordinated beating of thousands of flagella, that also cause the colonies to rotate about their vertical axes. When they are swimming in a chamber of finite depth, pairs (or more) of Volvox carteri colonies were observed by Drescher et al. (Phys. Rev. Lett., vol. 102, 2009, 168101) to exhibit hydrodynamic bound states when they are close to a rigid horizontal boundary. When the boundary is above, the colonies are attracted to each other and orbit around each other in a ‘waltz’; when the boundary is below they perform more complex ‘minuet’ motions. These dances are simulated in the present paper, using a novel ‘spherical squirmer’ model of a colony in which, instead of a time-independent but $\theta$-dependent tangential velocity being imposed on the spherical surface (radius $a$; $\theta$ is the polar angle), a time-independent and uniform tangential shear stress is applied to the fluid on a sphere of radius $(1+\epsilon )a, \epsilon \ll 1$, where $\epsilon a$ represents the length of the flagella. The fluid must satisfy the no-slip condition on the sphere at radius $a$. In addition to the shear stress, the motions depend on two dimensionless parameters that describe the effect of gravity on a colony: $F_g$, proportional to the ratio of the sedimentation speed of a non-swimming colony to its swimming speed, and $G_{bh}$, that represents the fact that colonies are bottom heavy; $G_{bh}$ is the ratio of the time scale to swim a distance equal to the radius, to the time scale for gravitational reorientation of the colony's axis to the vertical when it is disturbed. In addition to reproducing both of the dancing modes, the simulations are able to determine values of $F_g$ and $G_{bh}$ for which they are stable (or not); there is reasonable agreement with the experiments. A far-field model for the minuet motions is also shown to have qualitative agreement, but does not describe some features that are reproduced in the full simulations.

1. Introduction

Volvox is a genus of algae, several species of which form spherical, free-swimming colonies consisting of up to 50 000 somatic cells embedded in an extracellular matrix on the surface, with interior germ cells that develop into a small number of colonies of the next generation. The colony has an anterior–posterior axis of symmetry and each somatic cell bears a pair of beating flagella that enable the colony to swim approximately parallel to this axis. Each cell's flagella beat in approximately the same direction (relative to the colony), i.e. in a plane that is offset from a purely meridional plane by an angle of $10^\circ \text {--}20^\circ$. This offset causes the colony to rotate about its axis as it swims (Hoops 1993); the rotation is always clockwise when viewed from its anterior pole. In still water colonies tend to swim vertically upwards, on average, because they are bottom heavy (daughter colonies being clustered towards the rear), although they are slightly (approximately 0.3 %) denser than water and therefore sediment downwards when the flagella are inactivated. The beating of the flagella of cells at different polar angles, $\theta$, has been observed, in colonies held stationary on a micro-pipette, to be coordinated in the form of a symplectic metachronal wave, which propagates from anterior to posterior in the same direction as the power stroke of the flagellar beat (Brumley et al. 2012). Modelling suggests that hydrodynamic interactions between the flagella of different cells, coupled with flagellar flexibility, provide the mechanism for the coordination (Niedermayer, Eckhardt & Lenz 2008; Brumley et al. 2012, 2015). A detailed survey of the physics and fluid dynamics of green algae such as Volvox has been given by Goldstein (2015).

The radius $a$ of a Volvox colony increases with age (the lifetime of a V. carteri colony is approximately 48 h) although the number and size of cells do not. Drescher et al. (2009) measured the free-swimming properties of many colonies of V. carteri of different radii. Results for the mean upswimming speed $W$, sedimentation speed $V_g$, angular velocity $\varOmega$, mean density difference between a colony and the surrounding fluid (inferred from $V_g$) and time scale $\tau$ for reorientation by gravity when the axis is disturbed from the vertical (a balance between viscous and gravitational torques) are shown in figure 1(a–e). Note that, if the colony were neutrally buoyant but otherwise identical, driven by the same flagellar beating, its swimming speed would be $U = W + V_g$. The largest colonies cannot make upwards progress ($W<V_g$); they naturally sink towards the bottom of the swimming chamber, even when their flagella continue to perform normal upswimming motions. The colony Reynolds number is always less than approximately $0.15$ so that the hydrodynamics is dominated by viscous forces.

Figure 1. Swimming properties of V. carteri as a function of radius: (a) upswimming speed, (b) rotational frequency, (c) sedimentation speed, (d) bottom-heaviness reorientation time, (e) density offset and (f) components of average flagellar force density. (From Drescher et al. (2009), figure 3, with permission.)

Drescher et al. (2009) also observed the behaviour of V. carteri colonies as they swim up towards a horizontal glass plane above or sink towards a horizontal plane below. In the former case, the flagella on a colony that is close to the upper surface continue beating and applying tangential thrust to the nearby fluid. Since the fluid is prevented from flowing from above, the flagellar beating pulls in fluid horizontally from all round and thrusts it downwards. This was observed by seeding the fluid with $0.5\ \mathrm {\mu }\textrm {m}$ polystyrene microspheres and the velocity field measured in horizontal and vertical planes using particle image velocimetry (Drescher et al. 2009).

The simplest model of a swimming colony (Short et al. 2006) ascribes the total mean force exerted by the flagella to a uniform tangential shear stress exerted on the spherical surface, with components $f_\theta$ and $f_\phi$ in the directions of the polar angle $\theta$ and the azimuthal angle $\phi$. Drescher et al. (2009) estimated $f_\theta$ and $f_\phi$, as functions of colony radius $a$, from the measured values of $W$, $V_g$ and $\varOmega$ and low-Reynolds-number hydrodynamics (the Stokes law and the equivalent for rotation)

(1.1a,b)\begin{equation} f_\theta = 6\mu (W+V_g)/{\rm \pi} a, \quad f_\phi = 8\mu \varOmega/{\rm \pi}, \end{equation}

where $\mu$ is the fluid viscosity. The estimated values corresponded to a few pN per flagellar pair, as also found by Solari et al. (2006). It can be inferred from these results that there is a critical colony radius, $a_c$, at which a colony far from any boundaries will hover at rest ($W=-V_g$). For the experiments shown in figure 1, $a_c \approx 300\ \mathrm {\mu }\textrm {m}$.

When two Volvox colonies of approximately the same size, with $a < a_c$, were introduced into the chamber, and when they were both spinning near the upper surface, they were observed to attract each other and to orbit around each other in a bound state, termed a ‘waltz’, by Drescher et al. (2009) (figure 2a and supplementary movie M1 available at https://doi.org/10.1017/jfm.2020.613). When the individual angular velocity $\varOmega$ was approximately $1\ \textrm {rad}\,\textrm {s}^{-1}$, the orbiting frequency $\omega$ was approximately $0.1\ \textrm {rad}\,\textrm {s}^{-1}$. The mutual attraction is consistent with the radial inflow of small particles to a single colony, and the rate of approach of two nearby colonies close to the top wall can be well approximated by treating each colony as a point Stokeslet at the sphere's centre, together with its image system in the plane (Drescher et al. 2009). These results provided the first quantitative experimental verification of the prediction by Squires (2001) of a wall-mediated attraction between downward-pointing Stokeslets near an upper no-slip surface.

Figure 2. (a) Waltzing of V. carteri: top view. Superimposed images taken 4 s apart, graded in intensity. Scale bar is 1 mm; (b) ‘minuet’ bound state: side views 3 s apart of two colonies near the chamber bottom. Arrows indicate the anterior–posterior axes ${\boldsymbol {p}}_m$ at angles $\theta _m$ to the vertical. Scale bar is $600\ \mathrm {\mu }\textrm {m}$. (From Drescher et al. (2009), figures 1a and 5a, with permission.)

However, the orbiting is not a direct consequence of colony 2 translating in the swirling velocity field generated by colony 1, for example, because an isolated colony does not generate a swirl velocity field. The overall torque on a colony is zero; therefore the azimuthal ($\phi$-direction) torque generated by the beating flagella is balanced by an equal and opposite viscous torque on the colony as a whole, as if the flagella were trying to crawl along the inside surface of a shell of fluid, but succeed only in pushing the spherical colony surface in the opposite direction. The orbiting comes about because of near-field effects as the colonies approach each other, which can be seen to be important because the rotation rate of colony 1 is reduced by viscous forces in the gap between the two colonies. To predict the rate of orbiting, Drescher et al. (2009) added a vertical rotlet at the centre of each sphere, together with its image in the plane and, assuming that the surface of each sphere was rigid, used lubrication theory to calculate the force and torque exerted by one sphere on the other for a given rotation rate $\varOmega$. The torque provides the rotlet strength and this, together with the force, determines the orbiting frequency,

(1.2)\begin{equation} \omega \approx 0.069 \log {(d/2a)}\varOmega, \end{equation}

where $d$ is the separation between the two spherical surfaces (Drescher et al. 2009). Equation (1.2) is close to the average of measurements on 60 different waltzing pairs.

Some pairs of colonies with $a\approx {a_c}$, which individually hover, form time-dependent bound states near the bottom of the chamber, with one colony above the other, both colonies oscillating horizontally back and forth. This motion was called a ‘minuet’ by Drescher et al. (2009). In this regime the state of perfectly aligned colony axes is unstable, the flow generated by the swimming of one colony tilting and moving the other one away, while the latter's bottom heaviness and swimming tend to bring it back (see figure 2b and movie M2). The distance between two minuetting colonies is large enough for lubrication effects to be negligible, so Drescher (2010) modelled them as vertical gravitational Stokeslets of equal strength, the resulting sedimentation of each being balanced by steady swimming with speed $U$, directed at a small angle $\theta ^{(m)} (m=1,2)$ to the vertical. This angle is determined by a gyrotactic balance between viscous and gravitational torques, the latter arising from bottom heaviness. The height of each colony above the chamber bottom was taken to be fixed.

We may note that a single sedimenting colony, with $a > a_c$, will hover with its centre at height $H$ above the chamber bottom, where the sedimenting velocity, $V_g = \tilde {F}_g/(6{\rm \pi} \mu a)$ is balanced by the upward velocity of the fluid at that location due to both the squirming ($U = {\rm \pi} a f_\theta /(6\mu )$), and the image Stokeslet in the plane, $V_I = 3\tilde {F}_g / (8{\rm \pi} \mu H)$. Hence

(1.3)\begin{equation} \tilde{F}_g \left(1 - \frac{9}{8} \frac {a}{H}\right) = {\rm \pi}^2 a^2 f_\theta,\end{equation}

which defines the hovering height $H$. When $\tilde {F}_g$ is large compared with the squirming force ${\rm \pi} ^2 a^2 f_\theta$ the hovering height asymptotes to $9a/8$; thus the minimum gap width between squirmer and wall is $a/8$. Hence lubrication theory would not be very accurate, and the far-field Stokeslet model will no longer be very accurate either.

The model for the minuet by Drescher et al. (2009) consisted of two vertical Stokeslets located at the centres of the spheres, ${\boldsymbol {x}}^{(m)}$, plus their image systems in the horizontal plane below (Blake 1971; see figure 3). For small displacements from the vertically aligned state the model of Drescher et al. (2009) led to a two-dimensional dynamical system for the difference in horizontal displacement of the two squirmer centres, $\xi$, and the difference in the angles their two axes make with the vertical, $\varTheta$. This system shows that the equilibrium, $\xi = 0$ and $\varTheta = 0$, goes unstable to a Hopf bifurcation if the gyrotactic time $B$ exceeds a critical value $B_c$, and the instability is a Hopf bifurcation if $U/hB_c$ is large enough. These results were qualitatively consistent with the experimental data.

Figure 3. Model for the minuet bound state: the centres of the two colonies $\boldsymbol {1}$ and $\boldsymbol {2}$ are at ${\boldsymbol {x}}^{(1)}$ and ${\boldsymbol {x}}^{(2)}$, with their images in the plane $e_3 = 0$ at ${\boldsymbol {x}}^{(1')}$ and ${\boldsymbol {x}}^{(2')}$; ${\boldsymbol {r}} = {\boldsymbol {x}}^{(1)}\text {--}{\boldsymbol {x}}^{(2)}$, ${\boldsymbol {R}} = {\boldsymbol {x}}^{(1)}\text {--}{\boldsymbol {x}}^{(2')}$. In the model analysed by Drescher (2010), the angle $\theta ^{(m)}$ between the orientation vector of colony $m$ and the vertical is taken to be small, as is the angle $\psi$ between ${\boldsymbol {r}}$ and the vertical.

It should be noted that the above model is only a coarse approximation, as it neglects vertical motions of the two colonies, as well as their spinning motion about the vertical, which can give rise to orbiting motions when the colony axes are not vertical. Moreover, the main weakness of the model is that it assumes that the Stokeslet strengths of the two squirmers are the same but their hovering heights above the bottom are different, which is not consistent with (1.3). A far-field model that does not make this assumption is outlined in § 4 and appendix B below.

It is also worth pointing out that orbiting motions of two squirmers close together, such as occur near the top wall, are not expected near the bottom, because the horizontal component of the fluid motion generated by the squirming will tend to push fluid out radially and in from above, not suck it in radially and push it out axially, as at the top.

The purpose of this paper is to provide a more detailed fluid mechanical understanding of the pairwise interactions of Volvox by means of an improved model of the above phenomena, which confirms and extends the modelling results of Drescher et al. (2009). We will simulate the flow due to two identical, spinning squirmers in a semi-infinite fluid with a rigid horizontal plane either above or below, for a range of realistic values of the relevant parameters. In § 2 the problem is specified precisely and the numerical method (using the boundary-element method, or BEM) described. The results are presented in § 3 for the waltz and § 4 for the minuet (in § 4 the two squirmers may have different masses). They will consist of representative movies of both the waltz and the minuet (in supplementary material) with careful comparison with the experiments of Drescher et al. (2009). In particular we seek to delineate regions of parameter space in which the dancing modes are stable and investigate what happens when they are not.

2. Basic equations and numerical methods

2.1. A Volvox model

A single colony is modelled as a steady ‘spherical squirmer’, modified from that used previously to study the hydrodynamic interactions between two such model organisms and their behaviour in suspensions (Ishikawa, Simmonds & Pedley 2006; Ishikawa & Pedley 2007a,b; Ishikawa, Locsei & Pedley 2008; Pedley 2016). In those studies the velocity on the spherical surface of the squirmer was taken to be purely tangential and prescribed as a function of polar angle $\theta$, while remaining symmetric about the orientational axis, represented by unit vector $\boldsymbol {p}$. Moreover, the azimuthal, $\phi$, component of velocity was taken to be zero. Thus an isolated squirmer of uniform density would ‘swim’ in the direction of ${\boldsymbol {p}}$, at a constant speed, $U$, but would not rotate. Here, instead of the surface velocity, we prescribe the mean shear stress ${\boldsymbol {f}}_s$ generated by the beating flagella of Volvox as acting tangentially at a radius $\alpha a$, where $\alpha = 1 + \epsilon$ and $\epsilon a$ is proportional to the length of a flagellum (figure 4a); there is no slip on the colony surface $r = a$. The shear stress has components $f_\theta$ and $f_\phi$ in the $\theta$- and $\phi$-directions, i.e. ${\boldsymbol {f}}_s = (\,f_r, f_\theta , f_\phi )$ with $f_r = 0$. Prescribing stresses not velocities is probably more realistic, especially when colonies come close to each other or to a fixed boundary, and especially because it permits no slip on the surface $r=a$. Non-zero $f_\phi$ means that colony rotation is automatically included. The model is still greatly oversimplified because the stresses are taken to be constant, independent of both time and position (i.e. $\theta$ and $\phi$). A similar ‘stress and no-slip’ squirmer model was used, but not fully analysed, in the computations of Ohmura et al. (2018). We also note that the model bears some relation to the ‘traction-layer’ model for ciliary propulsion proposed by Keller, Wu & Brennen (1975), and to the model studied by Short et al. (2006).

Figure 4. Fluid mechanical model of Volvox. (a) The colony is modelled as a rigid sphere, and forces generated by flagella are expressed by a shell of shear stress ${\boldsymbol {f}}_s$ at the distance $\epsilon$ above the spherical surface. (b) Cartesian coordinate system used in the study, in which the gravity ${\boldsymbol {g}}$ acts in the ${\boldsymbol {e}}_3$ direction. A plane wall exists at $e_3 = 0$. The orientation vector of colony $m$ is ${\boldsymbol {p}}^{(m)}$ that has the angle $\theta _p^{(m)}$ from the ${\boldsymbol {e}}_3$ axis.

Solution of the Stokes equations shows that the swimming speed and rotation rate of a neutrally buoyant squirmer in an infinite fluid are given by

(2.1a,b)\begin{equation} U = \frac{a f_\theta}{\mu} \frac{{\rm \pi}}{6} \frac{4\alpha^3 - 3\alpha^2 - 1}{4\alpha} ,\quad \varOmega = -\frac{f_\phi}{\mu} \frac{{\rm \pi}}{8} (\alpha^3 - 1) \end{equation}

(see appendix A for details). Thus, for small $\epsilon = \alpha - 1$, the dimensionless shear stresses are given by $af_\theta / (\mu U) = 4/({\rm \pi} \epsilon )$ and $f_\phi /(\mu \varOmega ) = -8/(3 {\rm \pi} \epsilon )$, which can be approximately inferred from the experimental measurements of figure 1(a–c), as long as a value of $\epsilon$ is assumed (this is discussed further in § 5). Moreover, the stresslet strength, which is important in determining the effect of micro-organisms on the fluid flow around them (Simha & Ramaswamy 2002; Saintillan & Shelley 2008), is identically zero, so according to this model Volvox carteri is approximately a neutral squirmer (Michelin & Lauga 2010), consistent with the observations of Drescher et al. (2010).

As for previous models (Ishikawa & Pedley 2007a,b) we can incorporate bottom heaviness by supposing that the centre of mass of the sphere is displaced from the geometric centre by the vector $-l\boldsymbol {p}$, so when $\boldsymbol {p}$ is not vertical the sphere experiences a torque $-l\boldsymbol {p}\wedge \boldsymbol {g}$ that tends to rotate it back to vertical (${\boldsymbol {g}} = - g{\boldsymbol {k}}$ is the gravitational acceleration). The relevant dimensionless quantity representing the effect of bottom heaviness relative to that of swimming is

(2.2)\begin{equation} G_{bh} = \frac {4{\rm \pi} \rho g a l}{3 \mu U} = \frac {8 {\rm \pi} a}{B U}, \end{equation}

where $\rho$ is the density of the fluid. When $G_{bh} = 8 {\rm \pi}$, the angular velocity of a neutrally buoyant colony that is oriented horizontally in an infinite fluid becomes $U/a$. We also add a point Stokeslet at the centre of the sphere to represent the negative buoyancy of a Volvox colony. The dimensionless quantity representing the effect of sedimentation relative to that of swimming is

(2.3)\begin{equation} F_g = \frac {4{\rm \pi} {\rm \Delta} \rho g a^2}{3 \mu U}, \end{equation}

where ${\rm \Delta} \rho$ is the density difference between a colony and the fluid. When $F_g = 6 {\rm \pi}$, the sedimentation velocity in an infinite fluid is $U$.

2.2. Basic equations

Since the colony Reynolds number is always less than approximately $0.15$, we neglect inertia. In the Stokes flow regime, the velocity $\boldsymbol {u}$ is given by an integral equation over the colony surface $S_c$ and the shell of shear stress $S_f$ as (Pozrikidis 1992)

(2.4)\begin{equation} {\boldsymbol{u}}({\boldsymbol{x}}) = - \frac{1}{8{\rm \pi}\mu} \int_{S_c} {\boldsymbol{\mathsf{J}}}({\boldsymbol{x}}, {\boldsymbol{y}}) \boldsymbol{\cdot} {\boldsymbol{q}}({\boldsymbol{y}})\,\textrm{d}S_c({\boldsymbol{y}}) - \frac{1}{8{\rm \pi}\mu} \int_{S_f} {\boldsymbol{\mathsf{J}}}({\boldsymbol{x}}, {\boldsymbol{y}}) \boldsymbol{\cdot} {\boldsymbol{f}}_s({\boldsymbol{y}})\,\textrm{d}S_f({\boldsymbol{y}}), \end{equation}

where ${\boldsymbol{\mathsf{J}}}$ is the Green function for a flow bounded by an infinite plane wall (Blake 1971), and $\boldsymbol {q}$ is the traction force; ${\boldsymbol {q}}$ is defined as ${\boldsymbol {q}} = {\boldsymbol {\sigma }} \boldsymbol {\cdot } {\boldsymbol {n}} = (-p {\boldsymbol{\mathsf{I}}} + 2 \mu {\boldsymbol{\mathsf{E}}} ) \boldsymbol {\cdot } {\boldsymbol {n}}$, where ${\boldsymbol {\sigma }}$ is the stress tensor, $p$ is the pressure and ${\boldsymbol{\mathsf{E}}}$ is the rate of strain. As explained in § 2.1, the mean shear stress ${\boldsymbol {f}}_s$ generated by the beating flagella acts at a radius $(1 + \epsilon ) a$. On the surface of the rigid sphere $r = a$, the no-slip boundary condition is given by

(2.5)\begin{equation} {\boldsymbol{u}}({\boldsymbol{x}}) = {\boldsymbol{U}} + {\boldsymbol{\varOmega}} \wedge {\boldsymbol{r}},\quad {\boldsymbol{r}} \in S_c, \end{equation}

where ${\boldsymbol {U}}$ and ${\boldsymbol {\varOmega }}$ are the translational and rotational velocities of the colony.

The shear stress ${\boldsymbol {f}}_s$ expresses the thrust force per unit area generated by the flagellar beat. The thrust force should be balanced by the viscous drag force and the sedimentation force. Thus, the force condition for a colony can be given as

(2.6)\begin{equation} \int_{S_c} {\boldsymbol{q}} \,\textrm{d}S_c + \int_{S_f} {\boldsymbol{f}}_s \,\textrm{d}S_f + \frac{4 {\rm \pi} a^3 {\rm \Delta} \rho}{3} {\boldsymbol{g}} + {\boldsymbol{F}}_{rep}= 0 . \end{equation}

Here, ${\boldsymbol {F}}_{rep}$ is the non-hydrodynamic repulsive force between colonies and between a colony and a wall. Although lubrication flow can prevent a rigid sphere colliding with a plane wall, the shear-stress shell can easily collide with a plane wall or another shear-stress shell. In the case of a real Volvox, the collision tends to deform the flagella, and the repulsive force may be generated by the elasticity of flagella. Here, we do not model such a complex phenomenon, but follow Brady & Bossis (1985) and Ishikawa & Pedley (2007b) and use the following function

(2.7)\begin{equation} {\boldsymbol{F}}_{rep} =\alpha_1 \frac{\alpha_2 \exp(-\alpha_2 \lambda)\boldsymbol{s}} {(1-\exp(-\alpha_2 \lambda)) s}, \end{equation}

where ${\boldsymbol {s}}$ is the centre-to-centre vector between two colonies or the normal vector from the wall to the colony centre; $\alpha _1$, $\alpha _2$ are dimensionless coefficients and $\lambda$ is the minimum separation between two shear-stress shells or between a shear-stress shell and the wall, non-dimensionalized by $a$. The coefficients used in this study are $\alpha _1 = 10$ and $\alpha _2 = 10$ for colony–wall interactions, whereas $\alpha _1 = 1$ and $\alpha _2 = 10$ for colony–colony interactions. The parameters were chosen to avoid collision while keeping computational efficiency. Since the colony surfaces are at least $2 \epsilon a$ apart in the present study, the repulsive force remains much smaller than the lubrication forces, and is much less significant than in Ishikawa et al. (2006), in which the gap could become infinitely small. The minimum separation obtained with these parameters is of the order of $10^{-2}a$.

The torque condition is given by

(2.8)\begin{equation} \int_{S_c} {\boldsymbol{r}} \wedge {\boldsymbol{q}} \,\textrm{d}S_c + \int_{S_f} {\boldsymbol{r}} \wedge {\boldsymbol{f}}_s \,\textrm{d}S_f - \frac{4 {\rm \pi} a^3 \rho l}{3} {\boldsymbol{p}} \wedge {\boldsymbol{g}} = 0 . \end{equation}

The repulsive force does not contribute to the torque balance.

2.3. Numerical methods

The governing equations are discretized by the BEM (Pozrikidis 1992). By combining the governing equations and the boundary condition, a set of linear algebraic equations can be generated. Each spherical surface of a colony is discretized by 320 triangles, while each spherical shear-stress shell is discretized by 1280 triangles. The numerical integration is performed using 28-point Gaussian polynomials, and the singularity is solved analytically. Time marching is performed using a fourth-order Runge–Kutta method. The details of these numerical methods can be found in Ishikawa et al. (2006).

The coordinate axes are taken as shown in figure 4(b). Gravity acts in the $- {\boldsymbol {e}}_3$ direction, i.e. ${\boldsymbol {k}} = {\boldsymbol {e}}_3$, and an infinite plane wall exists at $e_3 = 0$. When we investigate a waltzing motion beneath the wall, colonies are placed in the negative $e_3$ half-space. When we investigate a minuet motion above the wall, on the other hand, colonies are placed in the positive $e_3$ half-space. ${\boldsymbol {p}}^{(m)}$ is the orientation vector of colony $m$. The angle of ${\boldsymbol {p}}^{(m)}$ from ${\boldsymbol {e}}_3$ is defined as $\theta _{p}^{(m)}$.

Parameter values are varied so as to cover experimental conditions. By assuming that the relaxation time, $B$, defined as $6 \mu / \rho g l$, is about 14 s (Drescher et al. 2009) and the colony swims about one body length per second, $G_{bh}$ is approximately 2. In the present study, $G_{bh}$ is varied in the range $0\text {--}100$. Small and young Volvox swim faster than the sedimentation speed, though large and old Volvox cannot swim upwards. In order to cover both conditions, $F_g$ is varied in the range $0 \text {--} 9 {\rm \pi}$. The tilt angle of the flagellar beating plane with respect to the colonial axis was approximately $15^\circ$ (Drescher et al. 2009). We thus set $\arctan (\,f_{\phi }/f_{\theta }) = 15^\circ$ throughout this study; $\epsilon$ is set as $0.05$ on the basis of experimental observations (Brumley et al. 2012).

3. Waltzing motion beneath a top wall

We first calculate the flow field around a single colony hovering beneath a plane wall. The colony is directed vertically upwards, and the hovering motion is stable when the colony is sufficiently bottom heavy. The results for velocity vectors in the $e_1\text {--}e_3$ plane are shown in figure 5(a) ($G_{bh}=25$ and $F_g = 3 {\rm \pi}$). We see that fluid is pulled in radially towards the colony and then goes downward. The white broken arrows in the figure schematically show the vortex structure.

Figure 5. A hovering colony beneath a top wall ($G_{bh} = 25$ and $F_g = 3 {\rm \pi}$). (a) Velocity vectors around a stably hovering colony beneath a top wall. The colony is directed vertically upwards. White broken arrows schematically show the vortex structure. (b) Time change of centre-to-centre distance s between two colonies, where $t_0$ is the time of collision. The broken line indicates experimental result Drescher et al. (2009), and the solid line indicates our simulation result. The simulation result is dimensionalized by assuming that the colony swims one body length per second in the absence of gravity.

It is not possible to predict precisely the gap width between the top of a single colony and the wall. Although lubrication forces prevent a rigid sphere from colliding with the wall, this is not the case for a squirmer, nor for a real Volvox colony, for which there will be contact between the flagella and the wall, and the flagellar beating will be modified. This is why we have incorporated the near-field repulsion force (2.7) between the squirmer and the wall. This guarantees that the minimum distance between the wall and the squirmer can be no less than $\epsilon a$.

The computed flow field is similar to that observed experimentally. When two colonies hover beneath a wall, they are attracted to each other due to the inward suction. The change with time of the centre-to-centre distance $s$ between two colonies with $G_{bh}=25$ and $F_g = 3 {\rm \pi}$ is shown in figure 5(b), where $t_0$ is the time of collision. The broken line indicates experimental results from Drescher et al. (2009) averaged over 60 colonies. The solid line indicates the simulation result, in which the time scale is dimensionalized by using the characteristic time of $a/U = 0.5$ s. The attraction velocity increases as the distance decreases, which is captured in the simulation. This result is not unexpected, as it was shown by Drescher et al. (2009) that the average experimental curve was almost identical with that predicted analytically by Squires (2001) for two equal, vertical Stokeslets close to a horizontal boundary, using the far-field, point-particle model.

Two nearby colonies beneath a wall orbit around each other in a ‘waltz’, as stated above. Figure 6 and supplementary movie 3 show the waltzing motion reproduced by the simulation under the condition of $G_{bh} = 25$ and $F_g = 3 {\rm \pi}$. We see that two colonies orbit around each other with a constant rotation rate. The radius of the orbiting is approximately 1.07, so the two shear-stress surfaces are very close to contact.

Figure 6. Waltzing motion of two colonies ($G_{bh} = 25$ and $F_g = 3 {\rm \pi}$). (a) Trajectories of two colonies. White or black circles indicate the centre positions of each colony, which are plotted with the time interval of $20a/U$. The colonies attracted each other and finally displayed waltzing motions. (b) Sample image of waltzing colonies, where two colonies are trapped just below the top wall and orbit around each other. Red and yellow arrows schematically show spin and orbit motions, respectively. (See supplementary movie 3.)

In order to discuss the stability of the waltzing motion, we calculated the change of orientations and distance between two nearby colonies. Figure 7(a) shows the definitions of parameters used in the analysis. Let ${\boldsymbol {x}}^{(m)} = ( x_1^{(m)}, x_2^{(m)}, x_3^{(m)})$ and ${\boldsymbol {p}}^{(m)} = ( p_1^{(m)}, p_2^{(m)}, p_3^{(m)})$ respectively be the position vector and the orientation vector of colony $m$. For simplicity, we assume that the two colonies align in the $e_2$-direction, i.e. $x_1^{(1)} = x_1^{(2)}$ and $x_3^{(1)} = x_3^{(2)}$. The orientation vectors were set as $p_1^{(1)} = -p_1^{(2)}, p_2^{(1)} = -p_2^{(2)}$ and $p_3^{(1)} = p_3^{(2)}$, so that a rotation of ${\rm \pi}$ around the $e_3$-axis leaves the configuration unchanged.

Figure 7. Stability of waltzing motion. White vectors indicate the angular velocity in spherical coordinates $\theta _p - \phi _p$. Colours indicate the separation velocity of two colonies. (a) Definition of ${\boldsymbol {s}}$ and $\phi _p$. (b) Stability in the case of $G_{bh} = 25$ ($F_g = 3 {\rm \pi}$). Stable waltzing motion is observed. Stable orientation ($\theta _p = 0.075, \phi _p = 0.092$) is shown by a black circle. Inset is the magnified image of the black rectangle. (c) Stability in the case of $G_{bh} = 5$ ($F_g = 3 {\rm \pi}$). Waltzing motion is unstable.

The length of the centre-to-centre vector is set as $2.14a$. The colour-coded values of $\textrm {d}s/\textrm {d}t$ indicate the separation velocity between the two colonies, i.e. $\textrm {d}s/\textrm {d}t < 0$ is attractive, whereas $\textrm {d}s/\textrm {d}t > 0$ is repelling. $\phi _p$ is the angle of the projection vector of ${\boldsymbol {p}}$ in the $e_1 e_2$ plane from the line connecting the two colony centres. Because of the condition $p_1^{(1)} = -p_1^{(2)}, p_2^{(1)} = -p_2^{(2)}$, $\phi _p$ is the same for each colony; $\theta _p$, defined in figure 4(b), is also the same for each colony.

Figure 7(b) shows the results of fluxes in phase space to understand the direction of motion with $G_{bh} = 25$ ($F_g = 3 {\rm \pi}$), in which stable waltzing motion was observed. The horizontal axis indicates $\phi _p$, the vertical axis indicates $\theta _p$, and the components of the white vectors are $\textrm {d} \phi _p / \textrm {d}t$ and $\textrm {d} \theta _p / \textrm {d}t$ at given $\phi _p$ and $\theta _p$. Moreover, the colour indicates the separation velocity $\textrm {d}s/\textrm {d}t$. By following the white vectors and considering the separation velocity, we can understand how the configurations of two colonies change with time. The black dot in figure 7(b) indicates the stable point, where a point sink of the white vector field exists with $\textrm {d}s / \textrm {d}t \le 0$. We can conclude that the waltzing motion with $G_{bh} = 25$ is stable with respect to small fluctuations in the colony configurations.

In the case of $G_{bh} = 5$ ($F_g = 3 {\rm \pi}$), on the other hand, there is no stable point (figure 7c). Thus, colonies with $G_{bh} = 5$ eventually repel each other and do not show the stable waltzing motion. Figure 8 shows the phase diagram for the stability of waltzing motion in $G_{bh} - F_g$ space. The waltzing becomes unstable in the bottom grey region, while it is stable in the top white region. The boundary lies between $G_{bh} = 5$ and 10, and $F_g$ has little influence on it. The mean $G_{bh}$ value in the experiments can be estimated as approximately 1.8 (Drescher et al. 2009), which is smaller than the stable limit in the simulation. There might be two possibilities to explain the discrepancy. First, the flagella beat might be disturbed in the experiments due to interaction with the top glass wall. If the flagella beat is disturbed, the torque generated by the flagella will be reduced, which effectively increases their bottom heaviness to stabilize the vertical orientation of the colony. Another possibility is that it was only colonies with large $G_{bh}$ that were observed in the experiment, because only they could stay near a top wall for a sufficiently long time.

Figure 8. Phase diagram on the stability of waltzing motion in $G_{bh} - F_g$ space. Circles indicate the simulation cases. The waltzing is unstable in the bottom grey region, and stable in the top white region

Next, we discuss the mechanism of the waltzing motion. For simplicity, we again assume the simple configuration shown in figure 7(a). The coordinate system, forces, torques and velocities are defined in figure 9. In the Stokes flow regime, the motions of two rigid spheres in the presence of a plane wall can be described by using the mobility tensor (Kim & Karrila 1992). Hence, the orbiting velocity of colony 1, $U_1^{(1)}$, which is equivalent to the orbiting rotation rate multiplied by the orbiting radius, can be given as follows

(3.1)\begin{equation} U_{x1} = M_{1,1} F_1^{(1)} + M_{1,5} T_2^{(1)} + M_{1,6} T_3^{(1)} + M_{1,7} F_1^{(2)} + M_{1,11} T_2^{(2)} + M_{1,12} T_3^{(2)} , \end{equation}

where $M_{i,j}$ is the $(i,j)$ component of the $12 \times 12$ mobility tensor; $F_2, F_3$ and $T_1$ do not contribute to $U_1$ due to the symmetry of the problem; $M_{i,j}$ can be calculated by BEM in the stable waltzing configurations with $G_{bh}=25$ and $F_g = 3 {\rm \pi}$, and the results are $(M_{1,1}, M_{1,5}, M_{1,6}, M_{1,7}, M_{1,11}, M_{1,12}) = 10^{-2} (2.39, -0.13, -0.8, 0.31, -0.01, -0.50)$. The forces and torques can also be calculated by BEM by fixing two colonies in space with the active shear stress ${\boldsymbol {f}}_s$. The results are $(F_1^{(1)}, T_2^{(1)}, T_3^{(1)}, F_1^{(2)}, T_2^{(2)}, T_3^{(2)}) = (-3.1, 1.5, -13.9, 3.1, -1.5, -13.9)$. The largest positive contribution comes from $M_{1,12} T_3^{(2)} = 0.069$, and other major positive contributions are $M_{1,7} F_1^{(2)} = 0.010$ and $M_{1,6} T_3^{(1)} = 0.011$. The largest negative contribution comes from $M_{1,1} F_1^{(1)} = -0.074$. Thus, one may roughly say that the orbiting velocity $U_1^{(1)}$ is mainly generated by $T_3^{(2)}$ and inhibited by $F_1^{(1)}$. $T_3^{(2)}$ is induced on the colony as a reaction torque from the flagellar beat. Negative $F_1^{(1)}$ is induced because the traction force ${\boldsymbol {q}}^{(1)}$ acting in regions A and $\textrm {A}^\prime$ in figure 9 are different. In region A, ${\boldsymbol {q}}^{(1)}$ is induced by the shear stress of colony 1, ${\boldsymbol {f}}_{s}^{(1)}$. In region $\textrm {A}^\prime$, on the other hand, ${\boldsymbol {q}}^{(1)}$ is induced by the shear stress of both colonies, ${\boldsymbol {f}}_{s}^{(1)}$ and ${\boldsymbol {f}}_{s}^{(2)}$, which tend to cancel each other out. Thus, smaller ${\boldsymbol {q}}^{(1)}$ is generated in region $\textrm {A}^\prime$ than in A.

Figure 9. Schematics of forces and torques exerted on two colonies fixed in space; ${\boldsymbol {f}}_s^{(m)}$ is the shear stress of colony $m$, and ${\boldsymbol {q}}_m$ is the traction generated on the surface of colony $m$; $F_1^{(m)}$ and $T_3^{(m)}$ are the $e_1$ component of the total force and the $e_3$ component of the total torque exerted on colony $m$, respectively. Magnified views of regions $\textrm{A}$ and $\textrm{A}'$ are indicated by the red broken lines.

The angular velocity of an individual spinning with $G_{bh} = 25$ and $F_g = 3 {\rm \pi}$ is $\varOmega _3^{(1)} \approx -0.41$. The angular velocity of orbiting, $\omega$ ($= -2U_1^{(1)}/s$) is approximately 0.013. The ratio of angular velocity of orbiting to that of spinning is about 0.03 in the simulation, which is considerably smaller than the experimental value of 0.19 from Drescher et al. (2009). The ratio, however, can be modified dramatically by reducing the value of $G_{bh}$, as shown in figure 10. When $G_{bh}$ becomes small, the colony orientations tend to tilt from the $e_3$-axis. Since the change in the colony orientations is towards the direction in which the colonies follow each other, the colonies tend to have large following velocities. The following velocity directly contributes to the angular velocity of orbiting. Thus, colonies with small $G_{bh}$ quickly follow each other and $\omega$ increases as $G_{bh}$ is decreased. We see from figure 10 that the effect of $G_{bh}$ on $\omega$ is significant, though that of $F_g$ is small.

Figure 10. Effect of $G_{bh}$ on the angular velocity of orbiting for various $F_g$ values.

4. Minuet motion above a bottom wall

4.1. Numerical results

When colonies become large as they age, so that $F_g$ exceeds approximately $6{\rm \pi}$, the sedimentation velocity exceeds the swimming velocity. Such colonies stay near a bottom wall and sometimes interact with each other, as discussed in § 1. Before going into the details of two-colony interactions, we first calculate the flow field around a solitary colony near the bottom wall. Figure 11(a) shows the simulated velocity vectors around a colony with $F_g = 9 {\rm \pi}$ hovering stably at a height of approximately 3.2 (non-dimensionalized with $a$) over a bottom wall ($G_{bh} = 5$). The wall is at $x_3 = 0$, and the $x_3$-axis is taken as shown in the figure. The colony is directed vertically upwards. We see that strong downward flow is generated around the colony. A toroidal vortex, shown by white arrows, is observed at the side of the colony. The height of stable hovering decreases as $F_g$ increases, as shown in figure 11(c) and predicted by (1.3). However, even for $F_g$ as large as $9{\rm \pi}$, a colony exhibits a positive upswimming velocity when its height above the wall is less than the height of stable hovering (figure 11b).

Figure 11. Hovering of a colony near a bottom wall ($G_{bh} = 5$). (a) Simulated velocity vectors around a stably hovering colony over a bottom wall ($F_g = 9 {\rm \pi}$). The wall exists at $e_3$ = 0, and the $e_3$-axis is taken as shown in the figure. The colony is directed vertically upwards. White arrows schematically show the vortex structure. (b) Upward velocity of a single colony for various $F_g$ values. (c) Stable height of a hovering colony for various $F_g$ values.

Next we examine the ‘minuetting’ bound state of two colonies, of different mass but otherwise identical, near a bottom wall. We show three examples; in each case colony 1 has $F_g = 7.5 {\rm \pi}$ and colony 2 has $F_g = 9 {\rm \pi}$. Both colonies are assumed to have the same $G_{bh}$ value, and $G_{bh}$ is varied from 2 to 6. Other parameters of the colonies, such as $a$ and $\epsilon$, are the same. Colonies 1 and 2 are initially placed at $(-1.5, 0, 5)$ and $(1.5, 0, 3)$, respectively. Trajectories of the two colonies near the bottom wall for time $t$ in the range $0-100$ are shown in figure 12. When $G_{bh} = 2$ (cf. figure 12a and supplementary movie 4), the two colonies attract each other when they are apart, but repel each other when they are close together. Attraction and repulsion are repeated, forming the ‘minuet’ bound state. In order to discuss the oscillation of trajectories in the horizontal direction, we calculate the distance between the two colonies projected onto the $x_1 x_2$ plane. The results are shown in figure 13. We see that the horizontal distance oscillates with amplitude up to 3.5 in the case of $G_{bh} = 2$.

Figure 12. Trajectories of two colonies near a bottom wall during $t = 0\text {--}100$. Trajectories of a colony with $F_g = 7.5 {\rm \pi}$ start from the black circles and end at the white circles. Trajectories of a colony with $F_g = 9 {\rm \pi}$ start from the black triangles and end at the white triangles. (a) Minuet motion with $G_{bh} = 2$ (supplementary movie 4). (b) Minuet motion with $G_{bh} = 3$ (supplementary movie 5). (c) Alignment of two colonies with $G_{bh} = 6$ (supplementary movie 6).

Figure 13. Time course of the changing distance between two colonies with $F_g = 7.5 {\rm \pi}$ and $9 {\rm \pi}$ projected in the $e_1\text {--}e_2$ plane. $G_{bh}$ is varied to 2, 3 and 6.

In the case of $G_{bh} = 3$ (cf. figure 12b and supplementary movie 5), the minuet motion is still observed, but the amplitude of the oscillation in the horizontal distance decreases to about 2 (cf. figure 13). This is because the orientation change induced by hydrodynamic interactions is suppressed by the bottom heaviness. We see that the centres of two colonies form almost two-dimensional trajectories up to $t = 30$, though the trajectories become gradually three-dimensional and the two colonies eventually orbit around each other in a bound state. We note that the direction of orbiting relative to the direction of spin, in this case, is opposite to the ‘waltzing motion’ observed near a top wall. Moreover, it seems that two-dimensional minuet motion can be unstable in the direction perpendicular to the plane. In the case of $G_{bh} = 6$, the two colonies eventually align vertically (cf. figure 12c and supplementary movie 6). Similar alignment was observed in the experiment (Drescher et al. (2009) and supplementary movie 2). The horizontal distance, shown in figure 13, gradually converges to zero in this case. We note that even when two colonies have the same $F_g$ values, such as $F^{(1)}_g = F^{(2)}_g = 7.5 {\rm \pi}$ or $9 {\rm \pi}$, we observed minuet motion, orbiting around each other or vertical alignment depending on the $G_{bh}$ values and the initial positions (not shown here). Figure 13 shows that the orbiting period decreases from nearly $20a/U$ when $G_{bh} = 2$ to approximately $5$ when $G_{bh} = 6$.

In figure 14, we show the phase diagram of two-colony interactions near a bottom wall ($F_g = 7.5 {\rm \pi}$ and $9 {\rm \pi}$). In this case the $G_{bh}$ values for the two colonies may be different. The black circle in the figure indicates unstable motion, in which the centre-to-centre distance between the two colonies exceeds $10a$. The white circles indicate the minuet motion or orbiting around each other in a bound state. The black triangles indicate vertical alignment, in which the distance in the $x_1\text {--}x_2$ plane is less than $0.3a$ for $t = 90\text {--}100$. We see that the colonies show the minuet motion when $G_{bh}$ is in the appropriate range, while they align vertically when $G_{bh}$ is large. The effects of the $G_{bh}$ values of colonies 1 and 2 on the stability are almost symmetric.

Figure 14. Phase diagram of two Volvox colonies interacting near a bottom wall ($F_g = 7.5 {\rm \pi}$ and $9 {\rm \pi}$). The black circle indicates ‘unstable motion’, in which the centre-to-centre distance between two colonies exceeds $10a$. The white circles indicate the ‘minuet motion’. The black triangles indicate ‘vertical alignment’, in which the distance in the $e_1\text {--}e_2$ plane is less than $0.3a$ during $t = 90\text {--}100$.

In further simulations we assume that colony 1 has $F_g = 6.5 {\rm \pi}$ and colony 2 has $F_g = 9 {\rm \pi}$, so that the two colonies have very different heights of stable hovering (cf. figure 11c) and may interact mainly in the far field. The reason why we changed the parameters is because we would like to focus on far-field fluid mechanics in comparing with the predictions of Drescher et al. (2009). Colonies 1 and 2 are initially placed at ($-1.5, 0, 7$) and ($1.5, 0, 3$), respectively. Trajectories of the two colonies near the bottom wall for time $t$ in the range 0–1000 or until centre-to-centre distance exceeds $10a$ are shown in figure 15. When $G_{bh} = 0.1$ (cf. figure 15a), the two colonies first show minuet motion, but eventually move apart from each other. The centre-to-centre distance between the two colonies, in this case, is shown in figure 15(d). We see that the distance oscillates due to the minuet motion, but gradually increases with time. In the range $t > 200$, the distance is larger than 4, so near-field hydrodynamics does not play a major role. Hence, we may say that the minuet motion in this case is unstable even in the far field.

Figure 15. Trajectories of two colonies near a bottom wall for time $t$ in the range 0–1000 or until centre-to-centre distance exceeds $10a$. Trajectories of a colony with $F_g = 6.5 {\rm \pi}$ start from the black circles and end at the white circles. Trajectories of a colony with $F_g = 9 {\rm \pi}$ start from the black triangles and end at the white triangles. (a) Unstable far-field interaction with $G_{bh}= 0.1$. (b) Unstable near-field interaction with $G_{bh} = 0.3$. (c) Stable bound state with $G_{bh} = 1$. Two colonies first show minuet motion, and then orbit around each other. (d) Time change of the centre-to-centre distance of the two colonies in (a–c).

When $G_{bh} = 0.3$ (cf. figure 15b), the two colonies first show minuet motion, then move apart from each other at approximately $t = 500$ (cf. figure 15d), then come close once again at approximately $t = 730$, and eventually separate fully. The second separation is induced by the near-field interactions at approximately $t = 730$, so the minuet motion becomes unstable due to the near-field hydrodynamics in this case. When $G_{bh} = 1$ (cf. figure 15c), the two colonies show a stable hydrodynamic bound state, in which they first show almost two-dimensional minuet motion, before the trajectories become gradually three-dimensional due to the instability of the two-dimensional minuet motion in the direction perpendicular to the plane. At around $t=600$, the two colonies orbit each other, and the motion continues until $t=1000$. These results illustrate that near-field hydrodynamics plays an important role in the hydrodynamic bound states of squirmers.

The graphs in figure 15(d) show that the period of the minuet oscillations decreases from approximately $33a/U$ to approximately $14a/U$ as $G_{bh}$ is increased from $0.1$ to $1.0$.

4.2. Far-field model

Since the pair of minuetting colonies in these simulations do not in general approach very close to each other, it is appropriate to investigate the extent to which a far-field model, such as that briefly outlined by Drescher et al. (2009), may provide a reasonably accurate description of the stability or instability of the system. We consider two vertical Stokeslets, of different strengths, $-8{\rm \pi} \mu F^{(1)}$ and $-8{\rm \pi} \mu F^{(2)}$, located at the centres of the spheres, and perturb them from an equilibrium configuration in which one is vertically above the other. They move, and are rotated from the vertical, because of their images in the plane boundary (see figure 3). This model, which does not include three-dimensional effects such as orbiting, is developed in appendix B, and conditions for instability of the aligned equilibrium are derived. In brief, the vertically aligned state is unstable if the reorientation time scale $B > B_c = P^{-1}$ and the instability represents a Hopf bifurcation if $4UQ > (P + 1/B)^2$, where

(4.1a,b)\begin{align} P = -\frac{1}{h^2}[F^{(2)}(1-\beta_1) - F^{(1)}(1+\beta_2)],\quad Q = \frac{1}{h^3}[F^{(2)}(1-\gamma_1) + F^{(1)}(1-\gamma_2)], \end{align}

and the quantities $\beta _{1,2}$ and $\gamma _{1,2}$ depend on the equilibrium heights of the two Stokeslets, $H$ and $H+h$, and are given in (B 6), (B 8), (B 11) and (B 13). The frequency of the resulting oscillations when $B = B_c$ is predicted to be

(4.2)\begin{equation} \lambda_1 = \sqrt{UQ - B^{-2}_c}. \end{equation}

We first consider the minuetting squirmers of figure 12, with $F^{(1)}_g = 7.5{\rm \pi}$ and $F^{(2)}_g = 9{\rm \pi}$, noting that $F^{(m)} = (Ua/8{\rm \pi} ) F^{(m)}_g$ where $U$ is the swimming speed in isolation, assumed the same for each squirmer (we use the values of $a$ and $U$ from Drescher et al. (2009): $a=300\ \mathrm {\mu } \textrm {m}$ and $U=300\ \mathrm {\mu }\textrm {m}\,\textrm {s}^{-1}$). If we take the values of $H$ and $h$ from figure 11(c), i.e. $H\approx 3a, h\approx 2a$, we find that $P<0$ so the equilibrium state is stable. However, if we take the apparent equilibrium values from figure 12, i.e. $H=h\approx 2a$, we find $P>1/B$, so the equilibrium is unstable if $B<B_c$; the critical value of $B=1/P\approx 32$ s. Moreover, the bifurcation is Hopf since $4UQ > (P+1/B)^2$, so for small amplitudes the squirmers’ positions will oscillate, as observed in both the current computations and the experiments, with a predicted oscillation period of $2{\rm \pi} /\lambda _i \approx 17\ \textrm {s}$; this is close to that shown for $G_{bh} = 2$ in figure 13. However, the predicted value of $B_c$ (about 32 s) corresponds to a value of $G_bh$ of approximately $0.79$; this is significantly lower than the range of critical values shown in figure 14. The discrepancy is more marked if we keep the same values of $h$ and $H$ but change the value of $F^{(1)}_g$ from $7.5{\rm \pi}$ to $6.5{\rm \pi}$, as for the computations leading to figure 15. Now $P$ takes a very small negative value, so the system is stable for all values of $G_{bh}$, however small, which is not consistent with the results shown in figure 15. These findings emphasize how sensitive the theory is to the precise parameter values.

Lastly, we reapply the far-field theory to the experiments of Drescher et al. (2009), using the parameter values quoted there: $h/a = 2, H/a = 1.5$ and $F^{(1)}_g = F^{(2)}_g = 6{\rm \pi}$ (in our notation). Then $P = 0.0756 >0$, so the equilibrium is unstable for large enough reorientation time $B$, with a critical value of $B_c = 13.2$ s ($G_{bh} \approx 2.0$; this is close to the value quoted by Drescher et al. (2009) and close to a typical value in our simulations. However, the assumption of equal $F_g$ values but different equilibrium heights casts doubt on the applicability of the model with these parameters.

We should mention here an interesting recent paper by Bolitho, Singh & Adhikari (2020), in which the dynamics of a pair of identical, self-propelling, self-spinning, active spheres between widely separated parallel planes is modelled as Stokeslets and rotlets, with built-in swimming speed and rotation rate, and their images in the plane boundaries. The analysis is thus also a far-field theory, but is more general than the above because the third dimension is not ignored. The nonlinear dynamical system that arises from the equations leads to limit-cycle oscillations, similar to those shown by Drescher et al. (2009), when the third dimension is ignored and only the bottom boundary is included. These authors also show that, far from either boundary, the system is Hamiltonian and also describes periodic orbits.

5. Discussion

The boundary-element computations in this paper, using the ‘shear-stress and no-slip’ spherical squirmer model for a swimming micro-organism, have succeeded in simulating the dancing motions performed by colonies of V. carteri in the experiments of Drescher et al. (2009). From a broader perspective, they further illustrate how non-trivial time-dependent motions can arise in systems closed to external flows, purely as a consequence of energy injection at the smallest scales. This is the essence of ‘active matter’, as seen in systems ranging from bacterial suspensions to cytoskeletal dynamics.

In the case of the waltzing bound state of pairs of Volvox colonies near the top wall of the chamber, the computations confirm the approximate analysis of Drescher et al. (2009), based on the earlier work of Squires (2001), which utilizes point Stokeslets and rotlets at the centres of the two colonies, and lubrication theory for the space between them when they are close together. In addition we examined the stability of the waltzing state, and found that it is stable if the bottom-heaviness parameter $G_{bh}$ exceeds a critical value between 5 and 10, more or less independently of the gravitational Stokeslet parameter $F_g$. A typical experimental value of $G_{bh}$ was estimated from the data in figure 1 to be approximately 1.8, which is below the computed critical value although the ‘waltzing’ appeared stable; the reason for this discrepancy has not been firmly established, though it seems likely that (a) the flagellar beating is reduced in the narrow gap between the colonies and the plane wall above, due to mechano-sensing and a reduction of the flagellar beat frequency, or due to the flagella sticking to the glass, as has been observed for Tetrahymena (Ohmura et al. 2018), and (b) only colonies with larger values of $G_{bh}$ would stay near the top surface for long enough to attract a neighbour into the waltz.

For the minuet bound state near the bottom boundary, the computational simulations of § 4.1 for two squirmers of different masses (i.e. different $F_g$) show that a bound state, in which the squirmers oscillate two- or three-dimensionally around each other, tends to arise when $G_{bh}$ is below one critical value (above which the vertically aligned state is stable) and above another, less well-defined one, below which the motion is totally unstable and the squirmers do not remain close to each other (figures 12–15). These findings are in qualitative agreement with both the observations and the far-field analysis of § 4.2, but not quantitative agreement since the critical values of $G_{bh}$ are larger in the simulations than in the analysis. A key feature of the minuet simulations, absent from the far-field model, is the three-dimensionality of the motion, as well as the significant tilt angles that the squirmers’ axes make with the vertical.

In the far-field model, of two vertical Stokeslets near the bottom wall, the vertically aligned equilibrium state remains stable if the two values of $F_g$ differ by too great an amount; when it does become unstable $G_{bh}$ falls below a critical value, but this value is lower than that found in the simulations. (If the far-field model is applied to two squirmers of equal $F_g$, as was done by Drescher et al. (2009), the results agree quite well with the simulations, although assuming two different vertical heights in equilibrium is not consistent.) A difference between the far-field model and the simulations is to be expected, because the former does not allow for significant tilt angles or the consequent three-dimensional motions. Further reasons for the discrepancy between the simulations and the model, as discussed in § 4, could be (a) that the heights of the two colonies in a three-dimensional minuet can vary with time, which is not accounted for in the model, and (b) that the far-field assumption is not accurate enough when the distance between the colonies is less than $10a$. Indeed, a major finding of § 4 is that near-field interactions are usually important at some stage during the minuet motions.

Finally, it is appropriate to give further discussion to the ‘shear-stress and no-slip’ squirmer model itself; here we neglect the density difference between the sphere and the fluid, so $F_g = 0$ and sedimentation is absent. The formulae (2.1a,b) relating the mean swimming speed $U$ and the mean angular velocity $\varOmega$ to the shear stresses $f_\theta$ and $f_\phi$ exerted at $r=(1+\epsilon )a$ give

(5.1)\begin{equation} f_\theta = \frac{4}{{\rm \pi}} \frac{\mu U}{\epsilon a} \quad \mathrm{and} \quad f_\phi = \frac{8}{3{\rm \pi}} \frac{\mu \varOmega}{\epsilon}, \end{equation}

to leading order in $\epsilon$ for $\epsilon \ll 1$. These are not the same as given by Drescher et al. (2009), quoted in (1.1a,b) above. Our model recognizes that the shear stress is effectively exerted by the beating flagella, at their tips in the power stroke, and lower down in the recovery stroke. The model requires that the resultant velocity field satisfies the no-slip condition on the (rigid) spherical surface of the Volvox colony, as well as the zero-Stokeslet condition for a self-propelled body. The earlier model balanced the total force exerted by the shear stress against the viscous (Stokes) drag on an inert sphere pulled through the fluid at the same speed. This ignores the fact that the force on the rigid sphere at $r=a$ consists of both the hydrodynamic (shear stress and pressure) force and the equal and opposite reaction force experienced by the flagella and transmitted by them to the rigid sphere. The force exerted by the flagella not only drives the outer flow, but also the high-shear flow in the flagella layer. Put another way, the previous model balanced the rate of viscous energy dissipation in the flow in $r > a_0$ driven by the shear stress against the rate of working of the Stokes drag on the inert sphere, but ignored the energy dissipation between the shear-stress shell and the no-slip spherical boundary, i.e. in the flagella layer. If we model the flow in this layer as a uniform shear flow, as in figure 16, the total rate of energy dissipation in the layer is $D_1 = \mu \times 4{\rm \pi} a^2\,h \times (\,f_\theta /\mu )^2$, where $h=\epsilon a$, which scales as $\epsilon a^3 f_\theta ^2 /\mu$. The dissipation in the outer flow, assumed to scale similarly to that for a translating rigid sphere, i.e. $6{\rm \pi} \mu a U^2$, which from (5.1) scales as $D_2 \sim \epsilon ^2 a^3 f_\theta ^2 /\mu$, is formally an order of magnitude smaller than that in the layer. In V. carteri $\epsilon$ is between $0.05$ and $0.1$ (Solari et al. 2011), which is not extremely small, so in view of numerical factors we cannot say that $D_1 \gg D_2$. But we can be confident that the dissipation rate in the layer is at least as important as that outside it. It follows that a greater shear stress is required to achieve the same swimming speed than in the previous model. Similar considerations apply to the angular velocity $\varOmega$ and the zero-torque condition.

Figure 16. Schematic of the flow field around the ‘shear-stress and no-slip’ squirmer model. There is a no-slip spherical boundary at $r=a$, and uniform tangential stresses $f_\theta$ and $f_\phi$ are applied to the fluid at radius $(1+\varepsilon )a$. Region 1 is defined as $a < r < (1+\varepsilon )a$, whereas region 2 is defined as $(1+\varepsilon )a < r$.

The consequence of the new formulae (5.1) is that the shear stresses for a sphere with $a=200\ \mathrm {\mu }\textrm {m}$, $U=380\ \mathrm {\mu }\textrm {m}\,\textrm {s}^{-1}$ and $\varOmega =1.3\ \textrm {rad}\,\textrm {s}^{-1}$ (figure 1), and flagella length $\epsilon a=15\ \mathrm {\mu }\textrm {m}$, are $f_\theta \approx 1.9\times 10^{-2}\ \textrm {N}\,\textrm {m}^{-2}$, $f_\phi \approx 1.6\times 10^{-2}\ \textrm {N}\,\textrm {m}^{-2}$. Noting that $1\ \textrm {N}\,\textrm {m}^{-2}= 10^3\ \textrm {fN}\,\mathrm {\mu }\textrm {m}^{-2}$, we see that these values are nearly a factor of 2 greater than the corresponding quantities in figure 1(f); this is mainly a consequence of the additional energy dissipation and the $\epsilon ^{-1}$ factors in (5.1).

The old steady-squirmer model, in which the tangential velocity is prescribed on the sphere surface, is very simple to apply and has therefore been used in numerous investigations of hydrodynamic interactions between micro-swimmers and bounding surfaces or each other. In particular, accurate computations of near-field interactions between pairs of such squirmers has been the essential first step in studies of suspensions of squirmers at non-dilute volume fractions (Ishikawa et al. 2006; Ishikawa & Pedley 2007a,b; Ishikawa et al. 2008). The prescribed shear-stress model of Short et al. (2006) has not been used so often, and the ‘shear-stress and no-slip’ model has not been used at all. A useful development would be to see how predictions of suspension rheology, coherent structures (clustering), and nutrient transport properties are changed by the explicit inclusion of no-slip on the inner sphere of radius $a$, given that the active shear stress is applied at radius $r = a(1+\epsilon )$. The new model focusses attention back onto the ciliary layer and hence on the boundary conditions to be applied, for example for nutrient uptake. Magar, Goto & Pedley (2003) considered two possible boundary conditions, one of constant solute concentration at $r = a$, and one of constant solute consumption rate in $0<r<a$. The time dependence of ciliary beating would interact nonlinearly with the mass transfer if the Péclet number $Pe$ were large enough. In this case the relevant Péclet number would be based on the length of a flagellum, $\epsilon a$, the velocity of the flagellar tip, $\sigma \epsilon a$, where $\sigma$ is the beating frequency, and the solute diffusivity, $D$. For Volvox, $\epsilon a \approx 15\ \mathrm {\mu }\textrm {m}$, $\sigma \approx 200\ \textrm {rad}\,\textrm {s}^{-1}$, and $D$ lies somewhat below $10^{3}\ \mathrm {\mu }\textrm {m}^2\,\textrm {s}^{-1}$, so the Péclet number is of the order of 50 or greater, which is large enough for the mass transfer problem to be interesting (cf. Magar & Pedley 2005).