Lattices of hydrodynamically interacting flapping swimmers

TL;DR

This paper introduces a novel discrete-time model to analyze hydrodynamic interactions in flapping swimmer lattices, revealing formation-dependent benefits and predicting bistability and instabilities in collective swimming behaviors.

Contribution

It develops a new dynamical system framework for modeling long-range hydrodynamic interactions in swimmer lattices, bridging a gap in theoretical understanding of collective fluid-mediated dynamics.

Findings

1D lattice models match experimental bistability data.

2D lattices, especially diamond formations, maximize hydrodynamic benefits.

Formation structure significantly influences collective swimming efficiency.

Abstract

Fish schools and bird flocks exhibit complex collective dynamics whose self-organization principles are largely unknown. The influence of hydrodynamics on such collectives has been relatively unexplored theoretically, in part due to the difficulty in modeling the temporally long-lived hydrodynamic interactions between many dynamic bodies. We address this through a novel discrete-time dynamical system (iterated map) that describes the hydrodynamic interactions between flapping swimmers arranged in one- and two-dimensional lattice formations. Our 1D results exhibit good agreement with previously published experimental data, in particular predicting the bistability of schooling states and new instabilities that can be probed in experimental settings. For 2D lattices, we determine the formations for which swimmers optimally benefit from hydrodynamic interactions. We thus obtain the following hierarchy: while a side-by-side single-row "phalanx" formation offers a small improvement over a solitary swimmer, 1D in-line and 2D rectangular lattice formations exhibit substantial improvements, with the 2D diamond lattice offering the largest hydrodynamic benefit. Generally, our self-consistent modeling framework may be broadly applicable to active systems in which the collective dynamics is primarily driven by a fluid-mediated memory.

Equations

$$u_{n+1}=u_n+\frac{T}{m_{\text{e}}}\left(F_D(u_n)+F_x\left[x_n,u_n,v_n,{\omega }_n(z)\right]\right),$$

$$x_{n+1}=x_n+u_{n+1}T,$$

$${\omega }_{n+1}(z)={\omega }_n(z){\mathrm{e}}^{-T/\tau }+(-1{)}^n\gamma \sum _{j=-1}^1\delta \left(z-(z_{n+1}+jL)\right),$$

$$\text{where }{z}_{n+1}=2a+{\hat{x}}_{n+1}+i(-1{)}^n{h}_0,$$

$$x_{n+1}=x_n+u_{n+1},$$

$$u_{n+1}=u_n+F_D(u_n)-F_0$$

$$-\sum _{j=-1}^1\sum _{k=-\infty }^{n}G\left(x_n-{\hat{x}}_k-jL,u_n,n-k\right),$$

$$F_D(U)=F_0+\sum _{j=-1}^1\sum _{m=0}^{\infty }G(U(m+1/2)-jL,U,m).$$

$$\mathcal{F}(z)$$

$$+g_{x}z,$$

$$F_x^{\text{v}}$$

$$\text{where}{V}_{\text{f}}$$

$$\pm 2U_1+F_D(U\pm U_1)=F_0$$

$$+\sum _{j=-2}^2\sum _{m=0}^{\infty }G\left(U\left(m+\frac{1}{2}\right)\pm \frac{U_1}{2}-\frac{jL}{2}\mp {\alpha }_{j}d,U\pm U_1,m\right),$$

Full text

††preprint: Preprint #

I Introduction

The complex collective dynamics of fish schools and bird flocks have long fascinated physicists, biologists and mathematicians Pavlov and Kasumyan (2000); Bajec and Heppner (2009). In addition to their biological relevance, they are living examples of active systems Cavagna and Giardina (2014); Marchetti et al. (2013); Ramaswamy (2010); Saintillan and Shelley (2013) in which energy input by the individual constituents gives rise to organized collective phenomena. While there has been considerable experimental and theoretical progress in characterizing “dry”  active systems (e.g. shaken granular rods Narayan et al. (2007); Kudrolli et al. (2008)) and the collective behavior of biological systems at the microscale (e.g. bacterial suspensions Sokolov and Aronson (2009); Wensink et al. (2012); Zhang et al. (2010); Saintillan and Shelley (2008)), significantly less is known about the role of hydrodynamic interactions in mediating schooling and flocking behavior in collectives of larger animals. More generally, the influence of inertial fluid flows and the consequent long-lived hydrodynamic interactions on collective behavior remains poorly understood. In contrast to the low Reynolds number (Stokes) regime in which microswimmers operate, fluid-mediated memory could significantly impact animal schools and also nonliving systems dominated by wave-particle interactions. As an example of the latter, oil droplets bouncing on a vertically vibrating fluid bath Bush (2015) are known to exhibit crystal-like bound states Lieber et al. (2007); Eddi et al. (2009) as a result of their coupling through surface waves. We present here a modeling framework for the long-lived hydrodynamic interactions between swimmers, with a view to understanding how their collective dynamics might be mediated by flow-induced forces.

The influence of hydrodynamics on schooling and flocking behavior in biological systems has been the subject of intense debate in the scientific literature. While some analyses of starling cluster flock data Bialek et al. (2012); Attanasi et al. (2014) focus on the behavioral mechanisms behind flocking, a recent analysis Portugal et al. (2014) of flying ibises in V-formations demonstrated coherence of the birds’ wing tip paths, which enables upwash capture from their neighbors to be maximized. Individual fish have been shown to sense and respond to environmental hydrodynamic cues Ristroph et al. (2015), and benefit from external flows by harnessing the energy from vortices Triantafyllou et al. (2000); Liao (2007); Liao et al. (2003). While some studies have argued against a hydrodynamic function for fish schools Partridge and Pitcher (1979) and instead focused on the social interactions between fish Jolles et al. (2017); Swain et al. (2015); Tunstrøm et al. (2013), a number of observations Partridge et al. (1983); Herskin and Steffensen (1998); Marras et al. (2015); Ashraf et al. (2017) have indicated that schooling fish could benefit from hydrodynamic interactions by realizing significant energy savings.

Theoretical models of flocks and schools have also largely ignored hydrodynamic interactions, and instead have shown that self-organized collective behavior may emerge from a relatively simple set of behavioral interaction rules Couzin et al. (2002); Lopez et al. (2012). Examples are the seminal works of Huth & Wissel Huth and Wissel (1992) and Viscek et al. Vicsek et al. (1995), which have been extended to other discrete-time models with more elaborate interaction rules Barberis and Peruani (2016). Far-field hydrodynamic interactions have been recently incorporated into self-propelled particle models of swimmers subject to phenomenological behavioral rules Gazzola et al. (2016); Filella et al. (2018), but such models neglect the vorticity-induced forces thought to be relevant for schooling fish Liao (2007).

