Synchronization in cilia carpets: multiple metachronal waves are stable, but one wave dominates

Carpets of actively bending cilia represent arrays of biological oscillators that can exhibit self-organized metachronal synchronization in the form of traveling waves of cilia phase. This metachronal coordination supposedly enhances fluid transport by cilia carpets. Using a multi-scale model calibrated by an experimental cilia beat pattern, we predict multi-stability of wave modes. Yet, a single mode, corresponding to a dexioplectic wave, has predominant basin-of-attraction. Similar to a"dynamic"Mermin-Wagner theorem, relaxation times diverge with system size, which rules out global order in infinite systems. In finite systems, we characterize the synchronization transition as function of quenched frequency disorder, using generalized Kuramoto order parameters. Our framework termed Lagrangian Mechanics of Active Systems allows to predict the direction and stability of metachronal synchronization for given beat patterns.

Motile cilia are slender cell appendages that bend rhythmically due to the activity of molecular dynein motors inside [1]. Collections of motile cilia can spontaneously synchronize their bending waves, e.g., in carpets of many cilia on airway epithelium [2], as well as on the surface of model organisms, e.g., green alga colonies or unicellular Paramecium [3,4]. Metachronal coordination manifests itself as a self-organized traveling wave of cilia phase (similar to a Mexican wave in a soccer stadium). Numerical models showed that this synchronization is important for efficient fluid transport [5,6]. Tissue-scale polarity systems align cilia bases [7], ensuring a common direction of the effective stroke of the cilia beat. In many species, cilia beat patterns are chiral, e.g., with counter-clockwise motion of cilia during their recovery stroke close to the surface [3]. The directions of metachronal waves enclose defined angles relative to the direction of the effective stroke [3,8], presumably set by the chirality of the cilia beat [9].
Already in 1952, Taylor proposed that hydrodynamic interactions between nearby cilia play a key role for their synchronization [10]. When a beating cilium performs its bending wave, it sets the surrounding fluid in motion, resulting in time-dependent hydrodynamic friction forces that act on nearby cilia. Recent experiments indeed demonstrated synchronization by hydrodynamic coupling in pairs of cilia [11], as well as phase-locking to external oscillatory flows with characteristic Arnold tongues [12,13]. Recent theoretical work predicts different synchronization modes between pairs of hydrodynamically coupled cilia, depending on their relative positions [14,15].
Yet, we still do not understand how hydrodynamic interactions and the shape of the cilia beat select the direction of metachronal waves in cilia carpets. Multiple wave directions are possible, yet these may not be stable to small perturbations (local stability) or be unlikely to be selected for random initial conditions (global stability). A key question thus concerns the local and global stability of different metachronal wave modes. The global stability of synchronization states in collections of coupled oscillators, not just interacting cilia, is still a field of active research [16][17][18].
The periodic sequence of shapes that a cilium assumes during its beat cycle represents a limit cycle in a high-dimensional shape space [19,20]. This limit cycle can be parameterized by a single phase variable such that phase speed is constant in the absence of perturbations and noise [21]. This allows to describe beating cilia as phase oscillators [19,22]. In the presence of external flows, which change the hydrodynamic load, the phase speed changes, i.e., cilia progress slower or faster along their beat cycle, while deviations from the limitcycle sequence of shapes remain small for moderate flows [13,23,24]. This load-response of cilia (reflected by the load-dependent speed of their phase variable) is a prerequisite for cilia synchronization by hydrodynamic interactions, and is implicit in previous minimal models [9,13,[25][26][27][28][29][30][31].
Here, we harness multi-scale simulations to combine the benefits of detailed hydrodynamic simulations based on experimentally measured cilia beat patterns, and those of minimal models amenable to local and global stability analysis. Our approach, termed Lagrangian mechanics of active systems [14], enables us to study global stability in arrays of hydrodynamically coupled cilia.
Beating cilia as coupled phase oscillators. We consider a carpet of N cilia positioned on a regular triangular lattice of base points x j in a rectangular domain with periodic boundary conditions, see Fig. 1(d). Each cilium is described as a phase oscillator whose phase ϕ j advances by 2π on each cycle, like a clock. This phase variable ϕ j parameterizes a periodic sequence of three-dimensional cilia shapes, previously measured for Paramecium [3,34], see Fig. 1(a). When the phase ϕ j increases, i.e., the cilium progresses along its beat cycle, the corresponding shape change of the cilium sets the surrounding fluid in motion, resulting in time-dependent hydrodynamic friction forces that act on the other cilia. For nearby cilia, the resultant hydrodynamic interactions can be computed from the Stokes equation valid at zero Reynolds number [14,35], see also Supplemental Material (SM). The plane containing the cilia base points is modeled as a non-slip boundary, thus hydrodynamic interactions decay as 1/d 3 as function of distance d [14,36].
We consider the dynamics of N cilia in a rectangular unit cell with periodic boundary conditions, which is characterized by a N -component vector Φ = (ϕ 1 , . . . , ϕ N ) ∈ R N of cilia phases. Because the Stokes equation is linear [37], the surface density of hydrodynamic friction forces f (x) at time t (defined on the combined surface S of all cilia and the boundary surface) is linear in the generalized velocityΦ. Thus, the power exerted by the moving cilia on the surrounding fluid R = S d 2 x f (x) ·ẋ becomes a quadratic form inΦ [14] R =Φ · Γ(Φ) ·Φ (1) with a symmetric N × N matrix of generalized hydrodynamic friction coefficients Γ = Γ(Φ). Here, Γ ii represents self-friction of cilium i, while Γ ij characterizes hydrodynamic interactions between cilia i and j. Below, we compute Γ(Φ) in a pairwise-interaction approximation.
For each cilium, we introduce the generalized hydrodynamic friction force P i as the friction force conjugate to the generalized coordinate ϕ i (following the formalism of Lagrangian mechanics of dissipative systems with R/2 as Rayleigh dissipation function [14,38]) Assuming low Reynolds numbers, there is at all times a force balance between the generalized friction force P i and an active driving force Q i that coarse-grains the active processes inside cilium i that drive the cilia beat The cilia driving force Q i is an intrinsic property of cilium i, hence only depends on ϕ i , and possibly load P i . We make the simplifying assumption that Q i is independent of load. Previous experiments in the green alga Chlamydomonas [24] as well as cilia bundles in external flow [13] showed that this assumption together with Eq. (3) quantitatively accounts for the load response of cilia [31,39], i.e., the experimental observation that cilia progress slower/faster along their beat cycle upon increase/decrease of hydrodynamic load. Next, we compute the generalized hydrodynamic friction forces P i = j Γ ijφj with friction coefficients Γ ij for a real cilia beat pattern, and calibrate the active driving forces Q i . Note that previous minimal models of hydrodynamically interacting spheres [25][26][27][28][29][30] can likewise be written in the form of Eq. (3), yet with simplified driving and friction forces.
Oscillator coupling calibrated from hydrodynamic simulations. Initial simulations showed that the friction coefficient Γ ij (Φ) is largely independent of the phases of the other cilia, ϕ k , k = i, j. This allows us to use an approximation of only pairwise-interactions for Γ(Φ) by averaging out all non-essential variables. In short, we set ϕ k = ϕ for k = i, j, and average over ϕ to obtain a function Γ ij (ϕ i , ϕ j ) of ϕ i and ϕ j only, see SM text for details. The active driving force Q i (ϕ i ) of each cilium is uniquely determined by a reference condition, namely that the phase speed of this cilium should be constant, ϕ i = ω 0 , if the other cilia do not beat. This condition yields Together, Eqs.
(2), (3) and (4) give an equation of motion in implicit forṁ (5) The normalized hydrodynamic interaction γ ij (ϕ i , ϕ j ) between cilium i and cilium j characterizes the relative amount by which the motion of cilium j changes the phase speed of cilium i. Fig. 1(c) shows γ ij (ϕ i , ϕ j ) as function of the respective phases ϕ i and ϕ j of the two cilia. In short, the effective stroke of cilium j (π ϕ j 2π) will speed up cilium i (γ ij <0, blue colors) if cilium i is also in its effective stroke (π ϕ i 2π), but will slow down cilium i (γ ij >0, red colors) if cilium i is in its recovery stroke (0 ϕ j π). When one of the two cilia transitions from effective stroke to recovery stroke, or vice versa (i.e., ϕ i ≈ 0, π or ϕ j ≈ 0, π), that cilium moves slowly and the hydrodynamic interaction between the two cilia is weak, γ ij ≈ 0. We emphasize that γ ij (ϕ i , ϕ j ) is not simply a function of the phase difference ϕ i − ϕ j as in a classical Kuramoto model, but is much richer.
Numerical computations further show that γ ij is very small except for close neighbors; we therefore set γ ij = 0 except for close neighbors i and j, see Fig. 1(d). We can now rewrite the equation of motion equivalently in explicit form asΦ = Γ −1 · Q. With pre-computed Γ ij (ϕ i , ϕ j ) and Q i (ϕ i ) at hand, this explicit ordinary differential equation can be efficiently integrated for tenthousands of cilia beat cycles.
Metachronal wave solutions. We are interested in dynamic steady-state solutions of the equation of motion, Eq. (5). As a reference, we first re-visit the classical  34], whose periodic shape sequence has been parameterized by a 2π-periodic phase variable ϕ (color code). (b) Computed flow field u for this beat pattern (colors: |u(x)|, arrows: projection of u on yz-plane). (c) Normalized hydrodynamic interaction γij(ϕi, ϕj) = Γij(ϕi, ϕj)/Γii(ϕi) between a pair of cilia as function of their phases ϕi and ϕj: Positive values γij cause cilium i to beat slower, see Eq. (5). Separation vector of cilia bases, xj − xi = a (cos ψ ex + sin ψ ey), ψ=π/3. (d) Triangular lattice of cilia base points xj (dots). The color-code represents the root-mean-square average γ 2 Kuramoto model with local sinusoidal coupling [40,41] specifically, we consider a Kuramoto model of coupled phase oscillators with phases ϕ i at respective lattice positions x i and equation of motionφ j (t) = ω 0 − i =j c ij (ϕ i , ϕ j ) with coupling function c ij = ε sin(ϕ j − ϕ i ) for all pairs (i, j) of neighbors and c ij = 0 else. For this Kuramoto model, the steady-state solution are perfect plane traveling waves with wave vector k Here, k is one of the N reciprocal lattice points in the Brillouin zone of the oscillator lattice (with unit cell of N oscillators and periodic boundary conditions), see also Fig. 2(a). Note ω k = ω 0 for this simple Kuramoto model. In our cilia carpet model, the hydrodynamic interaction coefficients γ ij are not perfect sinusoidal functions, but a superposition of many Fourier modes. As a consequence, periodic solutions of cilia carpet dynamics are not perfect plane traveling waves as in Eq. (6). Nonetheless, we numerically find N periodic wave solutions Φ * k (t) of cilia carpet dynamics, where each Φ * k (t) is close to one of the N plane traveling wave Φ k (t) of Eq. (6). We will refer to Φ * k (t) as metachronal wave solutions. The global frequency ω k of these periodic solutions decreases with inverse wavelength |k|, see Fig. 2(a). The numerical dispersion relation is well approximated by ω k /ω k=0 ≈ 1 + β[cos(π|k|/k max ) − 1] with β ≈ 0.04 and k max = 4π/(3a), inline with analytical results for a slightly more general Kuramoto model [42] with c ij = ε sin(ϕ j −ϕ i +δ) involving an additional phase shift δ in the coupling, see SM text for details.
Linear stability analysis of metachronal wave solutions. To analyze the stability of metachronal wave solutions with respect to small perturbations, we map periodic solutions onto fixed points of a suitable Poincaré map [43]. We can then analyze the local stability of these fixed points using standard linear stability analysis. We first define a continuous global phase as the mean ϕ(t) = j ϕ j (t)/N for a continuous trajectory Φ(t) ∈ R N in phase space. Note that the mean of angular values can only be defined modulo 2π/N ; yet this ambiguity is resolved if we define ϕ for an entire time- Traditionally, wave directions are classified as symplectic, antiplectic, dexioplectic, laeoplectic, depending on the direction of k relative to the direction ey of the cilia effective stroke [8]. (b) Linear stability: Linear stability analysis for each k reveals that multiple solutions are linearly stable (green colors: stable metachronal wave solution, color represents relaxation time τ relax of the slowest decaying perturbation mode, normalized by beat period T0 = 2π/ω0 of single cilium; red: unstable). For the computation, we define a global phase ϕ and analyze the stroboscopic dynamics of the cilia carpet given by ϕ = 0 modulo 2π: fixed points Φ * k of this Poincaré map correspond to metachronal wave solutions, see left inset. (c) The relaxation time of the slowest-decaying perturbation for the dominant wave solution increases with system length as ∼ L 2 , resembling a dynamic Mermin-Wagner theorem for cilia carpets, which rules out global order in infinite systems. Lattice of 16 × 16 cilia; other parameters as in Fig. 1. continuous trajectory. We now define a Poincaré plane H by setting this global phase to zero, ϕ = 0, and a Poincaré return map L : H → H, corresponding to an increase of the global phase ϕ by 2π [i.e., a trajectory Φ(t) starting at Φ(0) = Φ 0 ∈ H intersects the shifted Poincaré plane H + 2π 1 at Φ 1 = L(Φ 0 ) + 2π 1], see inset on the left in Fig. 2 To determine whether a metachronal wave solution is stable, we linearize the Poincaré map at the correspond- The eigenvalues λ 1 , . . . , λ N −1 of ln(L k ) represent dimensionless Lyapunov exponents (whose real parts are proportional to inverse relaxation times), while the corresponding eigenvectors ∆ 1 , . . . , ∆ N −1 represent funda-mental perturbation modes. The fixed point Φ * k is linearly stable if Re λ i < 0 for all i. In this case, all perturbation modes ∆ i decay with respective relaxation times τ i = 2π/|ω k Re λ i |. A non-zero imaginary part of the Lyapunov exponents implies that perturbations decay in a spiral-like fashion to the fixed point Φ * k with period (2π) 2 /|ω k Im λ i |. We observe that multiple metachronal wave solutions are simultaneously stable: Fig. 2(b) reports the relaxation time τ relax = max τ i of the slowest decaying perturbation mode for stable wave solutions. The multistability of wave solutions is inline with previous observations in minimal models [9].
Global stability: one wave dominates. Although many metachronal wave solutions with different wave vectors k are simultaneously stable to small perturbations, we find that trajectories with uniformly sampled random initial conditions will predominantly converge to just one wave solution. The fraction of trajectories converging to Φ * k , equals the volume fraction of the basinof-attraction of Φ * k , which yields 86% for the dominant wave solution with wave vector k I , see Fig. 3

