Introduction

Flow interactions are thought to allow flying and swimming animals to derive energetic benefits when moving in groups [1]. However, direct assessment of such benefits is challenging, chiefly because animal groups do not generally conform to regular patterns – individuals in these groups dynamically change their position relative to their neighbors [2–6]. Also, because direct energetic measurements in moving animals, flying or swimming, are notoriously difficult and often unreliable as proxy for hydrodynamic energy savings [4, 7–12]. These difficulties hinder a direct mapping from the spatial pattern of the group to the energetic benefits or costs incurred at each position within the group.

An understanding of how the spatial arrangement of individuals within a group influences their cost of locomotion can provide insights into the evolution of social structures, resource allocation, and overall fitness of each individual in cooperative activities such as foraging, mating, and evasion [13–19]. It could also guide the design of bio-inspired engineering systems and algorithms that steer groups of entities, such as swarms of autonomous robotic vehicles, underwater or in flight, that collaborate to achieve a desired task while minimizing energy consumption and improving the overall system efficiency [20–24].

To understand the potential energetic benefits of group movement, various direct and indirect approaches have been employed. Li et al. [9] associated energy savings in pairs of flapping robotic swimmers with a linear relationship between their flapping phase lag and relative distance. Based on this, a strategy, called Vortex Phase Matching, was extrapolated for how fish should behave to maximize hydrodynamic benefits: a follower fish should match its tailbeat phase with the local vorticity created by a leader fish. Pairs of freely swimming fish seemed to follow this linear phase-distance relationship even with impaired vision and lateral line sensing, that is, in the absence of sensory cues about their relative position and neighbor-generated flows. Interestingly, the same linear phase-distance relationship was uncovered independently in flapping hydrofoils and accredited solely to flow interactions [25–27]. It is therefore unclear whether vortex phase matching is an active strategy, mediated by sensing and feedback control, that fish employ to minimize energy expenditure, or if it arises passively through flow interactions between flapping swimmers. Importantly, active or passive, it is unclear if this strategy scales to groups larger than two.

In an effort to directly gauge the energetic benefits of schooling, metabolic energetic measurements were recently performed in solitary and groups of eight fish, and impressive energetic savings were attributed to schooling compared to solitary swimming when the fish were challenged to swim at high speeds [12]. Lamentably, the study made no mention of the spatial patterns assumed by these physically-thwarted individuals [12]. In an independent previous study [5], changes in spatial patterns and tailbeat frequencies were reported in similar experiments, albeit with no energetic measurements. Specifically, [5] showed that, when challenged to sustain higher swimming speeds, the fish in a group rearranged themselves in a side-by-side pattern as the speed increased, presumably to save energy.

Taking together the results of [5, 12], are we to conclude that side-by-side formations are more energetically beneficial than, say, inline or diagonal formations? The answer is not simple! The metabolic measurements of [4] in a school of eight fish report that side-by-side formations, though beneficial, produce the least energetic savings compared to diagonal formations [1]. In an experimental study of a single fish interacting with a robotic flapping foil, the freely-swimming fish positioned itself in an inline formation in the wake of the flapping foil, supporting the hypothesis that swimming in the wake of one another is an advantageous strategy to save energy in a school [11]. Why then the fish in the experiments of [5] self-organized in a side-by-side formation when challenged to swim at higher swimming speeds?

The answer is not simple because ample hydrodynamic mechanisms for energy savings in fish schools have been stipulated for each possible configuration – side-by-side, inline and diagonal (see, e.g., Fig. 1 of [12]) – but no assessment is provided of the relative advantage of these configurations. For example, side-by-side formations, where fish mirror each other by flapping antiphase, are thought to create a wall-like effect that reduces swimming cost [12, 28]. A fish swimming in the wake between two leading fish encounters a reduced oncoming velocity, leading to reduced drag and thrust production [1]. Inline formations, where fish swim in tandem, are thought to provide benefits to both leader and follower, by an added mass push from follower to leader [28, 29]and a reduced pressure on the follower [30]. All of these mechanisms can in principle be exploited by schooling fish as they dynamically change their relative spacing in the group. But are these mechanisms equally advantageous? Or is there a hierarchy of hydrodynamic benefits depending on the relative position within the school? The literature offers no comparative analysis of the energetic savings afforded by each of these configurations.

Figure 1:

The study of [4] is arguably the closest to addressing this question, but, to map the energetic benefits for pair-wise configurations, the authors employed statistical averages in a school of eight fish, thus inevitably combining the various hydrodynamic mechanisms at play and cross-polluting the estimated benefits of each configuration. A cleaner analysis in pairs of flapping foils shows that these relative positions – side-by-side, inline, and diagonal – all emerge spontaneously and stably due to flow interactions [27], but provides no method for estimating the energetic requirements of these formations, let alone comparing them energetically. Even vortex phase matching makes no distinction between side-by-side, inline, or diagonal pairs of fish [9]. It simply postulates that an unknown amount of energetic benefit is acquired when the linear phase-distance relationship is satisfied. Thus, to date, despite the widespread notion that group movement saves energy, a direct comparison of the energetic savings afforded by different spatial formations remains lacking. Importantly, it is unknown whether and how the postulated benefits scale with increasing group size.

Here, to circumvent the challenges of addressing these questions in biological systems, we formulate computational models that capture the salient hydrodynamic features of single and pairs of swimming fish. Namely, we represent each fish as a freely-swimming hydrofoil undergoing pitching oscillations about its leading edge. A single flapping hydrofoil shares many hydrodynamic aspects with its biological counterpart, including an alternating, long-lived pattern of vorticity in its wake [31–35]. These similarities have been demonstrated repeatedly, within biologically relevant ranges of flapping parameters [31, 36], for different geometries [37–40], material properties [40–42], and flapping kinematics [43–45]. In this study, we show, based on our own simulations and by conducting a thorough literature survey, that flow interactions, with no sensing and feedback control, lead to emergent formations that preserve the linear phase-distance relationship uncovered independently in live and robotic fish [9, 11] and in flapping hydrofoils [25–27]. This relationship is preserved irrespective of geometry [26, 27, 46], material properties [25, 47–49], and flapping kinematics [50, 51]. The universality of this relationship serves as strong validation of our models and anchors our subsequent exploration of the opportunities for hydrodynamic benefits available in a given flow field.