Fish schools and bird flocks can exhibit orderly lattice-like formations, although field observations are diverse and sometimes contradictory Pavlov and Kasumyan (2000); Bajec and Heppner (2009). Prior observations revealed that fish schools may adopt lattice configurations in a statistical sense Cullen et al. (1965). A number of fish species (e.g. minnow, bream, saithe, herring) adopt schooling formations reminiscent of 3D tetrahedral and cubic lattices, but others (e.g. cod) adopt less ordered configurations Partridge et al. (1980); Pitcher (1973). Relatively rigid school structures, which are treated in this paper, have also been observed. For example, there are a number of accounts of fish swimming in linear chains Gudger (1944). In their studies of red nose tetra fish, Ashraf et al. Ashraf et al. (2017, 2016) noted a prevalence of 2D diamond lattice formations at low swimming speeds and “phalanx”  formations at higher speeds, the latter being a side-by-side arrangement of swimmers roughly equispaced in a single line perpendicular to the swimming direction. These two schooling formations have also been observed in bluefin tuna Newlands and Porcelli (2008), for which small schools tend to adopt a phalanx formation while larger schools adopt a diamond lattice. With respect to bird flocks, ibises have been observed to obtain aerodynamic and energetic benefits from in-line formation flight Portugal et al. (2014). Certain species (e.g. pelicans Weimerskirch et al. (2001), Canada geese Badgerow and Hainsworth (1981), ibises Portugal et al. (2014)) are thought to benefit from their observed V-formation flight, while it is claimed that others (e.g. pigeons Usherwood et al. (2011), pink-footed geese Cutts and Speakman (1994)) do not. Corcoran & Hedrick Corcoran and Hedrick (2019) have recently identified a new type of ordered configuration in shorebird flocks, the “compound V-formation,”  wherein a given bird commonly flies roughly one wingspan to the side and to the back from the bird in front of it. Such an ordered structure is observed at all spatial scales within an extended flock, unlike the purely local structure exhibited by starling cluster flocks Bialek et al. (2012); Attanasi et al. (2014).

In an attempt to explain such ordered structures, Weihs’ seminal papers Weihs (1973, 1975) considered vortex-induced hydrodynamic interactions in a 2D lattice of swimmers. By positing that fish seek to minimize their hydrodynamic drag while avoiding large flow velocity gradients, Weihs argued that a diamond lattice is the energetically optimal arrangement. While this model has been highly influential, it has not been developed further, due to the lack of experimental confirmation and various modeling assumptions. Specifically, the school’s swimming dynamics was not accounted for; hence, the speed and efficiency were not self-consistently calculated, and there was no consideration of the stability of the most efficient state. Moreover, the influence of the streamwise spacing between swimmers was largely neglected, as the swimmers were assumed to be separated by at least five flapping wavelengths. These deficits are addressed through the model we present herein.

Various groups have conducted numerical simulations of the Navier-Stokes equations coupled to an immersed body’s dynamics, and have thus studied flapping wings Streitlien and Triantafyllou (1995); Anderson et al. (1998); Alben et al. (2012); Moore (2017); Akhtar et al. (2007) and deformable bodies with more realistic fish-like kinematics Zhu et al. (2014); Neveln et al. (2014); Tytell et al. (2010); Borazjani and Sotiropoulos (2008); Liu et al. (2017); Gazzola et al. (2011); Hieber and Koumoutsakos (2008); Maertens et al. (2017). Hemelrijk et al. Hemelrijk et al. (2014) and Daghooghi et al. Daghooghi and Borazjani (2015) modeled a fish school by numerically simulating a swimmer with doubly-periodic boundary conditions in 2D and 3D, respectively, and found that swimmers move faster in schools than in isolation. However, neither study examined the dependence of the speed on the streamwise distance between swimmers. While these studies allowed for the complex flow structures around flapping bodies to be quantitatively studied, simulations of fish schools are computationally challenging because of the large Reynolds number of the associated flow and the number of interacting bodies, prohibiting a detailed parametric study of lattice formations.

Understanding the role of hydrodynamic interactions in fish schools may thus benefit from a simplified physical system amenable to theoretical analysis. An example is the recent experimental work of Becker et al. Becker et al. (2015), who realized an in-line formation of swimmers using freely-translating, periodically heaving wings in a cylindrical water tank. They observed that the system exhibited a bistability of “schooling states”  and spontaneously locked into either a slow mode or fast mode, the latter of which exhibited a significant speed increase relative to an isolated swimmer. The experiments were extended by Ramananarivo et al. Ramananarivo et al. (2016) to allow tandem swimmers to dynamically select both their speeds and relative positions.

We present here a modeling framework for understanding the hydrodynamic interactions among flapping swimmers at high Reynolds number. Our conceptually simple model is rich enough to incorporate the essential features of the swimmers’ hydrodynamic interactions, while allowing for analytical determination of the model’s exact solutions and their stability properties. Crucially, our model differs from other self-propelled particle models in that the swimmers’ shed vortices are accounted for, so the thrust on a swimmer depends on the system’s history. The swimmers thus interact through a fluid-mediated memory, stored through the collective shed vorticity, and we self-consistently solve for the formation’s emergent speed. The model yields insight into experiments on interacting wings Becker et al. (2015) and predicts instabilities that can be explored in the laboratory. We also extend our model to 2D and determine the lattice formations that allow for the greatest speed and efficiency, thus addressing the questions first posed by Weihs Weihs (1973, 1975). Specifically, we obtain the following hierarchy: while a phalanx formation (§IV.2) offers a small improvement over a solitary swimmer, 1D in-line (§IV.1) and 2D rectangular lattice formations (§IV.3) exhibit substantial improvements, with the 2D diamond lattice (§IV.4) offering the largest hydrodynamic benefit.

II Simple model of an in-line formation

We first consider the simplest model for a school: an infinite line of swimmers modeled as heaving rigid wings driven periodically with prescribed vertical position $y(t)=h_{0}\sin(2\pi ft)$, and coupled by nearest-neighbor interactions, as shown in Fig. 1a. We assume the swimmers to be separated by a fixed distance $L$, so the only degree of freedom is the formation’s speed $U$. Our goal is to construct the governing equations for this system, and determine the dependence of $U$ on the kinematic parameters $h_{0}$ and $f$.

