Stability of dancing Volvox

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. 102, 168101 (2009)] 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...

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. 102, 168101 (2009)] 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 θ-dependent tangential velocity being imposed on the spherical surface (radius a; θ is the polar angle), a time-independent and uniform tangential shear stress is applied to the fluid on a sphere of radius (1 + )a, 1, where 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 quantitative agreement with the experiments.

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 • −20 • . 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 (about 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, θ, 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 et al. 2008;Brumley et al. 2012Brumley et al. , 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 about 48 hours) though 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 Ω, mean density difference between a colony and the surrounding fluid (inferred from V g ) and timescale τ for reorientation by gravity when the axis is disturbed from the vertical (a balance between viscous and gravitational torques) are shown in Fig. 1(a-e). Note that, if the colony were neutrally buoyant, the swimming speed would be U = W +V g , and that the colony Reynolds number is always less than about 0.15 so that the hydrodynamics is dominated by viscous forces. Note too that the largest colonies cannot make upwards progress; they naturally sink towards the bottom of the swimming chamber, even when their flagella continue to perform normal upswimming motions. 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 µm polystyrene microspheres and the velocity field measured in horizontal and vertical planes using Particle Image Velocimetry (PIV) (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 θ and f φ in the directions of the polar angle θ and the azimuthal angle φ. Drescher et al. (2009) estimated f θ and f φ , as functions of colony radius a, from the measured values of U , V and Ω and low-Reynolds-number hydrodynamics (the Stokes Law and the equivalent for rotation): where µ 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. For the experiments shown in Fig. 1, a c ≈ 300µ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) (Fig. 2a and movie M1). When the individual angular velocity Ω was about 1 rad s −1 , the orbiting frequency ω was about 0.1 rad 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.
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 (φ-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 succeeds only in pushing the spherical colony surface in the opposite direction. The orbiting could come about because of nearfield effects as the colonies approach each other: 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 Ω. The torque provides the rotlet strength and this, together with the force, determines the orbiting frequency, where d is the separation between the two spherical surfaces. Equation (1.2) is close to the average of measurements on 60 different waltzing pairs. Some pairs of colonies with a ≈ 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 Fig. 2(b) and movie M2). The distance between two minuetting colonies is large enough for lubrication effects to be negligible, so Drescher (2010) modelled each one as a vertical gravitational Stokeslet, the resulting sedimentation being balanced by steady swimming with speed U , directed at a small angle θ (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.
Thus the model consisted of two vertical Stokeslets located at the centres of the spheres, x (m) , plus their image systems in the horizontal plane below (Blake (1971); see Fig. 3). The motion of sphere m is given by where u (m) is the velocity field at sphere m generated by the Stokeslet of the other sphere and its image system in the plane, and is the unit orientation vector of sphere m (note that axes have been taken such that p (m) lies in the x 1 x 3 plane, x 3 being vertically upwards -see Fig. 3). If we consider m = 1, take the Stokeslet strength of each sphere to be (0, 0, −8πµF ), and consider only the x 1 -component of (1.3), then the results of Blake (1971) give where H is the height of the centre of sphere 2 above the plane, assumed constant. Here is the image of x (2) ; r and R are the magnitudes of r and R, respectively. For constant height r 3 = h, small displacement r 1 , and small angles θ (m) , equations (1.5) and (1.3) with m = 1 reduce to (1.7) The corresponding expression for (1.9) In addition, where k is a vertical unit vector, and ω (1) is the vorticity at x (1) due to the sphere at x (2) and the image system, so (1.9) becomes where γ = 1 − h 3 /(h + 2H) 3 , and similarly for dθ (2) /dt. Here B = 6µ/(lρg), where l is the distance from a colony's centre of buoyancy to its centre of mass, is the timescale for gyrotactic reorientation.
1 , Θ = θ (1) − θ (2) and β = (β 1 + β 2 )/2, the system of linearised equations (1.6), (1.8), (1.11) reduces to (1.12) Assuming that ξ and Θ are proportional to e λt , (1.12) gives a quadratic equation for the eigenvalues λ whose roots are (1.13) Thus the equilibrium steady state (ξ = Θ = 0) is unstable if B > B c = h 2 2F β . Moreover the bifurcation at B = B c is a Hopf bifurcation if the quantity in the square root is negative. The parameter values of the experiments by Drescher et al. (2009) (B = 14s, h = 600µm, H = 450µm, F = 6.75 × 10 4 µm 2 /s, U = 300µm/s), which give B c = 13.2s, satisfy these inequalities, which is consistent with the observed oscillations. Drescher et al. (2009) computed the solution to the nonlinear system (1.3) and (1.10) and indeed found that it exhibited limit-cycle oscillations for these parameters over a limited range of values of B: 12s < B < 20s. It should be noted that the nonlinear version of this model is still only a coarse approximation, as it neglects vertical motions of the two colonies, as well as their rotation about the vertical, which can give rise to orbiting motions when the colony axes are not vertical.
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 section 2 the problem is specified precisely and the numerical method (using the Boundary Element Method, or BEM) described. The results are presented in section 3 for the waltz and section 4 for the minuet. 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.

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 et al. (2006) (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 θ, while remaining symmetric about the orientational axis, represented by unit vector p. Moreover, the azimuthal, φ, component of velocity was taken to be zero. Thus an isolated squirmer of uniform density would 'swim' in the direction of p, at a constant speed, U , but would not rotate. Here, instead of the surface velocity, we prescribe the mean shear stress f s generated by the beating flagella of Volvox as acting tangentially at a radius αa, where α = 1 + and a is proportional to the length of a flagellum (Fig.  4a); there is no slip on the colony surface r = a. The shear stress has components f θ and f φ in the θ-and φ-directions, i.e. f s = (f r , f θ , f φ ) 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. Nonzero f φ 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. θ and φ). 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 et al. (1975), and to the model studied by Short et al. (2006).
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 Thus, for small = α − 1, the dimensionless shear stresses are given by (af θ /(µU )) = 4/(π ), which can be approximately inferred from the experimental measurements of Fig. 1(a-c), as long as a value of is assumed (this is discussed further in section 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 bottomheaviness by supposing that the centre of mass of the sphere is displaced from the geometric centre by the vector −lp, so when p is not vertical the sphere experiences a torque −lp ∧ g that tends to rotate it back to vertical (g = −gk is the gravitational acceleration). The relevant dimensionless quantity representing the effect of bottomheaviness relative to that of swimming is where ρ is the density of the fluid. When G bh = 8π, 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 where ∆ρ is the density difference between a colony and the fluid. When F g = 6π, the sedimentation velocity in an infinite fluid is U .

Basic equations
Since the colony Reynolds number is always less than about 0.15, we neglect inertia. In the Stokes flow regime, the velocity u is given by an integral equation over the colony surface S c and the shell of shear stress S f as (Pozrikidis 1992) where J is the Green function for a flow bounded by an infinite plane wall (Blake (1971)), and q is the traction force. q is defined as q = σ · n = (−pI + 2µE) · n, where σ is the stress tensor, p is the pressure and E is the rate of strain. On the surface of the rigid sphere, the no-slip boundary condition is given by where U and Ω are the translational and rotational velocities of the colony.
The shear stress 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 Here 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 where s is the centre-to-centre vector between two colonies or the normal vector from the wall to the colony centre; α 1 , α 2 are dimensionless coefficients and λ 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 α 1 = 10 and α 2 = 10 for colony-wall interactions, whereas α 1 = 1 and α 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 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) The repulsive force does not contribute to the torque balance.

Numerical methods
The governing equations are discretized by a boundary element method (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 fourthorder 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 Fig. 4b. Gravity acts in the −e 3 direction, i.e. k = 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. p (m) is the orientation vector of colony m. The angle of p (m) from e 3 is defined as θ Parameter values are varied so as to cover experimental conditions. By assuming that the relaxation time, B, defined as 6µ/ρgl, is about 14 seconds (Drescher et al. 2009) and the colony swims about one body length per second, G bh is about 2. In the present study, G bh is varied in the range 0 − 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 − 9π. The tilt angle of the flagellar beating plane with respect to the colonial axis was about 15 • (Drescher et al. 2009). We thus set arctan(f φ /f θ ) = 15 • throughout this study; is set as 0.05 on the basis of experimental observations (Brumley et al. 2012).

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 − e 3 plane are shown in Fig. 5a (G bh = 25 and F g = 3π). 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. 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π is shown in Fig. 5b Two nearby colonies beneath a wall orbit around each other in a 'waltz', as stated above. Fig. 6 and Supplementary Movie 3 show the waltzing motion reproduced by the simulation under the condition of G bh = 25 and F g = 3π. 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.
In order to discuss the stability of the waltzing motion, we calculated the change of orientations and distance between two nearby colonies. Fig. 7a shows the definitions of parameters used in the analysis. Let x (m) = (x 3 , so that a rotation of π around the e 3 -axis leaves the configuration unchanged. The length of the centre-to-centre vector is set as 2.14a. The colour-coded values of ds/dt indicate the separation velocity between the two colonies, i.e. ds/dt < 0 is attractive, whereas ds/dt > 0 is repelling. φ p is the angle of the projection vector of p in the e 1 e 2 plane from the line connecting the two colony centres. Because of the condition p (1) 2 , φ p is the same for each colony; θ p , defined in Fig. 4b, is also the same for each colony. Fig. 7b shows the results of the stability analysis with G bh = 25 (F g = 3π), in which stable waltzing motion was observed. The horizontal axis indicates φ p , the vertical axis indicates θ p , and the components of the white vectors are dφ p /dt and dθ p /dt at given φ p and θ p . Moreover, the colour indicates the separation velocity ds/dt. 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 Fig. 7b indicates the stable point, where a point sink of the white vector field exists with ds/dt 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π), on the other hand, there is no stable point (Fig.  7c). Thus, colonies with G bh = 5 eventually repel each other and do not show the stable waltzing motion. Fig. 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 about 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.
Next, we discuss the mechanism of the waltzing motion. For simplicity, we again assume the simple configuration shown in Fig. 7a. The coordinate system, forces, torques and velocities are defined in Fig. 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 , which is equivalent to the orbiting rotation rate multiplied by the orbiting radius, can be given as follows where M i,j is the (i, j) component of the 12 × 12 mobility tensor. F 2 , F 3 and T 1 do not contribute to U 1 due to the symmetry of the problem.
The angular velocity of an individual spinning with G bh = 25 and F g = 3π is Ω 3 ≈ −0.41. The angular velocity of orbiting, ω (= −2U (1) 1 /s) is about 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 Fig. 10. When G bh becomes small, the colony orientations tend to tilt from the e 3 -axis and appear to follow each other. Such inclination dramatically reduces F (1) 1 in Eq.(3.1), and therefore increases ω. We see from Fig. 10 that the effect of G bh on ω is significant, though that of F g is small.

Minuet motion above a bottom wall
When colonies become large as they age, so that F g exceeds approximately 6π, the sedimentation velocity exceeds the swimming velocity. Such colonies stay near a bottom wall and sometimes interact with each other, as discussed in section 1. Before going into the details of two-colony interactions, we first calculate the flow field around a solitary colony. Fig. 11a shows the simulated velocity vectors around a colony with F g = 9π hovering stably at a height of approximately 3.2 (non-dimensionalised 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 Fig. 11c. However, even for F g as large as 9π, a colony exhibits a positive upswimming velocity when its height above the wall is less than the height of stable hovering (Fig.  11b).
Next we examine the 'minuetting' bound state of two colonies near a bottom wall. We show three examples; in each case colony 1 has F g = 7.5π and colony 2 has F g = 9π. 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 , 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 Fig. 12. When G bh = 2 (cf. Fig. 12a and Supplementary Movie 4), the two colonies attract each other when they are apart, but repel each other when they are close to contact. 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 e 1 e 2 plane. The results are shown in Fig. 13. We see that the horizontal distance oscillates with amplitude up to 3.5 in the case of G bh = 2.
In the case of G bh = 3 (cf. Fig. 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. Fig. 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 3-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. Fig.  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 Fig. 13, gradually converges to zero in this case. We note that even when two colonies have the same F g values, such as F g,1 = F g,2 = 7.5π or 9π, we observed minuet motion, orbiting around each other or vertical alignment depending on the G bh values and the initial positions.
In Fig. 14, we show the phase diagram of two-colony interactions near a bottom wall (F g = 7.5π and 9π). 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 e 1 − e 2 plane is less than 0.3a for t = 90 − 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.
Last, we compare the present numerical results with the theory of Drescher et al. (2009) in which two-colony interactions were analyzed by assuming far-field hydrodynamics (see section 1 above). Each colony was assumed to have the same Stokeslet strength, i.e. the same F g , but different equlibrium heights above the bottom boundary. This is not fully compatible with our results in Fig. 11, where different equilibrium heights follow from different values of F g . In our simulations we assume that colony 1 has F g = 6.5π and colony 2 has F g = 9π, so that the two colonies have very different heights of stable hovering (cf. Fig. 11c) and may interact mainly in the far-field. 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 Fig. 15. When G bh = 0.1 (cf. Fig. 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 Fig. 15d. 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. In Eq. (1.12), the G bh value for unstable interactions can be estimated as G bh < 0.017 by assuming H = h = 2a as in Fig. 15c. The present results illustrate that the minuet motion can be unstable even with G bh = 0.1.
When G bh = 0.3 (cf. Fig. 15b), the two colonies first show minuet motion, then move apart from each other at around t = 500 (cf. Fig. 15d), then come close once again at around t = 730, and eventually separate fully. The second separation is induced by the near-field interactions at around t = 730, so the minuet motion becomes unstable due to the near-field hydrodynamics in this case. When G bh = 1 (cf. Fig. 15c), the two colonies show a stable hydrodynamic bound state, in which they first show almost twodimensional 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 also plays an important role in the hydrodynamic bound states of squirmers.

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 Volvox carteri in the experiments of Drescher et al. (2009).
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 Fig. 1 to be around 1.8, which is below the 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 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, Drescher et al. (2009) gave calculated results for (in our notation) a = 300µm, H/a = 2, h/a = 1.5 and F g = 0.75. Their Fig. 5c shows stable vertical alignment of two colonies for B < 12s (G bh > 2.1) and limit cycle oscillations for 12s < B < 20s (2.1 > G bh > 1.3, becoming unstable for larger B (smaller G bh ). Our simulations have the same qualitative features, but the threshold values of G bh are significantly smaller. As discussed in section 4, this discrepancy could be a consequence of the assumption of the same value of F g for both colonies, but in addition it may also have one of the following two causes: (a) the far-field assumption may not be accurate enough when the distance between the colonies is less than 10a, and (b) the minuet motion can become three-dimensional and the heights of the two colonies vary with time, though in Eq.(1.12) two-dimensional trajectories with constant heights were assumed.
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.1) relating the mean swimming speed U and the mean angular velocity Ω to the shear stresses f θ and f φ exerted at r = (1 + )a give f θ = 4 π µU a and f φ = 8 3π µω , (5.1) to leading order in for 1. These are not the same as given by Drescher et al. (2009), quoted in (1.1) above. Our model recognises 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 highshear 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 Fig.16, the total rate of energy dissipation in the layer is D 1 = µ × 4πa 2 h × (f θ /µ) 2 , where h = a, which scales as a 3 f 2 θ /µ. The dissipation in the outer flow, assumed to scale similarly to that for a translating rigid sphere, i.e. 6πµaU 2 , which from (5.1) scales as D 2 ∼ 2 a 3 f 2 θ /µ, is formally an order of magnitude smaller than that in the layer. In V. carteri is between 0.05 and 0.1 (Solari et al. 2011), which is not very small, so in view of numerical factors we cannot say that D 1 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 Ω and the zero-torque condition.
The consequence of the new formulae (5.1) is that the shear stresses for a sphere with a = 200µm, U = 380µm/s and Ω = 1.3rad/s (Fig. 1), and flagella length a = 15µm, are f θ ≈ 1.9 × 10 −2 N/m 2 , f φ ≈ 1.6 × 10 −2 N/m 2 . Noting that 1N/m 2 = 10 3 fN/µm 2 , we see that these values are nearly a factor of 2 greater than the corresponding quantities in Fig. 1(f ); this is mainly a consequence of the additional energy dissipation and the −1 factors in (5.1).

Appendix A.
Here we derive equations (2.1), with reference to Fig. 16. There is a no-slip spherical boundary at r = a and uniform tangential stresses f θ and f φ are applied to the fluid at radius a 0 . The squirmer is taken to swim at speed U in the θ = 0 direction so, relative to the sphere, the velocity at infinity is −U (in the θ = π direction). The sphere rotates with angular velocity Ω about the axis of symmetry; there is no azimuthal velocity at infinity. The squirmer swims freely, so the force and torque exerted on it by the fluid are zero. We solve the axisymmetric Stokes equations separately for the radial and meridional velocity components, and for the swirl velocity component.
We consider the flow in two regions, a < r < a 0 (region 1) and a 0 < r < ∞ (region2), and represent it by stream functions ψ (i) (r, θ), i = 1, 2. In region 1, the solution of the Stokes equation can be written where V n (θ) = 2 n(n+1) sin θP n (cos θ), the P n being Legendre polynomials, and A (1) n are constants to be determined. In region 2 the stream function is The first term incorporates the (unknown) uniform stream at infinity, and C (2) n = D (2) n = 0 for all n so the corresponding contributions to the velocity tend to zero at infinity. Moreover, A (2) 1 is also zero, because this is the Stokeslet term, proportional to the net force on the sphere, which is zero.
The velocity components, pressure and tangential shear stress in region 1 are , with corresponding equations for region 2. The boundary conditions at r = a are u r = u θ = 0 and at r = a 0 are continuity of u r , u θ and the normal stress −p + 2µ∂u r /∂r, and the jump in σ rθ from 1 to 2 is f θ . The constant f θ can also be expanded in a series of the V n : where F 2l = 0, F 2l+1 = (4l + 3)Γ (l + 1 2 )Γ (l + 3 2 ) 4Γ (l + 1)Γ (l + 2) .
Now, the object of this analysis is to calculate U , which appears only in the cos θ and sin θ terms in the above equations. Hence we need to use only the n = 1 terms in the equations; for example the relevant contribution to f θ is F 1 = 3π/8. A simple calculation gives the result where α = a 0 /a. It remains to perform a similar analysis for the swirl velocity u φ . The solution of the azimuthal component of the Stokes equation in region 1 is (A 11) Similar equations apply to region 2, except the the H (2) n terms are all zero because the swirl velocity must tend to zero at infinity. Moreover, the torque on the body is proportional to G (1) 1 , so this too must be zero. The boundary conditions are that u φ is a sin θΩ at r = a and continuous at r = a 0 , while the jump in azimuthal shear stress at r = a 0 is f φ F 1 . Hence we deduce that This completes the derivation of equations (2.1).
In the model analysed by Drescher (2010), the angle θ (m) between the orientation vector of colony m and the vertical is taken to be small, as is the angle ψ between r and the vertical.       projected in the e1 − e2 plane. G bh is varied to 2, 3 and 6. Figure 14. Phase diagram of two Volvox colonies interacting near a bottom wall (Fg = 7.5π and 9π). The black circle indicates 'unstable motion', in which the center to center 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 e1 − e2 plane is less than 0.3a during t = 90 − 100. Figure 15. Trajectories of two colonies near a bottom wall for time t in the range 0 − 1000 or until center-to-center distance exceeds 10a. Trajectories of a colony with Fg = 6.5π start from the black circles and end at the white circles. Trajectories of a colony with Fg = 9π 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 center-to-center distance of two colonies in (a), (b) and (c). 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 θ and f φ are applied to the fluid at radius (1 + ε)a. Region 1 is defined as a < r < (1 + ε)a, whereas region 2 is defined as (1 + ε)a < r.