Importantly, we go beyond two swimmers to investigate flow interactions in larger groups and find that inline formations differentially distribute hydrodynamic savings to members within the school, favoring trailing swimmers, but only up to a certain school size, while side-by-side formations equally distribute hydrodynamic savings and scale to arbitrary number of swimmers without loss of cohesion. Our findings provide a direct mapping from the school’s spatial pattern to the energetic savings experienced by its members; inversely, our methods allow us to identify beneficial positions in the flow field created by single or multiple swimmers in terms of energetic savings and overall school cohesion. Importantly, our results raise an interesting hypothesis that the dynamic repositioning of members within a fish school could be driven by greed and competition to occupy hydrodynamically advantageous positions and open up opportunities for analyzing and comparing the role flow physics in the evolution of cooperative versus greedy behavior in animal groups.

Results

Mathematical models of flow-coupled flapping swimmers

Inspired by the experiments of [9, 26], we study self-organization in the context of flapping swimmers, coupled passively via the fluid medium, with no mechanisms for visual [52–55], flow sensing [56–59], or feedback control [23, 60, 61] (Fig. 1). The swimmers are rigid, of finite body length L and mass per unit depth m, and undergo pitching oscillations of identical amplitude A and frequency f in the (x, y)-plane of motion, such that the pitching angle for swimmer j is given by θ_j = A sin(2πft + ϕ_j), j = 1, 2,…, N, where N is the total number of swimmers. In pairwise interaction, we set ϕ₁ = 0 and ϕ₂ = ϕ, with ϕ being the phase lag between the oscillating pair. We fixed the lateral distance ℓ between the swimmers to lie in the range ℓ [L, L], and allowed the swimmers to move freely in the x-direction in an unbounded two-dimensional fluid domain of density ρ and viscosity µ. Hereafter, all parameters are scaled using the body length L as the characteristic length scale, flapping period T = 1/f as the characteristic time scale, and ρL² as the characteristic mass per unit depth. Accordingly, velocities are scaled by Lf, forces by ρf ²L³, moments by ρf ²L⁴, and power by ρf ³L⁴.

The equations governing the free motion x_j(t) of swimmer j are given by Newton’s second law (here, the downstream direction is positive),

The hydrodynamic forces on swimmer j are decomposed into a pressure force F_j acting in the direction normal to the swimmer and a viscous drag force D_j acting tangentially to the swimmer. These forces depend on the fluid motion, which, in turn, depends on the time history of the states of the swimmers.

To maintain their pitching motions, swimmers exert an active moment M_a about the leading edge, whose value is obtained from the balance of angular momentum. The hydrodynamic power P expended by a flapping swimmer is given by P = M_aθ.

To compute the hydrodynamic forces and swimmers motion, we used two fluid models (Figs. 1, S.1, S.2 and S.3). First, we employed a computational fluid dynamics (CFD) solver of the Navier-Stokes equations tailored to resolving fluid-structure interactions (FSI) based on an adaptive mesh implementation of the immersed boundary method [62–64]. Then, we solved the same FSI problem, in the limit of thin swimmers, using the more computationally-efficient inviscid vortex sheet (VS) model [50, 65–68]. To emulate the effect of viscosity in the VS model, we allowed shed vorticity to decay after a dissipation time τ_diss; larger τ_diss correlates with larger Reynolds number Re in the Navier-Stokes model; see SI for a brief overview of the numerical implementation and validation of both methods.

Flow coupling leads to stable emergent formations

We found, in both CFD and VS models, that pairs of swimmers self-organize into relative equilibria at a streamwise separation distance d that is constant on average, and swim together as a single formation at an average free-swimming speed U (Figs. 1 and 2). We distinguished four types of relative equilibria: inline, diagonal, side-by-side inphase and side-by-side antiphase (Fig. 1).

Figure 2:

Inline formations at ℓ = 0 arise when the follower positions itself, depending on its initial distance from the leader, at one of many inline equilibria, each with its own basin of attraction (Fig. 2A). These inline equilibria occur at average spacing d that is approximately an integer multiple of UT, consistent with previous experimental [26, 46, 69] and numerical [25, 48–50, 70, 71] findings.

When offsetting the swimmers laterally at ℓ = 0 (Fig. 2B), the leader-follower equilibria that arise at ℓ = 0 shift slightly but persist, giving rise to diagonal leader-follower equilibria [27]. Importantly, at a lateral offset ℓ, inphase swimmers (ϕ = 0) that are initially placed side-by-side reach a relative equilibrium where they travel together at a close, but non-zero, average spacing d L. That is, a perfect side-by-side configuration of inphase flapping swimmers is unstable but the more commonly-observed configuration [9] where the two swimmers are slightly shifted relative to each other is stable. This configuration is fundamentally distinct in terms of cost of transport from the mirror-symmetric side-by-side configuration that arises when flapping antiphase at ϕ = π (Fig.3A). Both side-by-side equilibria were observed experimentally in heaving hydrofoils [27], albeit with no assessment of the associated hydrodynamic power and cost of transport.

Figure 3:

We next examined the effect of varying the phase ϕ on the emergent traveling formations. Starting from initial conditions so as to settle on the first equilibrium d/UT 1 when ϕ = 0, and increasing ϕ, we found, in both CFD and VS simulations, that the spacing d/UT at equilibrium increased with increasing ϕ (Fig. 2C). This increase is linear, as evident when plotting d/UT as a function of ϕ (Fig. 3B). Indeed, in Fig. 3B, we plotted the emergent average separation distance d/UT as a function of ϕ for various values of ℓ. Except for the antiphase side-by-side formation, the linear phase-distance relationship ϕ/2π d/UT persisted for ℓ = 0.

The key observation, that pairs of flapping swimmers passively self-organize into equilibrium formations, is independent of both scale and fluid model. In our CFD simulations (Fig. S.1 and S.3), we tested a range of Reynolds number Re = ρUL/µ from 200 to 2000, which covers the entire range of existing CFD simulations [49, 69], where Re O(10²), and experiments [26, 27, 69], where Re O(10³). In our VS simulations, we varied τ_diss from 2.45T to (Figs. S.2 and S.3). Note that the separation distance d is scale-specific and increases with Re; at low Re, a compact inline formation is reached where the two swimmers “tailgate” each other, as observed in [48]. However, the scaled separation distance d/UT remains nearly constant for all Re and τ_diss (Fig. S.3).

The fact that these equilibria emerge in time-forward simulations is indicative of stability [72]. A more quantitative measure of linear stability can be obtained numerically by perturbing each equilibrium, either by applying a small impulsive or step force after steady state is reached [27] or by directly applying a small perturbation to the relative equilibrium distance between the two swimmers and examining the time evolution of d and F to quantify variations in hydrodynamic force δF as a function of signed variations in distance δd from the equilibrium [50]. In either case, we found that the force-displacement response to small perturbations at each equilibrium exhibited the basic features of a linear spring-mass system, where δF/δd is negative, indicating that the hydrodynamic force acts locally as a restoring spring force that causes the initial perturbation to decay and that stabilizes the two swimmers together at their equilibrium relative position. Larger values of δF/δd imply faster linear convergence to the stable equilibrium and thus stronger cohesion of the pairwise formation. Results of this quantitative stability analysis are discussed in subsequent sections.

Emergent formations save energy compared to solitary swimming

We evaluated the hydrodynamic advantages associated with these emergent formations by computing the hydro-dynamic power P_single of a solitary swimmer and P_j of swimmer j in a formation of N swimmers. We calculated the cost of transport COT_j = P_j/mU, of swimmer j and the change in COT compared to solitary swimming ΔCOT_j = (COT_single − COT_j)/COT_single (Fig. 3A). We also calculated the average change in cost of transport for each formation (Fig. 3B). In all cases, except for the antiphase side-by-side formation, in both CFD and VS simulations, the swimmers traveling in equilibrium formations save power and cost of transport compared to solitary swimming. The savings are larger at tighter lateral spacing ℓ.

For inline and diagonal formations, these hydrodynamic benefits are granted entirely to the follower, whose hydrodynamic savings can be as high as 60% compared to solitary swimming (Fig. 3A) [50]. Intuitively, because in 2D flows, vortex-induced forces decay with the inverse of the square of the distance from the vortex location, flow coupling between the two inline or diagonal swimmers is non-reciprocal; the follower positioned in or close to the leader’s wake interacts more strongly with that wake than the leader interaction with the follower’s wake (Figs. 1A, B and S.4A, B).

In side-by-side formations, by symmetry, flow coupling between the two swimmers is reciprocal, or nearly reciprocal in inphase flapping (Figs. 1C, D and S.4C, D). Thus, hydrodynamics benefits or costs are expected to be distributed equally between the two swimmers. Indeed, for inphase flapping, the hydrodynamic benefits are shared equally between both swimmers. For antiphase flapping the cost is also shared equally (Fig. 3A).

The biased distribution of benefits in favor of the follower in inline and diagonal formations could be a contributing factor to the dynamic nature of fish schools [3, 6]. The egalitarian distribution of benefits in the inphase side-by-side formation could explain the abundance of this pairwise configuration in natural fish populations [9] and why groups of fish favor this configuration when challenged to swim at higher speeds [5].

Linear phase-distance relationship in emergent formations is universal

To probe the universality of the linear phase-distance relationship, we compiled, in addition to our CFD and VS results, a set of experimental [9, 26, 46] and numerical [47–51] data from the literature (Table S.1 and Suppl. Excel Table). Data including CFD simulations of deformable flapping flags (☆) [47], (□) [25] and flexible airfoil with low aspect ratio (▷) [49], physical experiments with heaving (○) [26, 46] and pitching (▽) [51] rigid hydrofoils, fish-foil interactions (∗) [11], and fish-fish interactions (△) measured in pairs of both intact and visually- and/or lateral line-impaired live fish [9] are superimposed on Fig. 3B. All data collapsed onto the linear phase-distance relationship ϕ/2π d/UT, with the largest variability exhibited by live fish with close streamwise distance, where the interaction between fish bodies may play a role. The side-by-side inphase formations trivially satisfy this linearity because d/UT ϕ/2π = 0, but the side-by-side antiphase formations don’t satisfy; in the latter, d/UT = 0 while ϕ/2π = 1.

These findings strongly indicate that flow-coupled flapping swimmers passively organize into stable traveling equilibrium formations with linear phase-distance relationship. This relationship is independent of the geometric layout (inline versus laterally-offset swimmers), flapping kinematics (heaving versus pitching), material properties (rigid versus flexible), tank geometry (rotational versus translational), fidelity of the fluid model (CFD versus VS versus particle model), and system (biological versus robotic, 2D versus 3D). Observations that are robust across such a broad range of systems are expected to have common physical and mechanistic roots that transcend the particular set-up or system realization.