The Reynolds number of the flow around the wings used in experiments Becker et al. (2015) is typically large, $\text{Re}=Uc/\nu\approx 10^{2}-10^{5}$, where $c$ is the chord length and $\nu$ the fluid’s kinematic viscosity. Such a flow is complex and difficult to quantitatively characterize, both experimentally and numerically. We thus make the simplifying assumption that the flow is two-dimensional, inviscid and incompressible, and that the flow structures may be approximated by point vortices shed from the swimmers’ trailing edges at the extrema of their trajectories (Fig. 1b). These vortices mediate the interactions between swimmers. We fix the vortex strength as $\gamma=(C_{v}/2\pi)\int_{0}^{T}\dot{y}(t)^{2}\,\mathrm{d}t=C_{v}\pi h_{0}^{2}f/2$, where $C_{v}$ is a free parameter Schnipper et al. (2009); Buchholz et al. (2011). We also assume that the vortex strength decays exponentially over a timescale $\tau$, an assumption that accounts for turbulent breakdown of vortical structures at high Reynolds number Daghooghi and Borazjani (2015); Ramananarivo et al. (2016); Higdon and Corrsin (1978). Both of these assumptions are discussed in detail in Supplemental Material §II B. We are primarily concerned with trajectories for which the swimming speed is much larger than the characteristic advection speed of vortices Milne-Thomson (1968), $U\gg\gamma/\lambda$, $\lambda$ being the distance between vortices of the same sign, and thus assume the vortices to be stationary.

We now construct evolution equations for the swimmers’ position and velocity, with the relevant variables listed in Supplemental Material Table I. Since the swimmers interact through vortices shed by their nearest neighbors, and the distance between them is fixed, the formation’s trajectory is determined by that of a single swimmer. Instead of modeling the swimmer’s continuous-time motion, we evolve its horizontal position $x_{n}$ and velocity $u_{n}$ on the midplane $y=0$ at the discrete times $t_{n}=nT$, where $T=1/(2f)$ is the flapping half-period. At each time step, the swimmer sheds a point vortex of positive (negative) circulation at the peak (trough) of its trajectory for odd (even) $n$, which generates the characteristic reverse von Kármán wake Triantafyllou et al. (2000) of a self-propelling swimmer (Fig. 1b, Supplemental Movie 1). The swimmer moves under the influence of two forces: a drag $F_{D}(u)$, and a propulsive thrust $F_{x}[x_{n},u_{n},v_{n},\omega_{n}(z)]$, where $\omega_{n}(z)$ is the fluid vorticity in the complex plane $z=x+iy$ and $v_{n}=(-1)^{n}2\pi fh_{0}$ is the swimmer’s vertical velocity. The equations of motion are thus

[TABLE]

the trailing edge of a swimmer centered at the origin is at $z=2a$, $\hat{x}_{n+1}=(x_{n+1}+x_{n})/2$ and $m_{\text{e}}$ is the swimmer’s effective mass per unit span, defined in Supplemental Material §II. We impose the boundary-layer drag law $F_{D}(u)=-C_{D}\rho\sqrt{c\nu}|u|^{1/2}u$, where $\rho$ is the fluid density and $C_{D}$ the drag coefficient. Crucially, the propulsive force $F_{x}$ depends on the swimmers’ dynamically generated vorticity field $\omega_{n}(z)$.

We compute the propulsive force $F_{x}$ using the method detailed in Supplemental Material §I. In summary, we model the swimmers as up-down symmetric wings in the complex plane. Such a wing is represented through the action of the so-called Joukowski map $J(\zeta)=\zeta+\zeta_{c}+a^{2}/(\zeta+\zeta_{c})$ on a circle of radius $r=|a-\zeta_{c}|$ in the $\zeta$-plane, where $\zeta_{c}\in\mathbb{R}$ sets the swimmers’ vertical thickness $d_{\text{f}}=\max_{|\zeta|=r}[\text{Im}(J(\zeta))]$ and $a\approx c/4$ roughly sets the chord length. The central vortex strength $\gamma_{c}$ is obtained by imposing the Kutta condition at the swimmer’s sharp trailing edge $z=2a$, which ensures that the flow remains smooth there. To make the calculation tractable, we assume that the swimmers influence each other only through their shed vortices and not through their body dynamics, an assumption we expect to be valid in the regime $L\gg c$. We thus obtain an explicit form for the iterated map (1) that evolves the swimmer’s position and velocity, given in Supplemental Material Eq. (12).

We note that a similar method was used by Ramananarivo et al. Ramananarivo et al. (2016) to calculate the hydrodynamic force on a wing due to shed point vortices. However, that work neglected consideration of the wings’ dynamics, which is explicitly accounted for by our iterated map (1). For this reason, our work goes beyond theirs in two significant ways. First, we are able to assess the stability of steady schooling states, which yields important information about which states are observed experimentally Becker et al. (2015) and thus which lattice configurations are hydrodynamically optimal. Second, we are able to capture time-dependent schooling states, which we find to emerge naturally in 1D flocks (§III). Moreover, while the prior work was limited to a 1D configuration, modeling 2D diamond lattices requires consideration of time-dependent schooling states, as shown in §IV.4.

Our theoretical model for the flow field generated by a flapping wing exhibits satisfactory agreement with Becker et al.’s measurements Becker et al. (2015) using particle image velocimetry, with upstrokes and downstrokes producing upward and downward fluid flows (Supplemental Fig. 1). Our model thus exhibits qualitative similarity with their ad hoc 1D kinematic model (i.e. swimmers’ inertia is neglected and forces directly determine the swimming speed), which posits an interaction force between swimmers that oscillates sinusoidally on the flapping period and also decays in time. However, our model is derived from a detailed physical description of vortex-induced fluid forces, thus permitting quantitative comparison to experimental data (§III) and treatment of 2D formations (§IV). We also explicitly incorporate the swimmers’ inertia, allowing for an assessment of the stability of steady schooling states (§III).

Our considerable simplification of the flow structures allows the formation’s dynamics to be analyzed mathematically. Moreover, simulation of the governing equations is inexpensive, which allows us to assess the dependence of the dynamics on the kinematic parameters $h_{0}$ and $f$. As shown in Supplemental Material §II, the governing equations can be written in the dimensionless form

[TABLE]

where $G$ is specified in Supplemental Material Eq. (17).

III Comparison with experiment

We seek steadily translating solutions to Eq. (2) and assess their stability, first with the goal of rationalizing the experimental observations of Becker et al. Becker et al. (2015). The analysis of Eq. (2) is nontrivial due to the temporal nonlocality of the formation’s dynamics: updating the velocity $u_{n}$ requires knowledge of the swimmer’s history, which is a generic feature of flow-induced interactions at high Reynolds number.

Substituting the steady state $x_{n}=Un,u_{n}=U$ into Eq. (2), we find that $U$ satisfies the algebraic equation

[TABLE]

This equation is solved numerically using a bisection method. The linear stability analysis of such steady state solutions is given in Supplemental Material §III. In summary, we linearize Eq. (2) around the steady state solution, and find the eigenvalues of the linear stability problem using the discrete Laplace transform. We show that the eigenvalues are given by the zeros of the function

[TABLE]

