Emergence of metachronal waves in a chain of symmetrically beating filaments

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):

$$\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

$$\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

$$\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

$$\partial_{t}h^{1}(x,t) =\mathcal{N}(h^{1})+\epsilon^{\rm eff}\partial_{t}h^{2}(x-\Delta l,t)$$

$$\partial_{t}h^{i}(x,t) =\mathcal{N}(h^{i})$$

$$+\epsilon^{\rm eff}\left[\partial_{t}h^{i-1}(x+\Delta l,t)+\partial_{t}h^{i+1}(x-\Delta l,t)\right]$$

$$\partial_{t}h^{N}(x,t) =\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,

$$\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.

$$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.

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

$$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}$,

$$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]

$$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:

$$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

$$\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]:

$$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

$$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

$$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:

$$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

$$f^{{\rm w}\rightarrow i}_{rep}(x^{i},t)/\zeta_{\perp} =\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}}$$

$$=\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

$$\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$mm² (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,

$$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$,

$$\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,

$$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]:

$$\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.

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

$$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

$$f_{x} =-f_{\perp}\sin\psi+f_{\parallel}\cos\psi$$

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

$$\approx\frac{1}{2}\xi_{\perp}\frac{\partial h^{i}}{\partial x}\frac{\partial h^{i}}{\partial t}$$

$$f_{y} =f_{\perp}\cos\psi+f_{\parallel}\sin\psi$$

$$=-(\xi_{\perp}\cos^{2}\psi+\xi_{\parallel}\sin^{2}\psi)\frac{\partial h^{i}}{\partial t}$$

$$\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

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

$$F^{i}_{drag,y} =-\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

$$F_{prop,x^{\prime}}=F_{prop,x}\cos\theta-F_{prop,y}\sin\theta$$

$$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

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

$$\bar{F}_{prop,y^{\prime}} =\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

$$\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

$$u^{i}_{x^{\prime}} =u^{i}_{x}\cos\theta-u^{i}_{y}\sin\theta$$

$$u^{i}_{y^{\prime}} =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:

$$\langle\dot{Q}\rangle =-\frac{1}{N}\sum_{i}^{N}\int_{0}^{L}f_{y}\frac{\partial h^{i}}{\partial t}dx^{i}$$

$$\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:

$$\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)