Importantly, this universal relationship indicates that flow physics passively positions a swimmer at locations d where the swimmer’s flapping phase ϕ matches the local phase of the wake ϕ_wake = 2πd/UT, such that the effective phase ϕ_eff = ϕ ϕ_wake is zero. Because the quantity UT is nearly equal to the wavelength of the wake of a solitary swimmer, ϕ_wake = 2πd/UT is practically equal to the phase of a solitary leader. These observations have two major implications. First, they are consistent with the vortex phase matching introduced in [9] as a strategy by which fish maximize hydrodynamic benefits. However, they proffer that vortex phase matching is an outcome of passive flow interactions among flapping swimmers, and not necessarily an active strategy implemented by fish via sensing and feedback mechanisms. blueSecond, they led us to hypothesize that emergent side-by-side formations can be predicted from symmetry arguments, while emergent inline and diagonal formations can be predicted entirely from kinematic considerations of the leader’s wake without considering two-way flow coupling between the two swimmers, as discussed next.

Leader’s wake unveils opportunities for stable emergent formations

To challenge our hypothesis that the leader’s wake contains information about the emergent pairwise equilibria, we examined the wake of a solitary swimmer in CFD and VS simulations (Fig. 4A, B), and treated the potential second swimmer as a particle located at a point (x_o, y_o) and undergoing prescribed oscillations y(t) = y_o + A sin(2πft ϕ), with velocity vector v(x_o, y_o, ϕ, t) = 2πAf cos(2πft + ϕ)e_y, where e_y is a unit vector in the transverse y-direction. The wake is thus blind to the existence of the second swimmer. This model differs from the minimal particle model used in [26, 69], which treats both swimmers as particles with minimal ‘wakes’ and considered two-way coupling between them. In our analysis, the leader and its wake can be described to any desired degree of fidelity of the fluid model, including using experimentally constructed flows when available.

Figure 4:

Inspired by [49] and following [50], we defined a flow agreement parameter , where t is chosen after steady state is reached, that describes how well the oscillatory motion v(x_o, y_o, t) of a potential second swimmer matches the local fluid velocity u(x_o, y_o, t) generated by the main solitary swimmer [50]. We normalized . Positive (negative) values of V indicate that the flow at (x, y) is favorable (unfavorable) to the follower’s flapping motion.

In Fig. 4A and B, we show V as a field over the physical space (x, y) for ϕ = 0. Blue regions indicate where the local flow favors the follower’s flapping motion. In both CFD and VS simulations, the locations with the maximum flow agreement parameter closely coincide with the stable equilibria (black circles) obtained from solving pairwise interactions. These findings suggest a simple rule for identifying the locations of stable equilibria in pairs of swimmers: they correspond to locations (x_o, y_o)^∗ of maximum flow agreement V(x_o, y_o, ϕ) between the follower’s flapping motion and the wake of a solitary swimmer. Importantly, these findings underscore that hydrodynamic coupling in pairs of flapping swimmers is primarily non-reciprocal – captured solely by consideration of the effects of the leader’s wake on the follower. This non-reciprocity allows one, in principle, to efficiently and quickly identify opportunities for hydrodynamic benefits in any wake, without the need to perform costly two-way coupled simulations and experiments.

To verify this proposition, we show in Fig. 5A, as a function of phase ϕ, the streamwise locations of the local maxima of V(x_o, y_o, ϕ) computed based on the CFD and VS models, and scaled by UT, where U is the speed of the solitary swimmer. We superimpose onto these results the equilibrium configurations obtained from pairwise interactions in the context of the CFD (), VS (■), and time-delay particle () models, where we modified the latter to account for non-zero lateral offset ℓ (SI, Sec. S.3 and Fig. S.5). Predictions of the equilibrium configurations based on maximal flow agreement parameter agree remarkably well with actual equilibria based on pairwise interactions, and they all follow the universal linear phase-distance relationship shown in Fig. 2B.

Figure 5:

The wake of a solitary leader contains additional information that allows us to evaluate the relative power savings of a potential follower and relative stability of the pairwise formation directly from the leader’s wake, without accounting for pairwise interactions. Assessment of the relative power savings follows directly from the maximal value of the flow agreement parameter: larger values imply more power savings and reduced cost of transport. To verify this, we calculated the maximal V in the wake of the solitary swimmer, where we expect the follower to position itself in pairwise interactions. In Fig. 5B, we plotted these V values as a function of lateral distance ℓ for ϕ = 0. We superimposed the power savings ΔP based on pairwise interactions of inphase swimmers using the CFD and VS simulations and normalized all quantities by the maximal value of the corresponding model to highlight variations in these quantities as opposed to absolute values. Power savings are almost constant for ℓ < 0.25L, but decrease sharply as ℓ increases. This trend is consistent across all models, with the most pronounced drop in the CFD-based simulations because the corresponding velocity field u decays more sharply when moving laterally away from the swimmer.

To assess the cohesion of the pairwise formation from the wake of the solitary leader, we estimated the thrust force acting on the virtual follower based on the fact that the thrust magnitude scales with the square of the swimmer’s lateral velocity relative to the surrounding fluid’s velocity [26, 31, 73]. We defined the thrust parameter . We normalized . At the locations of the maxima of V(x_o, y_o, ϕ), a negative slope ∂T/∂d of the thrust parameter is an indicator of linear stability or cohesion of the potential pairwise equilibria; that is, pairwise formations are expected to be stable if a small perturbation in distance about the maximum V-locations is accompanied by an opposite, restorative change in T. Indeed, in both CFD and VS wakes, ∂T/∂d at the maximum V-locations is negative (Fig. 4).