where the constants $g_{x},g_{u}$ and function $\mathcal{G}(z)$ are specified in Supplemental Material §III. The steady state solution is stable if all of the roots of $\mathcal{F}(z)$ lie inside the unit disc in the complex plane, $|z|<1$, and is unstable otherwise. We use a numerical contour integration method Delves and Lyness (1966) to find the roots of $\mathcal{F}(z)$ in the annular region $1<z<R$, where $R$ is a sufficiently large number. The stability properties of the steady state solution $u_{n}=U$ are dictated by the location of the root $z^{*}$ of $\mathcal{F}(z)$ that is largest in magnitude.

Figure 2 shows the comparison between theory (curves) and experiment Becker et al. (2015) (triangles) for an in-line formation. In the experiments of Becker et al. Becker et al. (2015), the distance $L$ between swimmers is fixed, while the flapping frequency $f$ and amplitude $h_{0}$ are varied. The measured swimming speed $|U|$ is shown in Fig. 2a. The data in Fig. 2b are plotted in terms of the schooling number $S=Lf/|U|$, which denotes the number of wavelengths separating the swimmers Becker et al. (2015). That is, integer values of $S$ indicate trajectories for which the swimmers traverse identical paths, and half-integer values indicate trajectories for which neighboring swimmers traverse paths that are perfectly out-of-phase. The steady state solutions predicted by Eq. (3) are color-coded according to their stability: specifically, blue denotes stable states; red denotes unstable states for which $\text{Re}(z^{*})>0$ and $\text{Im}(z^{*})=0$; and green denotes oscillatory states, which may destabilize via either a flip bifurcation ($\text{Re}(z^{*})<0$ and $\text{Im}(z^{*})=0$) or a Neimark-Sacker bifurcation ($\text{Im}(z^{*})\neq 0$). The physical significance of these instabilities is explained at the end of this section.

The theoretical predictions in Fig. 2 exhibit generally excellent agreement with the experimental data of Becker et al. Becker et al. (2015). At the lowest flapping amplitude considered, $h_{0}=0.75$ cm, the agreement between theory and experiment is less good, presumably owing to the breakdown of the point-vortex approximation at low flapping amplitudes. We note that three free parameters, namely, $C_{D}$ (drag coefficient), $\tau$ (vortex time decay) and $C_{v}$ (initial vortex strength), are chosen once to best fit all of the experimental data across flapping frequency $f$ and amplitude $h_{0}$. The numerical values of these parameters exhibit good agreement with results in the existing literature, as detailed in Supplemental Material §II B. Indeed, our fit value $\tau\approx 2$ seconds compares well with the value $\tau\approx 5$ seconds inferred from the data of Newbolt et al. (Newbolt et al. (2019), Supporting Information), who measured the temporal decay of the flow generated by a flapping wing in a water tank. Generally, the correspondence between theory and experiment makes clear that our simplified approach to modeling the flow and the swimmers’ dynamics captures the key features of the hydrodynamic interactions between swimmers at high Reynolds number.

The experiments of Becker et al. Becker et al. (2015) provided evidence for the bistability of steady states, and the emergence of coexisting “slow modes”  and “fast modes”  for the same flapping frequency $f$ and amplitude $h_{0}$. Our theoretical predictions explain these observations in terms of the stability properties of the steady state solutions. Indeed, branches of stable (blue) solutions are separated by unstable (red) branches, as shown in Fig. 2. Specifically, schooling numbers $S\approx s+1/4$ for $s\in\mathbb{N}$ are favored by the system when the wings are moving fast due to their large flapping frequency, a regime in which hydrodynamic interactions are the strongest. The unstable branches typically run through schooling numbers $S\approx s+3/4$, which are thus avoided by the system. The emergence of these schooling numbers is explained analytically in Supplemental Material §II A. We note that the oscillatory structure of the predicted solutions is a consequence of the hydrodynamic coupling between swimmers, as the speed is a monotonically increasing function of the flapping frequency $f$ for an isolated swimmer.

In addition to explaining the bistability of steady states, the linear stability analysis in Fig. 2 hints that in-line formations may show unsteady behavior, which is confirmed by numerical simulations of Eq. (2). Specifically, Supplemental Fig. 2a shows a simulation conducted in a parameter regime in which the steady state solution $u_{n}=U$ goes unstable via a flip bifurcation, so the swimmer’s velocity oscillates on the flapping period. In experiments, one would thus expect the swimming speed to change appreciably during a single flap, but to be roughly the same at the start of each flap. More interestingly, Supplemental Fig. 2c shows that a steady schooling state may also undergo a Neimark-Sacker bifurcation, in which the swimming speed oscillates over a period long relative to the flapping period, $T_{\text{osc}}\approx\pi\left[\text{Im}\left(\log z^{*}\right)\right]^{-1}\approx 6T$ for the parameter regime explored here. Due to its simplicity, our model furnishes testable predictions for the manner in which steady schooling states destabilize, and the parameter regimes in which they do so.

IV Comparative analysis of different formations

Having benchmarked our theoretical model against experimental data, we analyze how hydrodynamic interactions impact the performance of different lattice formations. We consider two performance measures: the formation’s speed $U$ and cost of transport $C$, and compare these with the corresponding speed $U_{0}$ and cost of transport $C_{0}$ of a single isolated swimmer. The cost of transport $C=\langle P_{n}/|u_{n}|\rangle$ is a “gallons-per-mile”  measure of efficiency that quantifies the formation’s energy consumption per unit distance Tokić and Yue (2012), where $P_{n}$ is the instantaneous mechanical power output of the swimmer at time $t_{n}=nT$ and $\langle\cdot\rangle$ denotes a time average. A formula for $P_{n}$ is given in Supplemental Material §IV.

We compute the speed $U$ and cost of transport $C$ of in-line (§IV.1), phalanx (§IV.2), rectangular lattice (§IV.3), and diamond lattice (§IV.4) formations as a function of the distance between swimmers. Our goal is to identify the lattice formations that maximize the speed and minimize the cost of transport relative to that of a single swimmer, values $U/U_{0}>1$ and $C/C_{0}<1$ indicating a benefit due to collective hydrodynamic interactions. In this section, we restrict our attention to a single representative set of flapping kinematics, $f=1.5$ Hz and $h_{0}=1.5$ cm, for which $h_{0}/c=0.25$ and $\text{St}\equiv 2h_{0}f/U_{0}\approx 0.1$, the low Strouhal number regime $\text{St}\ll 1$ being biologically relevant for fish schools Triantafyllou et al. (2000). All distances are reported in units of the swimmer’s approximate body length $4a$.

IV.1 In-line formation

