Lattices of Hydrodynamically Interacting Flapping Swimmers

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 oneand twodimensional 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-byside 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.


I. INTRODUCTION
The complex collective dynamics of fish schools and bird flocks have long fascinated physicists, biologists, and mathematicians [1,2]. In addition to their biological relevance, they are living examples of active systems [3][4][5][6] 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 [7,8]) and the collective behavior of biological systems at the microscale (e.g., bacterial suspensions [9][10][11][12]), 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 [13] are known to exhibit crystal-like bound states [14,15] as a result of their coupling through surface waves. We present here a modeling framework for the longlived 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 [16,17] focus on the behavioral mechanisms behind flocking, a recent analysis [18] 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 [19], and benefit from external flows by harnessing the energy from vortices [20][21][22]. While some studies have argued against a hydrodynamic function for fish schools [23] and instead focused on the social interactions between fish [24][25][26], a number of observations [27][28][29][30] 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 [31,32]. Examples are the seminal works of Huth and Wissel [33] and Viscek et al. [34], which have been extended to other discrete-time models with more elaborate interaction rules [35]. Far-field hydrodynamic interactions have been recently incorporated into selfpropelled particle models of swimmers subject to phenomenological behavioral rules [36,37], but such models neglect the vorticity-induced forces thought to be relevant for schooling fish [21].
Fish schools and bird flocks can exhibit orderly latticelike formations, although field observations are diverse and sometimes contradictory [1,2]. Prior observations revealed that fish schools may adopt lattice configurations in a statistical sense [38]. 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 [39,40]. 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 [41]. In their studies of red nose tetra fish, Ashraf and coworkers [30,42] 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 [43], 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 [18]. Certain species (e.g., pelicans [44], Canada geese [45], ibises [18]) are thought to benefit from their observed V-formation flight, while it is claimed that others (e.g., pigeons [46], pink-footed geese [47]) do not. Corcoran and Hedrick [48] 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 [16,17].
In an attempt to explain such ordered structures, the seminal papers of Weihs [49,50] 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 selfconsistently 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 [51][52][53][54][55] and deformable bodies with more realistic fishlike kinematics [56][57][58][59][60][61][62][63]. Hemelrijk et al. [64] and Daghooghi et al. [65] 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. [66], 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 a fast mode, the latter of which exhibited a significant speed increase relative to an isolated swimmer. The experiments were extended by Ramananarivo et al. [67] 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 selfconsistently solve for the formation's emergent speed. The model yields insight into experiments on interacting wings [66] 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 [49,50]. Specifically, we obtain the following hierarchy: while a phalanx formation (Sec. IV B) offers a small improvement over a solitary swimmer, 1D in-line (Sec. IVA) and 2D rectangular lattice formations (Sec. IV C) exhibit substantial improvements, with the 2D diamond lattice (Sec. IV D) 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πftÞ and coupled by nearest-neighbor interactions, as shown in Fig. 1(a). 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 [66] is typically large, Re ¼ Uc=ν ≈ 10 2 -10 5 , where c is the chord length and ν 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. 1(b)]. These vortices mediate the interactions between swimmers. We fix the vortex strength as γ ¼ where C v is a free parameter and T ¼ 1=ð2fÞ is the flapping half-period [68,69]. We also assume that the vortex strength decays exponentially over a timescale τ, an assumption that accounts for turbulent breakdown of vortical structures at high Reynolds number [65,67,70]. Both of these assumptions are discussed in detail in Supplemental Material Sec. II B [71]. We are primarily concerned with trajectories for which the swimming speed is much larger than the characteristic advection speed of vortices [72], U ≫ γ=λ, λ 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 [71]. 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. 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 [20] of a selfpropelling swimmer [ Fig. 1(b), and Movie 1 in Supplemental Material [71] ]. 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 ; ω n ðzÞ, where ω n ðzÞ is the fluid vorticity in the complex plane z ¼ x þ iy, and v n ¼ ð−1Þ n 2πfh 0 is the swimmer's vertical velocity. The equations of motion are thus T m e fF D ðu n Þ þ F x ½x n ; u n ; v n ; ω n ðzÞg; the trailing edge of a swimmer centered at the origin is at , and m e is the swimmer's effective mass per unit span, defined in Supplemental Material Sec. II [71]. We impose the boundary-layer drag law (1), presented in the complex plane z ¼ x þ iy. The formation's dynamics is determined by the swimmer at the center, and we track its position x n , horizontal velocity u n , and vertical velocity v n ¼ ð−1Þ n 2πfh 0 at the midplane y ¼ 0. Vortices of positive (negative) circulation AEγ are shed from the swimmers' trailing edges at the peaks (troughs) y ¼ AEh 0 of their trajectories, indicated by the dashed curves.
where ρ is the fluid density and C D the drag coefficient. Crucially, the propulsive force F x depends on the swimmers' dynamically generated vorticity field ω n ðzÞ.
We compute the propulsive force F x using the method detailed in Supplemental Material Sec. I [71]. 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ðζÞ ¼ ζ þ ζ c þ a 2 =ðζ þ ζ c Þ on a circle of radius r ¼ ja − ζ c j in the ζ plane, where ζ c ∈ R sets the swimmers' vertical thickness d f ¼ max jζj¼r fIm½JðζÞg and a ≈ c=4 roughly sets the chord length. The central vortex strength γ 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 ≫ c. We thus obtain an explicit form for the iterated map Eq. (1) that evolves the swimmer's position and velocity, given in Supplemental Material Eq. (12) [71].
We note that a similar method was used by Ramananarivo et al. [67] 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 Eq. (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 [66] and thus which lattice configurations are hydrodynamically optimal. Second, we are able to capture timedependent schooling states, which we find to emerge naturally in 1D flocks (Sec. 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 Sec. IV D.
Our theoretical model for the flow field generated by a flapping wing exhibits satisfactory agreement with the measurements of Becker et al. [66] using particle image velocimetry, with upstrokes and downstrokes producing upward and downward fluid flows (Supplemental Material Fig. 1 [71]). 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 (Sec. III) and treatment of 2D formations (Sec. IV). We also explicitly incorporate the swimmers' inertia, allowing for an assessment of the stability of steady schooling states (Sec. 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 Sec. II [71], the governing equations can be written in the dimensionless form where G is specified in Supplemental Material Eq. (17) [71].

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. [66]. 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 This equation is solved numerically using a bisection method. The linear stability analysis of such steady-state solutions is given in Supplemental Material Sec. III [71]. 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 where the constants g x , g u and function GðzÞ are specified in Supplemental Material Sec. III [71]. The steady-state solution is stable if all of the roots of F ðzÞ lie inside the unit disk in the complex plane, jzj < 1, and is unstable otherwise. We use a numerical contour integration method [73] to find the roots of 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 F ðzÞ that is largest in magnitude. Figure 2 shows the comparison between theory (curves) and experiment [66] (triangles) for an in-line formation.
In the experiments of Becker et al. [66], the distance L between swimmers is fixed, while the flapping frequency f and amplitude h 0 are varied. The measured swimming speed jUj is shown in Fig. 2(a). The data in Fig. 2(b) are plotted in terms of the schooling number S ¼ Lf=jUj, which denotes the number of wavelengths separating the swimmers [66]. 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 Reðz Ã Þ > 0 and Imðz Ã Þ ¼ 0, and green denotes oscillatory states, which may destabilize via either a flip bifurcation [Reðz Ã Þ < 0 and Imðz Ã Þ ¼ 0] or a Neimark-Sacker bifurcation [Imðz Ã Þ ≠ 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. [66]. 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), τ (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 Sec. II B [71]. Indeed, our fit value τ ≈ 2 s compares well with the value τ ≈ 5 s inferred from the data of Newbolt et al. ([74], 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. [66] 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 ≈ s þ 1=4 for s ∈ 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 ≈ s þ 3=4, which are thus avoided by the system. The emergence of these schooling numbers is explained analytically in Supplemental Material Sec. II A [71]. 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. 2(a) 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 [71]. 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. 2(c) 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 osc ≈ π½Imðlog z Ã Þ −1 ≈ 6T for the parameter regime explored here [71]. Because of 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 ¼ hP n =ju n ji is a "gallons-per-mile" measure of efficiency that quantifies the formation's energy consumption per unit distance [75], where P n is the instantaneous mechanical power output of the swimmer at time t n ¼ nT and h·i denotes a time average. A formula for P n is given in Supplemental Material Sec. IV [71].
We compute the speed U and cost of transport C of inline (Sec. IVA), phalanx (Sec. IV B), rectangular lattice (Sec. IV C), and diamond lattice (Sec. IV D) 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 St ≡ 2h 0 f=U 0 ≈ 0.1, the low Strouhal number regime St ≪ 1 being biologically relevant for fish schools [20]. All distances are reported in units of the swimmer's approximate body length 4a.

A. 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 Sec. 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 Sec. III, a slow mode (U=U 0 < 1) and fast mode (U=U 0 > 1) may coexist for a given value of L [ Fig. 3(a)]. The maximum speed-up 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. 3(b)]. 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. 3(c)]. Plotting U=U 0 as a function of S [ Fig. 3(b)] shows that states with S ≳ s þ 1=4 and S ≲ s þ 3=4 typically have the highest speeds, whereas those with S ≲ s and S ≲ s þ 1=4 have the lowest. Conversely, Fig. 3(c) shows that states with S ≲ s þ 3=4 typically have the highest cost of transport, whereas those with S ≲ s and S ≲ s þ 1=4 have the lowest. Comparing Figs. 3(b) and 3(c), 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 trade-off 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 v x and lift F v y on a swimmer due to the vortices shed by its neighbors, assuming that the vortices (with positions z k ∈ C and strengths γ k ) are far from the body, jz k j ≫ r. In Supplemental Material Sec. I [71], we show that is the vertical velocity of the flow induced by the neighboring swimmers' shed vortices and ðU; VÞ is the swimmer's velocity. Figure 3(d) 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 v x < 0 and thus a speed-up (U=U 0 > 1) for s þ 1=4 ≲ S ≲ s þ 3=4, and the opposite for s þ 3=4 ≲ S ≲ s þ 5=4, according to Eq. (5). We note that these conclusions may be justified mathematically by adapting the argument in Supplemental Material Sec. II A [71], in which we derive an approximate form for the interaction force Eq. (5) in the biologically relevant low-Strouhal number limit St ≪ 1 [20].
A similar argument may be used to qualitatively understand the trade-off 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 f < 0), which implies F v y < 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 St ≪ 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 Sec. IVA [71].
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 formulas in Eq. (5) for F v x and F v y would simply have their signs reversed. We thus expect a speed-up and larger cost of transport for schooling numbers s þ 3=4 ≲ S ≲ s þ 5=4, and a slowdown with lower cost of transport for s þ 1=4≲ S ≲ s þ 3=4. That is, Figs. 3(b) and 3(c) would be qualitatively unchanged, apart from a shift of the horizontal axis, S → S þ 1=2.

B. Phalanx
Motivated by the experimental observations of Ashraf et al. [30], we now consider a phalanx of swimmers: infinitely many swimmers equally spaced by a distance d in the lateral (y) direction and, following Weihs [50], flapping in antiphase with respect to their neighbors, as shown in Fig. 4 and Supplemental Material Movie 2 [71]. As discussed in Supplemental Material Sec. VA, Weihs [50] and Stöcker [76] 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 [71]. Such an argument is in agreement with the experimental observations of Ashraf et al. [42], 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 min ≡ 2ðh 0 þ d 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 Sec. II permits consideration of this formation, as shown in Supplemental Material Sec. V [71]. Following the procedure detailed in Sec. III, we find that the swimming speed U satisfies an algebraic equation of the form Eq. (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 multistability of steady states seen for in-line formations (Sec. IVA). By generalizing the linear stability analysis of steady-state solutions presented in Sec. III, we find that the steady state is stable for all values of d. As shown in Fig. 4, the formation exhibits a slight speed-up for all values of d, with the maximum speed-up of roughly 5% occurring when the swimmers are most tightly packed, d ¼ d min . However, such a formation also increases the cost of transport by roughly 4%, with the cost of transport decreasing as d → ∞. Similar to the 1D formations discussed in Sec. IVA, the phalanx formation exhibits a trade-off between speed and cost of transport, although both measures exhibit relatively small variations over a range of values d.

C. Rectangular lattice
We now consider the rectangular lattice of swimmers shown in Fig. 5(a) and Supplemental Material Movie 3 [71]: swimmers separated by a streamwise distance L and a vertical distance d, starting at the positions z ¼ jL þ ikd for j; k ∈ Z. Swimmers flap in phase (antiphase) with respect to their streamwise (lateral) neighbors. As in Sec. II, a swimmer at the origin [red box in Fig. 5(a)] interacts with the swimmers in the neighboring columns at x ¼ 0; AEL. Note that the lattice is effectively an in-line formation in the limit d → ∞. In Supplemental Material Sec. 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 av . Figure 5(b) shows the normalized velocity U av =U 0 as a function of L and d. As with the in-line formation, there may be multiple steadystate solutions for a given set of parameters; in such cases, we perform multiple simulations and plot the timeaveraged 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. 5(b). As the lateral distance d between swimmers is decreased progressively, these regions of multistability typically shrink, but new regions of multistability may also emerge.
We find that the formation experiences a maximum speed-up 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 av ðL; dÞ=U 0 , the optimal lateral spacing for a given streamwise spacing, and the corresponding speed U max ðLÞ ¼ jU av (L; d U ðLÞ)j. Note that, unlike the phalanx (Fig. 4), the greatest speed-up is not necessarily achieved by packing the swimmers tightly in the y direction, so d U ðLÞ is not identically equal to d min . A comparison between the rectangular lattice and in-line formation is shown in Fig. 5(d). The most salient feature is that the slow modes (U=U 0 < 1) for an in-line formation all exhibit speed-up in the corresponding rectangular lattice with the minimum lateral spacing, d U ðLÞ ¼ d 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 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. 6(b) 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 → ∞), for which C=C 0 ¼ 0.75 and S ¼ 1.26 [ Fig. 3(c)]. 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 min ðLÞ ¼ C(L; d C ðLÞ). Unlike d U ðLÞ, which exhibits a nontrivial dependence on the streamwise spacing L [ Fig. 5(b)], d C ðLÞ typically assumes the values d C ≈ d min and d C ¼ ∞ [ Fig. 6(b)]. Figure 6(d) shows that the states for which d C ≈ d 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 ¼ ∞ 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 Figs. 5(d) and 6(d) 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 Sec. IVA 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 speed-up (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 ≲ S ≲ s þ 3=4, the upstream swimmers at z ¼ −L AE id in a rectangular lattice contribute an upflow [ Fig. 5(a)], 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 ≲ S ≲ s þ 5=4, the parameter regime in which the upstream swimmers at z ¼ −L AE 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).

D. Diamond lattice
We now consider the diamond lattice shown in Fig. 7(a) and Supplemental Material Movie 4 [71]: 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 ∈ Z and l ∈ f0; 1g. 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. 7(a)] interacts with the swimmers in the neighboring columns at x ¼ 0; AEL=2; AEL.
The complete equations are detailed in Supplemental Material Eq. (45) [71]. 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 where α j ¼ mod ðj; 2Þ=2 accounts for the vertical shift in the columns at x ¼ AEL=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 Sec. III to assess the stability of these period-2 solutions, in Supplemental Material Secs. 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 Sec. IV C, the simulations reveal the existence of multiple coexisting states, indicated by the regions of parameter space bounded by the gray curves in Fig. 7(b). In such cases, Fig. 7(b) shows the normalized average speed-up U av =U 0 of the fastest state. We find that the formation experiences a maximum speed-up of 22% for L ¼ 4.4 and d ¼ d min , which corresponds to S ¼ 0.67. As with the rectangular lattice considered in Sec. IV C, d U ðLÞ is not identically equal to d 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. 7(d). As with the rectangular lattice, the slow modes (U=U 0 < 1) for an in-line formation exhibit speed-up in a diamond lattice with the minimum lateral spacing, d U ðLÞ ¼ d min . Most fast modes for an in-line formation (U=U 0 > 1) benefit minimally from the diamond geometry; however, those with S ≲ s þ 3=4 experience a speed-up of 3%-5% relative to the in-line formation, also with d ¼ d min . This is due to the existence of an entirely new branch of fast modes that emerges for d ≳ d min , as shown in Supplemental Material Fig. 3 [71]. These fast modes may be attributed to the beneficial dipolar structure generated by the swimmer's upstream neighbors at z ¼ −L=2 AE id=2, as shown in Fig. 7(a). 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 min , indicating that a tightly packed diamond lattice formation realizes a substantial increase in efficiency due to hydrodynamic interactions. Figure 8(b) shows that, for a given streamwise spacing L, the optimal lateral spacing is typically d C ðLÞ ≈ d min or  d C ðLÞ ¼ ∞. This effect is similar to that observed for the rectangular lattice, and may be rationalized by the physical argument presented in Sec. IV C. However, the diamond lattices for which S ≈ 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 AE id=2 [ Fig. 8(a)], which provides an additional lift force as the swimmer accelerates upward. While the distinct optimal diamond lattice formations shown in Figs. 7(a) and 8(a) 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 ≳ d 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 present 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 [66], while using only three fitting parameters (C D , τ, 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 multistability 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 Material Fig. 2 [71]) 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. 3(d). 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 ≲ 0.75 [ Fig. 3(b)]. By contrast, a phalanx formation of swimmers flapping in antiphase moves roughly 5% faster than a single swimmer (Fig. 4), with the optimal speed-up 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. 7(a)].
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 ≳ 1.25 [ Fig. 3(c)]. While our finding that the phalanx formation affords a speed-up is in agreement with the experimental observations of Ashraf et al. [30], we also find such formations to be slightly less efficient than an isolated swimmer (Fig. 4). Generally, we observe that inline and phalanx formations exhibit a trade-off 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. 8(a)]. 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 Sec. IVA, 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 → S þ 1=2, we still expect the diamond lattices with S ¼ s þ 1=2 to be particularly efficient [ Fig. 8(d)], 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 min is consistent with the results of Weihs [49,50]. However, our results differ in that the curves d U ðLÞ and d C ðLÞ [black curves in Figs. 7(d) and 8(d), 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's 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 (Sec. IVA). More sophisticated computational models [63,77] of fishlike 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 ∼Uτ between swimmers. Third, we explicitly account for the formation's dynamics and thus the schooling modes' stability, which was not possible in Weihs's 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. [64] conducted 2D numerical simulations using a multiparticle collision dynamics model and also found that the diamond lattice provides the largest speed-up 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. [64] did not consider the influence of the streamwise spacing L, which we find to play a dominant role. Daghooghi and Borazjani [65] 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 [Figs. 5(b) and 6(b)]. Tsang and Kanso [78] proposed a farfield 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 [67]. 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 multistability 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. [30,42] reported that highspeed 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 ≈ 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 [27,43]. When in relatively large schools, this fish species has also been observed to adopt diamond-lattice-like formations with a mean firstand second-nearest neighbor separation angle between 14°a nd 17° (Fig. 5 in Ref. [27]). 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°[ Fig. 2(b) in Ref. [30] ]. Moreover, this fish species tends to adopt the phalanx over the diamond lattice formation at high swimming speeds [30], 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 [48], ibis flocks were observed to preferentially assume V formations with median schooling number S ≈ s þ 1=4, and in-line formations with S ≈ s þ 1=2 (Supplemental Fig. 2 in Ref. [18]). Extensions of our model might shed light on these phenomena, for which lift generation is an important consideration [70]. 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 [67,74].
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 [13], and may be extended to probe the droplet's complex collective dynamics [14,15]. Generally, we expect models of the type described herein to be broadly applicable to systems of active particles interacting via their collective histories.