††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Emergence of metachronal waves in a chain of symmetrically beating filaments

Narina Jung  Korea Institute for Advanced Study, Seoul 02455, Korea  Won Kyu
Kim  Korea Institute for Advanced Study, Seoul 02455, Korea  Changbong Hyeon
Korea Institute for Advanced Study, Seoul 02455, Korea

###### Abstract

Recent experiments have shown that metachronal waves (MCWs) can emerge from a
chain of symmetrically beating nematodes aligned at the edge of sessile
droplets. Our study, employing a coupled elastohydrodynamic model of active
filaments, elucidates that a misalignment caused by a tilt against the
bounding wall disrupts the synchronization and generates a constant time lag
between adjacent filaments, leading to MCWs. The MCWs, enhancing the fluid
circulation, achieve their maximum thermodynamic efficiency over the same
range of tilt angles observed in the nematode experiments.

Introduction. Collective dynamics observed at the scale of animal population
are an outcome of local dynamics and coupling between individual agents [1, 2,
3, 4, 5, 6]. Among them, long wavelength modulations of collective beating
behavior called _metachronal waves_ (MCWs), characterized by a constant phase
or time lag between the dynamics of adjacent agents, enhances the fluid flow
and transport of nutrients essential for a host of microorganisms. As a
result, MCW has long been a topic of abiding interest [7, 8, 9, 10, 11, 12,
13].

Since the Taylor’s seminal works [7, 8], long-range hydrodynamic coupling has
widely been appreciated as the key interaction for the coordination and
generation of MCWs among filamentous micro-swimmers in low Reynolds number
environment [14, 15, 16, 17, 9, 10, 18]. Significance of short-range steric
interactions was also highlighted in the gait synchronization and generation
of collective waves with constant phase lags among _C. elegans_ [19, 20].
Applying transverse external flow to symmetrically undulating filaments,
Guirao and Joanny have suggested that breaking the left-right (LR) symmetry in
the filament dynamics as well as including the hydrodynamic interactions (HIs)
is necessary to generate MCWs [21].

Peshkov _et al._ have recently demonstrated formation of MCWs in a chain of
_Turbatrix aceti_ (_T. aceti_) aligned along the contact line of a water
droplet [22, 23]. In more broadly discussed MCWs in cilia, individual cilium
exhibits two-phase asymmetric beating consisting of power and recovery
strokes. In contrast, despite its functional importance, the MCW emerging from
symmetrically beating nematodes is physically not obvious and its origin
remains unclear. To elucidate how MCWs arise from the 1D array of nematodes,
we extend a phenomenological elastohydrodynamic model of a single active
filament [24] to many filaments, and study the physics underlying the
formation of MCWs and its effect on the flow of surrounding fluid.

Theoretical model. We consider that a 1D array of $N$ slender filaments with
radius $a$ and length $L$ are aligned along the $y^{\prime}$-axis, separated
by $d$ ($a\ll d\ll L$), with a tilt angle $\theta$ against $x^{\prime}$-axis
(see Fig. 1). The force-balance equation for $i$th filament coupled with its
neighbors and the bounding wall is written in a non-dimensionalized form (see
Supplemental Materials for details):

$\displaystyle\partial_{t}h^{i}(x^{i},t)=\mathcal{N}(h^{i})+\mathcal{F}_{\rm
HI}(h^{i})+\mathcal{F}_{\rm SI}(h^{i})+\mathcal{F}_{\rm WI}(h^{i})$ (1)

with the boundary conditions, $h(0,t)=\partial_{x}^{2}h(x,t)|_{x=0}=0$ and
$\partial_{x}^{2}h(x,t)|_{x=L}=\partial_{x}^{3}h(x,t)|_{x=L}=0$. Here
$h^{i}(x^{i},t)$ is the transversal displacement (parallel to the $y$-axis) at
position $x^{i}$ along the centerline of $i$th filament and at time $t$. In
$\mathcal{N}(h^{i})=[-c\partial_{x}h^{i}-2\partial^{2}_{x}h^{i}+\left(\partial^{2}_{x}h^{i}\right)^{3}]-\partial^{4}_{x}h^{i}$,
the three terms in the parenthesis model the effect of an active force
originating from a mechanical feedback of molecular motors [25, 26]. The
parameter $c$ in the first term, denoting a dimensionless speed of wave along
the centerline, reflects the filament activity.

To make the hydrodynamic forces, $\mathcal{F}_{\rm HI}(h^{i})$, acting on the
$i$th filament more tractable and interpretable, we consider the resistive
force theory [27, 24], and simplify $\mathcal{F}_{\rm HI}(h^{i})$ to

$\displaystyle\mathcal{F}_{\rm HI}(h^{i})=\sum\limits_{j\neq i}f^{j\rightarrow
i}_{\rm HI}/\zeta_{\perp}=\epsilon_{h}\sum\limits_{j\neq
i}\partial_{t}h^{j}(x^{j},t)$ (2)

where $x^{j}$ denotes the position along the $j$th filament. The coefficient
for the strength of HI is given by

$\displaystyle\epsilon_{h}=\ln{(d/L)}/\ln{(a/L)}+(1/2)(a/d)^{2}/\ln{(a/L)}$
(3)

with the first and the second terms from the Stokeslet and doublet,
respectively [28, 29]. In our system with $L=100a$ and $d=10a$, the
coefficient is $\epsilon_{h}=0.4989$. The terms $\mathcal{F}_{\rm SI}(h^{i})$
and $\mathcal{F}_{\rm WI}(h^{i})$ model the repulsion between the filaments
and the bounding wall effect on the filaments. Inclusion of SIs in Eq. 1 is
necessary in stabilizing the in-phase synchronized state against otherwise the
equally stable anti-phase synchronized state under hydrodynamic coupling,
especially for an array of filaments satisfying the condition $d<\langle\delta
h\rangle$ where $\langle\delta
h\rangle\equiv\left[(1/T)\int_{0}^{T}dt(1/L)\int_{0}^{L}dx(h(x,t)-\langle
h\rangle)^{2}\right]^{1/2}$ (see SM Text and Fig. S1 for details). We study
the collective dynamics of $\\{h^{i}(x^{i},t)|i=1,2,\ldots,N\\}$ by
numerically solving Eq.1 with varying $N$ and tilt angle $\theta$ [30].

Figure 1: Coordinate system in two-dimension used to study the beating
dynamics of $N$ filaments aligned in parallel to the $x$-axis with a finite
tilt angle $\theta(=\angle x^{\prime}Ox=\angle y^{\prime}Oy)$ towards the
boundary line ($\overline{Oy^{\prime}}$). The coordinate of the $i$th filament
in the reference frame $(x,y)$ is associated with the local coordinates
$(x^{i},h^{i})$ with $x^{i}\in[0,L]$ as $(x^{i}+(i-1)\Delta
l,h^{i}(x^{i},t)+(i-1)d)$ where $\Delta l=d\tan{\theta}$ is the size of
mismatch. The head position of the first filament ($i=1$) is set to the origin
$(0,0)$.

In what follows, we first characterize the single filament dynamics and next
explore how the inter-filament coupling modulates the symmetrically beating
dynamics of each filament and contributes to the formation of MCWs.

Single filament dynamics. For a single filament ($N=1$) without tilt
($\theta=0^{\circ}$), the wall interaction (WI) is minimal, which reduces Eq.
1 to $\partial_{t}h(x,t)=\mathcal{N}(h)$, and yields a periodic solution of
$h(x,t)\equiv h_{o}(x,t)$ for $t\gg 1$, where the subscript $``o"$ denotes the
waveform solution for a single filament ($N=1$) at $\theta=0^{\circ}$ which we
use as a reference waveform throughout the paper (Fig. 2(b), see also SM Movie
1). One of the key features of $h_{o}(x,t)$ is that it preserves LR symmetry
such that both $h_{o}$ and $-h_{o}$ are the solution of Eq.1, while displaying
a growing amplitude from the head to tail (Fig. 2(a)). The $x$-dependent beat
period, $T(x)$, can either be read off from the kymograph (Fig. 2(b)) that
satisfies $h(x,t)\approx h(x,t+T(x))$, or be obtained as
$T(x)=\operatorname*{argmax}_{T}S(x,T)$ from the power spectrum
$S(x,T)\sim\int h(x,t)e^{-2\pi it/T}dt$. The spectrum $S(x,T)$ (Fig. 2(c))
shows a stripe at $T(x)\approx 10\equiv T_{o}$, effectively independent of
$x$. Similarly, the wavelength at time $t$ is either estimated through
$h(x,t)\approx h(x+\lambda(t),t)$ or is obtained as
$\lambda(t)=\operatorname*{argmax}_{\lambda}F(\lambda,t)$ from the power
spectrum $F(\lambda,t)\sim\int_{0}^{L}h(x,t)e^{2\pi ix/\lambda}dx$ (Fig.
2(d)). The mean wavelength of a filament for $c=1$ is $\lambda_{o}=0.5L$ with
$L=20$. We confirm that the wave speed along the filament body calculated by
using $\lambda_{o}=10$ and $T_{o}=10$ is identical to the input wave speed
($c=1$), i.e., $v_{c}=\lambda_{o}/T_{o}=c=1$.