We first solve Eq. (3) to find the dependence of the swimming speed $U$ on the streamwise distance $L$ between swimmers in a line. As in §III, we assume the formation’s dynamics to be dominated by nearest-neighbor interactions in the streamwise direction. The results are shown in Fig. 3. As expected from the discussion in §III, a slow mode ($U/U_{0}<1$) and fast mode ($U/U_{0}>1$) may coexist for a given value of $L$ (Fig. 3a). The maximum speedup of 17% corresponds to a state with $S=0.73$, while the largest slowdown by 19% corresponds to a state with $S=0.96$ (Fig. 3b). The lowest cost of transport is $C/C_{0}=0.75$ and corresponds to a state with $S=1.26$, indicating a maximum energy savings of 25%, while the highest cost of transport $C/C_{0}=1.42$ corresponds to a state with $S=0.75$ (Fig. 3c). Plotting $U/U_{0}$ as a function of $S$ (Fig. 3b) shows that states with $S\gtrsim s+1/4$ and $S\lesssim s+3/4$ typically have the highest speeds, whereas those with $S\lesssim s$ and $S\lesssim s+1/4$ have the lowest. Conversely, Fig. 3c shows that states with $S\lesssim s+3/4$ typically have the highest cost of transport, whereas those with $S\lesssim s$ and $S\lesssim s+1/4$ have the lowest. Comparing Fig. 3b and 3c, we observe that high-speed states ($U/U_{0}>1$) are typically associated with an increased cost of transport ($C/C_{0}>1$), indicating a tradeoff between speed and energy consumption.

To understand the oscillatory dependence of $U/U_{0}$ on $L$, we derive an approximate form for the thrust $F^{\text{v}}_{x}$ and lift $F_{y}^{\text{v}}$ on a swimmer due to the vortices shed by its neighbors, assuming that the vortices (with positions $z_{k}\in\mathbb{C}$ and strengths $\gamma_{k}$) are far from the body, $|z_{k}|\gg r$. In Supplemental Material §I, we show that

[TABLE]

is the vertical velocity of the flow induced by the neighboring swimmers’ shed vortices and $(U,V)$ is the swimmer’s velocity. Figure 3d shows a schematic of the reverse von Kármán (thrust) wake and associated fluid flow generated by a flapping swimmer. A swimmer moving upward ($V>0$) in its neighbor’s wake would experience a constructive interaction force $F^{\text{v}}_{x}<0$ and thus a speedup ($U/U_{0}>1$) for $s+1/4\lesssim S\lesssim s+3/4$, and the opposite for $s+3/4\lesssim S\lesssim s+5/4$, according to Eq. (5). We note that these conclusions may be justified mathematically by adapting the argument in Supplemental Material §II A, in which we derive an approximate form for the interaction force (5) in the biologically-relevant low Strouhal number limit $\text{St}\ll 1$ Triantafyllou et al. (2000).

A similar argument may be used to qualitatively understand the tradeoff between speed and cost of transport. A wing moving up ($V>0$) in a high-speed state ($U/U_{0}>1$) experiences a downward flow from its neighbor’s wake ($V_{\text{f}}<0$), which implies $F_{y}^{\text{v}}<0$ by Eq. (5). High-speed states are thus typically associated with an increased power consumption, as the wing’s vertical motion is opposed by its neighbor’s induced flow. In the low Strouhal number limit $\text{St}\ll 1$, the increase in cost of transport $C$ due to the increased power consumption dominates the decrease in $C$ associated with the higher speed, as shown in Supplemental Material §IV A.

We may use the foregoing arguments to draw conclusions about an in-line formation in which the nearest neighbors flap perfectly out-of-phase with respect to each other. For such a configuration, the right-hand sides of the formulae in Eq. (5) for $F_{x}^{\text{v}}$ and $F_{y}^{\text{v}}$ would simply have their signs reversed. We thus expect a speedup and larger cost of transport for schooling numbers $s+3/4\lesssim S\lesssim s+5/4$, and a slowdown with lower cost of transport for $s+1/4\lesssim S\lesssim s+3/4$. That is, Figs. 3b and 3c would be qualitatively unchanged, apart from a shift of the horizontal axis, $S\rightarrow S+1/2$.

IV.2 Phalanx

Motivated by the experimental observations of Ashraf et al. Ashraf et al. (2017), we now consider a phalanx of swimmers: infinitely many swimmers equally spaced by a distance $d$ in the lateral ($y$) direction and, following Weihs Weihs (1975), flapping in antiphase with respect to their neighbors, as shown in Fig. 4 and Supplemental Movie 2. As discussed in Supplemental Material §V A, Weihs Weihs (1975) and Stöcker Stöcker (1999) argued that in-phase flapping would result in an increased drag force on the downstream fish, due to the anomalously large induced velocity in the $y$-direction. Such an argument is in agreement with the experimental observations of Ashraf et al. Ashraf et al. (2016) who found that pairs and triplets of red nose tetra fish preferentially flap in antiphase with respect to their lateral neighbors. We also restrict our attention to the parameter regime $d>d_{\text{min}}\equiv 2(h_{0}+d_{\text{f}})$ to ensure that the swimmers do not collide with each other during a flapping cycle.

A straightforward extension of the iterated map presented in §II permits consideration of this formation, as detailed in Supplemental Material §V. Following the procedure detailed in §III, we find that the swimming speed $U$ satisfies an algebraic equation of the form (3) with $L=0$. The interaction function $G$ describes the hydrodynamic thrust due to a side-by-side arrangement of swimmers flapping in antiphase, and is defined in Supplemental Material Eq. (41). The solutions to this equation are shown in Fig. 4. The phalanx formation evidently does not exhibit the multi-stability of steady states seen for in-line formations (§IV.1). By generalizing the linear stability analysis of steady state solutions presented in §III, we find that the steady state is stable for all values of $d$. As shown in Fig. 4, the formation exhibits a slight speedup for all values of $d$, with the maximum speedup of roughly 5% occurring when the swimmers are most tightly packed, $d=d_{\text{min}}$. However, such a formation also increases the cost of transport by roughly 4%, with the cost of transport decreasing as $d\rightarrow\infty$. Similar to the 1D formations discussed in §IV.1, the phalanx formation exhibits a tradeoff between speed and cost of transport, although both measures exhibit relatively small variations over a range of values $d$.

IV.3 Rectangular lattice

We now consider the rectangular lattice of swimmers shown in Fig. 5a and Supplemental Movie 3: swimmers separated by a streamwise distance $L$ and a vertical distance $d$, starting at the positions $z=jL+ikd$ for $j,k\in\mathbb{Z}$. Swimmers flap in phase (antiphase) with respect to their streamwise (lateral) neighbors. As in §II, a swimmer at the origin (red box in Fig. 5a) interacts with the swimmers in the neighboring columns at $x=0,\pm L$. Note that the lattice is effectively an in-line formation in the limit $d\rightarrow\infty$. In Supplemental Material §V, we show that the steady speed $U$ satisfies the algebraic equation (3), with the function $G$ defined in Supplemental Material Eq. (41).