In Fig. 5C, we plotted ∂T/∂d as a function of lateral distance ℓ for ϕ = 0. We superimposed the magnitude of the eigenvalues δF/δd obtained from the linear stability analysis of pairwise interactions in inphase swimmers using the VS and time-delay particle models. As in Fig. 5B, all quantities are normalized by the maximal value of the corresponding model to highlight variations in these quantities as opposed to absolute values. Also, as in Fig. 5B, all models produce consistent results: pairwise cohesion is strongest for ℓ < 0.25L, but weakens sharply as ℓ increases, with the most pronounced drop in the CFD-based simulations.

These findings show that, in the emergent inline and diagonal formations of flow-coupled flapping swimmers, the dominant dynamics is non-reciprocal; that is, the leader affects the follower, but the follower has a negligible (potentially higher-order) effects on the leader. Our diagnostics tools (flow agreement and thrust parameters) are agnostic to how the flow field is constructed (CFD vs VS simulations), and can be equally applicable to experimental data. That is, the approach we developed here could be applied broadly to analyze, predict, and test opportunities for schooling and hydrodynamic benefits in live and robotic fish when measurements of the flow field are available.

Parametric analysis over the entire space of phase lags and lateral offsets

Having demonstrated consistency in the emergence of flow-mediated equilibria in both CFD and VS simulations, we next exploited the computational efficiency of the VS model to systematically investigate emergent pairwise formations over the entire space of phase lag ϕ [0, 2π) and lateral offset ℓ [L, L], excluding side-by-side antiphase formations.

Equilibrium configurations are dense over the entire range of parameters: for any combination of phase lag ϕ and lateral offset ℓ, there exists an emergent equilibrium configuration where the pair of swimmers travel together at a separation distance d/UT (Fig. 6A). Perturbing one or both parameters, beyond the limits of linear stability, causes the swimmers to stably and smoothly transition to another equilibrium at different spacing d/UT. Importantly, increasing the phase lag ϕ shifts the equilibrium positions in the streamwise direction such that d/UT depends linearly on ϕ, but the effect of lateral distance for ℓ L is nonlinear and nearly negligible for small ℓ: increasing the lateral offset ℓ by an entire bodylength L changes the pairwise distance d/UT by about 15%. Our results explored emergent equilibria up to d/UT 2.5 and are consistent with the experimental findings in [27], which explored up to nine downstream equilibria.

Figure 6:

To assess the hydrodynamic advantages of these emergent formations, we calculated the average change in hydrodynamic power per swimmer. The pair saves power compared to solitary swimming (Fig. 6B). Power savings vary depending on phase lag ϕ and lateral distance ℓ: for the entire range of ϕ from 0 to 2π, the school consistently achieves over 20% power reduction, as long as the lateral offset is ℓ 0.25L. However, increasing ℓ from 0.25L to L reduces significantly the hydrodynamic benefit. That is, swimmers can take great liberty in changing their phase without compromising much the average energy savings of the school, as long as they maintain close lateral distance to their neighbor.

A calculation of the linear stability of each equilibrium in Fig. 6A shows that these emergent formations are linearly stable (Fig. 6C), and the degree of stability is largely insensitive to phase lag, with strongest cohesion achieved at lateral offset ℓ 0.25L. The results in Fig. 6A-C are constructed using pairwise interactions in VS simulations, but can be inferred directly from the wake of a solitary leader, as discussed in the previous section and shown in Fig. 6D-F.

Analysis of larger groups of inline and side-by-side swimmers

How do these insights scale to larger groups? To address this question, we systematically increased the number of swimmers and computed the emergent behavior in larger groups based on flow-coupled VS simulations.

In a group of six swimmers, all free to move in the streamwise x-direction, we found that the last three swimmers split and form a separate subgroup (Fig. 7A). In each subgroup, swimmer 3 experiences the largest hydrodynamic advantage (up to 120% power saving!), swimmer 2 receives benefits comparable to those it received in pairwise formation (65% power saving), and swimmer 1 no benefit at all (Fig. 7C).

Figure 7:

We asked if loss of cohesion is dependent on the number of inline swimmers. To address this question, we gradually increased the number of swimmers from two to six (Fig. S.6). We found that in a school of three inline swimmers, flow interactions led to a stable emergent formation with hydrodynamic benefits similar to those experienced by the three swimmers in each subgroup of Fig. 7A and C. When computing the motion of four inline swimmers (Fig. 8A), we found that the leading three swimmers maintained cohesion, at hydrodynamic benefits similar to a formation of three, but swimmer 4 separated and lagged behind, receiving no advantage in terms of power savings because it split from the formation (Figs. 8D and S.6). In a group of five, the last two swimmers split and formed their own subgroup. That is, in all examples, swimmer 4 consistently lost hydrodynamic advantage and served as local leader of the trailing subgroup. These observations are consistent with [48] and demonstrate that flow interactions alone are insufficient to maintain inline formations as the group size increases.

Figure 8:

We next explored the robustness of the side-by-side pattern to larger number of swimmers starting from side-by-side initial conditions (Fig. 7B). The swimmers reached stable side-by-side formations reminiscent of the configurations observed experimentally when fish were challenged to swim at higher swimming speeds [5]. The swimmers in this configuration saved power compared to solitary swimming (Fig. 7C): swimmers gained equally in terms of hydrodynamic advantage (up to 55% power saving for the middle swimmers in a school of six), except the two edge swimmers which benefited less. We tested these results by gradually increasing the number of swimmers from two to six (Fig. S.7). The robustness and overall trend of power saving among group members is robust to the total number of swimmers in these side-by-side formations.

Mechanisms leading to loss of cohesion in larger inline formations

