Collective synchronization states in arrays of driven colloidal oscillators

The phenomenon of metachronal waves in cilia carpets has been well known for decades; these waves are widespread in biology, and have fundamental physiological importance. While it is accepted that in many cases cilia are mainly coupled together by the hydrodynamic velocity field, a clear understanding of which aspects determine the collective wave properties is lacking. It is a difficult problem, because both the behavior of the individual cilia and their coupling together are nonlinear. In this work, we coarse-grain the degrees of freedom of each cilium into a minimal description in terms of a configuration-based phase oscillator. Driving colloidal particles with optical tweezers, we then experimentally investigate the coupling through hydrodynamics in systems of many oscillators, showing that a collective dynamics emerges. This work generalizes to a wider class of systems our recent finding that the non-equilibrium steady state can be understood based on the equilibrium properties of the system, i.e. the positions and orientations of the active oscillators. In this model system, it is possible to design configurations of oscillators with the desired collective dynamics. The other face of this problem is to relate the collective patterns found in biology to the architecture and behavior of individual active elements.


Introduction
Most eukaryotic cells are ciliated; in many cell types the cilia act as chemical or mechanical sensors, while some cell types have motile cilia which undergo periodic motion [1]. This filament motion can enable cell motility, as in the case of unicellular algae such as Chlamydomonas reinhardii [2]; in the case of many ciliated cells forming cilia carpets, they can drive flows [3]. These surface flows are essential in various physiological processes including clearance of mucus in the airways [4], circulation of cerebrospinal fluid in the brain [5] and transport of ovules in the fallopian tube [6]. Understanding of cilia-driven flows requires a combination of two different levels of physical description: the behavior of the single cilium, including how it manages to beat periodically and to transfer a net momentum into the fluid; the behavior of two or more cilia, in particular how they couple together to form a collective system that drives fluid efficiently [7]. The single cilium behavior is itself the result of various physical factors coupled together: the internal structure of eukaryotic cilia (usually nine microtubule doublets, plus two central single tubules, connected together in a well-defined architecture by molecular motors); the drive at the scale of the molecular motors, which slide pairs of microtubules along each other; the microtubule mechanics (bending modulus and compressive constraints); and the hydrodynamics of a thin filament driven through a viscous liquid. All these are aspects that are still being studied, although there are already very promising models aiming to represent the key physics [8,9]. From the description at this level of detail, one hopes to 3 conclude on what is the force during a cilium cycle, what effect a time-dependent perturbation has on the cilium dynamics and what effect is expected from changes in physical parameters such as the length of the cilium or the viscosity (or even viscoelasticity) of the surrounding fluid. Single cilium experiments are being carried out on various biological organisms, in particular the alga Chlamydomonas that is experimentally accessible [10,11] and sperm [12]. There is also growing effort at designing synthetic driven filaments that would be able to act as micro-pumps to drive fluids in confined environments [13][14][15][16].
The outcome of all the activities from the molecular up to the single cilium level is an oscillating filament. At a higher level of description, an open question is whether and how these oscillators will be coupled together. This is conceptually a distinct question from the behavior of the single filament, but it is clear that the two issues are coupled: two filaments will or will not exhibit coupled motion, in the first case having a certain well-defined phase difference, depending on the details of the force versus time in each oscillation. Indeed, a key question in the field is how to link the single cilium behavior to the collective dynamical properties. Collective dynamics has been observed in many systems with two or many cilia [2,17,18]. Understanding this link is very important in both directions: from the bottom up, it means that it will be possible to design active elements that will function cooperatively to induce a certain desired flow profile; from the top down, it will mean that from the observation of collective dynamics, which is in most cases much simpler than single cilium observations (think for example of the challenge posed by visualizing single cilia embedded in the mucus on the surface of an airway tissue), it will be possible to infer the properties of the individual active elements, the cilia. The latter aspect opens up the way for diagnostic tools, which could be developed to discriminate between various cilia pathologies [19].
In order to address with a tractable model the question of what determines the state of synchronization between filaments, it is necessary to perform many approximations. These may at first appear quite brutal, but they retain what we believe are the important physical aspects. It is increasingly accepted that the most important coupling term is due to hydrodynamic flow [3,4,20,21]. Firstly, the hydrodynamic flow induced by a slender semiflexible filament is very complex, but in the far field (i.e. at inter-cilia distances greater than the cilium length) it seems reasonable to consider instead the flow induced by a moving sphere; see figure 1(a). This greatly simplifies the calculation of drag forces, both those acting on the individual object and the force induced by one object on another. Indeed, in this far field approximation for spheres in low Reynolds number, there are classical analytical results due to Oseen [22], and the presence of a solid boundary can also be accounted for [23]. Secondly, the behavior of the individual cilium should be simplified. Given that the cilium has been replaced by a solid sphere, all the internal degrees of freedom have been removed in terms of the filament's configuration. There have been two classes of models proposed to account for this. The first is to consider the coupling of two or more objects driven by a force field over closed two-dimensional differentiable orbits [24]. Within this model for cilia, the phase of each 'rotor' is free, and there can be a synchronization of different rotors under certain conditions, studied very generally in [25]. A second class of models consists of driving an oscillation based on the configuration of the system. This 'geometric switch' was first proposed in [26]; here the force is discontinuous, and this is not described within the formalism of [25]. It is fair to say that it is still an open question which of these models represents better the internal degrees of freedom of a cilium, and work is under way to elucidate commonalities and differences in the collective dynamics. Both models describe phase oscillators coupled via hydrodynamic flow. The sketch shows parameters of the model, further described in the text. (c) Optical microscopy image obtained in the optical tweezers setup (a movie of this experiment is available as movie 1 in the online supplementary data, available at stacks.iop.org/NJP/14/105023/mmedia); the center of each bead is determined through image analysis on each frame, and by using that datum a decision on where to place the trapping potential is taken independently for each bead, on each frame. The average distance between nearest neighbors is d = 8 µm. Arrows indicate the orientation of driven oscillations, for each bead. This configuration of nine beads is the largest of the ones investigated experimentally in this work, and is presented in section 5.3.2.

