Energetics of synchronisation for model flagella and cilia
EEnergetics of synchronisation for model flagella and cilia
Weida Liao and Eric Lauga ∗ Department of Applied Mathematics and Theoretical Physics,University of Cambridge, Wilberforce Road, CB3 0WA, Cambridge, UK (Dated: January 18, 2021)Synchronisation is often observed in the swimming of flagellated cells, either for multiple ap-pendages on the same organism or between the flagella of nearby cells. Beating cilia are also seento synchronise their dynamics. In 1951, Taylor showed that the observed in-phase beating of theflagella of co-swimming spermatozoa was consistent with minimisation of the energy dissipated inthe surrounding fluid. Here we revisit Taylor’s hypothesis for three models of flagella and cilia: (1)Taylor’s waving sheets with both longitudinal and transverse modes, as relevant for flexible flagella;(2) spheres orbiting above a no-slip surface to model interacting flexible cilia; and (3) whirling rodsabove a no-slip surface to address the interaction of nodal cilia. By calculating the flow fields explic-itly, we show that the rate of working of the model flagella or cilia is minimised in our three modelsfor (1) a phase difference depending on the separation of the sheets and precise waving kinematics;(2) in-phase or opposite-phase motion depending on the relative position and orientation of thespheres; and (3) in-phase whirling of the rods. These results will be useful in future models probingthe dynamics of synchronisation in these setups. ∗ [email protected] a r X i v : . [ phy s i c s . b i o - ph ] J a n I. INTRODUCTION
The majority of cells able to move in fluids do so with the aid of slender cellular appendages called flagella or cilia [1],whose time-varying motion generates hydrodynamic stresses that propel the cells forward [2–6]. Although differentorganisms employ different methods, the same general physical principles of swimming at low Reynolds number applyto the locomotion of all cells, from simple bacteria [7] to spermatozoa [8, 9] and higher aquatic microorganisms [10, 11].One of the peculiar features of swimming cells is their ability to synchronise the motion of their appendages. Indeed,synchronisation is observed to take place either for organisms equipped with multiple flagella or cilia, or between theflagella of nearby cells [12]. Such synchronisation occurs, for example, for bacteria equipped with multiple helicalflagella [13–18]. Spermatozoa swimming close to one another are also known to cooperate [19] and synchronise theirwaving motion [20, 21], even though each cell is independently actuated. A unicellular alga can synchronise its twoflagella [22], while hydrodynamic interactions also lead to synchronised dynamics for appendages on different cells [23].For larger ciliated cells [24], or for multicellular aquatic organisms equipped with many short flagella [25], the synchro-nisation dynamics takes the form of metachronal waves, akin to spectator waves in sport stadiums. Synchronisationalso extends beyond swimming to the dynamics of cilia on biological tissues. Nodal cilia, which rotate rigidly aroundtilted conical orbits and produce left–right asymmetry in embryos, are a famous example of this [26, 27].While the synchronised dynamics can often be described mathematically within the framework of coupled oscil-lators [28], a key question lies in identifying the physical ingredients leading to synchronisation. In order for twooscillators to synchronise, it is intuitive that two features are required. First, there needs to be a physical means forthe oscillators to communicate. It has long been suspected that this is generically achieved for swimming cells throughhydrodynamic interactions [29], though “dry” synchronisation is also possible via coupling to the cell body [30, 31] orintracellular coupling [32, 33].The second required ingredient is a physical mechanism that allows the phase of the oscillators to evolve. Inthe context of cellular synchronisation via hydrodynamic interactions, two such mechanisms have been identifiedtheoretically. In the first one, it is the elastic compliance of the orbits that allows oscillators to speed up or slow downin response to hydrodynamic flows [34]. The second generic mechanism is phase-dependent cellular forcing, so thatflagella or cilia are able to change their phase dynamics in response to external hydrodynamic stresses [35].One interesting aspect of cellular synchronisation concerns its energetic consequences. The degrees of freedomavailable to interacting flagella and cilia could be used to minimise energetic costs, for example, the dissipation ofmechanical energy in the surrounding viscous fluid. This was the argument originally examined by Taylor to explainthe in-phase synchronisation of nearby spermatozoa [36]. While this energetic view has proven to be popular inaddressing the collective dynamics of cilia [37] and metachronal waves [38], theoretical calculations have shown thatthe dynamics of swimmers can cause cells to synchronise into a state where the energy dissipation is maximum [39].The states of minimum energy also depend critically on the relative position and orientation of the flagella [40].In this paper, we study the relationship between synchronised dynamics and the dissipation in the surrounding fluidusing three classical models for flagella and cilia. We start in Sec. II by revisiting Taylor’s two-dimensional swimmingsheet model, relevant to the synchronisation of spermatozoa [36, 39, 41–43], and extending it to the general case offlexible swimmers undergoing both longitudinal and transverse waving. In this case, we show that the minimum ofenergy dissipation does not necessarily occur at the in-phase configuration but at an optimal phase difference betweenthe two sheets that depends on the detailed kinematics of the sheets. We further use a long-wavelength calculationto provide physical intuition in the case of swimmers with purely longitudinal waves. In Sec. III, we then considerorbiting spheres, a minimal model used to address the dynamics of interacting cilia in past theoretical [34, 44–47] andexperimental work [48, 49]. We find that depending on the orientations and relative positions of the orbits, either in-phase or opposite-phase motion of the spheres minimises energy dissipated. In Sec. IV, we finally consider interactingnodal cilia [26] modelled by rods whirling above a no-slip surface [27, 50]. In this case, the energy is seen to be alwaysminimised at the in-phase configuration. We finish in Sec. V with a discussion of our results in the context of thedynamics of biological synchronisation. Our work, able to characterise analytically the states of minimum dissipatedenergy in these models, will be useful for future work investigating the complex synchronised dynamics of realisticflagella and cilia.
II. MODEL FOR INTERACTING FLAGELLA: TWO WAVING SHEETSA. Setup
In this first section, we consider the hydrodynamics of two waving sheets as models for two identical flagellatedspermatozoa, which interact through a viscous fluid. The setup is illustrated in Fig. 1 in dimensionless form. Thesheets are infinite and two-dimensional, as in Taylor’s original article [36]. The two sheets are separated by a mean
FIG. 1. Two identical waving sheets as models for two flagellated spermatozoa, which interact hydrodynamically. The sheets,of dimensionless period 2 π , are separated by a mean dimensionless distance d . The phase difference between the sheets isdenoted by ∆. The waving motion includes both longitudinal and transverse modes, so that a material point traces an ellipse. distance h , the first at the unperturbed position y = 0 and the second at unperturbed position y = h . The twosheets undergo identical prescribed waving motion but with an imposed constant phase difference ∆. The wavingmotion includes both longitudinal and transverse deformation modes with amplitudes A and B respectively. Foreach swimmer, we denote by φ a phase difference between these two modes, where φ = 0 means that the materialpoints move in ellipses with semi-axes parallel to the x and y axes. The shape of each sheet is periodic in space withwavenumber k and in time with angular frequency ω . The material points on the first sheet are given by (cid:16) x (1)s , y (1)s (cid:17) ,with similar notation for the second sheet.To proceed with the calculation, we will consider the limit of small-amplitude waving and solve for the leading-orderrate at which work is done by the sheets on the fluid between them, as a function of the phase difference ∆. (Thework done below sheet 1 and above sheet 2 is independent of the phase difference and therefore we do not need tocompute it.) We will compute the optimal phase difference ∆ ∗ , which minimises this rate of working with respectto ∆. Taylor [36] previously found that if the sheets undergo only transverse deformation ( A = 0), then the rate ofworking at leading order is always minimised by in-phase waving of the sheets (∆ ∗ = 0). In our paper, we allow thesheets to be flexible and thus undergo both longitudinal and transverse deformation; that is, both A and B may benonzero. For this setup, we then find ∆ ∗ as a function of the amplitudes, mean sheet separation and phase differencebetween the modes.We impose the kinematics of material points with unperturbed position x on each sheet at time t as x (1)s = x + A cos( kx − ωt − φ ) , (1) y (1)s = B sin( kx − ωt ) , (2) x (2)s = x + A cos( kx − ωt + ∆ − φ ) , (3) y (2)s = h + B sin( kx − ωt + ∆) . (4)Thus, a material point traces out an ellipse in time. In the case of longitudinal oscillation only or transverse oscillationonly of each material point, this ellipse becomes a line segment. If there is no longitudinal oscillation, then the materialpoint has position x s = x .We note that the phase difference ∆ is defined only up to a multiple of 2 π . In our later discussion, we will refer tovalues of ∆ that satisfy − π ≤ ∆ < π . Thus, if ∆ is small and positive, then the phase of the second sheet is slightlyahead of that of the first sheet; that is, the peaks of the second sheet are slightly to the left of those of the first sheet.Conversely, if ∆ is small and negative, then the phase of the second sheet is slightly behind that of the first sheet,and the peaks of the second sheet are slightly to the right of those of the first sheet.The fluid between the sheets is assumed to be Newtonian and the Reynolds number sufficiently small that thegoverning equations are the incompressible Stokes equations ∇ · σ = µ ∇ u − ∇ p = , (5) ∇ · u = 0 , (6)where σ is the stress tensor, u is the fluid velocity field, p is the dynamic pressure and µ is the dynamic viscosity [51].The stress tensor σ is given by σ = − p + µ (cid:2) ∇ u + ( ∇ u ) T (cid:3) , where is the identity tensor.We nondimensionalise the problem using k − and ω − as characteristic length and time scales, so that the pressurescale is µω . The sheets are then separated by a mean dimensionless distance d ≡ kh (see Fig. 1). In preparationfor calculation of the flow in the small-amplitude limit, we define the small parameter (cid:15) ≡ (cid:2) ( Ak ) + ( Bk ) (cid:3) , so thatwe have dimensionless longitudinal and transverse amplitudes a ≡ Ak/(cid:15) and b ≡ Bk/(cid:15) respectively. We now view a and b as independent of (cid:15) and let x, y, t, u , p and σ denote their corresponding dimensionless quantities. Then thedimensionless kinematics for the sheets become x (1)s = x + (cid:15)a cos( x − t − φ ) , (7) y (1)s = (cid:15)b sin( x − t ) , (8) x (2)s = x + (cid:15)a cos( x − t + ∆ − φ ) , (9) y (2)s = d + (cid:15)b sin( x − t + ∆) . (10)Note that we use the shared swimming frame of the two sheets, which has coordinates ( x, y ). This is possible at order (cid:15) by the following symmetry argument. A replacement (cid:15) (cid:55)→ − (cid:15) is equivalent to a translation of the setup by halfa wavelength in the x direction. Such a translation does not affect the relative translational velocity of the sheets,because the sheets are infinite and periodic in space. Hence, the relative translational velocity is even in (cid:15) . Therefore,the order (cid:15) relative translational velocity of the sheets is zero.The no-slip velocity boundary conditions for flow on the sheets are given by u x | (cid:16) x (1)s ,y (1)s (cid:17) = (cid:15)a sin( x − t − φ ) , (11) u y | (cid:16) x (1)s ,y (1)s (cid:17) = − (cid:15)b cos( x − t ) , (12) u x | (cid:16) x (2)s ,y (2)s (cid:17) = (cid:15)a sin( x − t + ∆ − φ ) , (13) u y | (cid:16) x (2)s ,y (2)s (cid:17) = − (cid:15)b cos( x − t + ∆) . (14)Since the problem is two-dimensional, we may introduce a stream function ψ , which enforces incompressibility ofthe flow, with u x = ∂ψ∂y , (15) u y = − ∂ψ∂x . (16)The Stokes equations imply classically that the stream function ψ is a solution to the biharmonic equation [52] ∇ ψ = 0 . (17)Written using this stream function, the velocity boundary conditions for flow on the sheets then become ∂ψ∂x (cid:12)(cid:12)(cid:12)(cid:12) (cid:16) x (1)s ,y (1)s (cid:17) = (cid:15)b cos( x − t ) , (18) ∂ψ∂y (cid:12)(cid:12)(cid:12)(cid:12) (cid:16) x (1)s ,y (1)s (cid:17) = (cid:15)a sin( x − t − φ ) , (19) ∂ψ∂x (cid:12)(cid:12)(cid:12)(cid:12) (cid:16) x (2)s ,y (2)s (cid:17) = (cid:15)b cos( x − t + ∆) , (20) ∂ψ∂y (cid:12)(cid:12)(cid:12)(cid:12) (cid:16) x (2)s ,y (2)s (cid:17) = (cid:15)a sin( x − t + ∆ − φ ) . (21)These boundary conditions are nonlinear in (cid:15) . We therefore use an asymptotic approach and consider in what followsthe limit of small amplitude. Since there is no flow if (cid:15) = 0, we may use the following regular perturbation expansionsin (cid:15) of the stream function ψ , fluid velocity field u , pressure p and stress tensor σ ψ = (cid:15)ψ + O (cid:0) (cid:15) (cid:1) , (22) u = (cid:15) u + O (cid:0) (cid:15) (cid:1) , (23) p = (cid:15)p + O (cid:0) (cid:15) (cid:1) , (24) σ = (cid:15) σ + O (cid:0) (cid:15) (cid:1) . (25)Using a Taylor expansion of the velocity boundary conditions, we thus obtain the boundary conditions for u on bothswimmers as u ,x | y =0 = a sin( x − t − φ ) , (26) u ,y | y =0 = − b cos( x − t ) , (27) u ,x | y = d = a sin( x − t + ∆ − φ ) , (28) u ,y | y = d = − b cos( x − t + ∆) . (29)Note that these order (cid:15) boundary conditions are now applied at the unperturbed positions of the two sheets, whichno longer involve the value of (cid:15) . B. Flow at order (cid:15)
The first-order stream function ψ in the fluid between the two sheets satisfies the biharmonic equation ∇ ψ = 0 . (30)Since the Stokes equations are linear, we may use complex notation for ψ and implicitly take real parts. Then theboundary conditions at order (cid:15) in Eqs. (26)–(29) become ∂ψ ∂x (cid:12)(cid:12)(cid:12)(cid:12) ( x, = b exp[ i ( x − t )] , (31) ∂ψ ∂y (cid:12)(cid:12)(cid:12)(cid:12) ( x, = − ai exp[ i ( x − t − φ )] , (32) ∂ψ ∂x (cid:12)(cid:12)(cid:12)(cid:12) ( x,d ) = b exp[ i ( x − t + ∆)] , (33) ∂ψ ∂y (cid:12)(cid:12)(cid:12)(cid:12) ( x,d ) = − ai exp[ i ( x − t + ∆ − φ )] . (34)The general unit-speed, 2 π -periodic solution to the biharmonic equation is ψ = Ky + Ly + M y + (cid:88) n ≥ { ( G n + yH n ) cosh[ n ( y − d )] + ( I n + yJ n ) sinh[ n ( y − d )] } exp[ in ( x − t )] , (35)where the capital letters denote constants to be determined and we have omitted an arbitrary, physically irrelevantadditive constant [6]. Since we impose zero mean pressure gradient in the x direction, it is necessary to have M = 0.Similarly, asking for zero mean shear acting on the swimmers means that we also need L = 0. Now, by comparingwith the boundary conditions in Eqs. (31)–(34), we see that K = 0 and that only the n = 1 mode in the infinite sumis nonzero. Dropping the subscript 1 from the n = 1 mode coefficients, we write down the relevant solution for ψ as ψ = [( G + yH ) cosh( y − d ) + ( I + yJ ) sinh( y − d )] exp[ i ( x − t )] . (36)Substituting this into the boundary conditions, we obtain a system of simultaneous equations i ( G cosh d − I sinh d ) = b, (37)( H + I ) cosh d − ( G + J ) sinh d = − ai exp( − iφ ) , (38) i ( G + dH ) = b exp( i ∆) , (39) H + I + dJ = − ai exp[ i (∆ − φ )] . (40) C. Leading-order rate of working
We now consider the energetics of the problem and calculate the rate at which work is done ˙ W by the sheets, atleading order in the waving amplitude, on the fluid between them (per wavelength). Since the rate of working variesquadratically with the flow velocity, its leading-order value occurs at order (cid:15) and we may use the regular perturbationexpansion ˙ W = (cid:15) ˙ W + O (cid:0) (cid:15) (cid:1) . (41)Thus, we only need the stream function to order (cid:15) . Note that work is also done on the fluid not between the sheets [4, 6].However, as noted earlier, this does not depend on the phase difference ∆, so is irrelevant to minimisation with respectto ∆ of the total rate of working.The leading-order rate of working ˙ W per wavelength is given by˙ W = − (cid:90) π u ,x σ ,xy | y =0 d x − (cid:90) π u ,y σ ,yy | y =0 d x + (cid:90) π u ,x σ ,xy | y = d d x + (cid:90) π u ,y σ ,yy | y = d d x, (42)where the relevant components of the stress tensor σ are σ ,xy = ∂u ,x ∂y + ∂u ,y ∂x , (43) σ ,yy = − p + 2 ∂u ,y ∂y . (44)We find the leading-order pressure by integrating the x component of the Stokes equations at order (cid:15)∂p ∂x = ∇ ∂ψ ∂y , (45)which gives p = − i [ H cosh( y − d ) + J sinh( y − d )] exp[ i ( x − t )] . (46)We may then use the boundary conditions on ψ in Eqs. (31)–(34) to obtain components of the stress tensor incomplex notation σ ,xy | y =0 =[( G + 2 J ) cosh d − (2 H + I ) sinh d − ib ] exp[ i ( x − t )] , (47) σ ,yy | y =0 = [2 i ( H cosh d − J sinh d ) − a exp( − iφ )] exp[ i ( x − t )] , (48) σ ,xy | y = d = [ G + dH + 2 J − ib exp( i ∆)] exp[ i ( x − t )] , (49) σ ,yy | y = d = { iH − a exp[ i (∆ − φ )] } exp[ i ( x − t )] . (50)Using the real and imaginary parts of the simultaneous equations in Eqs. (37)–(40), we may eliminate G and I toobtain in real notation σ ,xy | y =0 =[Re( G + 2 J ) cosh d − Re(2 H + I ) sinh d ] cos( x − t )+ [Im(2 H + I ) sinh d − Im( G + 2 J ) cosh d + b ] sin( x − t )=2[Re( J ) cosh d − Re( H ) sinh d ] cos( x − t )+ 2[ b − Im( J ) cosh d + Im( H ) sinh d ] sin( x − t ) , (51) σ ,yy | y =0 =2[Re( J ) sinh d − Re( H ) cosh d ] sin( x − t )+ 2[Im( J ) sinh d − Im( H ) cosh d ] cos( x − t ) − a cos( x − t − φ ) , (52) σ ,xy | y = d = Re( G + dH + 2 J ) cos( x − t ) − Im( G + dH + 2 J ) sin( x − t ) + b sin( x − t + ∆)=[ b sin ∆ + 2 Re( J )] cos( x − t ) + [ b cos ∆ − J )] sin( x − t ) + b sin( x − t + ∆) , (53) σ ,yy | y = d = − H ) sin( x − t ) − H ) cos( x − t ) − a cos( x − t + ∆ − φ ) . (54)Solving Eqs. (37)–(40) for H and J gives H = − i sinh d − d (cid:8) a sinh d exp[ i (∆ − φ )] + ad sinh d exp( − iφ ) − bd exp( i ∆) − b sinh d cosh d exp( i ∆) + b sinh d + bd cosh d (cid:9) , (55) J = − i sinh d − d (cid:8) − ad exp[ i (∆ − φ )] + a sinh d cosh d exp[ i (∆ − φ )] − a sinh d exp( − iφ )+ ad cosh d exp( − iφ ) − b sinh d exp( i ∆) + bd sinh d (cid:9) . (56)Then we find the four contributions to the rate of working per wavelength as (cid:90) π u ,x σ ,xy | y =0 d x = 2 π (cid:8) − a [sinh d cosh d − d + ( d cosh d − sinh d ) cos ∆] − abd cos φ + abd sinh d cos(∆ + φ ) (cid:9) sinh d − d , (57) (cid:90) π u ,y σ ,yy | y =0 d x = 2 π (cid:8) − b [sinh d cosh d + d − ( d cosh d + sinh d ) cos ∆] − abd cos φ − abd sinh d cos(∆ − φ ) (cid:9) sinh d − d , (58) (cid:90) π u ,x σ ,xy | y = d d x = 2 π (cid:8) a [sinh d cosh d − d + ( d cosh d − sinh d ) cos ∆] − abd cos φ + abd sinh d cos(∆ − φ ) (cid:9) sinh d − d , (59) (cid:90) π u ,y σ ,yy | y = d d x = 2 π (cid:8) b [sinh d cosh d + d − ( d cosh d + sinh d ) cos ∆] − abd cos φ − abd sinh d cos(∆ + φ ) (cid:9) sinh d − d . (60)Note that, due to the general lack of 1 ↔ W from the first sheet is not equal to the contribution from the second. (If one of thewaving modes disappears, then both sheets give equal contributions.) The total leading-order rate at which work isdone by the sheets on the fluid between them per wavelength is finally obtained as˙ W = 4 π sinh d − d (cid:8) a [sinh d cosh d − d + ( d cosh d − sinh d ) cos ∆]+ b [sinh d cosh d + d − ( d cosh d + sinh d ) cos ∆]+ 2 abd sinh d sin ∆ sin φ (cid:9) . (61) D. Minimisation of leading-order rate of working
We now consider the minimisation of the leading-order rate of working ˙ W in Eq. (61) with respect to the phasedifference ∆ between the waving sheets. Classically, the rate at which work is done by the sheets is equal to the rateof dissipation of mechanical energy in the fluid between the sheets. Thus, minimising ˙ W is equivalent to findingthe optimal phase difference of synchronised sheets, which we denote by ∆ ∗ , leading to minimum dissipation (in ourdiscussion we choose − π ≤ ∆ ∗ < π ).First, we may find by inspection the value of ∆ ∗ for the two special cases a = 0 (transverse waving only), whichwas considered by Taylor, and b = 0 (longitudinal waving only). To do this, we note that for d >
0, the quantitiessinh d cosh d ± d , d cosh d ± sinh d and sinh d − d , which appear in Eq. (61), are all strictly positive. Thus, we recoverTaylor’s result: if the sheets undergo only transverse deformation ( a = 0), then the in-phase waving of the two sheets(∆ ∗ = 0) minimises the rate of working ˙ W . Conversely, if the sheets undergo only longitudinal deformation ( b = 0),then, perhaps surprisingly, the opposite is true. The opposite-phase waving of the two sheets (∆ ∗ = − π ) minimises˙ W , while in-phase waving maximises it. This case will be interpreted further in Sec. II E using a long-wavelengthlimit of this calculation. FIG. 2. Variation of the optimal phase difference ∆ ∗ with the phase difference between the longitudinal and transversedeformation modes φ (left panel) and with the mean sheet separation d (right panel). Both panels are in the case a/b = 0 . d (left) and φ (right). Next we derive the value of the optimal phase difference ∆ ∗ for general values of a and b . We consider the part ofthe expression for the rate of working ˙ W that varies with ∆ and define S ≡ (cid:0) a − b (cid:1) d cosh d − (cid:0) a + b (cid:1) sinh d, (62) T ≡ abd sinh d sin φ. (63)These do not depend on the phase difference ∆ between the two sheets, but are functions of the waving amplitudes a and b , the mean sheet separation d and the phase difference between deformation modes φ . Then ˙ W as a functionof ∆ is minimised when the sum S cos ∆ + T sin ∆ is minimised. This is equivalent to minimising cos(∆ − Θ), wherecos Θ = S ( S + T ) , (64)sin Θ = T ( S + T ) . (65)The optimal phase difference ∆ ∗ is therefore given by∆ ∗ = Θ − π. (66)Conversely, we see that ˙ W is maximised when ∆ = Θ, which differs from ∆ ∗ by π .We illustrate the value of the optimal phase difference ∆ ∗ in Figs. 2, 3 and 4 for three different values of theamplitude ratio a/b : a/b = 0 . ∗ as a function of the wavingphase φ , for different values of the separation d . In these plots, we restrict values of φ so that 0 ≤ φ ≤ π/
2. Thisavoids redundancy, since T is proportional to sin φ , which has symmetries. Note that if φ = 0, then T = 0. Thismeans that ∆ ∗ = − π or 0 depending on whether S > < φ increases from 0 to π/ d fixed, T increases, while S remains fixed. Thus, the optimal phase difference ∆ ∗ also varies monotonicallyas φ increases from 0 to π/
2, and whether it increases from ∆ ∗ = − π or decreases from ∆ ∗ = 0 depends on the signof S . For example, for a/b = 10 with d = 0 .
2, we have
S <
0, which results in ∆ ∗ = 0 for φ = 0. This reflects thelimit of small d , which we discuss below. However, for a/b = 10 with the larger values of d = 0 . , . S >
0, which results in ∆ ∗ = − π for φ = 0. This reflects the fact that in these cases the waving is mostly longitudinal;as found earlier, the special case b = 0 has ∆ ∗ = − π .On the right panel of each figure, we next plot how ∆ ∗ varies with the mean sheet separation d , for selected valuesof the phase difference φ between longitudinal and transverse waving modes. In that case, we can examine analytically FIG. 3. Same as in Fig. 2 in the case a/b = 1, i.e., equal longitudinal and transverse waving amplitudes.
FIG. 4. Same as in Fig. 2 in the case a/b = 10, i.e., mostly longitudinal waving. the two limits d → d → ∞ . In the limitas d → W in powers of d and find˙ W = 4 π d (cid:2) a d (2 + cos ∆) + 6 b (cid:0) d (cid:1) (1 − cos ∆) + 30 abd sin ∆ sin φ + O (cid:0) d (cid:1)(cid:3) . (67)In particular, the leading-order result for b (cid:54) = 0 is that in-phase waving always minimises the rate of working. This isclearly reflected in the dependence of ∆ ∗ on d in Figs. 2, 3 and 4 (right panels). We also see that in general, if thesheets only deform longitudinally ( b = 0), then the rate of working is order 1 /d , but otherwise ( b (cid:54) = 0), the rate ofworking is order 1 /d . Thus, contributions to the rate of working due to longitudinal deformation are much smallerthan those due to transverse deformation.In the opposite limit where d → ∞ (i.e., the limit where the sheets are widely separated), we need to distinguishthe cases a (cid:54) = b and a = b . If a (cid:54) = b , then S ∼ (cid:0) a − b (cid:1) d cosh d and T ∼ abd sinh d sin φ as d → ∞ , so that tan ∆ ∗ → ab sin φ/ (cid:0) a − b (cid:1) as d → ∞ . In the case a = b , we instead have S = − (cid:0) a + b (cid:1) sinh d and T = 2 abd sinh d sin φ .Then we have cot ∆ ∗ = − (cid:0) a + b (cid:1) / (2 abd sin φ ) → d → ∞ , unless we have φ = 0, a case that was treated earlier.0These results are reflected in the plots.We further note that the two plots for a/b = 0 . ∗ is close to zero. This is therefore the relevant limit for inextensible flagella [36]. In otherwords, the leading-order rate of working is minimised when the sheets wave with a small phase difference betweenthem. This result is reminiscent of the small, but nonzero, phase differences observed in cilia arrays deforming asmetachronal waves. E. Longitudinal waving: long-wavelength limit and interpretation
We calculated above the rate at which work was expended by the swimmers with both longitudinal and transversemodes in the limit of small-amplitude waving. In particular, we found in Sec. II D that in the case of only longitudinalwaving, the optimal phase difference between the two sheets was ∆ ∗ = − π , i.e., opposite-phase synchronisation. Thisresult is the opposite of what was obtained by Taylor in the case of transverse waving [36]. To investigate this further,we consider here the long-wavelength limit, and allow only small-amplitude, longitudinal deformation of the sheets( b = 0). In the long-wavelength limit, the mean separation of the swimmers is much smaller than their wavelength, sothat the dimensionless mean sheet separation satisfies d (cid:28)
1. Therefore, we may use the classical lubrication theoryof hydrodynamics to find the rate at which work is done, ˙ W L , by the swimmers on the fluid between them. We willthen compare our result to the limit as d → (cid:15) ˙ W previously found for small-amplitude waving, in Eq. (67), whichwill allow us to gain physical intuition on the case of longitudinal waving.As in the earlier setup of Sec. II A, we impose the dimensionless kinematics of material points on the sheets, givenby Eqs. (7)–(10) with b = 0 in the shared swimming frame with coordinates ( x, y ). This results in velocity boundaryconditions correct to order (cid:15) given by Eqs. (26)–(29). Since for this lubrication theory calculation we consider onlylongitudinal deformation of the sheets, we may set the phase difference between the longitudinal and transverse wavingmodes φ to zero without loss of generality. Then the kinematics become x (1) s = x + (cid:15)a cos( x − t ) , (68) y (1) s = 0 , (69) x (2) s = x + (cid:15)a cos( x − t + ∆) , (70) y (2) s = d. (71)We now move from the swimming frame to the wave frame, which travels at unit speed in the x direction relative tothe swimming frame; that is, the wave frame has coordinates (˜ x, y ), where ˜ x ≡ x − t . In the wave frame, the velocityat each ˜ x on the sheets is constant in time; physically, the wave frame follows the compressions and extensions asthey travel along the sheet.In the long-wavelength limit, the incompressible Stokes equations become the lubrication equations [53] ∂p∂ ˜ x = ∂ u ˜ x ∂y , (72) ∂p∂y = 0 , (73) ∇ · u = 0 . (74)The no-slip velocity boundary conditions correct to order (cid:15) are given by u ˜ x | y =0 = (cid:15)a sin ˜ x − , (75) u y | y =0 = 0 , (76) u ˜ x | y = d = (cid:15)a sin(˜ x + ∆) − , (77) u y | y = d = 0 . (78)Since the pressure p is a function of ˜ x only, we may directly integrate the ˜ x component of the lubrication equations.Then the horizontal component of the velocity, correct to order (cid:15) , is given by u ˜ x = 12 ∂p∂ ˜ x y ( y − d ) + (cid:15)ayd [sin(˜ x + ∆) − sin ˜ x ] + (cid:15)a sin ˜ x − . (79)1This flow is clearly a sum of a pressure-driven flow (quadratic in y ), a shear flow (linear) and uniform flow (constant).The flow rate between the sheets Q is then a function of ˜ x given by Q ≡ (cid:90) d u ˜ x d y = − d ∂p∂ ˜ x + (cid:15)ad x + ∆) + sin ˜ x ] − d. (80)Using the incompressibility condition and the boundary conditions, we find the derivative of the flow rate ∂Q∂ ˜ x = (cid:90) d ∂u ˜ x ∂ ˜ x d y = − (cid:90) d ∂u y ∂y d y = u y | y =0 − u y | y = d = 0 . (81)Therefore, the flow rate Q is constant along the sheet. Rearranging the expression for Q gives the pressure gradient ∂p∂ ˜ x = 6 d { (cid:15)ad [sin(˜ x + ∆) + sin ˜ x ] − d − Q } . (82)To eliminate the constant Q , we impose the condition that pressure is periodic via (cid:90) π ∂p∂ ˜ x d˜ x = 0 . (83)This gives Q = − d , so that the pressure gradient is obtained explicitly as ∂p∂ ˜ x = 6 (cid:15)ad [sin(˜ x + ∆) + sin ˜ x ] . (84)In the long-wavelength limit, the rate of viscous dissipation in the fluid between the sheets, which is equal to the rateat which work is done by the sheets on the fluid between them, ˙ W L , is given by˙ W L = (cid:90) π (cid:90) d (cid:18) ∂u ˜ x ∂y (cid:19) d y d˜ x. (85)The shear rate is given by ∂u ˜ x ∂y = 12 ∂p∂ ˜ x (2 y − d ) + (cid:15)ad [sin(˜ x + ∆) − sin ˜ x ] , (86)so that the rate of working, correct to order (cid:15) , is finally obtained as˙ W L = 4 π(cid:15) a d (2 + cos ∆) . (87)Earlier, we calculated the order (cid:15) rate of working for small-amplitude waving, and found the expression in Eq. (67)for ˙ W , valid for small mean sheet separation d . The expression for ˙ W L in Eq. (87) agrees with (cid:15) ˙ W from Eq. (67)for b = 0 (longitudinal waving only) to leading order. Therefore, in the lubrication limit, we recover the result for d (cid:28) b = 0 that opposite-phase waving minimises the rate at which work is done by the sheets.To help explain intuitively why the value of ∆ ∗ = − π leads to minimum energy dissipation for synchronisedlongitudinal waving, we sketch in Fig. 5 the flow fields for in-phase (top) and opposite-phase (bottom) waving of thesheets in the swimming frame, so that we may see the order (cid:15) flow clearly. From the results above, we note that thepressure gradient given by Eq. (84) may be written as ∂p∂ ˜ x = 12 (cid:15)ad cos (cid:18) ∆2 (cid:19) sin (cid:18) ˜ x + ∆2 (cid:19) . (88)2 FIG. 5. Sketch of flow fields for in-phase (top) and opposite-phase (bottom) small-amplitude, longitudinal deformation ofsheets in the swimming frame ( x, y ). Arrows indicate direction of motion of the sheets (dark arrows) and the fluid (lightarrows). In-phase waving of the sheets induces backflows halfway between the swimmers and increases the total rate of viscousdissipation.
The periodic pressure gradient has an amplitude proportional to | cos(∆ / | . In particular, the pressure gradient is zeroat the optimal phase difference ∆ ∗ = − π between the sheets, for which there is no quadratic part of the flow in Eq. (79).In that case, in the swimming frame, the fluid circulates in cells of width π and height d between the two sheets ina linear, shearlike manner. If instead the sheets wave in phase (∆ = 0), then in the swimming frame, backflows areinduced (relative to the waving of the swimmers at each x ) halfway between the sheets due to incompressibility. Thefluid circulates in smaller cells of width π and smaller height d/ III. MODEL FOR INTERACTING CILIA: TWO SPHERES IN ELLIPTICAL ORBITS
After considering a model for two interacting flexible flagella in the previous section, we now investigate theenergetics of synchronised cilia interacting through the fluid.Following a classical modelling approach, we represent the anchored cilia as rigid spheres immersed in a fluid andorbiting periodically above a no-slip surface. We consider two such identical spheres that interact hydrodynamicallyin the far field. The flow outside the spheres is governed by the incompressible Stokes equations for a fluid of viscosity µ , Eqs. (5)–(6), and the far-field flow generated by the motion of a sphere will be approximated as that due to apoint force (stokeslet) [51]. Each sphere undergoes periodic motion in an elliptical orbit with a constant imposedphase difference ∆. An elliptical orbit is a minimal model for the two-stroke motion of a flexible cilium, capturing theessential time-irreversibility that allows it to generate net forces and flows. Here, similarly to the previous section,we impose the kinematics of each sphere and calculate their time-averaged rate of working as a function of the phasedifference ∆. We will consider two different cases. In Sec. III A, the orbits lie in a plane parallel to the no-slip surfaceand are circular, similar to the setup in Ref. [34]. Next, in Sec. III B, each sphere’s orbit lies in a plane perpendicularto the no-slip surface and is elliptical, a motion akin to that of the centres of mass of flexible cilia during their effectiveand recovery strokes.3 FIG. 6. Two identical spheres interacting hydrodynamically are models for two anchored cilia. Each sphere is of radius a andmoves in a circular orbit of radius R . The orbits lie in a plane a distance h above a parallel no-slip surface and their centresare separated by distance (cid:96) . FIG. 7. View from above of each orbit in Fig. 6. A. Circular orbits parallel to no-slip surface
1. Setup
The setup is illustrated in Fig. 6, with views of each orbit from above shown in Fig. 7. Each sphere, of radius a ,moves at constant angular speed ω in a circular orbit of radius R , clockwise when viewed from above, with a constantimposed phase difference of ∆. Each orbit lies in the plane a distance h above the no-slip surface. We use Cartesiancoordinates ( x, y, z ) with corresponding unit vectors e x , e y , and e z , where the plane z = 0 is the no-slip surface, andthe centres of the orbits of the spheres are separated by a distance (cid:96) in the x direction. Thus, we may choose thecentres of the orbits of the first and second spheres to have coordinates (0 , , h ) and ( (cid:96), , h ) respectively. The spheresare far from the no-slip surface ( a (cid:28) h ). In order to derive our results analytically, and as done in other work, welet the separation of orbit centres (cid:96) be much larger than all other imposed length scales in the setup ( a, h, R (cid:28) (cid:96) ).This means that the spheres are widely separated at all times and interact hydrodynamically only in the far field.We proceed to calculate the rate at which work is done by each sphere averaged over one period, including the firstcorrection due to interaction between the two spheres, which depends on the phase difference ∆.We impose the kinematics as follows. The position of the first sphere r is a function of time t given by r = R (cos ωt e x − sin ωt e y ) + h e z , (89)so that its velocity is U = d r d t = − ωR (sin ωt e x + cos ωt e y ) . (90)4The position of the second sphere r is given by r = R [cos( ωt + ∆) e x − sin( ωt + ∆) e y ] + (cid:96) e x + h e z , (91)and its velocity is U = d r d t = − ωR [sin( ωt + ∆) e x + cos( ωt + ∆) e y ] . (92)
2. Time-averaged rate of working
From the setup above, we can calculate the time-averaged rate at which work is done by the first sphere in the far-field limit, by using Fax´en’s first law. We denote by u → the flow induced by the motion of the second sphere at theposition of the first sphere and let r → be the position of the first sphere relative to the second (i.e., r → ≡ r − r ).We also denote by F i the force exerted by the i th sphere on the fluid. The flow u → may be approximated as thatdue to a point force F above a no-slip surface, since the spheres are widely separated. Using the classical imagesystem for flow singularities above rigid walls [54], and since the spheres are the same distance h above the no-slipsurface, the flow u → is given to leading order by u → (cid:39) h πµ ( F · r → ) r → | r → | . (93)Approximating the second sphere generating the force F above as moving in an unbounded fluid otherwise at rest,that force is given by Stokes’s law as F (cid:39) πµa U . (94)In other words, we neglect at leading order the impact of the flow due to the first sphere at the position of the secondon the hydrodynamic interactions, since this decays in the limit as the spheres become infinitely separated. We alsoneglect order a/h corrections to the viscous resistance coefficient for the sphere ( (cid:39) πµa ), since the sphere is assumedto be far from the no-slip surface. To leading order we have approximately r → (cid:39) − (cid:96) e x , since the radius of orbit issmall, R (cid:28) (cid:96) . Then the flow u → is given to leading order by u → (cid:39) − ah (cid:96) ωR sin( ωt + ∆) e x . (95)This then allows us to obtain the leading-order correction to the Stokes’s law value for the force F exerted by thefirst sphere on the fluid. This perturbation to the value given by Stokes’s law is due to the flow created by secondsphere in the far field. By Fax´en’s first law, the force F exerted by the first sphere on the fluid is approximatelygiven by F (cid:39) πµa ( U − u → ) (cid:39) πµaωR (cid:18) − (sin ωt e x + cos ωt e y ) + 9 ah (cid:96) sin( ωt + ∆) e x (cid:19) . (96)This gives the instantaneous rate of working F · U (cid:39) πµaω R (cid:18) − ah (cid:96) sin ωt sin( ωt + ∆) (cid:19) . (97)Then the rate at which work is done by the first sphere averaged over a period is obtained as ω π (cid:90) π/ω F · U d t (cid:39) πµaω R (cid:18) − ah (cid:96) cos ∆ (cid:19) . (98)To find the corresponding result for the second sphere, we replace the phase difference ∆ with − ∆ and find thatthe rate of working is unchanged; the total rate at which energy is dissipated in the fluid is thus twice the resultfrom Eq. (98). Consequently, the in-phase synchronised motion of the spheres always minimises the average rate ofdissipation of energy in the fluid. We will compare these results to past predictions for the dynamic synchronisationin Sec. V.5 FIG. 8. Two identical spheres of radius a interacting hydrodynamically as models for two anchored cilia. Each orbit is ellipticaland has centre at height h above the no-slip surface z = 0. Each orbit lies in a plane perpendicular to the no-slip surface. Therelative angle between the planes of the two orbits is α (see top view in Fig. 9). The centre of the second sphere’s orbit isdisplaced from that of the first sphere’s orbit by the vector (cid:96) cos β e x + (cid:96) sin β e y .FIG. 9. View from above of setup in Fig. 8. Each elliptical orbit has one semi-axis of length b lying in the plane z = h . B. Elliptical orbits perpendicular to no-slip surface
1. Setup
In Sec. III A, we considered a model in which the orbits of the spheres are circular and lie in a plane parallel to ano-slip surface. This allowed us to consider the energetic properties of the setup addressed in the classical work ofRef. [34]. In this section, we now allow the orbits to be elliptical and oriented perpendicular to the no-slip surface.Each of the two identical rigid spheres of radius a may be viewed as representing the centre of mass of a flexiblecilium [23]. The elliptical orbit is therefore an approximation to the cilium’s two-stroke periodic motion, with aneffective stroke where the cilium is further from the surface and a recovery stroke where the cilium is closer to thesurface. We again will calculate the rate at which work is done by each sphere as a function of the phase difference∆ between their orbiting motions.The setup is illustrated in Fig. 8, with a view from above in Fig. 9. Each sphere moves in an elliptical orbit withsemi-axes b and c lying in a plane perpendicular to the no-slip surface. The centres of the orbits are separated by adistance (cid:96) and lie a distance h above the no-slip plane, while the relative angle between the two orbital planes is α .The imposed motion of each sphere is identical but with a constant imposed phase difference of ∆. The period of themotion is 2 π/ω , but importantly we note that the frequency ω is not the angular velocity of the spheres, since theorbits are not circular in general. In other words, unless the orbit is circular, the motion does not occur at constantspeed.We again use Cartesian coordinates ( x, y, z ) with corresponding unit vectors e x , e y and e z , where the no-slip surfaceis given by z = 0. We choose the orbit of the first sphere to have centre with coordinates (0 , , h ), and its semi-axis6of length b aligned with the y axis when viewed from above as in Fig. 9. The semi-axes of the orbits have length c parallel to the z axis. We let the centre of the second sphere’s orbit have coordinates ( (cid:96) cos β, (cid:96) sin β, h ). The secondsphere’s orbit, when viewed from above, makes an angle α with the y axis. This setup breaks the symmetry betweenthe two spheres. The spheres are at all times far from the no-slip surface ( a, c (cid:28) h ) and widely separated ( a, b (cid:28) (cid:96) ).The separation (cid:96) of the centres of the orbits is the largest imposed length scale ( (cid:96) (cid:29) h ), allowing us to treat thehydrodynamic interactions between the spheres in the far field.We impose the kinematics as follows. The position of the first sphere r is given by r = b cos ωt e y + c sin ωt e z + h e z , (99)so that its velocity U is U = d r d t = ω ( − b sin ωt e y + c cos ωt e z ) . (100)Similarly, the second sphere has position r given by r = b cos( ωt + ∆)(sin α e x + cos α e y ) + c sin( ωt + ∆) e z + (cid:96) cos β e x + (cid:96) sin β e y + h e z , (101)with velocity U given by U = d r d t = ω [ − b sin( ωt + ∆)(sin α e x + cos α e y ) + c cos( ωt + ∆) e z ] . (102)
2. Time-averaged rate of working
As in Sec. III A, we use Fax´en’s first law to find the time-averaged rate at which work is done by the first sphere inthe far-field limit. We also retain the notation u → for the flow induced by the motion of the second sphere at theposition of the first sphere and r → ≡ r − r for the relative position vector. As above we denote by F i the forceexerted by the i th sphere on the fluid. Since we are in the far field, we may approximate the flow u → as that dueto a point force F . Furthermore, we may use the leading-order expression for u → in Eq. (93), since the orbits aresmall compared with the distance h to the no-slip plane ( a, c (cid:28) h ). As in Eq. (94), the force F is then obtained toleading order as F (cid:39) πµa U , (103)since the spheres are widely separated and far from the no-slip surface. This time, however, we have r → (cid:39) − (cid:96) (cos β e x + sin β e y ) , (104)which gives the leading-order flow u → (cid:39) ah (cid:96) (cid:2)(cid:0) U ,x cos β + U ,y cos β sin β (cid:1) e x + (cid:0) U ,x cos β sin β + U ,y sin β (cid:1) e y (cid:3) = − abωh sin( ωt + ∆) sin( β + α ) (cid:96) (cos β e x + sin β e y ) . (105)As above, we can now find the leading-order correction to the Stokes’s law value for the force F exerted by thefirst sphere, due to the motion of the second sphere. Substituting quantities into Fax´en’s first law gives F (cid:39) πµa ( U − u → ) (cid:39) πµaω (cid:18) − b sin ωt e y + c cos ωt e z + 9 abh sin( ωt + ∆) sin( β + α ) (cid:96) (cos β e x + sin β e y ) (cid:19) (106)and F · U (cid:39) πµaω (cid:18) b sin ωt + c cos ωt − ab h sin ωt sin( ωt + ∆) sin β sin( β + α ) (cid:96) (cid:19) . (107)7 FIG. 10. Diagram illustrating the required changes in angles in the setup of Fig. 9 to derive the rate at which work is done bythe second sphere from that by the first sphere. We replace the angle β with − [( π/ − α ) + ( π/ − β )] = α + β − π , the relativeangle α with − α and the phase difference ∆ with − ∆. Then the rate of work done by the first sphere averaged over a period is obtained as ω π (cid:90) π/ω F · U d t (cid:39) πµaω (cid:18) b + c − ab h sin β sin( β + α ) cos ∆ (cid:96) (cid:19) . (108)To find the corresponding result for the second sphere, we may use a symmetry argument illustrated in Fig. 10 (whichwas not evident a priori since the geometrical setup is not symmetric). If we replace the angle β with α + β − π , therelative angle α with − α and the phase difference ∆ with − ∆, we see that the setup has been switched 1 ↔
2. As aresult, we obtain that the rate of working of the second sphere is equal to the rate of working of the first sphere atthis order.With this, we can now consider the optimisation of the viscous dissipation with respect to the phase difference ∆between the synchronised motion of the two spheres. The sign of the coefficient of cos ∆ determines whether rate ofworking is minimised by in-phase or opposite-phase motion. This depends on the quantity sin β sin( β + α ), whichis unchanged by replacing β with π + β and whose sign is changed by replacing α with π + α . If in-phase motionminimises the rate of working, then opposite-phase motion maximises it, and vice versa. We see that if sin β sin( β + α )is positive, then in-phase motion of the spheres minimises the rate of working; that is, motion where the two spheresare at the same height at all times leads to minimum power dissipated in the fluid. This includes the case of preciselyaligned elliptical orbits ( α = 0). On the other hand, if sin β sin( β + α ) is negative, then opposite-phase motion of thespheres is the optimal motion. These results are illustrated in a plot of the ( α, β ) plane in Fig. 11. We see that if thetwo orbital planes are almost but not precisely aligned ( α close to 0 or 2 π ), then in-phase motion minimises energydissipation for most values of β , while opposite-phase motion leads to the minimum for a narrow range of β near π .The same result (either in-phase or opposite-phase motion is optimal, depending on geometry) was obtained for the8 FIG. 11. Regions in the ( α, β ) plane for the setup in Fig. 8 where in-phase or opposite-phase motion of spheres minimisestime-averaged rate of working. The relative angle between the planes of the orbits is α . The angle β indicates the relativeposition of the centres of the orbits.FIG. 12. Two identical rigid rods of length L as models for two nodal cilia interacting hydrodynamically. Each rod rotatesclockwise when viewed from above, sweeping out a cone with tip fixed on a no-slip surface. The semi-angle of each cone is ψ and each cone has a posterior tilt of θ about the x -axis. The second cone is displaced from the first cone by (cid:96) cos β e x + (cid:96) sin β e y . planar beating of three-dimensional flagella [40]. In Sec. V, we will compare these results to related studies showingthat both in-phase and opposite-phase beating can also arise dynamically, depending on the relative conformation ofmodel cilia [47]. IV. MODEL FOR INTERACTING NODAL CILIA: TWO WHIRLING RODS
In this final model, we consider nodal cilia, which, instead of beating with effective and recovery strokes of differentshapes, rotate rigidly along conical surfaces. We use the whirling rod model for nodal cilia proposed in past studies [50,55] to compute the impact of phase difference on the rate of working of two such interacting, synchronised cilia.
A. Setup
The setup is illustrated in Figs. 12 and 13. We model each nodal cilium as a rigid rod of length L , a simplificationthat neglects slight bending of the cilia due to viscous resistance [26]. The flow outside the rods is governed bythe incompressible Stokes equations for a fluid of viscosity µ , Eqs. (5)–(6). We consider two such rods interactinghydrodynamically in the far field. We impose their kinematics so that each rod rotates with constant angular frequency ω , sweeping out a cone of constant semi-angle ψ above a no-slip surface. The rotation is clockwise when viewed fromabove and the tip of each cone has a fixed position on the no-slip plane. Each cone axis is tilted by a constant angle θ ,called the “posterior tilt”, away from the normal to the no-slip surface; this tilt of nodal cilia, which is a mechanism9 FIG. 13. View of each of the orbits traced by the whirling rods in the setup of Fig. 12, looking down each cone axis towardsthe no-slip surface. for producing unidirectional flow, was confirmed in experimental studies of mouse embryos [56]. So that the rodsdo not make contact with the no-slip surface, we necessarily have the geometrical constraint θ + ψ < π/
2. The twocones are separated by a distance (cid:96) , assumed to be much larger than the rod length L . We impose a constant phasedifference ∆ between the otherwise-identical, synchronised, periodic motion of the two rods, and we calculate belowthe rate of working averaged over one period as a function of the phase difference ∆.We use Cartesian coordinates ( x, y, z ) with the no-slip surface given by z = 0. We choose the first cone to haveits tip at the origin and the second cone to have its tip at ( (cid:96) cos β, (cid:96) sin β, r ( s, t ), where s is the arclength with 0 ≤ s ≤ L and t is time. If the cone had no posteriortilt ( θ = 0), then the cone axis would be aligned with the z axis. To find the cone with posterior tilt θ , we rotate thecone with zero posterior tilt, about the x axis and away from the positive y direction. Thus, the position r ( s, t ) ofthe material point at arclength s along the first rod at time t is given by r ( s, t ) ≡ x ( s, t ) y ( s, t ) z ( s, t ) = θ − sin θ θ cos θ s sin ψ cos ωt − s sin ψ sin ωts cos ψ = s sin ψ cos ωt − s sin ψ sin ωt cos θ − s cos ψ sin θ − s sin ψ sin ωt sin θ + s cos ψ cos θ . (109)The corresponding velocity of this material point U ( s, t ) is then given by a time derivative U ( s, t ) = ∂ r ( s, t ) ∂t = − ωs sin ψ sin ωt − ωs sin ψ cos ωt cos θ − ωs sin ψ cos ωt sin θ . (110)Similarly, the position r ( s, t ) of the material point on the second rod is r ( s, t ) ≡ x ( s, t ) y ( s, t ) z ( s, t ) = s sin ψ cos( ωt + ∆) − s sin ψ sin( ωt + ∆) cos θ − s cos ψ sin θ − s sin ψ sin( ωt + ∆) sin θ + s cos ψ cos θ + (cid:96) cos β(cid:96) sin β , (111)0with velocity U ( s, t ) found as U ( s, t ) = ∂ r ( s, t ) ∂t = − ωs sin ψ sin( ωt + ∆) − ωs sin ψ cos( ωt + ∆) cos θ − ωs sin ψ cos( ωt + ∆) sin θ . (112) B. Time-averaged rate of working
We next use resistive-force theory [57, 58] to calculate the rate at which work is done by the first rod on the fluid,averaged over one period. We denote by c (cid:107) and c ⊥ the resistance coefficients of the rods in the tangential and normaldirections respectively. We approximate these coefficients as constant, even though the presence of the no-slip wallmeans that the coefficients may vary as the orientations of the rods change during their periodic motion [3]. Neglectingthe wall effect, as in Refs. [50, 55], will allow us to make progress analytically.We denote by F ( t ) the total force exerted by the second rod on the fluid at time t . To leading order in theirseparation, the second rod is unaffected by the first, so that using resistive-force theory, F is given by F = (cid:90) L c ⊥ U d s = − c ⊥ L ω sin ψ sin( ωt + ∆)cos( ωt + ∆) cos θ cos( ωt + ∆) sin θ . (113)Similarly to what we did in Sec. III, we denote by u → ( s, t ) the flow induced by the motion of the second rod,evaluated at the material point r ( s, t ) at arclength s along the first rod, at time t . We introduce the centre of massof the second rod R ( t ) with components ( X ( t ) , Y ( t ) , Z ( t )) given by R ( t ) ≡ X ( t ) Y ( t ) Z ( t ) ≡ r (cid:18) L , t (cid:19) ≡ x (cid:0) L , t (cid:1) y (cid:0) L , t (cid:1) z (cid:0) L , t (cid:1) = L sin ψ cos( ωt + ∆) − L sin ψ sin( ωt + ∆) cos θ − L cos ψ sin θ − L sin ψ sin( ωt + ∆) sin θ + L cos ψ cos θ + (cid:96) cos β(cid:96) sin β . (114)Since the rods are widely separated, we may approximate the flow u → as that due to a point force F exerted onthe fluid at position R . Using the image system near no-slip surfaces [54], the flow u → is given to leading order by u → (cid:39) z Z πµ | r − R | F ,x ( x − X ) + F ,y ( x − X )( y − Y ) F ,x ( x − X )( y − Y ) + F ,y ( y − Y ) , (115)where we recall that z is given by Eq. (109). Using this result and the approximations x − X (cid:39) − (cid:96) cos β and y − Y (cid:39) − (cid:96) sin β (justified since L (cid:28) (cid:96) ), the flow u → is given to leading order by u → (cid:39) − c ⊥ L ω sin ψz Z [cos β sin( ωt + ∆) + sin β cos( ωt + ∆) cos θ ]4 πµ(cid:96) cos β sin β . (116)We next denote by f ( s, t ) the force per unit length exerted on the fluid by the material point on the first rod atarclength s and time t . Then, using resistive-force theory, we have f ( s, t ) (cid:39) { c (cid:107) t ( t ) t ( t ) + c ⊥ [ − t ( t ) t ( t )] } · [ U ( s, t ) − u → ( s, t )] , (117)1where we recall that c (cid:107) and c ⊥ are the resistance coefficients and is the identity tensor. Here t ( t ) is the unit tangentvector to the first rod, given explicitly by t ( t ) = sin ψ cos ωt − sin ψ sin ωt cos θ − cos ψ sin θ − sin ψ sin ωt sin θ + cos ψ cos θ . (118)The final term in Eq. (117) captures the interaction between the two model cilia. Since the velocity of each rod isperpendicular to the rod’s tangent, we have (omitting the s and t dependence for brevity) f (cid:39) c ⊥ ( U − u → ) − ( c (cid:107) − c ⊥ )( t · u → ) t , (119)so that the instantaneous rate of working per unit length at each material point is f · U (cid:39) c ⊥ ( U − u → ) · U . (120)We can then calculate explicitly U · U = ω s sin ψ (121)and u → · U (cid:39) s c ⊥ L ω sin ψ πµ(cid:96) × ( − sin ψ sin ωt sin θ + cos ψ cos θ )[ − sin ψ sin( ωt + ∆) sin θ + cos ψ cos θ ] × [cos β sin( ωt + ∆) + sin β cos( ωt + ∆) cos θ ] (cos β sin ωt + sin β cos ωt cos θ ) . (122)From this, we find the average rate at which work is done by the first rod in one period as ω π (cid:90) π/ω (cid:90) L f · U d s d t (cid:39) c ⊥ L ω sin ψ × (cid:18) c ⊥ L sin θ sin ψ (cid:0) sin β cos θ − cos β (cid:1) πµ(cid:96) − c ⊥ L cos ∆ (cid:0) ψ cos θ + sin ψ sin θ cos ∆ (cid:1) (cid:0) sin β cos θ + cos β (cid:1) πµ(cid:96) (cid:19) . (123)Note that for physically relevant setups, the semi-angle of the cone ψ is nonzero, and since the posterior tilt satisfies θ < π/
2, the quantity sin β cos θ + cos β is nonzero. To find the corresponding result for the second rod, we simplyreplace the angle β with π + β and the phase difference ∆ with − ∆, and see (again) that the rate of working isunchanged. Thus, the total rate of energy dissipation in the fluid is given by twice the final result in Eq. (123).With this, we can now consider the optimisation of the rate of working. By inspection, it is straightforward to seethat in-phase motion of the rods (i.e., the case where ∆ = 0) always minimises the rate of working, since the coefficientsof both cos ∆ and cos ∆ are negative. We may additionally compute the phase difference ∆ that maximises the rateof working. In the case sin θ = 0, it is clear that opposite-phase motion maximises the rate of working. Otherwise, therate of working may be written as A − B cos ∆(2 cot ψ cot θ + cos ∆), where both A and B are positive constants.If there are no solutions to the equation cos ∆ = − cot ψ cot θ, (124)then opposite-phase motion of the rods (cos ∆ = −
1) maximises the rate of working. Luckily, for all physically relevantvalues of the parameters ψ and θ , there are no solutions to Eq. (124). This is due to the fact that, for a fixed value of ψ , the quantity cot ψ cot θ decreases from + ∞ to 1 as θ increases from 0 to π/ − ψ (which we recall is the value of θ corresponding to a cone that is tangent to the no-slip surface). Hence, in all cases, the opposite-phase synchronisedmotion of the rods maximises the rate of working, and in-phase motion minimises it. As discussed in the next section,this compares favourably to past theoretical predictions for the dynamic synchronisation of cilia [34].2 V. DISCUSSION
In this paper, we investigated the energetics of synchronised motion for three different models for interacting flagellaand cilia. In each case, we imposed periodic motion for the two appendages with a constant phase difference ∆. Fortwo sheets waving with small amplitude and with both longitudinal and transverse modes, we found that the optimalphase difference ∆ ∗ , which minimises the rate of working of the sheets, can take all values between − π and π ina manner that depends on the amplitudes of and the phase difference between the waving modes, and the meansheet separation. Our calculations reproduced Taylor’s result that in the case of purely transverse deformation,the rate of working is minimised for in-phase waving. We saw also that if the longitudinal mode has nonzero butsmall amplitude compared to the transverse mode, then the energetic minimum occurs for small, but nonzero, phasedifferences, reminiscent of what is seen in the metachronal waves of cilia and in continuum models of cilia arrays [38].Furthermore, if the sheets deform only longitudinally, then in-phase waving results in maximum dissipation, a resultthat we rationalised physically.In contrast to this, in the cases of spheres moving along circular orbits parallel to a no-slip plane (modelling two-stroke cilia) and whirling rods with a posterior tilt also above a no-slip surface (modelling nodal cilia), the energeticminimum always occurs at the in-phase configuration, ∆ = 0. For spheres in elliptical orbits perpendicular to a no-slipplane, the minimum rate of dissipation of energy occurs for either in-phase or opposite-phase motion, depending onthe relative position and orientation of the orbits. For all these models considered, the dissipation is maximised whenthe phase difference ∆ differs from that giving minimum dissipation by exactly π .The original motivation of our work was Taylor’s energetic argument as explaining the synchronisation of spermato-zoa observed in nature [36]. In his paper, Taylor also remarked that the dynamics of interacting swimmers might notnecessarily result in the minimum-dissipation configuration, depending on how the flagella are actuated. Equippedwith our energetics results in the case where we imposed the kinematics, we can now examine previous studies on themechanisms that result in the synchronisation of model flagella or cilia.In an important paper in the field, elastic compliance was proposed as a generic mechanism allowing hydrodynamicsynchronisation [34]. In the model from that work, cilia are replaced by spheres, which orbit in a plane parallel toa no-slip surface and interact hydrodynamically in the far field. However, unlike our setup in Sec. III A, the orbitsof the spheres are only approximately circular. The intrinsic preferred waving motion of each cilium is modelled asa preferred radius of circular orbit and a preferred angular frequency. An elastic restoring force then acts radiallyto bring the radius of the orbit back to this preferred value, and the angular frequencies of the spheres remainclose to their preferred values since the spheres are widely separated. With this setup, in the case of equal intrinsicfrequencies, the phase difference between the two spheres is governed by the Adler equation. This showed remarkableagreement with experimental observation of synchronisation of the two flagella on the unicellular biflagellate alga Chlamydomonas [22]. The theoretical prediction in that case was therefore that the model cilia synchronise to anin-phase configuration regardless of initial phase difference. In Sec. III A, we found that for imposed constant angularvelocity on precisely circular orbits, dissipation of energy is also minimised by in-phase motion. This is also the casefor the whirling rods modelling nodal cilia found in Sec. IV. There is therefore full agreement in these cases betweenthe dynamics and the energetics points of view.Another study also modelled cilia as spheres above a no-slip surface, but constrained them to move on fixed, tilted,precisely elliptical trajectories [47]. This led to coupled differential equations for the evolution of the phases of thetwo spheres. We recall that, in our paper, the angle β describes the relative position of the two orbits in Sec. III Bfor spheres in elliptical orbits and in Sec. IV for whirling rods. In Ref. [47], far-field calculation showed that forgeneric tilted orbits, the stable synchronised states are in-phase motion for 0 < β < π/ π < β < π/
2, andopposite-phase motion for π/ < β < π and 3 π/ < β < π . These stability results were found to also be consistentwith full numerical solutions to the dynamic equations.We may relate the above far-field result to our result for energetics in Sec. IV on whirling rods. We found that inthe far field, for rotation at constant angular frequency, dissipation is minimised by in-phase motion, while in Ref. [47]the stable synchronised states are opposite-phase or in-phase motion, depending on the relative position of the twonodal cilia. Thus, it is possible for two cilia to synchronise to a state where dissipation is maximum.We also recall our far-field result in Sec. III B, for elliptical orbits perpendicular to the no-slip surface, allowedto have different orientations. We found that the rate of dissipation of energy is minimised either by in-phase oropposite-phase motion of the spheres, depending on the relative position and orientation of orbits. Although thesetup is not quite identical to that in Ref. [47], we note the appearance of in-phase and opposite-phase motion inboth cases. The same also occurs for interacting three-dimensional flagella, with in-phase or opposite-phase motionminimising viscous dissipation depending on the relative position and orientation of the flagella [40].The swimming sheet model in Sec. II is the only one that we considered where the configuration with minimumenergy dissipation can occur for synchronised motion that is neither in phase nor in opposite phase. Past researchon dynamic models of synchronisation for swimming sheets only considered the case of transverse waving [39, 41, 43].3Further work will therefore be needed to probe the relationship between energetics and dynamics in the case of flexiblewaving swimmers with the most general swimming kinematics. ACKNOWLEDGMENTS
This project has received funding from the Faculty of Mathematics, University of Cambridge (Cambridge SummerResearch in Mathematics Programme studentship to W.L.) and from the European Research Council (ERC) underthe European Union’s Horizon 2020 research and innovation programme (grant agreement 682754 to E.L.). [1] D. Bray,
Cell Movements (Garland Publishing, New York, NY, 2000).[2] J. Lighthill, “Flagellar hydrodynamics,” SIAM Rev. , 161–230 (1976).[3] C. Brennen and H. Winet, “Fluid mechanics of propulsion by cilia and flagella,” Annu. Rev. Fluid Mech. , 339–398 (1977).[4] S. Childress, Mechanics of Swimming and Flying (Cambridge University Press, Cambridge U.K., 1981).[5] E. Lauga and T. R. Powers, “The hydrodynamics of swimming microorganisms,” Rep. Prog. Phys. , 096601 (2009).[6] E. Lauga, The Fluid Dynamics of Cell Motility (Cambridge University Press, Cambridge U.K., 2020).[7] E. Lauga, “Bacterial hydrodynamics,” Annu. Rev. Fluid Mech. , 105–130 (2016).[8] L. J. Fauci and R. Dillon, “Biofluidmechanics of reproduction,” Annu. Rev. Fluid Mech. , 371–394 (2006).[9] E. A. Gaffney, H. Gadelha, D. J. Smith, J. R. Blake, and J. C. Kirkman-Brown, “Mammalian sperm motility: Observationand theory,” Annu. Rev. Fluid Mech. , 501–528 (2011).[10] T. J. Pedley and J. O. Kessler, “Hydrodynamic phenomena in suspensions of swimming microorganisms,” Annu. Rev.Fluid Mech. , 313–358 (1992).[11] J. S. Guasto, R. Rusconi, and R. Stocker, “Fluid mechanics of planktonic microorganisms,” Annu. Rev. Fluid Mech. ,373–400 (2012).[12] E. Lauga and R. E. Goldstein, “Dance of the microswimmers,” Phys. Today , 30 (2012).[13] M. Kim, J. C. Bird, A. J. Van Parys, K. S. Breuer, and T. R. Powers, “A macroscopic scale model of bacterial flagellarbundling,” Proc. Natl. Acad. Sci. U.S.A. , 15481–15485 (2003).[14] M. Kim, M. Kim, J. C. Bird, J. Park, T. R. Powers, and K. S. Breuer, “Particle image velocimetry experiments on amacro-scale model for bacterial flagellar bundling,” Exp. Fluids. , 782–788 (2004).[15] M. Reichert and H. Stark, “Synchronization of rotating helices by hydrodynamic interactions,” Eur. Phys. J. E , 493–500(2005).[16] P. J. A. Janssen and M. D. Graham, “Coexistence of tight and loose bundled states in a model of bacterial flagellardynamics,” Phys. Rev. E , 011910 (2011).[17] S. Y. Reigh, R. G. Winkler, and G. Gompper, “Synchronization and bundling of anchored bacterial flagella,” Soft Matter , 4363–4372 (2012).[18] S. Y. Reigh, R. G. Winkler, and G. Gompper, “Synchronization, slippage, and unbundling of driven helical flagella,”PLOS One , e70868 (2013).[19] F. Hayashi, “Sperm co-operation in the Fishfly, Parachauliodes japonicus ,” Funct. Ecol. , 347–350 (1998).[20] Y. Yang, J. Elgeti, and G. Gompper, “Cooperation of sperm in two dimensions: Synchronization, attraction, and aggre-gation through hydrodynamic interactions,” Phys. Rev. E , 061903 (2008).[21] D. M. Woolley, R. F. Crockett, W. D. I. Groom, and S. G. Revell, “A study of synchronisation between the flagella ofbull spermatozoa, with related observations,” J. Exp. Biol. , 2215–2223 (2009).[22] R. E. Goldstein, M. Polin, and I. Tuval, “Noise and synchronization in pairs of beating eukaryotic flagella,” Phys. Rev.Lett. , 168103 (2009).[23] D. R. Brumley, K. Y. Wan, M. Polin, and R. E. Goldstein, “Flagellar synchronization through direct hydrodynamicinteractions,” eLife , e02750 (2014).[24] J. R. Blake and M. A. Sleigh, “Mechanics of ciliary locomotion,” Biol. Rev. Camb. Phil. Soc. , 85–125 (1974).[25] D. R. Brumley, M. Polin, T. J. Pedley, and R. E. Goldstein, “Hydrodynamic synchronization and metachronal waves onthe surface of the colonial alga Volvox carteri ,” Phys. Rev. Lett. , 268102 (2012).[26] Y. Okada, S. Takeda, Y. Tanaka, J.-C. I. Belmonte, and N. Hirokawa, “Mechanism of nodal flow: a conserved symmetrybreaking event in left-right axis determination,” Cell , 633–644 (2005).[27] J. H. E. Cartwright, O. Piro, and I. Tuval, “Fluid-dynamical basis of the embryonic development of left-right asymmetryin vertebrates,” Proc. Natl. Acad. Sci. U.S.A. , 7234–7239 (2004).[28] S. Strogatz, “From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators,”Physica D , 1–20 (2000).[29] R. E. Goldstein, E. Lauga, A. I. Pesci, and M. R. E. Proctor, “Elastohydrodynamic synchronization of adjacent beatingflagella,” Phys. Rev. Fluids , 073201 (2016).[30] R. R. Bennett and R. Golestanian, “Phase-dependent forcing and synchronization in the three-sphere model of Chlamy-domonas ,” New J. Phys. , 075028 (2013). [31] V. F. Geyer, F. J¨ulicher, J. Howard, and B. M. Friedrich, “Cell-body rocking is a dominant mechanism for flagellarsynchronization in a swimming alga,” Proc. Natl. Acad. Sci. U.S.A. , 18058–18063 (2013).[32] G. Quaranta, M.-E. Aubin-Tam, and D. Tam, “Hydrodynamics versus intracellular coupling in the synchronization ofeukaryotic flagella,” Phys. Rev. Lett, , 238101 (2015).[33] K. Y. Wan and R. E. Goldstein, “Coordinated beating of algal flagella is mediated by basal coupling,” Proc. Natl. Acad.Sci. U.S.A. , E2784–E2793 (2016).[34] T. Niedermayer, B. Eckhardt, and P. Lenz, “Synchronization, phase locking, and metachronal wave formation in ciliarychains,” Chaos , 037128 (2008).[35] N. Uchida and R. Golestanian, “Generic conditions for hydrodynamic synchronization,” Phys. Rev. Lett. , 058104(2011).[36] G. I. Taylor, “Analysis of the swimming of microscopic organisms,” Proc. R. Soc. A , 447–461 (1951).[37] S. Gueron and K. Levit-Gurevich, “Energetic considerations of ciliary beating and the advantage of metachronal coordi-nation,” Proc. Natl. Acad. Sci. U.S.A. , 12240–12245 (1999).[38] S. Michelin and E. Lauga, “Efficiency optimization and symmetry-breaking in a model of ciliary locomotion,” Phys. Fluids , 111901 (2010).[39] G. J. Elfring and E. Lauga, “Hydrodynamic phase locking of swimming microorganisms,” Phys. Rev. Lett. , 088101(2009).[40] C. Mettot and E. Lauga, “Energetics of synchronized states in three-dimensional beating flagella,” Phys. Rev. E , 061905(2011).[41] G. J. Elfring and E. Lauga, “Passive hydrodynamic synchronization of two-dimensional swimming cells,” Phys. Fluids ,011902 (2011).[42] G. J. Elfring, O. S. Pak, and E. Lauga, “Two-dimensional flagellar synchronization in viscoelastic fluids,” J. Fluid Mech. , 505–515 (2010).[43] G. J. Elfring and E. Lauga, “Synchronization of flexible sheets,” J. Fluid Mech. , 163–173 (2011).[44] M. C. Lagomarsino, B. Bassetti, and P. Jona, “Rowers coupled hydrodynamically. Modeling possible mechanisms for thecooperation of cilia,” Eur. Phys. J. B , 81–88 (2002).[45] M. C. Lagomarsino, P. Jona, and B. Bassetti, “Metachronal waves for deterministic switching two-state oscillators withhydrodynamic interaction,” Phys. Rev. E , 021908 (2003).[46] P. Lenz and A. Ryskin, “Collective effects in ciliar arrays,” Phys. Biol. , 285–294 (2006).[47] A. Vilfan and F. J¨ulicher, “Hydrodynamic flow patterns and synchronization of beating cilia,” Phys. Rev. Lett. , 058102(2006).[48] J. Kotar, M. Leoni, B. Bassetti, M. C. Lagomarsino, and P. Cicuta, “Hydrodynamic synchronization of colloidal oscilla-tors,” Proc. Natl. Acad. Sci. U.S.A. , 7669–7673 (2010).[49] J. Kotar, L. Debono, N. Bruot, S. Box, D. Phillips, S. Simpson, S. Hanna, and P. Cicuta, “Optimal hydrodynamicsynchronization of colloidal rotors,” Phys. Rev. Lett. , 228103 (2013).[50] D. J. Smith, A. A. Smith, and J. R. Blake, “Mathematical embryology: the fluid mechanics of nodal cilia,” J. Eng. Math. , 255–279 (2011).[51] S. Kim and J. S. Karrila, Microhydrodynamics: Principles and Selected Applications. (Butterworth-Heinemann, Boston,MA, 1991).[52] J. Happel and H. Brenner,
Low Reynolds Number Hydrodynamics (Prentice Hall, Englewood Cliffs, NJ, 1965).[53] L. G. Leal,
Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes , Vol. 7 (CambridgeUniversity Press, 2007).[54] J. R. Blake and A. T. Chwang, “Fundamental singularities of viscous flow. Part 1: Image systems in vicinity of a stationaryno-slip boundary,” J. Eng. Math. , 23–29 (1974).[55] D. J. Smith, J. R. Blake, and E. A. Gaffney, “Fluid mechanics of nodal flow due to embryonic primary cilia,” J. R. Soc.Interface , 567–573 (2008).[56] S. Nonaka, S. Yoshiba, D. Watanabe, S. Ikeuchi, T. Goto, W. F. Marshall, and H. Hamada, “De novo formation ofleft–right asymmetry by posterior tilt of nodal cilia,” PLOS Biol. , e268 (2005).[57] G. K. Batchelor, “Slender-body theory for particles of arbitrary cross-section in Stokes flow,” J. Fluid Mech. , 419–440(1970).[58] R. G. Cox, “The motion of long slender bodies in a viscous fluid. Part 1. General theory,” J. Fluid Mech.44