To understand why three swimmers form a stable inline formation but four don’t, we extended the analysis in Fig. 4 to analyze the wake created behind two-swimmer (Fig. 9A) and three-swimmer (Fig. 9B) groups. Specifically, we computed pairwise interactions in a two-swimmer school and considered the combined wake of both swimmers after they had settled onto an equilibrium state. Similarly, we computed the behavior of a three-swimmer school and analyzed the combined wake at steady state. Compared to the single leader wake in Fig. 4B, in the wake of a two-swimmer school, positive flow agreement in the (blue) region is enlarged and enhanced, corresponding to swimmer 3 receiving the largest power savings. On the other hand, behind three inline swimmers, the region of positive flow agreement is weakened and shrunk, indicating weaker potential for energy saving by a fourth swimmer.

Figure 9:

Importantly, in the wake of the pairwise formation, the downstream jet is modest at the location of maximum V, where swimmer 3 is expected to position itself for hydrodynamic benefit, thus allowing swimmer 3 to reach this position and stay in formation (Fig. 9E). Also, at this location, the wake has a substantial transverse velocity u e_y (Fig. 9G), which aids thrust production at a diminished cost. In contrast, three inline swimmers generate a much stronger downstream jet at the location of maximum V where swimmer 4 is expected to position itself (Fig. 9F). This jet prevents swimmer 4 from stably staying in formation, and the transverse flow velocity u e_y is nearly zero for the entire flapping period (Fig. 9H) indicating little opportunity for exploiting the flow generated by the three upstream swimmers for thrust generation. This limitation is fundamental; it results from the flow physics that govern the wake generated by the upstream swimmers. There is not much that a trailing swimmer can do to extract hydrodynamic benefits from an oncoming flow field that does not offer any.

Critical size of inline formations beyond which cohesion is lost

We sought to understand what determines the critical group size, here three, beyond which inline formations lose cohesion and split into subgroups. Because we have established that the flow agreement parameter V plays an important role in predicting emergent formations, we first examined V in the wake of a pair of flapping swimmers in CFD (Fig. S.1) and VS (Fig. S.2) simulations. These results show that at lower Re and smaller dissipation time τ_diss, the flow agreement parameter V decays rapidly downstream of the flapping swimmers, thus diminishing the opportunities for downstream swimmers to passively stay in cohesive formation and achieve hydrodynamic benefits. We thus hypothesized that the number of swimmers that passively maintain a cohesive inline formation is not a universal property of the flow physics, but depends on the flow regime.

We tested this hypothesis in VS simulations with increasing number of swimmers and increasing τ_diss. As we increased τ_diss, the number of swimmers that stayed in cohesive inline formation increased (Fig. 8B,C). These findings confirm that this aspect of schooling – the maximal number of swimmers that passively maintain a cohesive inline formation – is indeed scale dependent. Interestingly, an analysis of the power savings in these formations shows that, although swimmers 4 and 5 stay in formation at increased τ_diss, swimmer 3 always receives the most hydrodynamic benefit (Fig. 8D).

Discussion

We analyzed how passive flow interactions mediate self-organization in groups of flapping swimmers, all heading in the same direction. Our approach relied on a hierarchy of fluid-structure interaction models and aimed to distill which aspects of self-organization are universal and those which are scale-dependent.

We found that a pair of flapping swimmers self-organize into inline, diagonal, or side-by-side formations (Fig. 1). The emergent formation depends on the swimmers’ flapping phase and initial conditions. In fact, the distinction between these types of formation is somewhat arbitrary because, as phase varies, the emergent equilibria are dense over the space of lateral offset and separation distance (Fig. 6). These findings are consistent with experimental observations [9, 11, 26, 27], but go beyond these observations to quantify the hydrodynamic benefits to each member in these formations. Two side-by-side swimmers flapping inphase share the hydrodynamic benefits nearly equally, whereas contrary to common misconception [12], to flap antiphase, both swimmers need to exert extra effort compared to solitary swimming. In leader-follower formations, whether inline or diagonal, hydrodynamic benefits are bestowed entirely on the follower (Fig. 3A).

Importantly, we showed that the wake of a solitary leader contains information that unveils opportunities for the emergence of stable and energetically-favorable formations in pairs of swimmers. Equilibrium locations, and trends in power savings and school cohesion, can all be predicted entirely from kinematic considerations of the leader’s wake with no consideration of the two-way coupling between the two swimmers (Fig. 5). These results are important because they highlight the non-reciprocal or asymmetric nature of flow coupling in leader-follower configuratoins, inline or diagonal, of flapping swimmers at finite Re and open new avenues for future studies of non-reciprocal flow-coupled oscillators. This is in contrast to classic mechanical and biological oscillators, such as Huygens pendula or viscosity-dominant oscillators, where the coupling between oscillators is reciprocal; see, e.g., [74–78].

Our analysis has practical importance in that it provides efficient diagnostics and predictive tools that are equally applicable to computational models and experimental data and could, therefore, be applied broadly to analyze, predict, and test opportunities for schooling and hydrodynamic benefits in live and robotic fish when flow measurements are available.

Case in point, we used these diagnostic tools to explain the mechanisms leading to scattering in larger groups of inline swimmers and to predict when the wake of a leading group of swimmers offers no opportunities for a follower to benefit from passive hydrodynamics (Fig. 9). At an increasing number of flow-coupled swimmers, side-by-side formations remain robust, but inline formations become unstable beyond a critical number of swimmers (Figs. 7 and 8). The critical number depends on the fluid properties and can be predicted by analyzing the wake of the leading group of swimmers. Future work will focus on testing these findings experimentally and in CFD simulations with increasing number of swimmers, together with accounting for body deformations [68], lateral dynamics [51], and varying flapping amplitude and frequency [26].