5
We have up to now focused on studying the 'geometric switch' model for cilia synchronization mainly for two reasons: it can be readily realized experimentally, by making use of optical tweezers with feedback, and it is amenable to theoretical analysis, including to first order the effect of thermal noise. We have shown remarkably rich behavior, and have discovered that there are (at least) three important physical considerations that affect the emergence of a collective dynamical state: (a) the strength of coupling, relative to stochastic thermal forces; (b) the mean profile of the driving force in a cycle; and (c) the spatial disposition and orientation of the active elements.
The coupling strength between moving spheres is given by the Oseen tensor, which describes a long-range interaction decaying inversely with the separation distance. Therefore, if two spheres are sufficiently far, the forces they induce on each other are drowned in the Brownian thermal fluctuations of the solvent. We previously investigated the geometric switch model for the simple case of two active elements, showing the extent of synchronization in the presence of noise, and the effect of detuning the periods of the two oscillators [27]. As well as imposing a threshold on the coupling strength for synchronization, thermal noise was shown to induce distinct cross-particle phase correlations [27]. In a simpler system of one particle driven to oscillate by geometric switch, but also subject to an externally clocked periodic force, we showed analytically the effect of noise on the stability of synchronization [28].
Within the geometric switch model, the shape of the driving potential affects the state of synchronization. Wollin and Stark [29] showed this in the geometry of a linear chain of geometric-switch-driven oscillators studied numerically in the absence of noise. For two driven particles, the curvature of the potential determines if nearest neighbors synchronize in phase or in antiphase [29][30][31]. We considered other regular arrangements of oscillators in [30], showing that the collective state in an active system can be predicted from the geometric arrangement and that the curvature of the potential selects for the collective dynamical state with the symmetry of either the slowest or the fastest relaxation modes of the passive system with the same geometry. Recently, we studied the effect of noise on the transition that occurs with changing curvature [31].
The geometric disposition of the oscillators is also essential for establishing a certain collective dynamical pattern. This will be explained in detail in the next section, since the present paper extends our previous work [30,32] to the more general scenario in which oscillators are placed at generic positions and orientations in the plane.
The main aim of this work is to prove that even in apparently complex arrangements of oscillators (disordered and lacking rotational or translational symmetries) it is possible to predict the collective dynamics from just a knowledge of the positions and orientations of the oscillators. In a biological tissue with a carpet of cilia, plenty of oscillators at various positions will be oscillating, quite possibly exerting forces in different principal directions, and their phases will be free to couple. The result obtained here with harmonic potential driving is therefore a key step, albeit in a simple model system, for establishing the important link between microscopic and macroscopic behavior in systems of hydrodynamically coupled oscillators.