We numerically solve this equation to find the dependence of the steady speed $U$ on the geometric parameters $L$ and $d$. We then perform numerical simulations of the evolution equations, with initial conditions determined by these steady-state solutions, and compute the trajectory’s time-averaged velocity $U_{\text{av}}$. Figure 5b shows the normalized velocity $U_{\text{av}}/U_{0}$ as a function of $L$ and $d$. As with the in-line formation, there may be multiple steady-state solutions for a given set of parameters; in such cases, we perform multiple simulations and plot the time-averaged speed of the fastest state. The simulations reveal the existence of multiple coexisting states within the regions of geometric parameter space bounded by the gray curves in Fig. 5b. As the lateral distance $d$ between swimmers is decreased progressively, these regions of multi-stability typically shrink, but new regions of multi-stability may also emerge.

We find that the formation experiences a maximum speedup of 18% for a roughly square geometry, $L=2.1$ ($S=0.33$) and $d=1.94$. The black curves in Fig. 5 indicate $d_{U}(L)=\arg\max_{d}(U_{\text{av}}(L,d)/U_{0})$, the optimal lateral spacing for a given streamwise spacing, and the corresponding speed $U_{\text{max}}(L)=|U_{\text{av}}(L,d_{U}(L))|$. Note that, unlike the phalanx (Fig. 4), the greatest speedup is not necessarily achieved by packing the swimmers tightly in the $y$-direction, so $d_{U}(L)$ is not identically equal to $d_{\text{min}}$. A comparison between the rectangular lattice and in-line formation is shown in Fig. 5d. The most salient feature is that the slow modes ($U/U_{0}<1$) for an in-line formation all exhibit speedup in the corresponding rectangular lattice with the minimum lateral spacing, $d_{U}(L)=d_{\text{min}}$. However, the fast modes ($U/U_{0}>1$) for an in-line formation benefit minimally from the rectangular geometry, and the corresponding optimal lateral spacing $d_{U}(L)$ is often much larger than $d_{\text{min}}$. Indeed, the fastest rectangular lattice formation is only marginally faster than the fastest in-line formation.

Figure 6 shows the corresponding results for the cost of transport of a rectangular lattice of swimmers. For the regions of parameter space in which multiple states coexist, Fig. 6b shows the cost of transport of the most efficient one. The cost of transport assumes the minimum value for an effectively in-line formation ($d\rightarrow\infty$), for which $C/C_{0}=0.75$ and $S=1.26$ (Fig. 3c). The black curves in Fig. 6 correspond to $d_{C}(L)=\arg\min_{d}C(L,d)$, the optimal lateral spacing for a given streamwise spacing, and the corresponding cost of transport $C_{\text{min}}(L)=C(L,d_{C}(L))$. Unlike $d_{U}(L)$, which exhibits a nontrivial dependence on the streamwise spacing $L$ (Fig. 5b), $d_{C}(L)$ typically assumes the values $d_{C}\approx d_{\text{min}}$ and $d_{C}=\infty$ (Fig. 6b). Figure 6d shows that the states for which $d_{C}\approx d_{\text{min}}$ typically correspond to inefficient states ($C/C_{0}>1$) for which the rectangular lattice formation affords a slight hydrodynamic advantage over the in-line formation. Conversely, the states for which $d_{C}=\infty$ typically correspond to efficient states ($C/C_{0}<1$) for which the lattice geometry actually increases the cost of transport, making the in-line formation the most efficient.

Taken together, the results in Fig. 5d and Fig. 6d make clear the dominance of streamwise interactions between swimmers in determining both the speed and cost of transport of a rectangular lattice of swimmers. The physical picture given in §IV.1 may thus be used to qualitatively understand how the lattice geometry influences both the speed and cost of transport. Recall that the in-line formation typically experiences a speedup ($U/U_{0}>1$) and a decrease in efficiency ($C/C_{0}>1$) when a swimmer at $z=0$ swims up into the downflow generated by its upstream neighbor at $z=-L$ (Fig. 3). For such values of $L$, corresponding to schooling number $s+1/4\lesssim S\lesssim s+3/4$, the upstream swimmers at $z=-L\pm id$ in a rectangular lattice contribute an upflow (Fig. 5a), decreasing the thrust but increasing the lift according to Eq. (5). These effects contribute to a decrease in speed but an increase in efficiency for a rectangular lattice as compared to an in-line formation, making it advantageous to increase (decrease) the lateral spacing $d$ when optimizing for speed (efficiency). Conversely, the in-line formation is relatively efficient ($C/C_{0}<1$) but slow ($U/U_{0}<1$) for values of $L$ corresponding to $s+3/4\lesssim S\lesssim s+5/4$, the parameter regime in which the upstream swimmers at $z=-L\pm id$ in the rectangular lattice generate a downflow. This contributes to an increase in speed but a decrease in efficiency, making it advantageous to decrease (increase) the lateral spacing $d$ when optimizing for speed (efficiency).

IV.4 Diamond lattice

We now consider the diamond lattice shown in Fig. 7a and Supplemental Movie 4: swimmers separated by a streamwise distance $L$ and lateral distance $d$, starting at the positions $z=(j+l/2)L+i(k+l/2)d$ for $j,k\in\mathbb{Z}$ and $l\in\{0,1\}$. As in the previous section, swimmers flap in phase (antiphase) with respect to their streamwise (lateral) neighbors, and a swimmer at the origin (red box in Fig. 7a) interacts with the swimmers in the neighboring columns at $x=0,\pm L/2,\pm L$.

The complete equations are detailed in Supplemental Material Eq. (45). Since the lattice is no longer symmetric about the midplane $y=0$, the steady state $u_{n}=U$ is not a solution to the governing equations, but the period-2 trajectory $u_{n}=U+(-1)^{n}U_{1}$ is. The unknowns $U$ and $U_{1}$ satisfy the pair of algebraic equations

[TABLE]

where $\alpha_{j}=\text{mod}(j,2)/2$ accounts for the vertical shift in the columns at $x=\pm L/2$. This equation describes the balance of drag and vorticity-induced thrust on a swimmer, both of which depend on the swimmer’s instantaneous flapping phase. We describe how to solve these equations numerically, and extend the methodology described in §III to assess the stability of these period–2 solutions, in Supplemental Material §V and §VI, respectively.

We then numerically simulate the governing equations for a diamond lattice with initial conditions determined by the period–2 base states. As with the rectangular lattice described in §IV.3, the simulations reveal the existence of multiple coexisting states, indicated by the regions of parameter space bounded by the gray curves in Fig. 7b. In such cases, Fig. 7b shows the normalized average speedup $U_{\text{av}}/U_{0}$ of the fastest state. We find that the formation experiences a maximum speedup of 22% for $L=4.4$ and $d=d_{\text{min}}$, which corresponds to $S=0.67$. As with the rectangular lattice considered in §IV.3, $d_{U}(L)$ is not identically equal to $d_{\text{min}}$, although the lateral spacing is indeed minimized for the fastest lattice.