Figure 2: Beat dynamics of a single filament at $\theta=0^{\circ}$ (see SM
Movie 1). (a) Filament configurations generated at $c=1$ in increasing time
with the interval of $\Delta t/\tau_{s}=2$ (pale to dark). (b) Kymograph
($h(x,t)$) depicting the spatiotemporal evolution of filament dynamics. (c),
(d) Power spectra of $h(x,t)$ in time and space: (c) $S(x,T)$, and (d)
$F(\lambda,t)$.  Figure 3: Dynamics of two interacting filaments with
increasing mismatch. (a) Configurations of the filaments (#1 (blue) and #2
(red)) and their power spectra, $S^{\\#1}(x,T)$ and $S^{\\#2}(x,T)$,
calculated for $\Delta l/L=0.15$, $0.3$, and $0.6$, and $t\gg 1$. (b) Mean
undulation period of each filament, $\langle
T\rangle=\frac{1}{L}\int_{0}^{L}T(x)dx$, as a function of $\Delta l/L$. See
also SM Movie 2.

Dynamics of interacting filaments.  When multiple filaments interact with each
other, additional force terms, $\mathcal{F}_{\rm HI}$, $\mathcal{F}_{\rm SI}$,
and $\mathcal{F}_{\rm WI}$, come into play. Based on the numerical solutions
of Eq.1 at steady states, we posit that the presence of SI and WI terms
renormalizes the coefficient for HIs is modified from $\epsilon_{h}$ to
$\epsilon^{\rm eff}$. When the array of filaments are tilted at a finite angle
$\theta$ towards the bounding wall (Fig. 1), the filaments do not fully
overlap, engendering a mismatch, whose size $\Delta l$ is related to $\theta$
as $\Delta l=d\tan\theta$ (Fig. 1). Using the resistive force theory that
enables one to replace the hydrodynamic coupling with the velocity of
transversal displacement (Eq. 2), the force-balance equation for each filament
in the array made of $N$ filaments can be written as

$\displaystyle\partial_{t}h^{1}(x,t)$
$\displaystyle=\mathcal{N}(h^{1})+\epsilon^{\rm eff}\partial_{t}h^{2}(x-\Delta
l,t)$ $\displaystyle\partial_{t}h^{i}(x,t)$ $\displaystyle=\mathcal{N}(h^{i})$
$\displaystyle+\epsilon^{\rm eff}\left[\partial_{t}h^{i-1}(x+\Delta
l,t)+\partial_{t}h^{i+1}(x-\Delta l,t)\right]$
$\displaystyle\partial_{t}h^{N}(x,t)$
$\displaystyle=\mathcal{N}(h^{N})+\epsilon^{\rm
eff}\partial_{t}h^{N-1}(x+\Delta l,t).$ (4)

The segment of $(i\mp 1)$th filament located at $x\pm\Delta l$ interacts with
the segment of $i$th filament at $x$. When the filaments beat collectively
with $\Delta l>0$, the transversal displacements of the interacting segments
satisfy the inequalities, $|h^{i-1}(x+\Delta
l,t)|>|h^{i}(x,t)|>|h^{i+1}(x-\Delta l,t)|$. Thus, we assume that the
transversal displacements of the neighboring filaments satisfy the relation
$h^{i\pm 1}(x\mp\Delta l,t)\approx b_{\pm}(x,\Delta l)h^{i}(x,t)$, where
$b_{\pm}(x,\Delta l)$ are the factors modulating the transversal displacement
at $x$.

In full synchrony without tilt ($\Delta l=0$), the filaments beat
collectively, satisfying $h^{1}\approx\cdots\approx h^{N}\approx h_{\rm
sync}$, and do not alter the waveform of adjacent filament
($b_{\pm}(x,0)\approx 1$). Furthermore, HIs become the most dominant among the
interactions, $\epsilon^{\rm eff}\approx\epsilon_{h}$. By summing the
equations ($i=1,\ldots,N$) in Eq. 4 and dividing it by $N$ we obtain an
evolution equation for $h_{\rm sync}(x,t)$,
$\left[1-2\left(1-1/N\right)\epsilon^{\rm eff}\right]\partial_{t}h_{\rm
sync}(x,t)=\mathcal{N}(h_{\rm sync})$. Thus, one can consider either that
$h_{\rm sync}(x,t)$ is mapped to $h_{o}(x,t)$ by rescaling the time as $h_{\rm
sync}\left(x,t\right)=h_{o}\left(x,\frac{t}{1-2(1-1/N)\epsilon^{\rm
eff}}\right)$, or that reduction of the drag coefficient occurs from $1$ to
$\left[1-2\left(1-1/N\right)\epsilon^{\rm eff}\right]$. This argument can be
generalized to the case with a finite tilt ($\Delta l\neq 0$), enabling us to
obtain the mean period of undulation for an array of $N$ filaments,

$\displaystyle\langle T_{\theta}\rangle(N)=\left[1-2\bar{b}(\Delta
l)\left(1-\frac{1}{N}\right)\epsilon^{\rm eff}\right]T_{o}$ (5)

with $\bar{b}(\Delta l)\equiv\frac{1}{2L}\int_{0}^{L}[b_{-}(x,\Delta
l)+b_{+}(x,\Delta l)]dx$.

For a pair of coupled filaments, we study the effect of misalignment between
the adjacent filaments on the synchronization. As the extent of mismatch
($\Delta l$) increases, three different dynamical behaviors emerge: (i) full
synchronization; (ii) partial synchronization; (iii) decoupled dynamics. Fig.
3(a) displays typical configurations of the filaments and their power spectra
$S(x,T)$ for three cases of $\Delta l$.

(i) Sync regime ($0\leq\Delta l/L\leq 0.15$): The beat dynamics of filaments
are fully synchronized in the region of two filaments being overlapped, and
there is a minor distinction in the pattern of power spectra between
$S^{\\#1}(x,T)$ and $S^{\\#2}(x,T)$. The corresponding mean period of
undulation increases from $\langle T\rangle=(1-\epsilon_{h})T_{o}\approx 5.0$
to $\approx 5.5$ (see Fig. 3(b)). This slight increase in the period is due to
the reduced size of overlap between the two filaments. In fact, for a filament
characterized with a fixed wave speed, both $\lambda_{o}$ and $T_{o}$ tend to
increase with decreasing filament size ($L$) [24].

(ii) Mixed regime ($0.15<\Delta l/L<0.45$): The tail part of the filament #1,
whose transversal fluctuations is greater than that of head part, interacts
with the head part of the filament #2, which engenders a mismatch in the
velocity profiles, $\partial_{t}h^{1}(x^{1},t)$ and
$\partial_{t}h^{2}(x^{2},t)$. Since the strength of HIs is proportional to the
amplitude of transversal displacement for a fixed oscillatory period, i.e.,
$|f_{\rm HI}^{j\rightarrow i}(x,t)|\propto|h^{j}(x,t)|$ (Eq. 2), the mismatch
leads to the left-right symmetry breaking of hydrodynamic forces, $f_{\rm
HI}^{2\rightarrow 1}\neq f_{\rm HI}^{1\rightarrow 2}$. In this regime, the
synchronization and decoupling of filament dynamics coexist, and the power
spectra $S^{\\#1}(x,T)$ and $S^{\\#2}(x,T)$, exhibiting two main peaks at
$T=(1-\epsilon_{h})T_{o}\approx 5$ and $T=T_{o}\approx 10$, clearly differ
from each other, suggestive of the broken LR symmetry in the beat dynamics of
the two.

Figure 4:  MCWs generated from the 1D array of filaments at varying $\theta$
with periodic boundary (see SM Movie 3). (a) Configurations of many filaments
producing MCWs with increasing $\theta$ and the corresponding power spectra,
$S(x,T)$ and $F(\lambda,t)$. (b) Wavelength, period, and wave speed for the
individual filaments in the array with increasing $\theta$. (c) The mean
transversal fluctuations of filaments in the array with increasing $\theta$:
$\langle\delta h^{+}\rangle$, and $\langle\delta h^{-}\rangle$, where
$\langle\delta h^{\pm}\rangle$ is the mean squared positive (or negative)
fluctuations. (left) The large amplitudes of $\langle\delta h^{-}\rangle$ at
$\theta=60^{\circ}$ and $72^{\circ}$ is highlighted by the red arrows. (right)
The ratio of left/right amplitudes, $|\langle\delta h^{-}\rangle/\langle\delta
h^{+}\rangle|$. (d) Time evolution of the mass densities of the filaments in
arrays projected on the $y^{\prime}$-axis. The speed of MCWs can be estimated
from the slope of stripes: $v_{\rm M}^{\ast,45^{\circ}}\approx 40$, $v_{\rm
M}^{\ast,60^{\circ}}\approx 10$ and $v_{\rm M}^{\ast,72^{\circ}}\approx 4$.

(iii) Decoupled regime ($0.45<\Delta l/L\leq 1$): The peak of $S(x,T)$ is
uniquely formed at $T\approx T_{o}$. Although 40 $\%$ of the filament body is
subject to HIs for $\Delta l/L=0.6$, the dynamics of two filaments are
effectively decoupled, and the beat dynamics of each filament is restored to
$h_{o}(x,t)$. An increasing mismatch ($\Delta l$) weakens the strength of HIs
(or reduces $\bar{b}(\Delta l)$), and it consequently tends to increase the
mean period of undulation $\langle T\rangle$ from the case with zero mismatch
to the case with full decoupling, i.e., $(1-\epsilon_{h})T_{o}\rightarrow
T_{o}$. As $\Delta l$ increases from the regime (i) to (iii), the relative
weight of the two main peaks in Fourier space ($S(x,T)$) shifts from
$T\approx(1-\epsilon_{h})T_{o}$ to $T\approx T_{o}$ and undergoes a first
order type discrete transition (Fig. 3(a)).

Figure 5:  (a) Work production and dissipation, and (b) thermodynamic
efficiency with varying $\theta$. The inset shows a snapshot of nematode array
forming $\theta\approx 70^{\circ}$ from the Peshkov _et al._ ’s experiment
[22]. (c) The mean velocity of fluid flow calculated for an array of filaments
($N=400$) at $\theta=30^{\circ}$.

We now analyze behaviors of 1D arrays of many filaments ($N=10$) at a fixed
separation ($d=10a$) by imposing a periodic boundary condition to eliminate
the boundary effect (see SM Text and Fig. S2 for the case with varying $N$
without the periodic boundary).

First, for $\theta=0^{\circ}$, The filaments beat in perfect in-phase
synchrony at steady states, yielding the waveform, $h_{\rm sync}(x,t)$, whose
wavelength and transversal fluctuation are identical to those of a single
filament in isolation ($\langle\lambda\rangle=\lambda_{o}=0.5L$ and
$\langle\delta h\rangle=\delta h_{o}=0.17L$), yet the time is rescaled as
$t\rightarrow t/(1-2\epsilon_{h})$. the oscillatory period is reduced to
$\langle T\rangle=2.16\times 10^{-3}T_{o}$, which substantially accelerates
the wave speed to $v_{c}=462c$. This is in excellent agreement with a
theoretical estimate, $\langle
T_{0^{\circ}}\rangle=(1-2\epsilon_{h})T_{o}\approx 2.20\times 10^{-3}T_{o}$
and $v_{c}=c/(1-2\epsilon_{h})\approx 460c$, where we have used Eq. 5 with
$\bar{b}(\Delta l=0)=1$, $\epsilon_{h}\approx 0.4989$ (Eq. 3) and
$N\rightarrow\infty$. Individual filaments beat without phase or time lag
between adjacent filaments and retain the LR symmetry.

Next, for $\theta>0^{\circ}$, the waveforms for all the filaments are modified
from $h_{0^{\circ}}(x,t)$ to $h_{\theta}(x,t)$, while the individual filaments
collectively undulate with a constant time lag $\tau_{\theta}$ against the
adjacent filaments, generating MCWs that travel across the array of filaments
($y^{\prime}$-axis), which differ from the body waves traveling along the
centerline ($x$-axis) of each filament at speed $v_{c}$. $\langle
T_{\theta}\rangle$, calculated by analyzing $S(x,T)$, displays a clear
increase with $\theta$ ($S(x,T)$ in Figs. 4(a) and 4(b)). Although the
wavelength ($\langle\lambda_{\theta}\rangle$) changes non-monotonically,
maximized at $\theta=60^{\circ}$, the speed of body wave ($v_{c}$) decreases
with $\theta$ (Fig. 4(b)).

In light of the two-filament dynamics (Fig. 3(b)), the cases with
$\theta=60^{\circ}$ and $72^{\circ}$ (Fig. 4(a)) belong to the mixed regime,
and the deviation of the wave characteristics from those of a single filament
is evident (Fig. 4(b)). Due to the increasing mismatch between the two
filaments with $\theta$, the average amplitude $\langle\delta h\rangle$ as
well as $h_{\theta}(x,t)$ changes non-monotonically from $\theta=0^{\circ}$ to
$72^{\circ}$, showing a tendency of recovering $h_{o}(x,t)$ at
$\theta=72^{\circ}$ (Fig. 4(c)). Importantly, for $\theta=60^{\circ}$, the
broken LR symmetry is apparent, such that the transversal displacement of
filaments is significantly skewed towards the left side of the centerline,
i.e., $|\langle\delta h^{-}\rangle|>|\langle\delta h^{+}\rangle|$ (see the red
arrows in Fig. 4(c)). Accordingly, the hydrodynamic coupling becomes stronger
on the left side than on the right, breaking the LR symmetry of HIs.

The collective waveform (MCW), composed of the individual body waves,
propagates along the array with a constant time lag $\tau_{\theta}$ between
adjacent filaments. The $\tau_{\theta}$ is the time required for a position of
a filament to travel from $x-\Delta l$ to $x$ and mimic the movement of the
adjacent filament, whose head positions are separated by $d_{h}=\Delta
l/\sin\theta$. Note that $\tau_{\theta}$ is related to $v_{c}(\theta)$ as
$\tau_{\theta}=\Delta l/v_{c}(\theta)$. For $\theta>0$, the speed of the MCW,
$v_{M}$, relates to $v_{c}$ as follows.

$\displaystyle v_{\rm
M}(\theta)=\frac{d_{h}}{\tau_{\theta}}=\frac{d_{h}}{\Delta
l}v_{c}(\theta)=\frac{v_{c}(\theta)}{\sin{\theta}},$ (6)

For different tilt angles, the speeds of MCWs are: $v_{\rm
M}^{45^{\circ}}=42.4$, $v_{\rm M}^{60^{\circ}}=11.5$, and $v_{\rm
M}^{72^{\circ}}=3.4$. Notably, $v_{\rm M}(\theta)>v_{c}(\theta)$ for all cases
of $\theta$. Alternatively, the speed of the MCW can be determined by
analyzing the spatiotemporal evolution of mass density projected on the
$y^{\prime}$ axis, $\rho(y^{\prime},t)$. As clearly shown in Fig. 4(d),
$\rho(y^{\prime},t)$ is periodic in space and time as
$\rho(y^{\prime},t)\approx\rho(y^{\prime}+\Delta y^{\prime\ast},t)$ and
$\rho(y^{\prime},t)\approx\rho(y^{\prime},t+\Delta t^{\ast})$, which enables
to estimate the speed of MCW by using $v_{\rm M}^{\ast}=\Delta
y^{\prime\ast}/\Delta t^{\ast}$. The estimated speeds, $v_{\rm
M}^{\ast,45^{\circ}}\approx 40$, $v_{\rm M}^{\ast,60^{\circ}}\approx 10$ and
$v_{\rm M}^{\ast,72^{\circ}}\approx 4$, agree well with the values of $v_{\rm
M}$ obtained using Eq. 6.

Viscous dissipation, work production, and thermodynamic efficiency.  We
investigate the thermodynamics associated with the generation of MCWs and the
fluid flow. First, given that the swimming velocity along the $x$\- and
$y$-directions resulting from the propulsion force $\bar{F}_{prop,x}$ of a
filament are $U_{x}=\bar{F}_{prop,x}/(\xi_{\parallel}L)$ and $U_{y}=0$ [31,
32], the velocities along the $x^{\prime}$\- and $y^{\prime}$-directions can
be expressed as ${U}_{x^{\prime}}={U}_{x}\cos\theta$ and
${U}_{y^{\prime}}={U}_{x}\sin\theta$. Then, the mean work production rate per
filament averaged over a period, which enhances the fluid flow along the
$y^{\prime}$-direction, is
$\langle\dot{W}\rangle=\bar{F}_{drag,y^{\prime}}\bar{U}_{y^{\prime}}=(\bar{F}_{prop,x})^{2}\sin^{2}\theta/(\xi_{\parallel}L)$.
Meanwhile, the mean viscous dissipation per filament is evaluated as
$\langle\dot{Q}\rangle\approx\frac{1}{N}\sum_{i}^{N}\int_{0}^{L}\xi_{\perp}\left(\frac{\partial
h^{i}}{\partial t}\right)^{2}dx^{i}$. Note that both $\langle\dot{W}\rangle$
and $\langle\dot{Q}\rangle$ are maximized at $\theta=30^{\circ}$ (Fig. 5(a)).
On the other hand, the extent of the mechanical work transmitted to the
surrounding fluid relative to the total energy consumption, quantified by
calculating $\varepsilon(\theta)=\langle\dot{W}\rangle/\langle\dot{Q}\rangle$,
is maximized in the range of $\theta\approx 60^{\circ}-70^{\circ}$ (Fig.
5(b)). Remarkably, we find that the array of _T. aceti_ in the experiment [22]
adopt effectively the same range of tilt angle observed in the experiment
($\theta\approx 70^{\circ}$, see the inset of Fig. 5(b)). This suggests that
the nematode population immersed in the fluid media is spontaneously adjusted
to a configuration that maximizes the work production for a given energy
resource, namely, the thermodynamic efficiency ($\varepsilon$).

The fluid flow generated by the collective motion of filaments is calculated
by using the linear response between the propulsion force and the fluid
velocity associated through a hydrodynamic tensor (Eq. S5, see SM Text for the
details). Fig. 5(c) shows the velocity fields of medium along
$y^{\prime}$-axis, generated by an array of beating filaments at
$\theta=30^{\circ}$, as a function of distance $x^{\prime}$ calculated by
using a viscosity value corresponding to 1 cP. It is noteworthy that the
velocity field decreases with the distance from the $y^{\prime}$-axis as
$u_{y^{\prime}}\sim 1/x^{\prime}$ for $x^{\prime}\gg 1$, reflecting a
$1/r$-decay characteristic of the fundamental solution of the hydrodynamic
tensor (Eq. S5).

Concluding Remarks.  One of the key features of MCWs observed in the
experiment on a chain of _T. aceti_ is the modulation of the waveform of
individual filaments [22]. To explore MCWs, previous studies adopted a phase
reduction approach, representing a filament dynamics by means of a frequency-
locked phase oscillator [33, 34, 35, 23, 36]. Such a strategy may only prove
effective for perfectly aligned filaments whose waveform remains unaltered in
synchrony (see Fig. S2) [37, 38, 39]; however, it overlooks the physical
reality that the wave characteristics of each filament is subject to change
upon interacting with its neighbors and generating MCWs. Our study makes it
clear that the full evolution of waveforms is essential to evaluate the
thermodynamics of a system accurately as well as to decipher the physics that
gives rise to MCWs.

Unlike the asymmetrically beating dynamics of cilium, the dynamics of a single
nematode is symmetric with respect to its centerline. Despite its biological
importance of facilitating the fluid circulation, generation of MCWs from the
nematode population appears unlikely at first sight. Our study clarifies that
a misalignment between adjacent filaments disrupts the synchronization and
generates a constant time lag, which is the essential mechanism for the
formation of MCWs. Importantly, a large tilt ($\theta=60^{\circ}-70^{\circ}$)
introduced to the array of filaments generates the energetically most
efficient MCWs enhancing the circulation of surrounding fluid.

###### Acknowledgements.

This study was supported by the KIAS individual Grants, CG085701 (NJ),
CG076002 (WKK), and CG035003 (CH) at Korea Institute for Advanced Study. We
thank the Center for Advanced Computation in KIAS for providing the computing
resources.

## References

  * Couzin [2007] I. Couzin, Collective minds, Nature 445, 715 (2007).
  * Elgeti  _et al._ [2015] J. Elgeti, R. G. Winkler, and G. Gompper, Physics of microswimmers—single particle motion and collective behavior: a review, Rep. Prog. Phys. 78, 056601 (2015).
  * Cavagna and Giardina [2014] A. Cavagna and I. Giardina, Bird flocks as condensed matter, Annu. Rev. Condens. Matter Phys. 5, 183 (2014).
  * Moussaïd  _et al._ [2011] M. Moussaïd, D. Helbing, and G. Theraulaz, How simple rules determine pedestrian behavior and crowd disasters, Proc. Natl. Acad. Sci. U. S. A. 108, 6884 (2011).
  * Lee  _et al._ [2018] S. Lee, C. Hyeon, and J. Jo, Thermodynamic uncertainty relation of interacting oscillators in synchrony, Phys. Rev. E 98, 032119 (2018).
  * Hong  _et al._ [2020] H. Hong, J. Jo, C. Hyeon, and H. Park, Thermodynamic cost of synchronizing a population of beating cilia, J. Stat. Mech. , 074001 (2020).
  * Taylor [1951] G. I. Taylor, Analysis of the swimming of microscopic organisms, Proc. R. Soc. A: Math. Phys. Eng. Sci. 209, 447 (1951).
  * Taylor [1952] G. I. Taylor, Analysis of the swimming of long and narrow animals, Proc. R. Soc. A: Math. Phys. Eng. Sci. 214, 158 (1952).
  * Gueron  _et al._ [1997] S. Gueron, K. Levit-Gurevich, N. Liron, and J. J. Blum, Cilia internal mechanism and metachronal coordination as the result of hydrodynamical coupling, Proc. Natl. Acad. Sci. 94, 6001 (1997).
  * Gueron and Levit-Gurevich [1999] S. Gueron and K. Levit-Gurevich, Energetic considerations of ciliary beating and the advantage of metachronal coordination, Proc. Natl. Acad. Sci. U. S. A. 96, 12240 (1999).
  * Niedermayer  _et al._ [2008] T. Niedermayer, B. Eckhardt, and P. Lenz, Synchronization, phase locking, and metachronal wave formation in ciliary chains, Chaos 18 (2008).
  * Gu  _et al._ [2020] H. Gu, Q. Boehler, H. Cui, E. Secchi, G. Savorana, C. De Marco, S. Gervasoni, Q. Peyron, T.-Y. Huang, S. Pane, _et al._ , Magnetic cilia carpets with programmable metachronal waves, Nature Commun. 11, 2637 (2020).
  * Wang  _et al._ [2024] T. Wang, T. ul Islam, E. Steur, T. Homan, I. Aggarwal, P. R. Onck, J. M. den Toonder, and Y. Wang, Programmable metachronal motion of closely packed magnetic artificial cilia, Lab on a Chip 24, 1573 (2024).
  * Hancock [1953] G. Hancock, The self-propulsion of microscopic organisms through liquids, Proc. R. Soc. A: Math. Phys. Eng. Sci. 217, 96 (1953).
  * Yang  _et al._ [2008] Y. Yang, J. Elgeti, and G. Gompper, Cooperation of sperm in two dimensions: synchronization, attraction, and aggregation through hydrodynamic interactions, Phys. Rev. E. 78, 061903 (2008).
  * Golestanian  _et al._ [2011] R. Golestanian, J. M. Yeomans, and N. Uchida, Hydrodynamic synchronization at low Reynolds number, Soft Matter 7, 3074 (2011).
  * Kim and Netz [2006] Y. W. Kim and R. R. Netz, Pumping fluids with periodically beating grafted elastic filaments, Phys. Rev. Lett. 96, 158101 (2006).
  * Zheng  _et al._ [2023] E. Zheng, M. Brandenbourger, L. Robinet, P. Schall, E. Lerner, and C. Coulais, Self-oscillation and synchronization transitions in elastoactive structures, Phys. Rev. Lett. 130, 178202 (2023).
  * Yuan  _et al._ [2014] J. Yuan, D. M. Raizen, and H. H. Bau, Gait synchronization in Caenorhabditis elegans, Proc. Natl. Acad. Sci. U. S. A. 111, 6865 (2014).
  * Chelakkot  _et al._ [2021] R. Chelakkot, M. F. Hagan, and A. Gopinath, Synchronized oscillations, traveling waves, and jammed clusters induced by steric interactions in active filament arrays, Soft Matter 17, 1091 (2021).
  * Guirao and Joanny [2007] B. Guirao and J.-F. Joanny, Spontaneous creation of macroscopic flow and metachronal waves in an array of cilia, Biophys. J. 92, 1900 (2007).
  * Peshkov  _et al._ [2022] A. Peshkov, S. McGaffigan, and A. C. Quillen, Synchronized oscillations in swarms of nematode _Turbatrix aceti_ , Soft Matter 18, 1174 (2022).
  * Quillen  _et al._ [2021] A. Quillen, A. Peshkov, E. Wright, and S. McGaffigan, Metachronal waves in concentrations of swimming _Turbatrix aceti_ nematodes and an oscillator chain model for their coordinated motions, Phys. Rev. E. 104, 014412 (2021).
  * Goldstein  _et al._ [2016] R. E. Goldstein, E. Lauga, A. I. Pesci, and M. R. Proctor, Elastohydrodynamic synchronization of adjacent beating flagella, Phys. Rev. Fluids 1, 073201 (2016).
  * Heuser  _et al._ [2009] T. Heuser, M. Raytchev, J. Krell, M. E. Porter, and D. Nicastro, The dynein regulatory complex is the nexin link and a major regulatory node in cilia and flagella, J. Cell Biol. 187, 921 (2009).
  * Jülicher and Prost [1997] F. Jülicher and J. Prost, Spontaneous oscillations of collective molecular motors, Phys. Rev. Lett. 78, 4510 (1997).
  * Gray and Hancock [1955] J. Gray and G. Hancock, The propulsion of sea-urchin spermatozoa, J. Exp. Biol. 32, 802 (1955).
  * Rotne and Prager [1969] J. Rotne and S. Prager, Variational treatment of hydrodynamic interaction in polymers, J. Chem. Phys. 50, 4831 (1969).
  * Yamakawa [1970] H. Yamakawa, Transport properties of polymer chains in dilute solution: hydrodynamic interaction, J. Chem. Phys. 53, 436 (1970).
  * Tornberg and Shelley [2004] A.-K. Tornberg and M. J. Shelley, Simulating the dynamics and interactions of flexible fibers in stokes flows, J. Comp. Phys. 196, 8 (2004).
  * Lauga and Powers [2009] E. Lauga and T. R. Powers, The hydrodynamics of swimming microorganisms, Rep. Prog. Phys. 72, 096601 (2009).
  * Wiggins and Goldstein [1998] C. H. Wiggins and R. E. Goldstein, Flexive and propulsive dynamics of elastica at low Reynolds number, Phys. Rev. Lett. 80, 3879 (1998).
  * Nakao [2016] H. Nakao, Phase reduction approach to synchronisation of nonlinear oscillators, Contemporary Physics 57, 188 (2016).
  * Kawamura and Tsubaki [2018] Y. Kawamura and R. Tsubaki, Phase reduction approach to elastohydrodynamic synchronization of beating flagella, Phys. Rev. E. 97, 022212 (2018).
  * Man and Kanso [2020] Y. Man and E. Kanso, Multisynchrony in active microfilaments, Phys. Rev. Lett. 125, 148101 (2020).
  * Quillen [2023] A. Quillen, Robust formation of metachronal waves in directional chains of phase oscillators, Phys. Rev. E 107, 034401 (2023).
  * Guo  _et al._ [2018] H. Guo, L. Fauci, M. Shelley, and E. Kanso, Bistability in the synchronization of actuated microfilaments, J. Fluid Mech. 836, 304 (2018).
  * Young [2009] Y.-N. Young, Hydrodynamic interactions between two semiflexible inextensible filaments in stokes flow, Phys. Rev. E. 79, 046317 (2009).
  * Vilfan and Jülicher [2006] A. Vilfan and F. Jülicher, Hydrodynamic flow patterns and synchronization of beating cilia, Phys. Rev. Lett. 96, 058102 (2006).
  * Lighthill [1975] S. J. Lighthill,  _Mathematical biofluid dynamics_ (SIAM, 1975).
  * Chakrabarti and Saintillan [2019] B. Chakrabarti and D. Saintillan, Spontaneous oscillations, beating patterns, and hydrodynamics of active microfilaments, Phys. Rev. Fluids 4, 043102 (2019).
  * Fang-Yen  _et al._ [2010] C. Fang-Yen, M. Wyart, J. Xie, R. Kawai, T. Kodger, S. Chen, Q. Wen, and A. D. Samuel, Biomechanical analysis of gait adaptation in the nematode Caenorhabditis elegans, Proc. Natl. Acad. Sci. U. S. A. 107, 20323 (2010).
  * Sznitman  _et al._ [2010] J. Sznitman, P. K. Purohit, P. Krajacic, T. Lamitina, and P. E. Arratia, Material properties of Caenorhabditis elegans swimming at low Reynolds number, Biophys. J. 98, 617 (2010).
  * Xu  _et al._ [2016] G. Xu, K. S. Wilson, R. J. Okamoto, J.-Y. Shao, S. K. Dutcher, and P. V. Bayly, Flexural rigidity and shear stiffness of flagella estimated from induced bends and counterbends, Biophys. J. 110, 2759 (2016).
  * Gilpin  _et al._ [2015] W. Gilpin, S. Uppaluri, and C. P. Brangwynne, Worms under pressure: bulk mechanical properties of C. elegans are independent of the cuticle, Biophys. J. 108, 1887 (2015).
  * Elfring and Lauga [2009] G. J. Elfring and E. Lauga, Hydrodynamic phase locking of swimming microorganisms, Phys. Rev. Lett. 103, 088101 (2009).
  * Leptos  _et al._ [2013] K. C. Leptos, K. Y. Wan, M. Polin, I. Tuval, A. I. Pesci, and R. E. Goldstein, Antiphase synchronization in a flagellar-dominance mutant of chlamydomonas, Phys. Rev. Lett. 111, 158101 (2013).

##  I Supplemental Materials

##  II Supplemental Movies

  * •

SM Movie 1: Beat dynamics of a single filament at $\theta=0^{\circ}$.

  * •

SM Movie 2: Dynamic of two interacting filaments at $\theta=0^{\circ}$ with
varying sizes of mismatch $\Delta l$.

  * •

SM Movie 3: Dynamics of $N=10$ filaments at $\theta=30^{\circ}$, $45^{\circ}$,
$60^{\circ}$, and $72^{\circ}$.

  * •

SM Movie 4: Demonstration of the effect of SI on two-filament dynamics
starting from a V-shape configuration.

  * •

SM Movie 5: Many filament dynamics at $\theta=0^{\circ}$ with varying $N$,
without periodic boundary condition.

  * •

SM Movie 6: Many filament dynamics at $\theta=60^{\circ}$ with varying $N$,
without periodic boundary condition.

##  III Elastohydrodynamic model

We consider $N$ filaments, each with a length of $L$ and a cross-sectional
radius of $a$, tilted with $\theta$ while maintaining a body-to-body
separation distance $d$ within a two-dimensional plane (Fig. 1). The head of
each filament is hinged along the $y^{\prime}$-axis. The head-to-head distance
depends on the tilt angle $\theta$ and is separated by $d_{h}=d/\cos{\theta}$
for the given separation $d$. An aspect ratio of the filament is $L/(2a)=50$
and the scale of $d$ used in the model is in the range of $a\ll d\ll L$. Under
an assumption of small undulation compared to the body length, the curvilinear
coordinate along the filament ($0\leq s\leq L$) is approximated to $s\approx
x$.

At low Reynolds number, the force balance on the $i$th filament should be
established among viscous drag ($f^{i}_{vis}$), bending ($f^{i}_{ben}$),
active forces ($f^{i}_{act}$), hydrodynamic couplings
($\sum\limits_{j=i-1,i+1}f^{j\rightarrow i}_{hyd}$) and steric interactions
($\sum\limits_{j=i-1,i+1}{f^{j\rightarrow i}_{rep}}$) from the neighboring
filaments, such that the sum of all the force components should vanish,
$\sum_{\alpha}f^{i}_{\alpha}=0$.

Each force term can be related to the waveform of filaments $h^{i}(x,t)$. The
force associated with viscous drag is given as

$\displaystyle f^{i}_{vis}=-\zeta_{\perp}\partial_{t}h^{i}$ (S1)

where $\zeta_{\perp}\approx 4\pi\mu/(\ln{L/a})\approx\pi\mu$ is the drag
coefficient for a slender body [40, 31]. The force arising from the bending
penalty can be obtained from the functional derivative of the bending energy
$(\kappa/2)\int_{0}^{L}dx(\partial^{2}_{x}h^{i})^{2}$,

$\displaystyle f^{i}_{ben}=-\kappa\partial_{x}^{4}h^{i}$ (S2)

where $\kappa$ is the coefficient of bending stiffness.

For the active force we adopt the phenomenological model by Goldstein et. al
[24]

$\displaystyle
f^{i}_{act}=-A\partial_{x}h^{i}-B\partial^{2}_{x}h^{i}+C\left(\partial^{2}_{x}h^{i}\right)^{3}$
(S3)

where the first term ($A>0$) generates propagating waves along the centerline
of filament ($x$-axis in Fig. 1). The second and the third term break head-
tail symmetry of the amplitude of transversal undulations, which is magnified
by the second term ($B>0$) and suppressed by the third term ($C>0$), saturates
to a finite value.

The hydrodynamic force from $j$th to $i$th filaments is defined via the
fundamental singular solutions of Stokes equations:

$\displaystyle f^{j\rightarrow
i}_{hyd}(x^{i},t)/\zeta_{\perp}=\textbf{e}_{y}\cdot\int_{0}^{L}dx^{j}\left(\textbf{G}(r)\cdot\textbf{f}^{j}_{vis}(x^{j},t)\right)$
(S4)

where $\textbf{e}_{y}=(0,1)$, $\textbf{G}(r)$ is the sum of Stokeslet and
doublet

$\displaystyle\textbf{G}(r)=\frac{1}{8\pi\mu
r}\left[\left(\textbf{I}+\frac{\textbf{r}\textbf{r}}{r^{2}}\right)+\frac{a^{2}}{2r^{2}}\left(\textbf{I}-3\frac{\textbf{r}\textbf{r}}{r^{2}}\right)\right]$
(S5)

with r denoting the relative position vector from $(x^{j},h^{j})$ to
$(x^{i},h^{i})$ and $r=|\textbf{r}|$. The second term, corresponding to the
doublet, corrects the boundary condition when the velocity is computed on the
surface of the slender body near the singularity distribution. In Eq.S4, the
term associated with viscous force of $j$-th filament,
$\textbf{f}^{j}_{vis}(x^{j},t)$, is also obtained from the resistive force
theory,
$f^{j}_{vis}(x^{j},t)=\textbf{e}_{y}\cdot\textbf{f}^{j}_{vis}(x^{j},t)=\zeta_{\perp}\partial_{t}h^{j}(x^{j},t)$.
It should be emphasized that for a non-zero tilt ($\Delta l\neq 0$), the
nearest point of adjacent chain to $x^{i}$ is $x^{j}=x^{i}\pm\Delta l$ for
$j=i\mp 1$. Taken together, Eq.S4 is simplified to the following form [24]:

$\displaystyle f^{j\rightarrow
i}_{hyd}(x^{i},t)/\zeta_{\perp}=\epsilon_{h}\partial_{t}h^{j}(x^{j},t).$ (S6)

For the steric repulsion between the adjacent filaments, we use the Weeks-
Chandler-Anderson (WCA) potential

$\displaystyle
U_{R}(r)=\begin{cases}4\epsilon_{r}\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]+\epsilon_{r}&\text{for
}r\leq 2^{1/6}\sigma\\\ 0&\text{otherwise}\end{cases}$ (S7)

where $\sigma$ is a length scale at which $U_{R}(\sigma)=0$ and $\epsilon_{r}$
is a scale of interaction energy. The repulsive force due to the steric
interaction is obtained from

$\displaystyle f^{j\rightarrow
i}_{rep}(x^{i},t)={\textbf{e}_{y}}\cdot\int_{0}^{L}dx^{j}\left(-\frac{\partial
U_{R}(r)}{\partial r}\right)\frac{\textbf{r}}{r}.$ (S8)

For a single active elastic filament demonstrating a beating dynamics, two
dimensionless parameters can characterize its dynamics: (i) activity number,
the active force relative to the bending force; (ii) sperm number, the beat
frequency relative to the rate of the bending wave relaxation in a viscous
fluid [32, 41]. A waveform generated from a collective dynamics of filaments
in a dense suspension is also characterized by the two parameters of a single
filament [15]. In our system, the activity and sperm numbers are given as
$\mu_{a}\equiv\frac{AL}{\kappa/(2aL)}=160$,
$Sp\equiv\left(L^{4}\zeta_{\perp}/\kappa\tau_{a}\right)^{1/4}=\left(L^{4}\tau_{b}/l_{s}^{4}\tau_{a}\right)^{1/4}=L/l_{s}=20$.
Here, we assume that the characteristic time scale of active motion,
$\tau_{a}$, is the same as that of bending relaxation time scale, $\tau_{b}$,
and use it as a characteristic time scale of the system,
$\tau_{s}=\tau_{b}=\tau_{a}$.

The following expression that considers the WCA potential is numerically
computed:

$\displaystyle f^{j\rightarrow
i}_{rep}(x^{i},t)/\zeta_{\perp}=\begin{cases}(4\epsilon_{r}/\zeta_{\perp})\int_{0}^{L}\left(\frac{12\sigma^{12}}{r^{14}}-\frac{6\sigma^{6}}{r^{8}}\right)[h^{i}-h^{j}+(i-j)d+2a]dx^{j}&\text{if
}r\leq 2^{1/6}\sigma,\\\ 0&\text{otherwise}.\end{cases}$ (S9)

for the steric repulsion between the neighboring filaments, and

$\displaystyle f^{{\rm w}\rightarrow i}_{rep}(x^{i},t)/\zeta_{\perp}$
$\displaystyle=\zeta_{\perp}^{-1}{\textbf{e}_{y}}\cdot\int_{0}^{L_{\rm
w}}dx^{\rm w}\left(-\frac{\partial U_{R}(r_{\rm w})}{\partial r_{\rm
w}}\right)\frac{\textbf{r}_{\rm w}}{r_{\rm w}}$
$\displaystyle=\begin{cases}(4\epsilon_{r}/\zeta_{\perp})\int_{0}^{L_{\rm
w}}\left(\frac{12\sigma^{12}}{r_{\rm w}^{14}}-\frac{6\sigma^{6}}{r_{\rm
w}^{8}}\right)[y^{i}-y^{\rm w}]dx^{\rm w}&\text{if }r_{\rm w}\leq
2^{1/6}\sigma,\\\ 0&\text{otherwise}.\end{cases}$ (S10)

for the wall repulsion, where $\mathbb{r}_{\rm w}=(x^{i},y^{i})-(x^{\rm
w},y^{\rm w})$, $||\mathbb{r}_{\rm w}||=r_{\rm w}=\sqrt{(x^{i}-x^{\rm
w})^{2}+(y^{i}-y^{\rm w})^{2}}$, $\mathbb{e}_{y}\cdot\mathbb{r}_{\rm
w}=y^{i}-y^{\rm w}$, and $\left(x^{i}-x^{\rm w},y^{i}-y^{\rm
w}\right)=\left(x^{i}+(i-1)\Delta l-x^{\rm w}\sin\theta+d_{\rm
w}\cos\theta,h^{i}(x^{i},t)+(i-1)d-x^{\rm w}\cos\theta-d_{\rm
w}\sin\theta\right)$. The boundary conditions were imposed to create a hinged
head ($h^{i}(0,t)=h^{i}_{xx}(0,t)=0$) and a free tail ($h^{i}_{xx}(L,t)=0$ and
$h^{i}_{xxx}(L,t)=0$). As initial conditions, we used $N$ straight filaments
in a parallel alignment, and studied the dynamics for various $N$. The
governing equation (Eq.S11) was integrated using the finite difference and the
semi-implicit scheme [30].

Using the relations, $(\tau_{s}\kappa/\zeta_{\perp}l_{s}^{4})=1$,
$(\tau_{s}A/\zeta_{\perp}l_{s})=c$, $(\tau_{s}B/\zeta_{\perp}l_{s}^{2})=2$,
and $(\tau_{s}C/\zeta_{\perp}l_{s}^{4})=1$, which determine the characteristic
length and time scales of the system as $l_{s}=\left(2\kappa/B\right)^{1/2}$
and $\tau_{s}=(4\kappa\zeta_{\perp}/B^{2})$, we non-dimensionalize the force
balance equation for each filament as

$\displaystyle\frac{\partial\hat{h}^{i}}{\partial\hat{t}}=\left[-c\frac{\partial\hat{h}^{i}}{\partial\hat{x}}-2\frac{\partial^{2}\hat{h}^{i}}{\partial\hat{x}^{2}}+\left(\frac{\partial^{2}\hat{h}^{i}}{\partial\hat{x}^{2}}\right)^{3}\right]-\frac{\partial^{4}\hat{h}^{i}}{\partial\hat{x}^{4}}+\frac{\tau_{s}}{l_{s}}\left[\sum_{j=i-1,i+1}\frac{f^{j\rightarrow
i}_{hyd}}{\zeta_{\perp}}+\sum_{j=i-1,i+1}\frac{f^{j\rightarrow
i}_{rep}}{\zeta_{\perp}}\right]$ (S11)

where $\hat{h}=h/l_{s}$, $\hat{x}=x/l_{s}$, $\hat{t}=t/\tau_{s}$, and
$c=A\kappa^{1/2}(2/B)^{3/2}$ is a dimensionless speed of wave. For the
parameter values of $\mu=1$ $\text{cP}=10^{-3}$ Pa$\cdot$sec, $B=0.8$ $\mu$N
[22], and $C=\kappa=10^{-3}$ $\mu$N$\cdot$mm2 (corresponding to the Young’s
modulus of $E=125$ kPa) [42, 43, 44, 45], we obtain $l_{s}=0.05$ mm and
$\tau_{s}=20$ ms. When $c=1$, $A=8$ $\mu$N/mm. We also set the values of
$\epsilon_{r}=A$ $(=5\times 10^{9}k_{B}T/l_{s}^{2})$ and $\sigma=2a$. In Eq.
1, the hat symbol for non-dimensional length ($\hat{h}$ and $\hat{x}$) and
time ($\hat{t}$) is dropped from the expression for simplicity.

The stability analysis on the linearized elastohydrodynamic equation,
$\partial_{t}\hat{h}=-c\partial_{x}\hat{h}-2\partial_{x}^{2}\hat{h}-\partial_{x}^{4}\hat{h}$
with $\hat{h}\sim e^{i(\omega t-kx)}$ yields a dispersion relation,

$\displaystyle i\omega=-ick+2k^{2}-k^{4}.$ (S12)

Thus, for the dynamics of the filament to be oscillatory, the condition of
$k^{2}(k^{2}-2)<0$ imposes the wavelength ($k=2\pi/\lambda<\sqrt{2}$) of a
single filament always greater than $\sqrt{2}\pi$,

$\displaystyle\lambda>\sqrt{2}\pi.$ (S13)

##  IV Effect of steric interactions on the synchronization of two
hydrodynamically coupled filaments

The effect of SI on the two hydrodynamically coupled filaments are
investigated for a system without tilt ($\theta=0^{\circ}$), separated by
$d=10a$. To quantify the process of synchronization between two active
filaments ($i,j=1,2$ and $i\neq j$), we calculate a correlation function,

$\displaystyle
C(t)=\frac{1}{L}\int_{0}^{L}\frac{1+\left(\partial_{x}h^{1}\right)\left(\partial_{x}h^{2}\right)}{\sqrt{1+\left(\partial_{x}h^{1}\right)^{2}}\sqrt{1+\left(\partial_{x}h^{2}\right)^{2}}}dx.$
(S14)

Here, $C(t)=1$ for perfect synchrony and $C(t)=0$ when their dynamics are
decoupled.

Figure S1: Dynamics of a pair of hydrodynamically-coupled filaments from
V-shape initial configuration with and without SI (see SM Movie 4). (a)
Evolution of correlation functions ($C(t)$) for the two filaments coupled
either only through HIs (black) or through HIs and SIs (red). Depicted in the
insets are the steady-state configurations in (top) IP and (bottom) AP
synchrony. (b) The kymographs of two filaments ($h^{1}(x,t)$ and $h^{2}(x,t)$)
coupled through HIs and SIs. (c) Phase diagrams of in-phase (IP) and anti-
phase (AP) synchronization between two active filaments from the V-shaped
initial condition for various values of $c$ and $d$. The steady state at given
$c$ and $d$ is determined by using $C(t)$ (Eq. S14), such that $C(t)>0.6$ for
the IP state, and $C(t)\leq 0.6$ for the AP state. The solid line in magenta
plots the $c$-dependent fluctuations of a single filament ($\delta h(c)$, Eq.
(S15)).

In the absence of SIs, depending on the initial condition, the two
hydrodynamically coupled filaments demonstrate either in-phase (IP) or anti-
phase (AP) synchronization after a transient time [35]. Adding the short-range
SI little affects the hydrodynamically coupled filaments already in IP
synchrony; however, the same SI alters the two-filament configuration in AP
synchrony into the one in IP synchrony, particularly, when the amplitude of
undulation is greater than the inter-filament separation ($\langle\delta
h\rangle\gtrsim d$).

Two-filament dynamics starting from a V-shape configuration demonstrates the
effect of SI dramatically (see SM Movie 4). The filaments, coupled only
through HI, converge to a configuration in AP synchrony (the black line in
Fig. S1(a)). However, an introduction of SI destabilizes the AP configuration,
nudging the system towards the one in IP synchrony (the red graph in Fig. S1A)
[46, 47, 35]. The kymographs of $h^{1}(x,t)$ and $h^{2}(x,t)$ depict the
transition from AP ($h^{2}=-h^{1}$) to IP ($h^{2}=h^{1}$) at $t\gtrsim 40$,
accompanied with the decrease of beat period from $(1+\epsilon_{h})T_{o}$ to
$(1-\epsilon_{h})T_{o}$ (Fig. S1(b)).

The phase diagram of the steady state configurations of two filaments (Fig.
S1(c)) shows the stability of AP and IP modes with respect to $c$ and $d$ [37,
24], along with their phase boundary which can be described using the
$c$-dependent amplitude of transversal displacement, $\delta h(c)$, obtained
by assuming that a single traveling wave for the tail part is given as
$h(x,t)=\delta h\cos(kx-\omega t)$. Since $h(x,t)$ should satisfy Eq. 1, we
obtain the amplitude of transversal displacement [24]:

$\displaystyle\langle\delta
h\rangle(c)=\frac{2\sqrt{2-\left[\left(c/8\right)^{2/3}+1/3\right]}}{\sqrt{3}\left[\left(c/8\right)^{2/3}+1/3\right]}.$
(S15)

IP state is stable for $c$ satisfying $\langle\delta h\rangle(c)>d$ [46, 35].

Figure S2: Dynamics among filaments at $\theta=0^{\circ}$ and $c=1$ with
varying $N$. (a) Configuration of $N=20$ filaments in in-phase synchrony. (b)
$\langle T\rangle$, $\langle\lambda\rangle/L$, $\langle\delta h\rangle/L$
calculated for $i=1,\ldots,N$ filament. For the demonstration, the filament in
the middle is centered. (c) Mean period of filament dynamics in synchrony,
$\langle T\rangle$ as a function of $N$. The red line is the fit using Eq. 5.

##  V Many filaments with non-periodic boundary conditions

We consider a case where a group of active filaments with non-periodic
boundary conditions to study the impact of the number of filaments $N$ on the
collective dynamics at two tilt angles, $\theta=0^{\circ}$ and $60^{\circ}$.

For a set of $N$-filaments with $\theta=0^{\circ}$ (see Fig. S2(a) for a
configuration of $N=20$ filaments), the filaments beat in perfect synchrony,
so that the mean periods ($\langle T\rangle$) and wavelengths
($\langle\lambda\rangle$) are identical for different filaments except for
those at the ends in the case of $N=15$ and $20$, whereas the amplitude of
transversal displacement ($\langle\delta h\rangle$) is the largest for the
filaments in the middle (Fig. S2(b)). Notably, the wavelength of the waveform
is effectively identical as $\langle\lambda\rangle/L=0.5$, irrespective of
$N$; however, the mean period ($\langle T\rangle$) decreases with $N$ (Fig.
S2(c)). Thus, the wave propagation speed ($v_{c}=\langle\lambda\rangle/\langle
T\rangle$) along the centerline ($x$) increases with $N$ (see SM Movie 4).

Figure S3: Dynamics among filaments at $\theta=60^{\circ}$ and $c=1$ with
varying $N$. (a) Configuration of N = 20 filaments in synchrony. (b) Same as
Fig. S2(b).

In fact, Eq. 5 with the fit parameters, $T_{o}\approx 9.98$, $\bar{b}=1$, and
$\epsilon^{\rm eff}\approx 0.52$, well describes the periods of collective
motion of filaments ($\theta=0^{\circ}$) for varying $N$. Note, however, that
in the collective wave generated from the $N$ filaments, there is neither the
phase lag between the adjacent filaments nor the significant change in the
waveform. The waveform of each filament bears the LR symmetry with no wave
packet on average being transmitted to left or right along the
$y^{\prime}$-axis (see SM Movie 5).

In contrast, for $\theta=60^{\circ}$, it is clearly seen that MCWs travel from
left to right along the $y^{\prime}$-axis, although the beat dynamics of
filaments are only partially synchronized (see SM Movie 6).

For $N=20$, the wave characteristics of the filaments for the two different
tilt angles averaged over all the filaments excluding those contributing at
the boundary are given as follows: (i) $\langle\lambda\rangle=0.5$ $L$ and
$\langle T\rangle=0.13$ for $\theta=0^{\circ}$; (ii)
$\langle\lambda\rangle=0.71$ $L$ and $\langle T\rangle=2.07$ for
$\theta=60^{\circ}$. Thus, the wave propagation speeds along the centerline
are $v_{c}^{0^{\circ}}=76.9$ and $v_{c}^{60^{\circ}}=6.76$. Note that
$v_{c}^{0^{\circ}}$ is greater than $v_{c}^{60^{\circ}}$ by 10 folds. The
transversal fluctuations, $\langle\delta h\rangle$, are also greater as well
for $\theta=60^{\circ}$. It is of note that more significant reduction is made
in $\langle T\rangle$ with increasing $N$ at $\theta=0^{\circ}$ than at
$\theta=60^{\circ}$, which is primarily due to the enhancement of HI in
multiple filaments in full synchrony at $\theta=0^{\circ}$ in comparison with
those in partial synchrony at $\theta=60^{\circ}$.

##  VI Propulsion Force

With the resistive force theory [31], we compute the anisotropic drag forces
on the $i$th filament in the perpendicular and parallel directions with
respect to the centerline of the filament body as follows.

$\displaystyle f_{\perp}(x^{i},t)=-\xi_{\perp}\cos\psi\frac{\partial
h^{i}}{\partial t}$ $\displaystyle
f_{\parallel}(x^{i},t)=-\xi_{\parallel}\sin\psi\frac{\partial h^{i}}{\partial
t}$ (S16)

where $\psi$ is the angle formed between the filament orientation at $x^{i}$
and the centerline ($x$-axis) and $\xi_{\perp}=2\xi_{\parallel}=4\pi\mu/(\ln
L/a)$. Then, the drag forces along the $x$\- and $y$-axis are calculated as

$\displaystyle f_{x}$ $\displaystyle=-f_{\perp}\sin\psi+f_{\parallel}\cos\psi$
$\displaystyle=(\xi_{\perp}-\xi_{\parallel})\cos\psi\sin\psi\frac{\partial
h^{i}}{\partial t}$ $\displaystyle\approx\frac{1}{2}\xi_{\perp}\frac{\partial
h^{i}}{\partial x}\frac{\partial h^{i}}{\partial t}$ $\displaystyle f_{y}$
$\displaystyle=f_{\perp}\cos\psi+f_{\parallel}\sin\psi$
$\displaystyle=-(\xi_{\perp}\cos^{2}\psi+\xi_{\parallel}\sin^{2}\psi)\frac{\partial
h^{i}}{\partial t}$ $\displaystyle\approx-\xi_{\perp}\frac{\partial
h^{i}}{\partial t}$ (S17)

where we have used the condition of small undulation in comparison with the
body length ($\sin\psi\approx\partial_{x}h^{i}$, $\cos\psi\approx 1$, and
$\sin^{2}\psi\approx 0$) to obtain the approximate expressions for the drag
forces.

Thus, the total drag forces generated by the $i$th filament along the $x$ and
$y$-axes are

$\displaystyle F^{i}_{drag,x}$
$\displaystyle=\frac{1}{2}\xi_{\perp}\int_{0}^{L}\frac{\partial
h^{i}}{\partial x}\frac{\partial h^{i}}{\partial t}dx^{i}$ $\displaystyle
F^{i}_{drag,y}$ $\displaystyle=-\xi_{\perp}\int_{0}^{L}\frac{\partial
h^{i}}{\partial t}dx^{i}$ (S18)

From the force balance, the propulsion force exerted to the fluid due to the
beat motion of a filament are $F^{i}_{prop,x}=-F^{i}_{drag,x}$ and
$F^{i}_{prop,y}=-F^{i}_{drag,y}$. Thus, the total propulsion forces are
obtained by summing over the contributions from $N$ filaments,
$F_{prop,x}=\sum_{i=1}^{N}F^{i}_{prop,x}$ and
$F_{prop,y}=\sum_{i=1}^{N}F^{i}_{prop,y}$, and the tilt angle ($\theta$)
dependent total propulsion forces along the $x^{\prime}$\- and
$y^{\prime}$-axes are

$\displaystyle F_{prop,x^{\prime}}=F_{prop,x}\cos\theta-F_{prop,y}\sin\theta$
$\displaystyle F_{prop,y^{\prime}}=F_{prop,x}\sin\theta+F_{prop,y}\cos\theta$
(S19)

The mean propulsion forces over the period of beat motion in the
$x^{\prime}$\- and $y^{\prime}$-directions are

$\displaystyle\bar{F}_{prop,x^{\prime}}$
$\displaystyle=\bar{F}_{prop,x}\cos{\theta}$
$\displaystyle\bar{F}_{prop,y^{\prime}}$
$\displaystyle=\bar{F}_{prop,x}\sin{\theta}$

Next, the velocity of the fluid flow at position $(x,y)$ at time $t$ generated
by the undulatory motion of the $i$th filament is

$\displaystyle\mathbb{u}^{i}(x,y,t)=-\int_{0}^{L}\mathbb{G}(r)\cdot\mathbb{f}(x^{i},t)dx^{i}$
(S21)

where $\mathbb{f}=(f_{x},f_{y})$ defined in Eq. (S17), $\mathbb{G}(r)$ is in
Eq. (S5), $\mathbb{r}=(x-x^{i},y-h^{i})$, and $r=|\mathbb{r}|$. The flow
velocity along the $x^{\prime}$ and $y^{\prime}$-axes are

$\displaystyle u^{i}_{x^{\prime}}$
$\displaystyle=u^{i}_{x}\cos\theta-u^{i}_{y}\sin\theta$ $\displaystyle
u^{i}_{y^{\prime}}$ $\displaystyle=u^{i}_{x}\sin\theta+u^{i}_{y}\cos\theta.$
(S22)

Therefore, the fluid flow generated by the collective motion of filaments is
$u_{x^{\prime}}=\sum_{i}^{N}u^{i}_{x^{\prime}}$ and
$u_{y^{\prime}}=\sum_{i}^{N}u^{i}_{y^{\prime}}$. Fig. 5 shows the mean values
of $u_{y^{\prime}}$ at location $x^{\prime}$ from the $y^{\prime}$-axis.

##  VII Viscous dissipation and work production

First, the mean viscous dissipation per filament (Fig. 5A) is evaluated as:

$\displaystyle\langle\dot{Q}\rangle$
$\displaystyle=-\frac{1}{N}\sum_{i}^{N}\int_{0}^{L}f_{y}\frac{\partial
h^{i}}{\partial t}dx^{i}$
$\displaystyle\approx\frac{1}{N}\sum_{i}^{N}\int_{0}^{L}\xi_{\perp}\left(\frac{\partial
h^{i}}{\partial t}\right)^{2}dx^{i}.$ (S23)

Next, since the swimming velocity along the $x$\- and $y$-directions [31, 32],
resulting from the propulsion force $\bar{F}_{prop,x}$ of a filament are
$U_{x}=\bar{F}_{prop,x}/\xi_{\parallel}L$ and $U_{y}=0$, the velocities along
the $x^{\prime}$\- and $y^{\prime}$-directions can be expressed as
${U}_{x^{\prime}}={U}_{x}\cos\theta$ and ${U}_{y^{\prime}}={U}_{x}\sin\theta$.
Consequently, the rate of work generated and transmitted to the fluid, which
enhances the fluid flow along the $y^{\prime}$-direction, is:

$\displaystyle\langle\dot{W}\rangle=\bar{F}_{drag,y^{\prime}}\bar{U}_{y^{\prime}}=\frac{(\bar{F}_{drag,x})^{2}\sin^{2}\theta}{\xi_{\parallel}L}$
(S24)