Geometric switch model of active oscillators
To understand from the theoretical and experimental points of view the hydrodynamic coupling of phase oscillators, we consider here spheres of radius a driven by harmonic potentials 6 U (x) = κ x 2 /2; for an isolated or single oscillator, the sphere falls towards the minimum of the potential (subject also to any thermal noise). When it reaches a boundary position which is set at a distance ξ from the minimum, the potential is replaced by another harmonic well, with the center a distance λ from the first well; see figure 1. The sphere thus inverts its mean direction and approaches the new minimum. Again, when it is within ξ of the new minimum, the potential is replaced by the first and so on; see figure 1(b). This gives rise to a periodic motion that is bounded in amplitude A = λ − 2ξ but free in phase and period.

Driven dynamics of a single oscillator
For a single oscillator, the equation of motion is a balance between the configuration-dependent driving force F(x, config), the drag and the thermal forces: where γ = 6π ηa, η is the viscosity of the solvent and f(t) is white noise with mean square amplitude set by the fluctuation dissipation theorem (see below). If the potentials are harmonic, and the motion is along the x-axis, then F x (x, config) = ±κ(x ± λ/2), with the signs depending on the configuration history; see figure 1(b). Up to a geometric switch boundary, equation (1) describes on average an exponentially decaying trajectory with relaxation time γ /κ.

Dynamics of many oscillators
Here, we consider N oscillators at positions r n = x nêx + y nêy . The Oseen tensor describes the relation between the set of all velocities {v i } and forces {F j } [22]: with whereê n,m = x n,m /r n,m . The steady-state dynamics of the geometric switch model can be obtained, in principle, by solving the dynamic system of equations where the second term describing the Oseen coupling forces scales with distance d as a/d, and the force F i acting on the ith particle is harmonic. Then [33][34][35], the motion of the jth particle originates a coupling force f c on the ith particle f c i, j = (H −1 ) , that depends on the whole configuration, where H(r) is the Oseen tensor [22]. The equation is valid in between any two switches, with F(r i , config) i driving the ith particle towards some fixed minimum; F(r i , config) i then changes form once the ith particle reaches a geometric switch position. The stochastic force f i (t) in equation (4) represents the thermal noise on the ith particle, and it can be assumed that . This stochastic force is present experimentally, and in previous work we showed that it can be accounted for in numerical Brownian dynamics simulations [27,32]. It complicates considerably the behavior of the nonlinear system, and we will neglect it to make analytical progress below.
The fact that equation (4) changes structure at every switch makes its solution difficult even in the absence of thermal noise. In the presence of noise, it then becomes an intractable manyparticle first passage problem. In a recent paper [30], we have shown how to build analytical solutions in the absence of noise, for some special configurations with high symmetry, starting from a knowledge of the hydrodynamic modes. Here below, we generalize the class of systems on which hydrodynamic modes of the actively driven system can be calculated and related to the steady-state solution.

An effective Oseen tensor accounts for the constraints
The oscillators are driven at fixed angles θ n , i.e. directionsê n = (cos(θ n ), sin(θ n )); see figure 2. In other words, the driving forces and to a good approximation also the displacements in equation (4) are oriented at a fixed angle θ i . We can impose this as a condition that the active driving of each bead along its fixed direction is also constraining each bead to move only along that fixed direction. As in [30,32], since for each bead most of its displacement takes place along the direction of its drive, this is a fair approximation of the main component of the motion. Based on this, we now show how the general coupling equations of equation (4) can be recast in a more transparent way for the system with driven motion. This will enable us to show that the main properties of the steady-state solutions can be predicted quite simply from the equilibrium coupling tensor, and hence depend on N and on the spatial distribution and arrangement of angular drives.
The direction of the forces is known: The non-diagonal terms of Oseen are perturbations to the motion, so as described above we approximate the velocities as v n ≈ v n ·ê n . 8 Substituting (3) and (5) into (6) n,m + x 2 n,m cos θ m cos θ n + r 2 n,m + y 2 n,m sin θ m sin θ n +x n,m y n,m (cos θ m sin θ n + sin θ m cos θ n ) F m .
Equation (7) is a linear relation between the scalar velocities v n and the scalar driving forces F m . By analogy with the general Oseen tensor coupling of equation (2), we can define an effective Oseen tensor that satisfies Confining the motion to one dimension for each bead, we have reduced the system from 2N (set of equation (2)) to N (set of equation (8)) coupled equations.
When working with systems where the amplitude of driving is much smaller than the distance between particles, i.e. A d, it is a good approximation to consider the average distance r n,m in the Oseen tensor, rather than the time-dependent quantity. This timeindependent tensor is then what we consider in the following.
From here onwards, to maintain light notation, we will denote as r n the displacement of bead n from its time average position.