We next reflect on the distribution of hydrodynamic benefits among group members, its dependence on the spatial pattern of the emergent formation, and the implications that these observations have on biological schools of fish. Field and laboratory experiments [2, 4, 5] have shown that actual fish schools do not generally conform to highly regularized patterns, and schooling fish dynamically change their position in the school. Neighboring fish vary from side-by-side to inline and diagonal configurations [11]. Importantly, in laboratory experiments that challenged groups of fish to sustain high swimming speeds, the fish rearranged themselves in a side-by-side pattern as the speed increased, much like the pattern in Fig. 7B, presumably to save energy [5]. These empirical observations, together with our findings that side-by-side formations provide the fairest distribution of efforts among school members (Fig. 7B and C), offer intriguing interpretations of the results in [5]: when the fish are not challenged by a strong background current to sustain high swimming speeds, they position themselves as they please spatially, without much consideration to equal sharing of hydrodynamic benefits. But when challenged to swim at much higher speeds than their average swimming speed, fish are forced to cooperate.

To take this notion a step further, our results hint at a connection between flow physics and what is traditionally thought of as social traits: greed versus cooperation. We posit that there is a connection between the resources that arise from flow physics – in the form of energetic content of the wake of other swimmers – and greedy versus cooperative group behavior. In cohesive inline formations, the leader is always disadvantaged and hydrodynamic benefits are accorded entirely to trailing swimmers (Figs. 3A and 7C). Importantly, flows generated by these inline formations present serious impediments for additional swimmers to join the line downstream (Figs. 7 and 8). Thus, we could call these formations greedy, leaving no resources in the environment for trailing swimmers. This thought, together with our interpretation of the observations in [5] that cooperation to achieve an egalitarian distribution of hydrodynamic benefits is forced, not innate, raise an interesting hypothesis. The dynamic repositioning of members within the school could be driven by greed and competition to occupy hydrodynamically advantageous positions, much like in peloton of racing cyclists [79]. On a behavioral time scale, these ideas, besides their relevance to schooling fish, open up opportunities for analyzing and comparing the collective flow physics in cooperative versus greedy behavior in animal groups from formations of swimming ducklings [80] and flying birds [81] to peloton of racing cyclists [79]. From an evolutionary perspective, it is particularly exciting to explore the thought that flow physics could have acted as a selective pressure in the evolution of social traits such as cooperation and greed in aquatic animal groups.

Acknowledgements

Many thanks to members of the McHenry Laboratory at UC Irvine and Kanso Laboratory at USC for numerous discussions about fish schooling behaviour and for comments on the manuscript.

Funding

Funding provided by the National Science Foundation grants RAISE IOS-2034043 and CBET-210020 (EK) and the Office of Naval Research grants N00014-22-1-2655 and N00014-19-1-2035 (EK).

Author contributions

E.K. conceptualized the study. S.H. and H.H. developed the CFD and VS codes and performed the numerical simulations. S.H., H.H. and E.K. analyzed the data. E.K. wrote the manuscript. S.H. and H.H. provided manuscript edits and comments and approved the final version.

Competing interests

Authors declare that they have no competing interests.

Data and materials availability

All data generated or analysed during this study are included in this published article (and its supplementary information files).

Supplementary Materials

Supplementary information; Figs. S1 to S7; Table S1; and Excel Sheet.

Supplementary materials

CFD model

In our CFD simulations, a swimmer is modeled as a symmetric 2D Joukowsky airfoil [82]. The chord length of the airfoil is the characteristic length L, and the maximum thickness is 0.12L. The airfoil undergoes pitching motion around its leading edge. Fluid-structure interactions are governed by the incompressible Navier-Stokes equations,

where u(x, t) and p(x, t) are the velocity and pressure field, respectively. We solved for these fields numerically using immersed boundary method (IBM) that handles the two way coupled fluid structure interaction [63, 64, 83–86].

The immersed boundary formulation involves an Eulerian descriptions of the flow field and a Lagrangian description of the immersed swimmers, modeled as Joukowsky airfoils. The boundary condition is mapped to a body force exerted on the fluid. The Lagrangian and Eulerian variables are correlated by the Dirac delta function, which is smoothed during discretization. Here, we used the implementation developed by the group of Professor Boyce Griffith, IBAMR [62], which has long been used to solve problems such as blood flow in heart [83, 87], water entry/exit problems [88], fish’s swimming [89–91], insect’s flight [92, 93], flexible propulsors [94–96], self propulsion of pitching/heaving airfoil [96, 97], and fish schooling [97, 98]. This implementation is based on an adaptive mesh, which enables us to accurately simulate self propulsion and reach steady state in a large computational domain with a reasonable computational cost.

The computational domain is a rectangle of dimensions 80L 20L, with periodic boundary conditions on the computational domain and no-slip boundary condition on the surface of airfoils. The initial location of the first swimmer is 12L away from the right boundary in streamwise direction. The initial distance between the two swimmers d(t = 0) ranges from 1.5L to 4L to ensure we access different equilibria the emerge in these pairwise formations. The coarsest Eulerian mesh is a uniform 500 125 Cartesian grid. The computational domain close to the airfoils and their wake are refined. There are 3 layers of refinement mesh, and the refinement ratio for each layer is 4. The simulation timestep is adaptive, with the maximum timestep Δt_max = 2.5 10⁻³.

The hydrodynamic forces F_x, F_y and moment M acting on each swimmer are calculated by the integrating, over the surface of that swimmer, the traction force σ n and moment x (σ n), where σ = pI +µ(u + u^T) is the fluid stress tensor, x denotes positions on the surface of that swimmer and n the unit normal to that airfoil into the fluid. The hydrodynamic forces and moment are used to solve the equations of motion (1) for each swimmer and obtain the active moment M_a in (7) needed to evaluate the hydrodynamic power in (8) that is expended by each swimmer.

Vortex sheet (VS) model