(a).
Slice-visualization of basins-of-attraction. To visualize basins-of-attractions of metachronal wave solutions, we additionally considered a specific set of initial conditions of the form ϕ j = −m · x j with "off-lattice" wave vectors m; these initial conditions correspond to a two-dimensional slice through the N -dimensional phase space, see Fig. 3(b). As expected, the majority of initial conditions converged to the dominant wave mode k I , while initial conditions m ≈ k in a small neighborhood of other stable modes k converged to the respective Φ * k . A magnification shows that the boundaries between basins-of-attraction are rough (and possibly fractal). Finally, a small number of initial conditions did not converge to any Φ * k within the simulation time [gray squares in Fig. 3(b)], but presumably converged to more exotic states, e.g., chimeras states consisting of multiple ordered domains [44], see SM text for examples.
Diverging relaxation time. We investigated cilia carpets of different size, and consistently found that the local stability patterns of metachronal waves remain similar to Fig. 2(b), see SM text. Similarly, we observe a single dominant wave solution for all system sizes tested, with corresponding wave vectors close to k I throughout. Nonetheless, in larger systems, perturbation modes with longer wavelengths and longer relaxation times appear. The relaxation time τ relax = max i τ i of the slowest-decaying perturbation mode for the respective dominant wave solution increases with system length L = max(L x , L y ) of the L x × L y -simulation domain approximately as see Fig. 2(c). While we demonstrate this power law only numerically for cilia carpets, one can in fact prove this power law analytically for a minimal Kuramoto model with local sinusoidal coupling, see SM text. This dynamic behavior parallels the Mermin-Wagner theorem from statistical mechanics for two-dimensional equilibrium systems with continuous symmetries [45]. For example, in the classical XY model of interacting spins in the plane with short-range interactions, so-called Goldstone modes appear; the energy-per-area of these these long-wavelength perturbation modes scales as 1/L 2 with system length L [46,47]. In a dynamic re-formulation, the relaxation times of these perturbation modes diverge as ∼ L 2 if we impose over-damped dynamics. In this sense, one may interpret Eq. (8) as a dynamic Mermin-Wagner theorem of a non-equilibrium system. [48] The analogy between synchronization and the XY model can be made more explicit for the classical Kuramoto model with local sinusoidal coupling [40]. Synchronization in presence of quenched frequency disorder. In real cilia carpets, the intrinsic beat frequencies of individual cilia will slightly differ. In a Kuramoto model with all-to-all coupling, a second-order phase transition occurs as function of a frequency disorder parameter, whereas in Kuramoto models with local coupling a synchronization transition can only be observed in finite systems [49,50].
We now investigate a cilia carpet, where each cilium has a slightly different intrinsic beat frequency ω i , with equation of motion given by Eq. (5), but with ω 0 replaced by ω i for cilium i, i.e.,φ i = ω i − j =i γ ijφj . Cilia beat frequencies are drawn from a normal distribution with mean ω 0 and standard deviation ∆ω > 0. [As a technical point, we rejected frequency sets whose sample standard deviation differed by more than ≈ 1% from ∆ω.] We are interested in the synchronization behavior of the cilia carpet as function of ∆ω, averaged over different frequency sets and initial conditions, see SM for details.
To characterize steady-state solutions, we introduce a generalized Kuramoto order parameter, see also [51] This order parameter r k is close to one, whenever the cilia phases approximately form a plane traveling wave Φ k (t) with wave vector k, i.e., ϕ j ≈ ϕ − k · x j . The inequality r k (Φ) > 2 −1/2 defines mutually disjoint neighborhoods for each k (each of which occupies only a tiny fraction < 10 −10 of the whole phase space). Fig. 3(c) shows the fraction of trajectories Φ(t) as function of ∆ω that both (i) converge to the neighborhood of a metachronal wave solution Φ * k (t) with r k [Φ(t)] > 2 −1/2 , and (ii) exhibit global frequency synchronization, i.e., phase differences between different cilia remain bounded. This definition for global metachronal coordination generalizes a previous definition for the case k = 0, which required both 'phase cohesiveness' and 'frequency synchronization' [41]. We find that the fraction of synchronized trajectories sharply decreases near a character-  [16]), by drawing 400 random initial conditions from a uniform distribution, of which 86% converged to one dominant wave mode kI, while 13% converged to the adjacent wave mode kII [introduced in Fig. 2(a)]. (b) Slice of sync basins. To visualize basins-of-attraction, we show limit points Φ * k for special initial conditions ϕj(t=0) = −xj · m with off-lattice wave vector m; this choice corresponds to a two-dimensional slice through N -dimensional phase space. Gray dots indicate initial conditions, for which trajectories did not converge to any Φ * k . Upon magnification, the boundaries of the basins-of-attraction appear rough, see inset to the left. (c) Frequency disorder. Relative size of basinsof-attraction for different metachronal wave modes as function of increasing quenched disorder ∆ω/ω0 of intrinsic cilia beat frequencies with ∆ω 2 ≈ ω 2 i − ωi 2 : synchronization is lost at a characteristic disorder threshold. For intermediate ∆ω, some realization display high order parameters r k ≥ 2 −1/2 for some k, but not all cilia adopt a common frequency, corresponding to a regime of partial synchronization (red). Parameters as in Fig. 2. istic value of frequency disorder, ∆ω c /ω 0 ≈ 2.5 × 10 −3 . This value likely depends on system size, as suggested by previous work on two-dimensional Kuramoto models with local coupling [49,50]. For intermediate values of ∆ω close to the transition point, ∆ω ≈ ∆ω c , we observe a fraction of trajectories that exhibit partial synchronization, i.e., trajectories satisfy condition (i) [large Kuramoto order parameter], but not condition (ii) [frequency synchronization], apparently because a few cilia did not synchronize and displayed phase drift instead.
Discussion. We analyzed global stability of metachronal synchronization in cilia carpets using a multi-scale model, and found that a single dominant wave solution has a basin-of-attraction that spans almost the entire phase space of initial conditions (generalizing early observations for oscillator rings [16]). The wave direction of this dominant metachronal wave solution encloses an angle of ≈ 60 • with the direction of the effective stroke of the cilia beat, which is close to the experimentally observed value ≈ 90 • , corresponding to a so-called dexioplectic wave [3]. The experimentally observed wavelength ≈ 11 µm is smaller than the wavelength of the dominant wave mode 2π/|k I | ≈ 34 µm predicted here; this discrepancy may simply be a consequence of the cilia density used in our model, which does not yet allow us to study smaller wavelengths.
Linear stability analysis showed that long-wavelength perturbations of the dominant synchronized state relax only slowly with relaxation time-scales that increase quadratically with system size. This dynamic behavior in a non-equilibrium system parallels the Mermin-Wagner theorem for two-dimensional equilibrium systems with continuous symmetries (such as the XY models of interacting spins in a plane) [45]. In these systems, longwavelength perturbations known as Goldstone modes appear in large systems, whose energy-per-area becomes arbitrarily small and hence their relaxation times diverge if we impose over-damped dynamics. Noise excites these Goldstone modes, which rules out global order in infinite systems. Based on the observed divergence of relaxation times, we expect a similar behavior for metachronal synchronization in cilia carpets [52]. The non-equilibrium dynamics in cilia carpets is thus different from other non-equilibrium dynamical models such as the Toner-Tu model of flocking birds [53]: in that two-dimensional model, global order is possible, because the active motion of agents results in a continuous exchange of neighbors. In contrast, the set of neighbors remains fixed in the cilia carpet model.
Our analysis became possible by a multi-scale simulation approach that describes beating cilia as phase oscillators [14,30,54]. We describe the cilia carpet as an array of phase oscillators, similar to a Kuramoto model with local coupling [55], yet where direction-dependent coupling functions are calibrated from detailed hydrodynamic simulations using a measured cilia beat pattern from Paramecium [3, 34]. Our approach tries to combine the mathematical elegance of popular minimal models that idealize beating cilia as orbiting spheres [9,13,[25][26][27][28][29][30][31], and the quantitative predictive power of full-scale numerical simulations that are computationally expensive [6,32,33].
For technical reasons, cilia spacing in our model (a=18 µm) is larger than in real cilia carpets (2 µm [3]), similar to the dilute limit considered in most theoretical studies. Therefore, we underestimate hydrodynamic interactions, which are expected to scale as inverse cubed distance of cilia distance in the far field [14,36]. In dense cilia carpets, near-field hydrodynamic interactions can change though and even steric repulsion can become important. As a consequence, we likely underestimate the characteristic value of disorder of intrinsic beat frequencies at which synchronization is lost.
Our model could be extended to systems consisting of separated cilia bundles found in airway epithelia [13]. Future refined models may include internal friction of cilia beating [13,24,56], and cilia waveform compliance [27,57], which are expected to reduce and increase synchronization strength, respectively. A putative role of basal coupling of cilia contributing to synchronization [12,57,58] remains open for cilia carpets, and has therefore not been included here. Real cilia carpets are characterized also by quenched disorder of cilia position, and non-perfect alignment of cilia [7], which should reduce the regularity of emergent metachronal waves. Intriguingly, some disorder of metachronal coordination might actually be beneficial for transport of suspended particles, e.g., virus clearance from ciliated airways [59].
AS and BMF are supported by the German National Science Foundation (DFG) through the Microswimmers priority program (DFG grant FR3429/1-1 and FR3429/1-2 to BMF), a Heisenberg grant (FR3429/4-1), as well as through the Excellence Initiative by the German Federal and State Governments (Clusters of Excellence cfaed EXC-1056 and PoL EXC-2068). We thank Christa Ringers and Nathalie Jurisch-Yaksi (NTNU), as well as all members of the 'Biological Algorithms' group for stimulating discussions.
Data availability. Python code used to generate results in this manuscript is available in public repositories [60] Data availability. We deposited code used to generate results in this manuscript as Python packages in three publicly accessible repositories, specifically: (i) digitalization of three-dimensional cilium beat from stereographic recordings, including coordinate files of the final cilium beat pattern (ii) routines for generating the triangulated mesh of cilia and boundary surfaces, and for solving the hydrodynamic Stokes equation and computing generalized hydrodynamic friction coefficients, (iii) routines for numerical integration of the equation of motion Eq. (5), as well as linear stability, global stability, and additional analyses [60].
Applicability of Stokes equation. In the presence of a no-slip boundary surface, the flow field generated by a static force monopole decays as 1/d 3 as function of distance d parallel to the plane in the limit of zero Reynolds number [36]. For an oscillating force monopole, whose amplitude oscillates with angular frequency ω 0 , the linearized Navier-Stokes equation predicts that the leading order singularity of the induced flow field becomes exponentially attenuated beyond a characteristic distance δ = [2µ/(ρω 0 )] 1/2 , where µ is the dynamic viscosity of the fluid, and ρ its density; for distances d δ, the flow field decays as 1/d 3 far from boundaries and as 1/d 5 close to a plane boundary [35,61,62]. Using a typical cilia beat frequency ω 0 /2π = 32 Hz and parameters for water at room temperature, we estimate δ ≈ 100 µm. Thus, hydrodynamic interactions from nearby cilia should contribute most to synchronization by hydrodynamic interactions.
Additionally, the flow induced by an oscillating force monopole exhibits a distance-dependent phase lag. For neighboring cilia, however, this phase lag is small. Correspondingly, we employ the approximation of zero Reynolds number and compute the interactions between nearby cilia using the Stokes equation.
Mesh generation. Cilia are modeled as slender curved rods with a radius of 0.125 µm with prescribed centerline, using a digitalization of cilia beat pattern from unicellular Paramecium recorded by [3] and represented by [34]. The simulation geometry representing a local region of a cilia carpet consists of a boundary surface modeled as a disk of radius 60 µm represented as a triangular mesh, whose upper face is coplanar with the xy plane containing the cilia base points x j . Triangulated meshes of the shape-changing cilia are anchored to the upper surface of this disk at the respective base points. For numerical accuracy, we performed local mesh refinement of the mesh in the vicinity of the base points, resulting in a mesh with a total of typically 8 · 10 3 node points, see Fig S1(a). This cilia carpet is immersed in an unbounded, Newtonian fluid with dynamic viscosity µ = 10 −3 Pa s (corresponding to viscosity of water at 20 • C). For details on mesh generation, see [14].
To solve for the surface density of hydrodynamic friction forces resulting from a shape change of the cilia, we employ fastBEM, a fast multipole solver for the Stokes equation [63].
Generalized hydrodynamic friction coefficients. We compute hydrodynamic interaction coefficients Γ ij = Γ ij (ϕ i , ϕ j ) in a series of numerical experiments, where only one cilium with index j beats at a constant frequency ω 0 , while other cilia are standing still, i.e.,φ k = 0 for k = j. Using the hydrodynamic solver, we obtain surface force densities f j (x) on the combined surface S of all cilia and the boundary surface. We compute the hydrodynamic friction coefficients Γ ij as where w i = ∂x/∂ϕ i is a rate of displacement of the surface S corresponding to a change of ϕ i , while all other ϕ k , k = i, do not change. Note that we can restrict the surface integral in Eq. (S1) to the surface S i of cilium i, since w i (x) = 0 on the rest of the surface S \ S i . For each relative orientation of cilia d = x j − x i , we computed generalized hydrodynamic friction coefficients Γ ij = Γ ij (ϕ i , ϕ j ) characterizing hydrodynamic interaction between cilia. Specifically, we sampled the respective phases ϕ i and ϕ j of the two cilia equidistantly with step size ∆ϕ = 2π/20, while the phases of all other neighboring cilia were set to a constant value of either 0, π/2, π or 3π/2, see Fig. S1(b). We then averaged over the constant phase the of other cilia, by fitting a truncated bi-variate Fourier series in ϕ i , ϕ j , of maximal order 4 (corresponding to (2 · 4 + 1) 2 = 81 Fourier terms for each Γ ij ). In rare cases (< 1%), the hydrodynamic solver would unexpectedly fail to converge to the prescribed tolerance (10 −3 ); these data points were excluded from the fit. The self-friction coefficients Γ ii (ϕ i ) are computed in a similar way, with one cilium phase sampled with step size ∆ϕ = 2π/20, and averaged over a constant phase of its 6 neighboring cilia (only 2 ·4 + 1 = 9 terms in Fourier series are kept), see Fig. S1(c). This provided 'look-up tables' for subsequent dynamic simulations of the equations of motion of the cilia carpet, Eq. (3).
While these hydrodynamic simulations consider only a finite cilia array, they are sufficient to calibrate relevant nearest-and next-to-nearest-neighbor hydrodynamic interactions, which are later used to simulate larger cilia carpets with periodic boundary conditions.

S2
Visualization of hydrodynamic interaction. For Fig. 1(c), we computed the pairwise normalized hydrodynamic interaction using the Fourier sum representation of Γ ij and Γ ii described above. For Fig. 1(d), we computed the root-mean-square average of γ ij (ϕ i , ϕ j ) as γ 2 ij 1/2 = (2π) −1 [ dϕ i dϕ j γ 2 ij (ϕ i , ϕ j )] 1/2 for nearest and next-tonearest neighbors. As a technical point, for some nextto-nearest neighbors (specifically, for distance d = √ 3a and direction angles ψ = ±π/6, ±5π/6 relative to x-axis, where γ ij is already very small), more than 1% but still less than 5% of the hydrodynamic computations did not converge to the prescribed tolerance. For the visualization of γ 2 ij 1/2 in Fig. 1(d), we included all data points in the fit of the Fourier sum, including those for which the hydrodynamic computation did not converge. Note that these problematic next-to-nearest neighbor interactions were not included in the final dynamic computations because the corresponding hydrodynamic interactions are already very small.
Approximation of pairwise interactions. We highlight the two simplifications underlying our effective multiscale simulation framework. (i) We introduced a minimal set of effective degrees of freedom, and constrained the full dynamics to these degrees of freedom. With these constraints imposed, the balance Eq. (3) is exact. (ii) We approximated the N -body hydrodynamic interaction as a superposition of pairwise interactions and introduced a distance cut-off. While the force balance is not exact anymore with these approximations, we numerically confirmed that it still holds to very good accuracy. Thus, the force balance equation with approximation of pairwise interactions reads Here, N i is the set of neighbors of cilium i, which includes all six nearest neighbors (at distance a = 18 µm) and two next-to-nearest neighbors located at ±d e y with d = √ 3a (corresponding to direction angle ψ = ±π/2), i.e., located along direction of the cilia effective stroke, where hydrodynamic interactions are the strongest, see Fig. 1(d). Next-to-nearest neighbor interactions along the other directions are much weaker, and were therefore not included in the final simulations for reasons of computational performance. Initial simulations showed that including these interactions with next-to-nearest neighbors only slightly changed quantitative results, and did not affect any of our qualitative conclusions. Active cilia driving forces. For our choice of reference condition, the active driving forces Q i (ϕ i ) are given by corresponding to a single cilium that beats at a constant frequency (while its neighbors are at rest and only act as obstacles for the fluid). Equation of motion. Numerically, we solve the equation of motion Eq. (5) in the forṁ The coupling functions Γ ij depend only on the phases ϕ i and ϕ j and the relative positions of cilia i and j, allowing for efficient storage. Alternatively, we could introduce the generalized mobility matrix M = Γ −1 , and the vector of active driving forces Q with components Q j (ϕ j ). The equation of mo-tionΦ = M · Q can then be written as a system of N coupled phase oscillatorṡ

S3
with coupling functions c ij = (M · Q) ij − ω 0 δ ij . Diagonal entries c ii characterize a modulation of beat frequency due to the presence of nearby cilia. As consequence of the no-slip boundary surface, hydrodynamic interactions decay with inverse cubed distance close to the surface [36]. Thus, in the limit of low cilia density with a where denotes cilia length, we have c ij ∼ ( /a) 3 for neighbor cilia with j ∈ N i . Yet, even for j / ∈ N i , c ij is in general non-zero albeit small, decaying at least as ( /a) 6 . Thus, although the generalized friction matrix Γ is sparse (given the approximation of including only nearest-neighbor interactions), the generalized mobility matrix M will be non-sparse in general.
Eq. (S5) represents a generalized Kuramoto model with local coupling. Indeed, if we set c ij = ε sin(ϕ i − ϕ j ) for nearest neighbors, and c ij = 0 else, we would obtain the classical Kuramoto model with local sinusoidal coupling in two space dimensions.
Numeric integration of equation of motion. We used a 4(5)-Runge-Kutta scheme with adaptive time-step (Python package scipy) to numerically integrate the deterministic equation of motion, Eq. (S5). We used numerical tolerance 10 −8 to determine fixed points and Lyapunov exponents from the linear stability analysis, and a numerical tolerance of 10 −6 for all other computations. Intersections with the Poincaré plane H defined by ϕ = 0 were detected using the integrated event handler. In each time-step, we compute the right side of the equation of motionΦ = Γ −1 · Q using a sparse linear solver.
Reciprocal lattice of metachronal wave vectors and Brillouin zone. We introduce basis vectors d x and d y of the reciprocal lattice defined by a tiling of the plane by copies of the unit cell of N cilia where L x = N x a and L y = √ 3N y a/2 denote the length of the rectangular unit cell in x and y direction, respectively. Any wave vector k in the reciprocal lattice can be written as with integers n x , n y ∈ Z, or, alternatively, with vector components k x = n x 2π/L x and k y = n y 2π/L y with respect to the normalized unit vectors e x = (1, 0) T and e y = (0, 1) T . The regular spacing of cilia at lattice positions x j inside the unit cell defines a Brillouin zone K: in the case of a triangular lattice, this Brillouin zone can be chosen as a hexagon with edge length k max = 4π/(3a), see Fig. 2(a). This Brillouin zone contains N = |K| unique wave vectors. Any other wave vector k of the reciprocal lattice can be mapped either inside or on the border of this hexagon using the equivalence relation exp(i k · x j ) = exp(i k · x j ) for all j. A visualization of the dominant wave mode k I is shown in Fig. S2. For a classical Kuramoto model with sinusoidal nearest-neighbor coupling, each wave vector k ∈ K defines a periodic solution Φ k with components ϕ i = ω 0 t − k · x i (also called k-twist [64,65] or splay states [66] in one-dimensional oscillator chains), see also section on the Kuramoto model below. For the cilia carpet model considered in the main text, we find periodic solutions that deviate slightly from these perfect traveling waves.
Numeric search for periodic solutions. To find periodic solutions Φ * k (t) of the generalized Kuramoto model given by Eq. (3), we numerically searched in the vicinity of the periodic solutions Φ k (t) of the classical Kuramoto model. Specifically, we searched for fixed points Φ * of the Poincaré map L for the Poincaré plane H given by ϕ = 0, where ϕ = j ϕ j /N denotes the global phase Here, Φ 0 = Φ(t 0 ) ∈ H is the start point of a trajectory Φ(t) that intersects the shifted Poincaré plane H + 2π 1 at Φ 1 = Φ(t 1 ), i.e., ϕ(t 0 ) = 0 and ϕ(t 1 ) = 2π. Numerically, it turned out to be easier to start also with initial phase vectors that had a non-zero global phase, i.e., ϕ(t 0 ) = ϕ 0 and ϕ(t 1 ) = 2π + ϕ 0 . We found fixed points Φ * by numerically searching for zeros of the following vector function, where the last term effectively restricts the search to the Poincaré plane H Note that the condition D(Φ 0 ) = 0 actually implies both L(Φ 0 ) − Φ 0 = 0 and ϕ(Φ 0 ) = 0. Hence, D(Φ * ) = 0 yields a fixed point Φ * ∈ H with zero global phase. By running the numerical search algorithm N times with start vectors Φ 0 given by plane waves ϕ i = −k · x i for each k ∈ K, we found N different fixed points Φ * k . The Kuramoto order parameters r k defined in Eq. (9) evaluated at the fixed points almost equal one with r k (Φ * k ) > 1 − 2 · 10 −3 . This confirms that these fixed points correspond to periodic solutions Φ * k (t) that are indeed close to perfect traveling waves.
FIG. S3: Linear stability analysis for systems of different size. We performed linear stability analyses for each wave vector k inside a Brillouin zone for systems of different sizes similar to Fig. 2(b) in the main text. In all three cases, stability patterns are similar: left: 8 × 8 carpet with N = 64 cilia, middle: 16×16 carpet with N = 256 cilia, right: 20×20 carpet with N = 400 cilia. Green colors represent max Re λj of respective Lyapunov exponents λj for linearly stable wave modes k; red dots represent modes that are linearly unstable.
The absolute values of eigenvalues tend to zero as system size increases, as discussed in Fig. 2(c) in the main text (which reports the relaxation time τ relax = 2π/|ω k max Re λj|).
Linear stability analysis We numerically find the linearized Poincaré map L k near a fixed point Φ * k [see Eq. (7)], by computing the Poincaré map L for small perturbations. Specifically, we apply small perturbations ∆  Fig. S3 shows results of a linear stability analysis for cilia carpets of different sizes.
Basins-of-attraction. To estimate the relative size of basins-of-attractions, we computed n = 400 trajectories with initial conditions given by uniformly sampled random phase vectors. For each trajectory, we integrated the equation of motion Eq. (5) for m = 4000 beat cycles (corresponding to an integration time t m ≈ 4000 T 0 ). All n trajectories converged to a neighborhood of a wave Φ k for a suitable wave vector k, as determined by a Kuramoto order parameter r k > 0.99. Additionally, we observed that each of the n trajectories apparently converged to a fixed point, by checking that the Euclidean norm of the change of the phase vector d l = N −1/2 L l [Φ(0)] − L l−1 [Φ(0)] 2 during one beat cycle was decreasing and sufficiently small after m = 4000 cycles, d m < 2 · 10 −4 .
In fact, all trajectories converged to just five waves (all of which are very close to each other in terms of both wave direction and wavelength); the majority of trajectories converged to either k I (86% ± 2%) or k II (13% ± 2%) (introduced in the main text Fig. 2(a)). The error e k was computed as the standard error of a Bernoulli trial [18] where B k ∈ [0, 1] is the relative size of the basin-ofattraction, and n is the total number of trajectories. Slice-visualization of basins-of-attraction. In Fig. 2(b), we additionally visualize convergence for a specific set of initial conditions of the form ϕ j = −m · x j , with "off-lattice" wave vectors m ∈ K (thus m does not necessarily respect the periodicity of the lattice). Note that these special trajectories were not used in calculating the relative size of the basins-of-attraction, as the initial conditions were not drawn randomly.
Each of the trajectories was integrated until it converged to one of the fixed points Φ * k (using the same numerical convergence criterion as detailed in the section Basins-of-attraction above). Intriguingly, for some of the initial conditions, trajectories did not converge do any of the fixed points Φ * k [gray squares in Fig. 2(c)]. Some of these special trajectories apparently became attracted to some other fixed point Φ * different from any of the Φ * k , k ∈ K (i.e., the distances d l introduced in the section above approached zero, but all Kuramoto order parameters remained below the threshold at the end of the integration time, r k < 0.99 for all k). Other special trajectories were not attracted to any fixed point Φ * even after a long integration time t ≈ 10 4 T 0 . Nonetheless, the dynamics of these later trajectories had become stationary in the sense that the Kuramoto order parameters r k , k ∈ K, either did not change in time anymore or oscillated in a regular way. Visual inspection revealed that these initial conditions had become attracted to more exotic states, such as chimera states (i.e., states with at least two ordered sub-domains) [44], see Fig. S4 for an example. However, the combined relative size of the basins-of-attraction of these exotic states is negligible; therefore, these states are not in focus of this study.
As expected, the basin of the dominant wave vector k I comprises a large portion of initial conditions in the slice of phase space shown in Fig. 3(b). Specifically, most initial conditions corresponding to unstable waves vectors k became attracted to the dominant wave k I . On the other hand, initial conditions in the vicinity of a stable wave vector k different from k I are likely to become attracted to this wave k, see the magnified region in Fig. 3(b). For t ≈ 5 · 10 3 T0, the order parameters reached steady state.

Quenched frequency disorder
In Fig. 3(c), we show the fraction of synchronized trajectories as a function of a frequency disorder parameter ∆ω. Specifically, we drew sets of random intrinsic beat frequencies Ω = (ω 1 , . . . , ω N ), where the intrinsic frequency ω i of cilium i was drawn from a normal distribution with mean ω 0 and standard deviation ∆ω. As a technical point, the (biased) sample variance where µ(Ω) = N −1 i ω i denotes the sample mean, may vary from its expectation value ∆ω 2 . We rejected frequency sets, where Var(Ω) 1/2 differed from ∆ω by more than 1%. Without this rejection (which amounts to about 82% of frequency sets), the synchronization transition in Fig. 3(c) would appear more gradual.
For each value of ∆ω > 0 considered, we first generated m = 25 valid frequency sets. For each frequency set, we then integrated n = 10 trajectories, starting from a fixed sub-sample of initial conditions. This sub-sample had been selected before from a larger sample of random initial conditions (uniformly distributed), such that the previously determined relative sizes of basins-of-attraction for the case without frequency disorder was faithfully reproduced (for nine initial conditions, the trajectories converged to wave k I and for one initial condition, the trajectory converged to k II for ∆ω = 0). Using only a small number of initial conditions reduced computation times considerably.
For different frequency set Ω, periodic solutions of the system (and corresponding fixed points of the Poincaré map) will slightly differ from the periodic solutions found for the case ∆ω = 0. Therefore, in order to compute the relative size of basins-of-attractions in Fig. 3(c), we employ a sufficiently large neighborhood of the plane wave solution Φ k (t). More precisely, we say that a trajectory Φ(t) synchronized to wave k if the following two conditions are satisfied (i) The respective Kuramoto order parameter was large at the end of the integration time, r k > √ 2/2.
(ii) Φ(t) converges to a fixed point of the Poincaré map.
Each of the trajectories was integrated until condition (ii) was met or a maximum integration time t ≈ 1.6 × 10 4 T 0 was reached. To check convergence to a fixed point, the same criterion as in section 'Basins-of-attraction' was used.
Without frequency disorder, ∆ω = 0, conditions (i) and (ii) are essentially equivalent, except for few rare cases, where initial conditions converged to exotic states, e.g., chimera states. However, in the case of frequency disorder with ∆ω > 0, the two conditions (i) and (ii) are no longer approximately equivalent, and we observe trajectories that satisfy condition (i) but not (ii), especially close to the synchronization transition. We refer to these trajectories with partial synchronization [red color in Fig. 3(c)].

Kuramoto models with local coupling
For the convenience of the reader, we review basic facts on the classical Kuramoto model with local coupling, part of which can be found in the standard literature [21].
One-dimensional chain of phase oscillators with nearest-neighbor sinusoidal coupling We consider a one-dimensional chain of N coupled phase oscillators with periodic boundary conditions. The oscillators in this ring topology are supposed to have equal angular frequency ω 0 and are coupled to their neighbors by a symmetric sinusoidal coupling with total coupling strength K For notational convenience, oscillator indices are considered modulo N (i.e., oscillator number N is coupled again to oscillator number 1). We assume a positive synchronization strength K > 0; correspondingly, the in-phase synchronized state is stable.
Traveling waves with angular wave number k define where k = 2π m/N for some integer m ∈ Z.
The fundamental perturbation modes of the Poincaré map for these periodic solutions Φ k are simply the Fourier modes for the chain with angular wave number ν where ν = 2π n/N for some n = 1, . . . , N − 1 (n = 0 would correspond to a trivial phase shift). The corresponding eigenvalues of the linearized Poincaré map ln L, which we call dimensionless Lyapunov exponents, read This can be proven by substituting the perturbation Eq. (S16) and keeping only terms to linear order. The periodic solution for wave number k is linearly stable if and only if the real parts of all eigenvalues λ kν are strictly negative; hence, according to Eq. (S17), exactly the solutions with |m| < N/2 are linearly stable.
We can now read off the dimensionless Lyapunov exponents of the slowest decaying mode for each stable periodic solution and find max ν =0 Here, we introduced a system length L = N a, where a is the spacing between oscillators. Thus, the long wavelength perturbations (|k| → 0) are indeed those that decay the slowest, with a decay rate that scales as the inverse square of system length L = N a.
In the main text, we describe a similar scaling for the relaxation time τ relax , which is inversely proportional to Lyapunov exponent max λ j of the slowest decaying perturbation mode, for the periodic solution Φ * kI (t) corresponding to the dominant wave mode k I , see Fig. 2(d).
In addition, we numerically checked that the largest dimension L = max (L x , L y ) dominates the scaling also if N x = N y (both for the Kuramoto and the cilia carpet models).
Dispersion relation for the one-dimensional Kuramoto model with local coupling As a generalization of Eq. (S14), we can consider the Sakaguchi-Kuramoto model with local coupling [42] c ij (ϕ i , ϕ j ) = ε sin(ϕ j − ϕ i + δ), i.e, with additional phase shift δ (as introduced in the main text). This generalized one-dimensional Kuramoto model can be written aṡ where K = 2ε cos δ and U = 2ε sin δ. We make an Ansatz of traveling waves with frequencies ω k and k = 2π m/N for some integer m ∈ Z. Substituting this Ansatz into Eq. (S19), yields periodic wave solutions with frequencies ω k with where ω k=0 = ω 0 + U and β = U/ω k=0 , Eq. (S21). Thus, an additional cosine term in the coupling function c ij causes a characteristic frequency dispersion relation. The stability of wave solutions, however, is not altered, as can be shown analogous to the previous section.
Kuramoto model with nearest-neighbor sinusoidal coupling in d dimensions More generally, we can consider a Kuramoto model of phase oscillators with identical frequencies on a cubic lattice with lattice spacing a and lattice positions x i in d-dimensional space and local sinusoidal coupling. Each oscillator with phase variable ϕ i is coupled to its 2d nearest neighbors (enumerated by an index set N i ) with total coupling strength K We assume periodic boundary conditions with system size N 1 × . . . × N d .
Let N i = max j N j be the number of oscillators along the longest direction of the N 1 × . . . ×N d -unit cell. The slowest decaying perturbation mode is then m max = 2π/L e i , where we introduce system length L = N i a. For the Lyapunov exponent of the slowest decaying perturbation mode m of the dominant wave solution k = 0, we thus find, analogous to the one-dimensional case treated above to leading order in a/L, where L = a max{N 1 , . . . , N d } denotes system length. This maximal Lyapunov exponent sets a relaxation time of the dominant wave solution, τ relax = T 0 |max λ 0m | −1 .

Relation to XY model
One can map the Kuramoto model with identical phase oscillators and sinusoidal coupling, Eq. (S22), to an equilibrium system by switching to a co-rotating frame with variables θ i = ϕ i − ω 0 t. Specifically, we consider the Hamilton of the classical XY model 1 For the calculation, notė and consider the over-damped dynamics Here, γ denotes an effective friction coefficient. Eq. (S27) is equivalent to Eq. (S22) for Fixed points θ * k of Eq. (S27) [over-damped XY model] correspond exactly to periodic solutions Φ * k (t) = ω 0 t 1 + θ * k of Eq. (S22) [Kuramoto model with local coupling]. For small perturbations ε∆ from a stable fixed point θ * k , we can approximate the Hamiltonian H as a harmonic potential where ∆H = ∆·∇ 2 H |θ=θ * k ·∆ † and † denotes the complex conjugate of a transposed vector. We can interpret ∆H either an effective spring stiffness along the direction of the perturbation ∆, or as a normalized energy penalty of the perturbation mode ∆. We have a direct relationship between the Lyapunov exponents λ 0m of the Kuramoto model for the dominant wave solution k = 0, as given in Eq. (S24), and the energy penalties ∆H m = ∆H(∆ m ) of the fundamental perturbation modes ∆ m defined in Eq. (S23). A short calculation shows 2 Here, T 0 = 2π/ω 0 is the period of the periodic solutions. The Hamiltonian H possesses O(2)-symmetry; any spontaneous "magnetization" with | e iθj | > 0 corresponds to spontaneous symmetry breaking. For d ≥ 3 space dimensions (i.e., Λ ⊂ R d ), the classical XY model is known to exhibit a conventional phase transition with spontaneous magnetization below a critical temperature T c . For d = 2 dimensions, there is no long-range order at any finite temperature, and thus no conventional phase transition. This is a consequence of the famous Mermin-Wagner theorem that rules out long-range order

S8
in two-dimensional systems with local coupling and continuous symmetries [45]. In these systems, the energy penalty for long-wavelength perturbations of the ordered ground state is independent of system size; hence these Goldstone bosons become thermally excited at any finite temperature. Nonetheless, for d = 2, the classical XY model exhibits a so-called Kosterlitz-Thouless transition, from a disordered high-temperature state with exponential decay of spatial correlations, to a quasi-ordered lowtemperature state with algebraic decay of spatial correlations [46], at a critical temperature k B T c /J ≈ 0.89 [47].