The steady-state dynamics is related to the hydrodynamic modes
The effective Oseen tensor describes the coupling of particles constrained on fixed directions of motion, and can be analyzed in terms of its normal modes. These normal modes would represent patterns of correlated motion that are independent of each other, are driven by harmonic potentials and therefore have each a well-defined relaxation timescale. For example, thermal fluctuations around equilibrium positions would be described by normal modes, as observed in [33][34][35]. However, obtaining the normal modes is only the first step in building a solution in the driven system with the geometric switch rule. In a recent paper, we showed how to proceed from the modes to solutions of the system, for the case of systems with high symmetry [30]. We summarize here the main points, although we do not attempt to build solutions for the less symmetric configurations considered here-it is not useful.
To build a solution of the dynamical system (deterministic case, without noise), it is necessary to propagate the initial conditions up to the first geometric switch event. This propagation is simple, since it can be decomposed into the normal modes. So one needs to work out the set of N mode amplitudes, and these decay exponentially up to the switch time. At that 9 time, the dynamical system changes (since the position of at least one driving harmonic potential changes). The solution needs to maintain continuity in the displacements of each bead, but a new set of normal mode amplitudes now describes the solution up to the next switch. Periodic solutions are those that, after all N beads had switched, returned to the initial displacement configuration. The point we made in [30] and here is that while there is no simple recipe for predicting the stable synchronized state (out of a great variety of solutions that can be found), its main properties are determined by the normal mode with the longest relaxation time.
We now show how the effective Oseen tensor allows us to describe the main properties of the synchronized state of any system of several oscillators (extending the case of driven colloids on circles [30,32]).

There is a dominant mode in the steady-state dynamics
H eff can be diagonalized (since H eff is real and symmetric) to obtain the different normal modes of the coupled system. This follows a similar logic to our previous work [27,30,32], and is summarized in the appendix. Projecting the velocities onto the eigenvectors gives the set of uncoupled equations: The relaxation time of the mode k is τ k = 1/(κλ k ). The N -bead system has N normal modes and, in general, N time constants (although certain modes may be degenerate and have the same τ ). The collective dynamics of the driven system is dominated by the mode (or one of the degenerate modes) with the longest relaxation time: the reason why this holds is currently understood only qualitatively, as discussed in [30,32], and the purpose of this work is to show experimentally that this observation is very general: several experimental configurations are considered in section 5.