These results may also be qualitatively understood using a simple argument based on fluid flows, and by comparing the diamond lattice and in-line formation in Fig. 7d. As with the rectangular lattice, the slow modes ($U/U_{0}<1$) for an in-line formation exhibit speedup in a diamond lattice with the minimum lateral spacing, $d_{U}(L)=d_{\text{min}}$. Most fast modes for an in-line formation ($U/U_{0}>1$) benefit minimally from the diamond geometry; however, those with $S\lesssim s+3/4$ experience a speedup of 3–5% relative to the in-line formation, also with $d=d_{\text{min}}$. This is due to the existence of an entirely new branch of fast modes that emerges for $d\gtrsim d_{\text{min}}$, as shown in Supplemental Fig. 3. These fast modes may be attributed to the beneficial dipolar structure generated by the swimmer’s upstream neighbors at $z=-L/2\pm id/2$, as shown in Fig. 7a. The fluid flow induced by this dipole at the midplane $y=0$ does not have a vertical component, and its horizontal component is negative, effectively imparting additional thrust to the swimmer and increasing its speed. The next flap generates a dipole of opposite sign, which imparts a drag; however, this contribution is weaker than the previous thrust contribution, since the associated vortices are farther from the body.

Figure 8 shows the corresponding results for the cost of transport of a diamond lattice of swimmers. The cost of transport assumes the minimum value $C/C_{0}=0.67$ for a state with $S=0.5$ and $d=d_{\text{min}}$, indicating that a tightly packed diamond lattice formation realizes a substantial increase in efficiency due to hydrodynamic interactions. Figure 8b shows that, for a given streamwise spacing $L$, the optimal lateral spacing is typically $d_{C}(L)\approx d_{\text{min}}$ or $d_{C}(L)=\infty$. This effect is similar to that observed for the rectangular lattice, and may be rationalized by the physical argument presented in §IV.3. However, the diamond lattices for which $S\approx s+1/2$ are noticeably more efficient than their in-line formation counterparts, an effect that is absent for rectangular lattices. This finding may also be explained by the beneficial dipolar structure generated by the swimmer’s upstream neighbors at $z=-L/2\pm id/2$ (Fig. 8a), which provides an additional lift force as the swimmer accelerates upward.

While the distinct optimal diamond lattice formations shown in Fig. 7a and Fig. 8a correspond to stable period-2 states, closer examination of the simulated solutions reveals that they often exhibit a complex nonlinear dynamics. This finding indicates that hydrodynamic interactions may significantly influence schooling behavior when the swimmers are interacting strongly, particularly in the regime $d\gtrsim d_{\text{min}}$. We note that we also observed chaotic solutions in our model for the in-line formation, and leave the complete characterization of the system’s nonlinear dynamics for future work.

V Discussion

We have presented a new model for the hydrodynamic interactions between swimmers in high-Reynolds number flows. While numerical simulations of such systems are computationally challenging, our model offers a simple framework by which to interpret the formation’s dynamics: the swimmers shed vortices during each stroke, and in turn are propelled due to the vorticity-induced flow field. Despite neglecting consideration of the details of the flapping kinematics and the associated flow structures, our model exhibits good agreement with experimental data on interacting wings Becker et al. (2015), while using only three fitting parameters ($C_{D}$, $\tau$ and $C_{v}$). As shown in Fig. 2, the observed bistability of slow and fast states is a consequence of overlapping branches of stable (blue) steady-state solutions, which are separated by unstable (red) branches. The multi-stability of states has not been observed in prior theoretical investigations, but is a generic feature of our model: indeed, multiple coexisting states are found in the in-line, rectangular and diamond lattice formations for appropriate values of the geometric parameters. Our model also predicts new schooling instabilities (Supplemental Fig. 2) through which the speed oscillates in time, an effect that can be probed experimentally. Animal schools might employ active control mechanisms in order to mitigate the effects of such instabilities.

Despite the apparent complexity of the hydrodynamic interactions, we show that the interaction thrust force between two swimmers is approximately proportional to the vertical velocity of a swimmer and the vertical velocity of the flow induced by its neighbor (Eq. (5)). This shows that the speed of an in-line formation is sensitive to the streamwise spacing, as a self-propelling flapping swimmer generates a spatially oscillatory flow field in its wake (Fig. 3d). For the set of flapping kinematics considered, we find that an in-line formation may move up to 17% faster than an isolated swimmer, provided the distance between swimmers $L$ is such that the schooling number $S\lesssim 0.75$ (Fig. 3b). By contrast, a phalanx formation of swimmers flapping in antiphase moves roughly 5% faster than a single swimmer (Fig. 4), with the optimal speedup occurring when the lateral distance between swimmers is minimized. While the fastest rectangular lattice provides a marginal advantage over an in-line formation (Fig. 5), the fastest diamond lattice is able to move 22% faster than a single swimmer (Fig. 7), provided that the lateral distance between swimmers is minimized. This effect may be attributed to the advantageous horizontal flow generated by a swimmer’s upstream neighbors (Fig. 7a).

By using the cost of transport to measure the energetic efficiency, we find that in-line formations may realize an energy savings of 25% over an isolated swimmer, provided the distance between swimmers $L$ is such that the schooling number $S\gtrsim 1.25$ (Fig. 3c). While our finding that the phalanx formation affords a speedup is in agreement with the experimental observations of Ashraf et al. Ashraf et al. (2017), we also find such formations to be slightly less efficient than an isolated swimmer (Fig. 4). Generally, we observe that in-line and phalanx formations exhibit a tradeoff between speed and efficiency. While all rectangular lattices have a higher cost of transport than the most efficient in-line formation (Fig. 6), the most efficient diamond lattice has an energy savings of 33% over an isolated swimmer, an effect that can also be attributed to the effective vortex dipole generated by the swimmer’s upstream neighbors (Fig. 8a). Taken together, our results show that the fastest and most efficient diamond lattice formations are distinct, and that they outperform the other geometries in speed and efficiency, respectively. Based on the argument presented at the end of §IV.1, we expect qualitatively similar results for a model in which the swimmers flap perfectly out-of-phase with respect to their streamwise neighbors in both in-line and 2D lattice formations. While most of the relevant schooling numbers will be shifted, $S\rightarrow S+1/2$, we still expect the diamond lattices with $S=s+1/2$ to be particularly efficient (Fig. 8d), owing to the vortex dipole shed by a swimmer’s upstream neighbors.