In the vortex sheet model, each swimmer is modeled as a rigid plate of length L, and it is approximated by a bound vortex sheet, denoted by l_b, whose strength ensures that no fluid flows through the rigid plate, and the separated shear layer is approximated by a free regularized vortex sheet l_w at the trailing edge of the swimmer. The total shed circulation Γ in the vortex sheet is determined so as to satisfy the Kutta condition at the trailing edge, which is given in terms of the tangential velocity components above and below the bound sheet and ensures that the pressure jump across the sheet vanishes at the trailing edge. To express th_√ese concepts mathematically, it is convenient to introduce the complex notation, such that z = x + iy, where i = −1. For more details on the vortex sheet model and our implementation of it, we refer the readers to appendix A of [50], as well as to [65–68, 99–103].

The pressure difference across the infinitely thin swimmers n[p]_∓ is given by integrating the balance of momentum equation for inviscid planar flow along a closed contour containing the vortex sheet and trailing edge.

Here, n = n_x + in_y is the unit normal, in complex notation, and [p]_∓ is the jump in pressure across the swimmer. The hydrodynamic force F = _lb n[p]_∓ds acting on the swimmer is given by,

The hydrodynamic moment M acting on the swimmer about its leading edge is given by

where z_l.e. is position of the leading edge s = 0.

We introduce a drag force D that emulates the effect of skin friction due to fluid viscosity in the context of the vortex sheet model [50, 104]. Namely, following [50, 104], we write the drag force for a swimming plate

where C_dis a drag coefficient and U_± are the spatially-averaged tangential fluid velocities on the upper and lower side of the plate, respectively, relative to the swimming velocity U,

where u_∓(s, t) denotes the tangential slip velocities on both sides of the plate. We estimate C_d to be approximately 0.05 in the experiments of [46].

Additionally, following [50, 66–68], we emulate the effect of viscosity on the wake itself by allowing the shed vortex sheets to decay gradually by dissipating each incremental point vortex after a finite time τ_diss from the time it is shed into the fluid. Larger τ_diss implies that the incrementally-shed vorticity along the vortex sheet stays in the fluid for longer times, mimicking the effect of lower fluid viscosity. We studied the effect of dissipation time in Figs S.2 and S.3, and used τ_diss = 2.45T in the rest of this paper for the sake of computational cost.

Pitching motions are produced by an active moment M_a imposed by the swimmer on the surrounding fluid about the leading edge. The value of M_a is obtained from the balance of angular momentum about the swimmer’s leading edge (l.e.),

In the VS model, I = mL²/3 is the swimmer’s moment of inertia about the leading edge, w_l.e. is the swimmer’s velocity at the leading edge (in complex form), and M is the hydrodynamic moment about the leading edge given in (4).

The hydrodynamic power P expended by the flapping swimmer to maintain its pitching motion is given by

Time-delay particle model

We supplemented these CFD and VS models by studying pairwise interactions in the context of the minimal particle model used in [26, 69]. This particle model was designed for inline swimmers. Here, we modify it slightly to account for lateral offset ℓ between the swimmers as shown in Fig S.5. In a nutshell, the model assumes that the leader leaves behind a vertical wake speed equal to the leader’s flapping speed at the tail. The speed of the ‘wake’ that is left behind decays exponentially in time, as an approximation for viscous dissipation. Details of the model can be found in the Supplementary Information of [26].

In this model, each swimmer is a particle of mass per unit depth m undergoing vertical oscillations such that y₁ = a sin(2πft) and y = a sin(2πft ϕ), where a = L sin A is the tailbeat amplitude of the pitching airfoil. In this model, each particle is assumed to experience a thrust force F_j that is proportional to the square of its vertical velocity relative to the surrounding fluid, and a drag force D_j relative to its relative horizontal speed. Namely,

Where

Here, u_x and u_y are the x- and y-components of the fluid velocity, C_T, C_D are constant thrust and drag coefficients. Since the leader swims into quiescent fluid, we assume that u_y(x₁, y₁) = 0. The follower swims into the wake of the leader, which we assume to have zero horizontal velocity (u_x = 0) and vertical velocity that decays exponentially both in time and in the lateral direction,

where Δt is the delay time between the leader and follower, that is, the time past since the leader passed by the follower’s current location: x₁(t Δt) = x₂(t) We added the last term in (11) to consider decay in the lateral direction ℓ; to estimate the parameters p and h, we used a best curve fit to the data of period-average velocity magnitude versus lateral distance in the wake of a single swimmer in the vortex sheet model (Fig. S.5). We numerically integrated (10) and solved for the motion of the follower. At steady state, we computed the separation distance d = x₂ x₁ between the pair acting on the follower for a range of ϕ [0, 2π] and ℓ [L, L] (Fig. S.5). The parameter values are chosen to be consistent with the experiments of Newbolt et al. (provided in Table S1 of SI in [26]). Specifically, ρ = 1 g/cm³, L = 4 cm, m = 5.3 g/cm, C_D = 0.25, C_T = 0.96, τ = 0.5 sec.

Data collection from published literature

To probe the relationship between the scaled separation distance d/UT in pairs of swimmers and their flapping phase lag ϕ, we collected data from published literature on pairs of interacting swimmers [9, 26, 47, 50, 51]; see Table S.1 and the Supplemental Excel Sheet. We excluded a few studies from this literature survey [70, 71], because we were unable to extract d/UT from the data they provided. We also excluded [69] because their study involved fixed inter-swimmer distances.

We found that the relationship between phase lag ϕ and separation distance d/UT is approximately linear ϕ/2π = d/UT + β. The value of β depends on how one defines the inter-swimmer separation distance. The two definitions that are commonly used in the literature are tail-to-tail or tail-to-head (gap) distance. These definitions do not change the linear scaling between d/UT and ϕ/2π; they only change the value of β. For the data provided in Fig. 3 of the main text, we used the definition of separation distance d that results in β = 0. Not all literature provided values for U or T. In the very few cases when this information were missing, we estimated UT from the linear relationship ϕ/2π = d/UT [47].