Experimental methods
Optical traps are used to confine colloidal beads within harmonic potentials, the system hardware is described in greater detail in [27,36], and a consideration of the errors that can affect this type of experiment is presented in [32]. In this work, a varying number of silica beads of radius a = 1.74 µm (Bangslabs, catalog no. SS05N) are trapped from below by a time-shared laser beam, focused by a water immersion objective (Zeiss, Achroplan IR 63×/0.90 W). A pair of acousto-optical deflectors (AOD) allows the positioning of the laser beam focus anywhere in one plane, with sub-nanometer precision; time sharing is at a rate of 2 × 10 4 Hz, corresponding to negligible diffusion of the beads in each laser cycle. The solvent is 31% glycerol (Fisher, Analysis Grade) and 69% water (Ultrapure grade, ELGA) by weight, corresponding to a nominal viscosity of η = 3 mPa s at 20 • C [37]. Experiments are performed in a temperaturecontrolled laboratory, T = 21 • C. The trapping plane is positioned several particle diameters above the flat bottom of the sample, in a sample volume that is about 150 µm thick.
To realize the geometric switch condition, active driving of each colloid is implemented here, similarly to [27,32], but for 3, 5 and 9 particles, driving the colloids on segments at predetermined orientations.
The geometric switch is implemented experimentally by monitoring the system configuration via image analysis, and positioning the laser beams as appropriate: images acquired at 100 Hz are analyzed by correlation filtering in the computer, giving particle positions on each camera frame. When the geometric switch is observed to have been reached by a bead, feedback is sent to the AOD to actuate laser deflection. Since the configuration is analyzed experimentally only at each frame, the corresponding time interval should be considered as a feedback time. Colloidal particles are always being driven, never reaching the minimum of the active trap.
The optical trap potential is harmonic to a very good approximation, with stiffness κ of 0.2 pN µm −1 when trapping two beads (stiffness was maintained fairly constant when changing the number of beads trapped and is calibrated from the distribution of displacements in static traps, typically with 10% precision). The relaxation time τ 0 = γ /κ is of the order of 0.44 s. The experiments have been performed with λ = 1 µm, ξ = 0.31 µm and d = 8 µm. The deterministic period of an isolated oscillator is T 0 = 2 0 log[(λ − ξ )/ξ ] [27], which under the experimental conditions is about 0.7 s. With image acquisition through an AVT Marlin F-131B CMOS camera, operating at a shutter aperture time of 1.5 ms and a frame rate of 100 fps (depending on the system size, hence the captured region of interest) there are multiple frames captured within the relevant timescales τ 0 , T 0 . Video is acquired for over 5 min, i.e. over 30 000 frames. There is typically a transient lasting around a few periods before the systems reach the steady state discussed below.

Two coupled oscillators
The simplest case we consider here is that of two oscillators with variable angles θ 1 and θ 2 ; see figure 2. Within this class, three particular cases can be solved without approximations: horizontal-parallel, vertical-parallel and horizontal-vertical, giving coupling terms 3a/2d, 3a/4d and 0, respectively.
In general, for this N = 2 case the effective Oseen tensor H eff is with C the hydrodynamic coupling term between the two beads, which depends on the directions of oscillation θ 1 and θ 2 : (2 cos θ 1 cos θ 2 + sin θ 1 sin θ 2 ) .
The diagonalization of H eff gives the eigenvalues 1 = (1 − C)/ and 2 = (1 + C)/ and the eigenvectors: The first mode (smaller λ, hence higher relaxation time) represents a motion of the beads in opposite directions, i.e. antiphase oscillations, while the second mode is oscillations in phase. To visualize a mode n, we represent its decomposition on (v i ) i as follows: for each oscillator i a In the simplest case of two beads, oscillating on axis (i.e. (θ 1 , θ 2 ) = (0, 0)) the two normal modes of oscillationsṽ 1 ,ṽ 2 are simply antiphase motion (a) and in-phase motion (b). The geometry of the modes is illustrated by the vector pairsṽ 1 ,ṽ 2 ; the eigenvectors are labeled from 1, . . . , N (in this case N = 2) in the order of decreasing relaxation timescale, and this convention is followed throughout this paper. The percentage indicates the fractional amplitude observed experimentally from the decomposition of the steady-state solution onto the normal modes. vector is drawn, centered on the average position r i , with a directionê n (angle θ n ) and a length (positive or negative) proportional to the coefficient of the component v i in the eigenvector describing the mode. As an example, the two modes for the (θ 1 , θ 2 ) = (0, 0) configuration are shown in figure 3. The trajectories of two beads in this horizontal parallel configuration are shown in figure 4(a).
In configurations (θ 1 , θ 2 ) for which the coupling C is non-zero, the mode with the longest relaxation time isṽ 1 . Therefore, we expect the system to synchronize in antiphase. This is evident in the trajectories of figure 4(a) and is confirmed by the projection of these experimental trajectories onto the eigenvectors, as shown in figure 4(b).
To estimate more quantitatively from experimental data the fraction of oscillations in each mode, we proceed in two steps. (a) We calculate the projections of the oscillations onto the modes. (b) We define the relative amplitude of the modes n from where std(r n ) denotes the standard deviation of ther n (t) mode. For example, in the case of (θ 1 , θ 2 ) = (0, 0) of figure 4, we found a high fraction f 1 = 87% in the 'antiphase' mode. This is as expected and consistent with [27], where the small existing fraction f 2 in the 'in phase' mode was shown to be due to the thermal fluctuations. In the particular configuration (θ 1 , θ 2 ) = (0, π/2), the fraction of the antiphase mode and in-phase mode are almost equal (52 and 48%). This is because in this case the coupling C between the beads is zero: the two modes are degenerate as they have the same relaxation time and there is no preference for the system to oscillate in either mode. The two particular cases above show two possible regimes: synchronization in AP and no synchronization. We have studied intermediate configurations by varying both θ 1 and θ 2 . Eight values of each of the angles were fixed in experiments, mapping out 64 points in the (θ 1 , θ 2 ) space. The magnitude of the theoretical value of coupling, given by |2 cos(θ 1 ) cos(θ 2 ) + sin(θ 1 ) sin(θ 2 )|, from equation (12), is shown in figures 5(a) and (b), while the experimentally measured fraction f 1 describing the decomposition into the antiphase mode is shown in figures 5(c) and (d). There is a remarkable similarity between these two surfaces. In this simple N = 2 case, this fraction is equivalent to the degree of correlation, falling to 1/2 when there is no correlation. In particular, the locus of points with zero coupling is known analytically from equation (12), and describes two arcs in (θ 1 , θ 2 ) space, in agreement with the experimental points where f 1 0.5. Clearly, it is not just the horizontal/vertical perpendicular pair of oscillators that is un-coupled.
To make a quantitative link between |C| and the fraction in antiphase requires not just the construction of deterministic solutions from the normal modes, but also a consideration of the effect of noise, which is beyond the scope of this work.