Our finding that the most efficient state is realized by a diamond lattice with the minimal lateral spacing $d=d_{\text{min}}$ is consistent with the results of Weihs Weihs (1973, 1975). However, our results differ in that the curves $d_{U}(L)$ and $d_{C}(L)$ (black curves in Fig. 7d and Fig. 8d, respectively) are not monotonic, implying that there is not an optimal lattice angle. While Weihs argued that it is disadvantageous for a swimmer to swim directly behind another, we find that this is not necessarily the case (Fig. 3). We note that our model differs from Weihs’ in some important ways. First, we model flapping swimmers, and find that the formation’s speed is influenced by the relationship between the swimmer’s kinematics and oncoming fluid flow (§IV.1). More sophisticated computational models Maertens et al. (2017); Verma et al. (2018) of fish-like swimming have also highlighted the importance of this relationship. Second, we allow the vortex strength to decay in our model, which sets an effective interaction distance $\sim U\tau$ between swimmers. Third, we explicitly account for the formation’s dynamics and thus the schooling modes’ stability, which was not possible in Weihs’ framework.

Our results may also be compared with recent theoretical and computational studies of a swimmer in a doubly periodic domain, which generates a lattice of swimmers flapping in phase. Recently, Hemelrijk et al. Hemelrijk et al. (2014) conducted 2D numerical simulations using a multi-particle collision dynamics model and also found that the diamond lattice provides the largest speedup relative to a single swimmer ($23\%$), but that the maximum occurred for $d=2$, the largest lateral spacing considered. Similarly, the diamond lattice for $d=1.75$ was found to have the highest Froude efficiency, or the ratio of the useful power to the total power input. However, Hemelrijk et al. Hemelrijk et al. (2014) did not consider the influence of the streamwise spacing $L$, which we find to play a dominant role. Daghooghi & Borazjani Daghooghi and Borazjani (2015) conducted 3D large eddy simulations of rectangular lattices of swimmers at high Reynolds number. They found that the swimming speed and efficiency increase as the lateral distance $d$ decreases, but also did not consider the influence of the streamwise spacing and instead fixed $L=1$. They attributed the observed hydrodynamic advantage to channeling or wake blockage; this effect is beyond the scope of our model, since we assume the vortices to remain in place once shed. Nevertheless, in our simulations of rectangular lattices of swimmers that flap in antiphase with respect to their lateral neighbors, we find that the swimming speed and efficiency do not necessarily exhibit a monotonic dependence on $d$, but instead are nontrivially influenced by both $L$ and $d$ (Fig. 5b and Fig. 6b). Tsang & Kanso Tsang and Kanso (2013) proposed a far-field hydrodynamic model of swimmers as finite-sized vortex dipoles, and found that swimmers in both rectangular and diamond lattices actually move slower than they would in isolation. They attributed this result to the absence of shed vorticity in their model, which has been shown to mediate the near-field interactions between swimmers Ramananarivo et al. (2016). Our dynamical model builds on these studies by explicitly modeling both the flapping kinematics and shedding of vortices. The model’s mathematical simplicity allows us to analytically show the multi-stability of states in in-line, rectangular and diamond lattice formations, a result absent from all of the studies described above. Moreover, the model’s computational tractability allows us to conclusively determine the dependence of the speed and cost of transport on the geometric parameters $L$ and $d$.

We now compare some of our results to schooling formations reported in the literature, despite the current sparsity of quantitative data. In their observations of red nose tetra fish, Ashraf et al. Ashraf et al. (2017, 2016) reported that high-speed schools typically adopt phalanx formations with nearest neighbors separated by 0.5–0.6 body lengths. This observation is consistent with our theoretical prediction that the fastest phalanx formation is realized for the minimum lateral distance considered, $d\approx 0.7$ body lengths (Fig. 4). Similarly, Atlantic bluefin tuna have been observed to adopt phalanx formations with an average lateral spacing between 1 and 1.5 body lengths Newlands and Porcelli (2008); Partridge et al. (1983). When in relatively large schools, this fish species has also been observed to adopt diamond lattice-like formations with a mean first- and second-nearest neighbor separation angle between $14^{\circ}$ and $17^{\circ}$ (Fig. 5 in Partridge et al. (1983)). This observation is roughly consistent with our theoretical prediction that the diamond lattice with the lowest cost of transport has a separation angle of 13*∘* (Fig. 8). However, the red nose tetra fish has been observed to adopt diamond lattice formations with a typical separation angle of approximately $37^{\circ}$ (Fig. 2b in Ashraf et al. (2017)). Moreover, this fish species tends to adopt the phalanx over the diamond lattice formation at high swimming speeds Ashraf et al. (2017), an observation that runs counter to our prediction that the fastest and most efficient states are realized by diamond lattice formations. We note that the spatial phase synchronization between neighboring fish, as measured through the schooling number $S$, is not typically reported, which prevents further quantitative comparison between our theoretical predictions and field observations. Such detailed comparison against biological swimmers would also benefit from more accurate modeling of the swimming kinematics and body shape. Considerations beyond hydrodynamics, such as social cues and predator avoidance, undoubtedly also impact the structure of schools observed in nature.

We expect the results presented herein to be most relevant for understanding fish schools, since we neglected the influence of lift forces on the dynamics. While a recent study of shorebird flocks found no evidence of temporal or spatial phase synchronization between birds Corcoran and Hedrick (2019), ibis flocks were observed to preferentially assume V-formations with median schooling number $S\approx s+1/4$, and in-line formations with $S\approx s+1/2$ (Supplemental Fig. 2 in Portugal et al. (2014)). Extensions of our model might shed light on these phenomena, for which lift generation is an important consideration Higdon and Corrsin (1978). We also note that a conceptually similar iterated map model may readily be applied to 3D flocks and schools; however, new techniques would be required to capture fully 3D dynamics, since the complex-variable techniques used herein cannot simply be extended to calculate the vortex-induced flow and associated forces on the bodies. Our work may also have limited application to “disordered”  schools and cluster flocks, which deviate significantly from ordered lattice formations. The model may be generalized to allow for more general body kinematics, including pitching, turning and adaptive spacing between swimmers Ramananarivo et al. (2016); Newbolt et al. (2019).

While the quantitative results presented in this paper depend on the model parameters used, our reduced modeling framework may be applied to understand temporally long-lived hydrodynamic interactions in active systems. Indeed, models for the pilot-wave dynamics of droplets bouncing on a vibrating fluid bath have a similar mathematical structure Bush (2015), and may be extended to probe the droplet’s complex collective dynamics Lieber et al. (2007); Eddi et al. (2009). Generally, we expect models of the type described herein to be broadly applicable to systems of active particles interacting via their collective histories.

Acknowledgements.

We thank Hassan Masoud, Fang Fang, Joel Newbolt, Sophie Ramananarivo and Stephen Childress for helpful discussions. A. U. O. acknowledges support from the Simons Foundation (Collaboration Grant for Mathematicians, award #587006), L. R. acknowledges support from the NSF (DMS-1847955), and M. J. S. thanks the Lilian and George Lyttle Chair of Applied Mathematics.