Three oscillators on a circle
To test the theoretical framework of the effective Oseen tensor for the driven system, and its predictive power for the symmetry and dominant properties of the steady-state dynamics, a number of three-bead configurations is explored.
In figure 6(a), radial oscillations are considered. The normal modes are shown, together with the experimental fractions f i , i = 1, 2, 3. The same in figure 6(b) for tangential oscillations. In both cases the normal mode with the highest relaxation time appears dominant in the solution.
In figure 7, two configurations that lack symmetry are considered. In figure 7(a), three beads on the vertices of an equilateral triangle are driven in a parallel direction. If one thinks of this arrangement by extension from the two bead system, the bead at the top vertex could 'choose' to synchronize with either of the bottom ones, with no preference; the longest-lived mode, shown in the figure, has zero amplitude for the top bead (hence it is clear that it cannot by itself represent the solution to the driven dynamics, which requires that once per period each bead should have a displacement of amplitude A). Aside from these considerations, which can help to build the solution from the normal modes, here too the dominant mode isṽ 1 . This is found also for the system in figure 7(b), where the bottom beads are moved vertically, and the top bead horizontally.

Synchronization through a 'master bead'
In our recent work [32], we introduced the concept of 'dynamical motifs' as those patterns of dynamics that one might expect to identify locally in a cluster of strongly correlated oscillators, coupled more weakly within a larger system. The most basic of all the motifs is the behavior of two beads, which as we have seen above can be tuned to be in antiphase all the way to noncoupled, by a choice of orientation. Developing this idea, it is particularly interesting to consider the behavior of two beads in a condition of no coupling (e.g. perpendicular) when a third bead is introduced into the system. The third bead is a free-phase oscillator, equivalent to the first two. Its role is only special due to its particular position and orientation. Through this third bead the first two become coupled, and their synchronization state is tuned by the position and orientation of the third bead, which acts as a control or 'master'. Understanding this simple structure might allow some design rules to be found, to make a system of oscillators at prescribed positions and orientations that will self-organize into a dynamical state with some desired properties such as inducing a fluid flow over a lengthscale larger than the distance between neighboring oscillators. Three beads oscillating (a) radially and (b) tangentially to a circle. Panels illustrate the configuration, and the three normal modes for each case. In both systems, the normal mode with the longest relaxation time,ṽ 1 , is found to be dominant in the experimentally observed steady-state solution. The percentages indicate the fraction of amplitude from the decomposition of the steady-state solution onto modes, as found in experiment.

Three beads.
In figure 8(a), the bottom beads 1 and 2 are oscillating along orthogonal axes and are not coupled directly. However, each of the first two beads are coupled to the top bead 3 and we can expect synchronization of the bottom beads through the top one.
The equations of motion given by the effective Oseen tensor are Figure 7. Even in configurations that lack the rotational symmetry considered in figure 6, the dominant mode in the steady-state solution is that with the longest relaxation time. Three beads are studied driven to oscillate in (a) a parallel configuration, and (b) in 'bridge' configuration on the vertices of an equilateral triangle. As before, the normal modes are shown, in the order of decreasing relaxation timescale, with the percentages indicating the decomposition of the steady-state solution onto modes, as found in experiment.  As bead 1 is coupled to bead 3 with a coefficient s and bead 2 is coupled to bead 3 with a coefficient u, the coupling term between 1 and 2 is expected to be of the form su: This second-order coupling term between beads 1 and 2 is not null in most cases. Between 0 and π, |C 1,2 | is maximum for θ 3 = 1.07 and 2.64 rad and it is null for θ 3 = 0.24 (u = 0) and 1.90 rad (s = 0). An order parameter Q 1,2 can be defined to describe the state of synchronization between beads 1 and 2, by taking the product of the instantaneous trap states σ 1 (t) = ±1 and σ 2 (t) = ±1. This is done as an average over the experiment duration t tot : In previous work [27,31], the directions were parallel, and we had chosen the convention that Q = 1 represented antiphase motion and Q = −1 in-phase motion. Here the directions are general, and we have picked an arbitrary choice of the sign. The theoretical coupling strength C 1,2 reproduces the same θ 3 dependence on the state of correlation as seen experimentally in the correlations between beads 1 and 2, shown by Q 1,2 in figure 8(b). In particular, the expected maxima and points of zero coupling occur exactly at the master angle θ 3 expected theoretically. The master bead can induce a (weak) P or AP synchronization, as well as no synchronization, depending on the choice of θ 3 . This can be 'exploited' to construct more complex collective dynamics, as explored below.

Nine beads.
Here a system is made, as shown in figures 1(c) and 9, in which the nearestneighbor bottom beads are not directly coupled with each other. This is a more complex arrangement, but the nine modes are readily calculated from the effective Oseen tensor; the corresponding eigenvalues are γ λ = {0.6320, 0.6541, 0.7969, 0.8347, 0.8566, 0.9314, 1.2393, 1.2820, 1.7731}. As expected, the modes with the lowest eigenvalue dominate the experimentally observed steady-state solution; see figure 9.

Large collective wavelengths in large arrays
The knowledge built with experimental systems gives us confidence that as a general principle the longest-lived normal mode will dominate the steady-state dynamics for hydrodynamically coupled systems driven by the 'geometric switch' rule. The effective Oseen tensor can be readily diagonalized for systems of thousands of oscillators. This should enable configurations to be tested to find patterns that will synchronize into large-scale collective dynamics. As a proof of principle, a 2 × 100 array is considered, testing the effect of choosing angles and distances to show a collective dominant mode with large lengthscale.
As for the case of nine beads, the lower set of beads is configured with the nearest neighbors vertical/horizontal, so it is not directly coupled. If beads are put on the top row with a fixed direction, as was the case for nine beads, this leads to collective synchronization but the emergent states have very short lengthscale: the nearest neighbors of the bottom row can oscillate in-phase or in antiphase. One way to design a system in which the dominant mode has large wavelength (which presumably could be a feature of biological metachronal waves) is to embed a spatial frequency in the structure of the oscillators. In the configuration sketched in figure 10(a), where only a nine bead out of 100 section on the bottom row is shown, the top row beads are set to oscillate at angles θ i that rotate by π radians over p beads. The period p was tested in a wide range; data in figures 10(b) and (c) are for p = 20. What emerges clearly is that the longest-lived modes (the first four are shown in the figure) have a large lengthscale, which was 'coded' in the geometric arrangement (the orientational order in this case). Comparing panels (b) and (c), the observation can be made that subtle adjustments in the positional arrangement influence the details of the dominant eigenmodes, but their large-scale structure remains unaffected. Such adjustments are key to designing arrays to optimally generate particular flow patterns-but as remarked often in this work this requires obtaining the deterministic solutions or running a Brownian dynamics simulation or performing the experiment: the first approach is difficult, and the last two become impractical above a few hundreds of particles (simulations) or a few tens (driving with optical traps). The set of systems for which we can get steady-state dynamics is therefore limited-but the normal modes remain easily computable for much larger arrays. This is what makes simple models powerful in terms of establishing a 'bottom up' prediction of emerging properties. The complementary  Only a part of the system, which is made up of 2 × 100 beads, is drawn here. The angles in the top row are chosen so that an angle of π is rotated over 20 beads. This spatial structure is observed in the four longest-lived modes, represented in (b) by showing the amplitude of the mode on the odd-index bottom row beads (blue, green, red, cyan solid lines, from modes 1 to 4, respectively). In the system studied in (b), the top row beads are staggered exactly at the mid point between beads on the bottom row, whereas in (c) the top row beads are shifted by 0.4d from the mid-point. In all cases d = 8 µm, and radius a = 1.74 µm, as in the rest of this work.
challenge-possibly more difficult-remains to determine the microscopic effective properties from the observation of collective dynamics in a large system, for example an array of biological cilia.

Conclusion
It has been shown that through the analysis of the mean configurations (position and angles) of driven oscillators we can predict the main properties of the steady state of the collective system undergoing geometric switch oscillations. The method is very general and can describe any planar configuration of oscillators. The mode that contributes the most (and hence sets the symmetry) of the steady-state solution is the one with the longest relaxation time. It is remarkable that in this (nonlinear, strongly coupled and actively driven) system it is therefore possible to make predictions on the out of equilibrium behavior, based purely on static equilibrium properties.
Periodic solutions of the dynamical system, deterministic i.e. without noise, are a basic ingredient for understanding and predicting the behavior of the physical system. In general, one such solution, we may call this fundamental, is the attractive trajectory for every initial condition. Other periodic trajectories are unstable. Different orientations of the line of oscillations of the beads and different average distances of the beads affect the hydrodynamic coupling and change the shape of the periodic solutions, still leading the system to synchronization. The effect of the thermal noise is twofold: (i) it distorts the analytic fundamental solution, leading to a measured fraction f i larger than expected from the theoretical fraction evaluated from the fundamental solution. (ii) It produces transitions from the fundamental solution to another, unstable, periodic solution. In these unstable solutions, often the normal modes associated with the longest relaxation time have smaller contribution than normal modes with shorter relaxation times. We cannot at the moment separate these two, and we have simply shown that in a variety of systems the measured fraction f 1 , associated with the normal mode with longest relaxation time, is always larger then other fractions.
We have shown that in certain geometrical arrangements, the position and orientation of one of the oscillators (having 'free' phase, the same as the others) is particularly important for the collective symmetry that arises; this oscillator has been called a 'master' bead. Of possible interest for future work are the ideas of having one or more oscillators with fixed phase and the related question of phase front propagation upon an external localized perturbation.
While much work remains to be done to explain in detail the properties of biological cilia, and biologically observed metachronal waves, this work shows that in order to understand a phenomenon that emerges at a certain level of complexity (here, the collective dynamics of many cilia) it can be valuable to coarse grain the degrees of freedom of the lower-level elements (here the cilia, which are themselves complex systems and are reduced to very simple phase oscillators). The physics of hydrodynamically coupled actively driven systems is very rich and many aspects remain to be investigated. Scalar coordinates of the bead, scalar velocities, scalar forces (each set being collected in an N -vector) may be decomposed on this basis: |r (t) = (x 1 (t), . . . , x n (t)), |v(t) = d dt (x 1 (t), . . . , x n (t)) = (v 1 (t), . . . , v n (t)), |F(t) = (F 1 (t), . . . , F n (t)) , By the use of this basis, the system of coupled differential equations (2) is translated into a system of un-coupled differential equations for the normal modes: d dtr j (t) = λ jFj (t), j = 1, . . . , N , which is equation (10). The experimental time evolution {r j (t)} is analyzed by computing the time averaged fractions f i , defined in (14).