1. Introduction
Fish schools are ubiquitous in aquatic life, with half of the known fish species thought to exhibit schooling behaviour during some phase of their life cycle (Shaw Reference Shaw1978). However, the role of the fluid medium as a mediator of the physical interactions between swimming fish remains unclear (Partridge & Pitcher Reference Partridge and Pitcher1979; Partridge Reference Partridge1982). Experimental evidence suggests that fish modify their motions and reduce muscular effort when swimming in vortex-laden flows (Liao et al. Reference Liao, Beal, Lauder and Triantafyllou2003). These findings support a long-standing but controversial hypothesis that schooling provides hydrodynamic benefits as fish move within the flows generated by others (Weihs Reference Weihs1973, Reference Weihs1975; Abrahams & Colgan Reference Abrahams and Colgan1985; Liao Reference Liao2007). A direct assessment of this hypothesis in biological and physical models remains a challenge because of the complexity in resolving the hydrodynamics of unsteady swimming at high Reynolds numbers with single (Wolfgang et al. Reference Wolfgang, Anderson, Grosenbaugh, Yue and Triantafyllou1999; Triantafyllou, Triantafyllou & Yue Reference Triantafyllou, Triantafyllou and Yue2000; Borazjani Reference Borazjani2008) and multiple interacting swimmers (Liao Reference Liao2007; Gazzola et al. Reference Gazzola, Tchieu, Alexeev, de Brauer and Koumoutsakos2016; Verma, Novati & Koumoutsakos Reference Verma, Novati and Koumoutsakos2018). Simplifications based on crystalline school arrangements and ideal flow models indicate that fish within a planar formation, with diamond-shaped unit cell, benefit energetically from near-field interactions with the wakes of upstream neighbours (Weihs Reference Weihs1973), whereas far-field interactions serve to passively stabilize the formation (Tsang & Kanso Reference Tsang and Kanso2013). These crystal lattice models do not capture that fish exhibit variable arrangements in field and laboratory experiments (Partridge & Pitcher Reference Partridge and Pitcher1979; Marras et al. Reference Marras, Killen, Lindström, McKenzie, Steffensen and Domenici2015), and the broader question of how flow interactions benefit schooling remains unresolved.
Physical models and numerical simulations of mechanically actuated foils found that, at the single swimmer level, flapping foils share with their biological counterparts many common aspects of the flows, forces and energetics (Blondeaux et al. Reference Blondeaux, Fornarelli, Guglielmini, Triantafyllou and Verzicco2005; Dong, Mittal & Najjar Reference Dong, Mittal and Najjar2006; Buchholz & Smits Reference Buchholz and Smits2008; Dabiri Reference Dabiri2009; Lauder et al. Reference Lauder, Lim, Shelton, Witt, Anderson and Tangorra2011; Wen & Lauder Reference Wen and Lauder2013). A key similarity is the reverse von Kármán wake left by both flapping foils and fish (Taneda Reference Taneda1965; Triantafyllou, Triantafyllou & Grosenbaugh Reference Triantafyllou, Triantafyllou and Grosenbaugh1993). Subsequently, several numerical and experimental studies used pairs of flapping foils to understand multi-swimmer interactions. Zhu, He & Zhang (Reference Zhu, He and Zhang2014) were first to examine, in the context of the immersed boundary method, the effects of pairwise hydrodynamic interactions on the self-propulsion of flapping flexible swimmers in tandem configuration. Flow-mediated interactions were found to stabilize the swimmers in particular spacings and to reduce the energetics cost of swimming in the follower. Experimental studies on heaving rigid foils confined to in-line positions and freely swimming in tandem were also found to assume one of several particular spacings, stabilized by the flow interactions (Becker et al. Reference Becker, Masoud, Newbolt, Shelley and Ristroph2015; Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016; Newbolt, Zhang & Ristroph Reference Newbolt, Zhang and Ristroph2019). These observations have since been confirmed in several numerical studies (Dai et al. Reference Dai, He, Zhang and Zhang2018; Park & Sung Reference Park and Sung2018; Peng, Huang & Xi-Yun Reference Peng, Huang and Xi-Yun2018; Lin et al. Reference Lin, Wu, Zhang and Yang2020). Here, we investigate the speed, energetics and stability of these planar formations using a mathematical model of self-propelling and interacting swimmers that flap by either heaving or pitching.
Existing mathematical models of flow interactions in fish schools vary in the degree of fidelity to the fluid dynamics and sensory-feedback control at the swimmer level. Ideal flow models – based on a dipolar far-field approximation (Tchieu, Kanso & Newton Reference Tchieu, Kanso and Newton2012) – with no feedback control have been used to assess the effect of passive flow interactions on the stability of pairwise (Kanso & Tsang Reference Kanso and Tsang2014, Reference Kanso and Tsang2015) and diamond lattice formations (Tsang & Kanso Reference Tsang and Kanso2013) and the advantages of flapping out of phase (Kanso & Newton Reference Kanso and Newton2009). This far-field flow model coupled to visual feedback control, either in the form of behavioural rules (Filella et al. Reference Filella, Nadal, Sire, Kanso and Eloy2018) or learning algorithms (Gazzola et al. Reference Gazzola, Tchieu, Alexeev, de Brauer and Koumoutsakos2016), was used to analyse the fish collective dynamics. Fish were shown to exhibit a novel collective turning mode and to swim faster thanks to the fluid (Filella et al. Reference Filella, Nadal, Sire, Kanso and Eloy2018). Near-field fish–wake interactions were also accounted for in ideal flow models with no feedback control, such as the vortex street model used by Weihs (Reference Weihs1973) or the phenomenological model derived in Oza, Ristroph & Shelley (Reference Oza, Ristroph and Shelley2019) to assess the efficiency of lattice formations. High-fidelity computational fluid dynamics coupled to reinforcement learning algorithms were recently implemented in pairwise interactions to optimize the flapping motion of the follower fish for harnessing the wake of the leader (Verma et al. Reference Verma, Novati and Koumoutsakos2018).
In this paper, we analyse pairwise interactions of heaving and pitching swimmers in the context of the vortex sheet model (see figure 1). The vortex sheet model has been used extensively to analyse problems of fluid–structure interactions, including ring formation at the edge of a circular tube (Nitsche & Krasny Reference Nitsche and Krasny1994) and wakes of oscillating plates (Jones Reference Jones2003; Sheng et al. Reference Sheng, Ysasi, Kolomenskiy, Kanso, Nitsche and Schneider2012), falling cards (Jones & Shelley Reference Jones and Shelley2005), flapping flexible flags (Alben & Shelley Reference Alben and Shelley2008; Alben Reference Alben2009), swimming plates (Wu Reference Wu1971) and hovering flyers (Huang, Nitsche & Kanso Reference Huang, Nitsche and Kanso2016; Huang et al. Reference Huang, Ristroph, Luhar and Kanso2018). Here, we use the implementation of Nitsche & Krasny (Reference Nitsche and Krasny1994), which we vetted in comparison to Navier–Stokes simulations and other implementations of the vortex sheet method in Sheng et al. (Reference Sheng, Ysasi, Kolomenskiy, Kanso, Nitsche and Schneider2012), Huang et al. (Reference Huang, Nitsche and Kanso2016). This study focuses on the effect of streamwise flow interactions on the swimming motion of heaving and pitching plates, and finds that ordered formations emerge spontaneously via these interactions, independent of the flapping mode, consistent with heaving foil experiments (Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016; Newbolt et al. Reference Newbolt, Zhang and Ristroph2019) and numerical simulations (Zhu et al. Reference Zhu, He and Zhang2014; Park & Sung Reference Park and Sung2018; Peng et al. Reference Peng, Huang and Xi-Yun2018; Lin et al. Reference Lin, Wu, Zhang and Yang2020). However, the flapping mode, heaving or pitching, affects the speed and energetics of these formations as well as their robustness to streamwise perturbations. We describe a specific hydrodynamic mechanism that explains the energetic and stability differences associated with each flapping mode.
Figure 1. A pair of swimmers undergoing (a) heaving motions at amplitude 
$A_{{h}} = 0.3$ and (b) pitching motions at amplitude 
$A_{{p}} = 15^{\circ }$. Snapshots of the velocity field (grey arrows) and free vortex sheet of the leader (blue) and follower (red) are taken after steady-state swimming is reached at a time instant where both swimmers are flapping downwards. Insets depict the pressure forces acting on each swimmer in the pairwise formation in comparison to a single swimmer undergoing the same prescribed motion.
2. Problem formulation
A swimmer is modelled as a rigid plate of length 
$2l$, small thickness 
$e\ll l$ and homogenous density 
$\rho$, submerged in an unbounded, planar, fluid domain of density 
$\rho _f$. The swimmer's mass per unit depth is given by 
$m = 2 \rho e l$. An inertial frame 
$(\boldsymbol {e}_x,\boldsymbol {e}_y,\boldsymbol {e}_z)$ is introduced, such that 
$(\boldsymbol {e}_x,\boldsymbol {e}_y)$ span the plane of motion. The vector 
$\boldsymbol {x} \equiv (x,y)$ denotes the position of the leading edge of the swimmer in the 
$(\boldsymbol {e}_x,\boldsymbol {e}_y)$ plane, and the angle 
$\theta$ its orientation relative to the 
$\boldsymbol {e}_x$-direction (see Appendix A and figure 7)
The swimmer is free to move in the 
$\boldsymbol {e}_x$-direction under periodic heaving or pitching motions. Heaving consists of periodic lateral motions in the 
$y$-direction, of amplitude 
$A_{{h}}$, at fixed angle 
$\theta = 0$. Pitching refers to angular oscillations 
$\theta$ of amplitude 
$A_{{p}}$, with zero lateral motion 
$y=0$ at the leading edge. The frequency of these heaving and pitching motions is denoted by 
$f$. Hereafter, we scale all parameter values using 
$l$ as the characteristic length scale, 
$1/f$ as the characteristic time scale and 
$\rho _f l^{2}$ as the characteristic mass per unit depth. Accordingly, velocities are scaled by 
$lf$, forces by 
$\rho _f f^{2} l^{3}$, moments by 
$\rho _f f^{2} l^{4}$ and power by 
$\rho _f f^{3} l^{4}$.
In dimensionless form, the heaving and pitching motions are given by
\begin{equation} \left.\begin{array}{c@{}}
\textrm{Heaving:} \ y(t) = A_{{h}} \sin(2 {\rm \pi}t), \quad
\theta(t) = 0, \\ \textrm{Pitching:} \ \theta(t) = A_{{p}}
\sin(2{\rm \pi} t),\quad y(t) = 0 . \end{array}\right\}
\end{equation}
The equation of motion governing the free swimming 
$x(t)$ is given by Newton's second law
\begin{equation} m \ddot{x}={-} F \sin\theta + S \cos\theta - D \cos\theta.\end{equation}
Here, the hydrodynamic forces acting on the swimmer consist of a leading edge suction force 
$S$, a pressure force 
$F$ acting in the direction normal to the swimmer and a skin drag force 
$D$ acting tangentially to the swimmer in the opposite direction to its motion. The drag force 
$D$ is introduced to emulate the effect of fluid viscosity, while the hydrodynamic pressure force 
$F$ is calculated in the context of the inviscid vortex sheet model. A detailed description of the method and its numerical implementation can be found in Nitsche & Krasny (Reference Nitsche and Krasny1994), Huang et al. (Reference Huang, Ristroph, Luhar and Kanso2018) and a brief overview is given in Appendix A. Detailed expressions of the fluid forces and moments acting on the swimmer are given in Appendix B.
To assess the swimming performance, we use four metrics: the period-averaged swimming speed 
$U = \int _t^{t+1} \dot {x} \, {\textrm {d}}t$ at steady state, the thrust force 
$T = S \cos \theta - F \sin \theta$, the input power 
$P$ required to maintain the prescribed heaving or pitching motions (see details in Appendix D) and the cost of transport defined as the input power 
$P$ divided by the swimming speed 
$U$.
3. Single swimmers: numerical results and scaling analysis
We solve (2.2) in the case of a single swimmer and compute the period-average swimming speed at steady state. In figure 2(a,b), we show the steady-state speed for heaving and pitching swimmers, respectively, as a function of the flapping amplitude. In both cases, the speed increases monotonically, albeit that, when pitching, the increase scales differently at small amplitudes. To get insight into how the swimming speed 
$U$ scales with the heaving and pitching amplitudes and frequency, it is instructive to use a simple scaling analysis.
Figure 2. Swimming speed vs. flapping amplitude for single swimmers. (a) Average swimming speed at steady state for a heaving swimmer. (b) Average swimming speed at steady state for a pitching swimmer. At small 
$A_{{p}}$, skin drag is dominant and the speed scales super-linearly with 
$A_{{p}}$. For 
$A_{{p}}>10^{\circ }$, pressure drag is dominant and speed scales linearly with 
$A_{{p}}$. (c) Experimental data (black symbols) of average swimming speed of a heaving foil (Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016, figure 2); the data collapse when scaled by the heaving frequency 
$f^{4/3}$ (yellow symbols). (d) Comparing the swimming speed of our heaving swimmer model (blue circles) to the frequency-scaled experimental data shown in (a), on a log–log scale. Both model and experimental results scale super-linearly with heaving amplitude.
At steady state, the sum of forces acting on the swimmer is zero on average. For heaving swimmers, the dominant forces are those due to leading edge suction and viscous skin drag (Garrick Reference Garrick1937). In dimensional form, the suction force scales as 
$\rho _{{f}} (2l)C_s^{2}U^{2}$, where the coefficient 
$C_s$ scales linearly with the effective angle of attack. In a heaving flat plate the effective angle of attack is given by 
$\dot {y}/ U \sim A_{{h}}f / U$ (Garrick Reference Garrick1937; Floryan et al. Reference Floryan, Van Buren, Rowley and Smits2017; Franck & Breuer Reference Franck and Breuer2017; Smits Reference Smits2019), and the suction force scales as 
$\rho _{{f}} (2l)(A_{{h}}f)^{2}$. Skin drag scales as 
$\rho _f (2l) C_{{f}}U^{2}$, where 
$C_{{f}} \sim \sqrt {\mu /\rho _f (2l) U}$ is the drag coefficient based on adapting Blasius theory to this inviscid fluid model (see Appendix C and White Reference White1979). Balancing suction and drag forces, we arrive at 
$(A_h f)^{2} \sim U^{3/2}$, which leads to
\begin{equation} \textrm{Heaving: } U\sim (A_{{h}}f)^{4/3}. \end{equation}
The swimming speed scales super-linearly with the heaving amplitude and frequency. We test this scaling law in light of the experimental results of (Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016, figure 2). The black data points in figure 2(c) represent the experimentally measured swimming speed as a function of heaving amplitude. The different marker shapes represent three different heaving frequencies used in the experiments (
$f=1, 2, 3$). We scaled the data by the heaving frequency according to our derived scaling law in (3.1). The scaled data (coloured symbols) collapse on a single curve, indicating that our scaling analysis is sound. In figure 2(d), we plot, using a log–log scale, the swimming speed obtained from our model in figure 2(a) (blue dots) and experimental data (coloured symbols) vs. the heaving amplitude. The slope of each line represents the power law that governs the relationship between the two quantities. In both the model and the experiment, the swimming speed depends super-linearly on the amplitude of heaving, however, the dependence is slightly stronger in the model.
The steady-state speed of the pitching swimmer scales differently depending on the flapping amplitude because the dominant drag forces acting on the swimmer differ. At small pitching amplitude 
$A_{{p}}$, the swimmer is almost parallel to the swimming direction, hence skin drag is dominant leading to the same scaling law as in the heaving case. At large amplitude 
$A_{{p}}$, pressure drag is dominant; it is well known that pressure drag scales as 
$U^{2}$; see, e.g. Moored & Quinn (Reference Moored and Quinn2019). Balancing inertia and pressure drag, we arrive at 
$U \sim A_p f$. Put together, we have
\begin{equation} \textrm{Pitching:} \begin{cases} \textrm{small} \ A_{{p}}: U \sim (A_{{h}}f)^{4/3},\\ \textrm{large} \ A_{{p}} : U \sim A_{{p}}f. \end{cases} \end{equation}These scaling laws fit remarkably well the numerical results in figure 2(b).
4. Pairwise formations: stability, speed and energetics
We examine the steady-state behaviour of a pair of swimmers undergoing heaving and pitching motions while freely interacting via the fluid medium. In figure 1, we show snapshots of the flow field (grey arrows) and free vortex sheets in the case when the leader (blue) and follower (red) are heaving at 
$A_{{h}} = 0.3$ (figure 1a) and pitching at 
$A_{{p}} = 15^{\circ }$ (figure 1b). The snapshots are taken after the pair has reached steady-state swimming in the positive 
$x$-direction, and passively locked into a constant separation distance. At these flapping amplitudes, the heaving swimmers experience longer transience and swim faster, whereas the pitching swimmers rapidly lock into a tighter formation (see supplementary movie available at https://doi.org/10.1017/jfm.2021.551).
An analysis of the hydrodynamic pressure forces 
$-F \sin \theta \boldsymbol {e}_x + F \cos \theta \boldsymbol {e}_y$, where 
$F$ is given in Appendix B, acting on each swimmer shows that compared to a single swimmer, the distribution on the leader remains relatively unchanged. However, the force distribution on the follower is affected by the wake of the leader, and the effect is more pronounced for pitching swimmers; see insets in figure 1(a,b). Specifically in the pitching case, the follower experiences less resistance from the fluid, and a favourable force distribution (in the same direction of flapping) at the swimmer's tail. At the instant shown in figure 1(b), the downward flow due to the vortex sheet created by the leader helps the follower in its downward pitching motion.
In figure 3, we vary the initial separation distance between the two swimmers for the examples shown in figure 1. We find that for both heaving and pitching, the follower tends to settle in one of several discrete locations behind the leader at nearly digital values of 
$d_{{h}}/\lambda$ and 
$d_{{p}}/\lambda$, respectively, where 
$d_{{h}}$ is the tail-to-head distance, 
$d_{{p}}$ the tail-to-tail distance and 
$\lambda = U/f$ the wavelength of the leader's swimming trajectory; see figure 3(a,b). Depending on initial conditions, the leader and follower reach one of these separation distances and swim together in ordered formation. These findings are consistent with the observations of Zhu et al. (Reference Zhu, He and Zhang2014), Ramananarivo et al. (Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016), Park & Sung (Reference Park and Sung2018), Peng et al. (Reference Peng, Huang and Xi-Yun2018) and Lin et al. (Reference Lin, Wu, Zhang and Yang2020).
Figure 3. Emergence of passive stable formations in a pair of heaving swimmers (
$A_{{h}} = 0.3$) and of pitching swimmers (
$A_{{p}} = 15^{\circ }$). (a) For heaving swimmers, the follower stabilizes at one of many discrete positions behind the leader where the gap (tail-to-head) distance 
$d_{{h}}$ is close to integer multiple of the wavelength 
$\lambda = U/f$ of the leader motion. (b) For pitching swimmers, the follower stabilizes at locations such that the tail-to-tail distance 
$d_{{p}}$ is close to integer multiples of 
$\lambda$. Basins of attraction of each the first three equilibria are depicted in gradually more faint shades of grey. (c,d) Linear stability analysis: we perturb the position of the follower about each of these equilibria and compute the total hydrodynamic force 
$F_x$. We simultaneously sample data from the change in 
$F_x$ and perturbation strength 
$\delta x$, and plot 
$\delta F_x$ vs. 
$\delta x$. Clearly, 
$\delta F_x$ acts as a restoring force. Taking the slope of 
$\delta F_x$, we construct the hydrodynamic potential 
$V$ on the follower. The potential well is deepest at the first equilibrium where the hydrodynamic interactions are strongest.
To examine the nonlinear basins of attraction of these equilibria, we vary the initial separation distance 
$d_{{h}}$ and 
$d_{{p}}$ between the two swimmers and keep track of the corresponding steady-state formation; The basin of attraction of each equilibrium is highlighted in a different shade of grey in figure 3(a,b). The pitching swimmers converge more rapidly to the corresponding equilibria, indicating that these equilibria are stronger attractors in pitching than in heaving. Further, the wavelength 
$\lambda = U/f$ is smaller in pitching, and so is the actual separation distance at equilibria (
$d_{{p}}< d_{{h}}$), indicating that pitching swimmers move in tighter formations.
To quantitatively assess the linear stability of these equilibria, we perturb the position of the follower about each equilibrium in the positive and negative 
$x$-direction with an initial perturbation of size 
$\delta x/l = 0.5$ and we calculate the corresponding change in 
$\delta x$ and change in the total hydrodynamic force 
$\delta F_x = \delta (- F \sin \theta + S \cos \theta - D \cos \theta )$ acting on the follower in the 
$x$-direction. We scale the change in total force by 
$U^{2}/l$ and the perturbation from equilibrium by 
$d_\cdot /\lambda$, where 
$d_\cdot$ is either 
$d_{{h}}$ or 
$d_{{p}}$. We sample simultaneously the scaled change in total force 
$\delta F_x$ and scaled perturbation strength 
$\delta x$ and we plot the results in the first row of figure 3(c,d). The results are depicted in red triangles for the first stable position, and in orange circles and yellow squares for the second and third positions, respectively. Straight line fit for each of these data sets results in straight lines with negative slopes, implying that, for each of these equilibrium positions, the hydrodynamic force acts as a restoring force 
$\delta F_x = - K\delta x$ that keeps the formation stable. Here, 
$K$ is obtained numerically from the straight line fit. The value of 
$K$ depends monotonically on the equilibrium position of the follower, with highest value at the first equilibrium (
$d_{{h}}/\lambda \approx 1$ and 
$d_{{p}}/\lambda \approx 1$). The first equilibrium is most stable because hydrodynamic interactions are strongest at closer distance. We write 
$\delta F_x = - \partial V/\partial (\delta x)$, where 
$V = K (\delta x)^{2}/2$ is the hydrodynamic potential function around the equilibrium 
$\delta x =0$. For both pitching and heaving, the formation is stable with weaker stability for larger inter-swimmer distance. In the pitching formation the potential well is deeper (by approximately 
$50\,\%$) for all equilibria, indicating faster convergence to the respective equilibrium; see Appendix E for detailed analysis of the hydrodynamic forces during transient and equilibrium states.
We next evaluate the advantages of these formations in terms of the speed and energetics of the pair of swimmers in comparison with swimming alone. Figure 4 shows details of the time evolution at steady state of a single and pair of swimmers for the first relative equilibrium 
$d_{{h}}/\lambda \approx 1$ and 
$d_{{p}}/\lambda \approx 1$ shown in figure 3, where hydrodynamic interactions are strongest. From top to bottom, we report the swimming speed, thrust force, input power and cost of transport vs. time. Instantaneous values are shown in solid lines and period-average values in dashed lines. For the heaving motion, the average speed of the pair is approximately 10 % higher than the speed of the single swimmer, consistent with experimental observations on heaving foils (Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016). However, the input power required to maintain these heaving motions in the presence of hydrodynamic interactions is also higher (approximately 30 %). Consequently, the cost of transport of the heaving pair is approximately 20 % higher than a single heaving swimmer. These results suggest that heaving swimmers can enhance their speed by swimming in a pair. However, this enhancement in swimming speed is achieved at an energetic cost.
Figure 4. Instantaneous swimming performance (time-dependent speed, thrust, input power and cost of transport vs. time) for a single and pair of swimmers undergoing (a) heaving at 
$A_h = 0.3$ and (b) pitching at 
$A_p = 15^{\circ }$, respectively. Results are shown after the swimmers have reached steady state. From top to bottom, the swimming speed, thrust force, input power and cost of transport are shown. Solid lines represent the instantaneous values and dashed lines represent time-period averages.
For pitching swimmers, the speed of the formation is comparable to that of the single swimmer (approximately 2 % slower). However, the follower's input power is significantly reduced (approximately 70 % less than the single pitching swimmer). This reduction in input power is due to the hydrodynamic benefits highlighted in figure 1(b). Correspondingly, the cost of transport of the pair of pithing swimmers drops by 30 % compared with swimming alone.
Figure 5 explores the effect of the flapping amplitude on the period-average values of the swimming speed, thrust force, input power and cost of transport, after the swimmers have reached steady state. Specifically, we examine the range 
$A_{{h}} \in [0, 0.7]$ and 
$A_p \in [0^{\circ }, 45^{\circ }]$ for single swimmers and 
$A_{{h}} \in [0.3, 0.7]$ and 
$A_p \in [10^{\circ }, 45^{\circ }]$ for pairs of swimmers, where small amplitudes are ignored to ensure that hydrodynamic interactions are sufficient for the spontaneous emergence of order formations. In pairwise interactions, we report all period-average values normalized by the corresponding values for a single swimmer.
Figure 5. Swimming performance (average speed, thrust, input power and cost of transport) vs. flapping amplitude for a single and pair of swimmers undergoing (a) heaving and (b) pitching motions, respectively. From top to bottom, average values of the swimming speed, thrust force, input power and cost of transport. Left columns (black symbols) in (a,b) show the results for single swimmers. For the pair of swimmers, all of the results are scaled by the corresponding quantity values for a single swimmer. The blue and red symbols represent the results for the follower and leader, respectively. The grey symbols are the school average.
When swimming alone, whether by heaving or pitching, an increase in the flapping amplitude monotonically increases the swimming speed, thrust, input power and cost of transport; see left columns of figure 5(a,b). Here, the swimming speed vs. flapping amplitude for single swimmers is a reproduction of the results in figure 2(a,b).
Across all heaving amplitudes, the pairwise formation is approximately 5 %–10 % faster than that of a single heaving swimmer. Both the leader and follower experience an increase in thrust compared with the single swimmer, but require more power to swim in formation compared with swimming alone, with extra power demand on the follower. The cost of transport of the heaving formation is thus slightly higher (approximately 
$15\,\%$) compared with swimming alone. Thus, heaving swimmers slightly enhance their swimming speed when in formation, albeit at a higher cost of transport.
The formation of pitching swimmers is approximately 5 % slower than swimming alone for almost all flapping amplitudes. The leader experiences consistently lower thrust and the follower consistently higher thrust compared with swimming alone. However, while the power demand on the leader is comparable to the single swimmer, the power demand on the follower is significantly reduced for all amplitudes. Taken together, these results lead to slightly higher cost of transport for the leader and significantly lower cost of transport for the follower compared with swimming alone. Indeed, the cost of transport of the follower is a fraction of the single swimmer (approximately 
$25\,\%$ at best), which in turn, causes the formation to save a significant amount of power (approximately 
$35\,\%$ at best) compared with swimming alone. These results imply that although the pairwise formation of pitching swimmers experiences no enhancement in swimming speed compared with swimming alone, it reduces the cost of transport by a significant amount.
5. Single swimmer wake informs pairwise formation
To gain additional insights into the information contained in the wake of the leader and the hydrodynamic mechanisms that mediate the power reduction and stability of the pairwise formation, we examine the flow field induced by a single swimmer. Namely, we compute the flow field generated behind a single heaving or pitching swimmer, and we consider a virtual ‘point’ follower placed at any location 
$(x_o,y_o)$ in the swimmer's wake and undergoing lateral oscillations 
$y(t) = y_o + A \sin (2{\rm \pi} t)$, where 
$A$ is the oscillation amplitude. We set 
$A$ to 
$A_{h}$ in the wake of the heaving swimmer and 
$A_{p}$ in the wake of the pitching swimmer. The wake is blind to the existence of the virtual follower. We ask whether there are particular locations in the swimmer's wake that are favourable to the follower's flapping motion. To address this question, we define a flow agreement parameter 
$\mathbb {Z}(x_o,y_o)$ that quantifies the agreement between the flow velocity in the wake of the single swimmer and the prescribed oscillations of the virtual follower,
\begin{equation} \textrm{flow agreement parameter:} \ \mathbb{Z} = \frac{1}{T} \int_{t_s} ^{t_s+T} \dot{y}(x_o,y_o,t) v(x_o,y_o,t) \, {\textrm{d}}t,\end{equation}
where 
$t_s$ is an arbitrary time after the single swimmer has reached steady state, 
$T$ is the flapping period, 
$\dot {y}(x_o,y_o,t)$ is the lateral velocity of the follower and 
$v(x_o,y_o,t)$ is the 
$y$-component of the flow velocity evaluated at the follower's location. Positive values of the flow agreement parameter imply a beneficial interaction between the flow and the follower's flapping motion, whereas negative values indicate a detrimental one.
The flow agreement parameter is shown in the top row of figure 6. The hypothetical follower is undergoing the same oscillatory motion irrespective of its location in the wake of the single swimmer at amplitude 
$A_{{h}} = 0.3$ (for heaving) and 
$A_{{p}} = 15^{\circ }$ (for pitching). Red regions indicate where the flow velocity in the swimmer's wake and the hypothetical follower's motion agree. Interestingly, regions of maximum flow agreement are located at almost integer multiples of the wavelength 
$\lambda = U/f$ of the single swimmer, similarly to the locations of the stable equilibria in fully coupled pairwise formations. Superimposed onto figure 6(a,c), we show a snapshot of the free vortex sheet of the swimmer, as well as the location of the actual follower at steady state obtained from our pairwise interacting swimmers. As noted previously, in the heaving case, the leading edge of the follower is located close to the integer multiples of 
$\lambda$, while in pitching, the follower's trailing edge is located at integer multiples of 
$\lambda$. In the heaving case, the leading edge of the follower is located at the intersection of the red and blue regions of the flow agreement parameter, that is at the location where the flow agreement parameter transitions from favourable to unfavourable. For the pitching swimmer, the follower is mostly located within the red region where the flow agreement parameter is favourable. This effectively means that a higher surface area of the pitching follower experiences a flow field favourable to its motion, whereas part of the heaving follower undergoes negative flow agreement. This mechanism could be responsible for the increased efficiency of the pitching formation in comparison with the heaving formation.
Figure 6. Flow agreement parameter and thrust parameter experienced by a hypothetical point follower undergoing prescribed oscillations in the wake of a single swimmer that does not see the follower. Top row shows flow agreement parameter field in the wake of (a) heaving and (b) pitching swimmers. The grey plates represent the steady-state position of the followers in the first, second and third stable spacings found from solving the system with pairwise interactions (figure 3). In both cases, the distance of the regions with maximum flow agreement from the leading plate is very close to integer multiples of the wavelength (
$d_{{{h}}, {{p}}}/\lambda = 1, 2,3$). Bottom row shows thrust parameter as a function of distance. The dashed lines represent the head and tail positions of the heaving and pitching follower, respectively. The negative slopes of plot at the steady-state distances imply linear stability of the follower to in-line perturbations. The prescribed amplitudes are 
$A_{{h}} = 0.3$ and 
$A_{{p}} = 15^{\circ }$.
We next examine the stability of pairwise formation in the context of the simpler model based on the wake of a single swimmer and a hypothetical follower. We specifically consider the case where the virtual follower is positioned in line behind the single swimmer. It is well established that the thrust of a self-propelled flapping swimmer scales with the square of the swimmer's lateral velocity relative to the surrounding fluid's velocity (Triantafyllou et al. Reference Triantafyllou, Triantafyllou and Grosenbaugh1993; Floryan et al. Reference Floryan, Van Buren, Rowley and Smits2017; Newbolt et al. Reference Newbolt, Zhang and Ristroph2019). We thus define the thrust parameter
\begin{equation} \textrm{Thrust parameter:} \ \mathbb{X}= \frac{1}{T} \int_{t_s} ^{t_s+T}(v - \dot{y} )^{2}\, {\textrm{d}}t,\end{equation}
which acts as a measure of the period-average thrust but not an exact value of thrust. We plot the thrust parameter as a function of the follower's downstream location of heaving and pitching swimmers in the bottom row of figure 6. The thrust parameter is minimum at three locations where the flow agreement parameter is maximum. This is due to the fact that higher agreement between the follower's oscillation and the flow implies smaller difference in the follower's lateral speed relative to the flow and therefore smaller thrust. Superimposed onto these plots are the three equilibria at steady state obtained from our pairwise simulations in figure 3 (vertical dashed lines). We next argue that the slope of the thrust parameter at these locations is an indicator of the stability of the pairwise formation. To this end, recall that at steady state, the thrust 
$F_x(x_o,t)$ and skin drag 
$D(x_o, t)$ balance each other on average, and the follower experiences zero net acceleration. Namely, 
$\langle F_x(x_o) \rangle - \langle D(x_o) \rangle = 0$, where the time-average notation 
$\langle (\cdot ) \rangle = (1/T) \int _{t_s}^{t_s+T} (\cdot )\, {\textrm {d}}t$ is introduced for brevity. If we perturb the horizontal position of the follower by 
$\delta x$, since skin drag depends only on the relative fluid's velocity tangential to the plate, it is reasonable to assume that its change due to in-line positional perturbations is negligible 
$\langle D(x_o + \delta x) \rangle \approx \langle D(x_o) \rangle$. We arrive at the period-average equation 
$\langle F_x(x_o + \delta x)\rangle - \langle F(x_o) \rangle = m \langle \delta \ddot {x}\rangle$. This equation provides a condition for the linear stability of the pairwise formation in the context of the (single swimmer/virtual follower) model: if the slope of the period-average thrust relative to the horizontal position is negative, the system is linearly stable to perturbations in the horizontal position. Otherwise, the perturbation grows and the pair leaves their relative spacing at steady state. Since the thrust parameter 
$\mathbb {X}$ is an approximation of period-average thrust, it suffices to obtain the slope of 
$\mathbb {X}$ with respect to 
$\delta x$ to gauge the stability of the formation. The slope is negative at the steady-state positions in both heaving and pitching swimmers (bottom row of figure 6). Further, the slope of these locations decreases as the distance between the two swimmers increases. This is consistent with figure 3 where the third stable position is less stable than the second and the second slightly less stable than the first. Finally, the significantly higher slope of the thrust parameter in pitching compared with heaving is consistent with the observations in figure 3, where pitching formations are more stable.
6. Conclusion
We analysed the locomotion dynamics of actively flapping swimmers interacting passively via the fluid medium in the context of the vortex sheet model. Within the two-swimmer model, we showed that hydrodynamic interactions lead to stable ordered formations, in which the follower falls into specific positions in the wake of the leader, and the pair travel together at the same speed. This well-ordered ‘schooling’ behaviour occurs for both heaving and pitching swimmers. Group cohesion is tighter and more stable for pitching swimmers. In heaving alone, the school swims slightly faster compared with swimming alone, approximately 5 %–10 % faster, albeit at a similar increase in cost of transport, especially for the follower (approximately 20 % higher cost for the follower and 15 % for the formation). When pitching, the school swims at a slightly (approximately 5 %) lower speed but has significant energetic benefits, with up to 35 % reduction in cost of transport for the formation and up to 75 % for the follower. Simultaneous heaving and pitching also leads to flow-mediated stable formations (see supplemental movie), indicating that this phenomenon is robust to the flapping mode.
Detailed comparison of our findings to previously known results are in order. Physical experiments and numerical simulations report stable pairwise formations in hydrodynamically interacting swimmers. The experiments of Ramananarivo et al. (Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016) using pairs of purely heaving rigid foils in tandem found that the foils stabilize at particular discrete gap distances, and that these formations were usually accompanied by an increase in the swimming speed of the pair (10 %–20 % compared with swimming alone). The increase in speed was observed up to three wavelengths away from the leader, however, its effect quickly diminished with distance. Numerical simulations of pairs of interacting flapping swimmers provided more details on swimming energetics. Zhu et al. (Reference Zhu, He and Zhang2014) used an immersed boundary method to study the dynamics of two flexible filaments undergoing heaving oscillations at their leading edges at Reynolds number 
$= 200$. They reported an increase in both the swimming speed and input power of the pair compared with swimming alone. These changes were only reported for pairs in compact configurations. In this configuration, the leading edge of the follower is almost touching the trailing edge of the leader and the narrow space between them causes the pair to behave like one long filament. The increased speed and power requirements seemed to completely disappear for pairs in regular configurations characterized by an increased distance between the swimmers and velocity and power equal to a single swimmer. Dai et al. (Reference Dai, He, Zhang and Zhang2018) studied the swimming dynamics of multiple flexible filaments under combined pitching and heaving motions at the leading edge. However, the heaving motion's amplitude was much smaller than the tail's displacement due to pitching. For two filaments swimming in tandem, they reported a decrease of about 
$18\,\%$ in the cost of transport when the swimmers were in compact configurations. The regular configurations was found to be energetically beneficial, but only by approximately 2 %–3 % compared with swimming alone. Park & Sung (Reference Park and Sung2018) also found a decrease of approximately 
$15\,\%$ in power for a pair of flexible filaments, when swimming close to one another. The increase in speed relative to swimming alone was found to be negligible.
We examined pairwise interactions of purely heaving and pitching rigid swimmers, thus isolating heaving from pitching as opposed to the studies of flexible heaving filament that combine both effects. We found that for each flapping mode, the swimmers reach stable steady-state formations with constant distances. The flapping mode had a significant impact on the stability and swimming energetics of the pair. We observed a slight increase in the swimming speed of the heaving pair (up to 
$10\,\%$) at the expense of higher cost of transport. For pitching swimmers, the swimming speed was not affected much by the pairwise interaction, but we found a significant decrease in the input power of the follower (up to 
$70\,\%$ for small amplitudes). In contrast to the findings of Zhu et al. (Reference Zhu, He and Zhang2014) and Park & Sung (Reference Park and Sung2018), where the effects of the pairwise interactions quickly vanished with increasing distance, our vortex sheet model observed these effects at longer distances, up to three swimming wavelengths, consistent the experiments of Ramananarivo et al. (Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016). The discrepancy is most likely due to the relatively small Reynolds number in Zhu et al. (Reference Zhu, He and Zhang2014) (
$Re = 200$), causing the wake-induced flow to diffuse faster due to higher viscous forces. Ramananarivo et al. (Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016) reported a much larger Reynolds number (
$Re = 10^{3} - 10^{4}$) in their experimental set-up. The higher Reynolds numbers in the experiments are consistent with our inviscid model. At this inviscid regime, the flow inertia is dominant, causing the wake of the leader to live longer in the fluid. In this regime, the effects of hydrodynamic interaction on stability and energetics decreased with distance, but much more gradually.
In sum, our results are consistent with numerical and experimental findings of heaving foils (Zhu et al. Reference Zhu, He and Zhang2014; Becker et al. Reference Becker, Masoud, Newbolt, Shelley and Ristroph2015; Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016; Park & Sung Reference Park and Sung2018; Peng et al. Reference Peng, Huang and Xi-Yun2018; Newbolt et al. Reference Newbolt, Zhang and Ristroph2019; Lin et al. Reference Lin, Wu, Zhang and Yang2020), but go beyond these results in two major ways. Firstly, we completely separated the flapping modes, heaving and pitching, probed the effect of each one on the stability, speed and energetic performance of the school, and we showed that the flapping mode affects the tightness and stability of the formation, as well as the cost of transport in school compared with swimming alone. Secondly, we analysed these formations in the context of a simpler model consisting of the wake of a single swimmer and a hypothetical point follower. We defined an empirical flow agreement parameter and showed that regions where the wake-induced flow and the follower's periodic motion agree are consistent with the stable formations observed in pairwise interactions of heaving and pitching swimmers. The reduced-order model also highlights that the heaving mode is less favourable energetically because, in steady-state formations of heaving swimmers, the follower is positioned such that it experiences negative agreement with the ambient flow. We also employed the simpler model to make predictions about the stability of the pairwise formation, consistent with our findings that the pitching mode leads to tighter and more stable formation. Indeed, an alternative interpretation of our results is that they reveal how active changes in the flapping mode can be used to control, via hydrodynamic interactions, the school emergent properties, including the school speed, energetics and cohesion. For example, to save energy or quickly overcome large perturbations, swimmers can adopt a pitching mode.
These findings could be instrumental for understanding the role of the fluid medium as a mediator of the physical interactions between swimming fish, and to assess the hydrodynamic benefits to fish schooling. Fish have more complex flapping motions than simple heaving and pitching (Lin, Wu & Zhang Reference Lin, Wu and Zhang2019; Van Buren, Floryan & Smits Reference Van Buren, Floryan and Smits2019; Ayancik, Fish & Moored Reference Ayancik, Fish and Moored2020), and the compliance of the fish body is believed to play an important role in the flapping efficiency and ability to extract energy from ambient flows (Beal et al. Reference Beal, Hover, Triantafyllou, Liao and Lauder2006; Lucas et al. Reference Lucas, Johnson, Beaulieu, Cathcart, Tirrell, Colin, Gemmell, Dabiri and Costello2014; Quinn, Lauder & Smits Reference Quinn, Lauder and Smits2014; Jusufi et al. Reference Jusufi, Vogt, Wood and Lauder2017). These considerations, as well as extensions to arrays of swimmers in tandem and side by side, potentially flapping at different amplitudes and phases as in Newbolt et al. (Reference Newbolt, Zhang and Ristroph2019), will be treated in future works.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2021.551.
Acknowledgements
The authors would like to thank M.J. Shelley and L. Ristroph for useful conversations.
Funding
This work is partially supported by the National Science Foundation grant CBET 21-00705, the Office of Naval Research grants 12707602 and N00014-17-1-2062 and the Army Research Office grant W911NF-16-1-0074.
Declaration of interest
The authors report no conflict of interest.
Appendix A. Vortex sheet model
The coupled fluid–structure interaction between the swimming plate and the surrounding fluid is simulated using an inviscid vortex sheet model. In viscous fluids, boundary layer vorticity is formed along the sides of the swimmer, and it is swept away at the swimmer's tail to form a shear layer that rolls up into vortices. In the vortex sheet model, the swimmer is approximated by a bound vortex sheet, denoted by 
$l_{{b}}$, whose strength ensures that no fluid flows through the rigid plate, and the separated shear layer is approximated by a free regularized vortex sheet 
$l_{{w}}$ at the trailing edge of the swimmer. The total shed circulation 
$\varGamma$ in the vortex sheet is determined so as to satisfy the Kutta condition at the trailing edge, which is given in terms of the tangential velocity components above and below the bound sheet and ensures that the pressure jump across the sheet vanishes at the trailing edge.
To express these concepts mathematically, it is convenient to use the complex notation 
$z = x + {\textrm {i}}y$, where 
${\textrm {i}} = \sqrt {-1}$ and 
$(x ,y )$ denote the components of an arbitrary point in the plane. The bound vortex sheet 
$l_{{b}}$ is described by its position 
$z_{{b}}(s,t)$ and strength 
$\gamma (s,t)$, where 
$s\in [-l,l]$ denotes the arc length along the sheet 
$l_{{b}}$. The separated sheet 
$l_{{w}}$ is described by its position 
$z_{{w}}(\varGamma ,t)$, 
$\varGamma \in [0,\varGamma _{{w}}]$ where 
$\varGamma$ is the Lagrangian circulation around the portion of the separated sheet between its free end in the spiral centre and the point 
$z_{{w}}(\varGamma ,t)$. The parameter 
$\varGamma$ defines the vortex sheet strength 
$\gamma ={\textrm {d}} \varGamma /{\textrm {d}}s$.
By linearity of the problem, the complex velocity 
$w(z,t) = u(z,t) - \mathrm {i}v(z,t)$ is a superposition of the contributions due to the bound and free vortex sheets
\begin{equation} w(z,t)= w_{{b}}(z,t) + w_{{w}}(z,t).\end{equation}
In practice, the free sheet 
$l_{{w}}$ is regularized using the vortex blob method to prevent the growth of the Kelvin–Helmholtz instability. The bound sheet 
$l_{{b}}$ is not regularized in order to preserve the invertibility of the map between the sheet strength and the normal velocity along the sheet. The velocity components 
$w_{{b}}(z,t)$ and 
$w_{{w}}(z,t)$ induced by the bound and free vortex sheets, respectively, are given by
\begin{equation} w_{{b}}(z,t) = \int_{{-}l}^{l} K_o(z-z_{b}(s,t))\gamma(s,t)\, {\textrm{d}}s, \quad w_{{w}}(z,t) = \int_0^{\varGamma_{{w}}}K_{\delta}(z-z_{{w}}(\varGamma,t))\, {\textrm{d}} \varGamma,\end{equation}
where 
$K_\delta$ is the vortex blob kernel, with regularization parameter 
$\delta$,
\begin{equation} K_\delta (z) = \frac{1}{2{\rm \pi} \mathrm{i}}\frac{\bar{z}}{|z|^{2}+\delta^{2}}, \quad \bar{z} = x - {\textrm{i}}y. \end{equation}
If 
$z$ is a point on the bound sheet for which 
$\delta = 0$, 
$w_{{b}}$ is to be computed in the principal value sense.
The position of the bound vortex sheet 
$z_{{b}}$ is determined from the plate's flapping 
$(y(t),\theta (t))$ and swimming 
$x(t)$ motions. The corresponding sheet strength 
$\gamma (s,t)$ is determined by imposing the no penetration boundary condition on the plate, together with conservation of total circulation. Let 
$n(s,t)=-\sin \theta +{\textrm {i}} \cos \theta$ be the upward normal to the plate, the no-penetration boundary condition is given by
\begin{equation} {\textrm{Re}} [ w n]_{z_b}= {\textrm{Re}} [ w_{swimmer} n], \end{equation}where
\begin{equation} w_{swimmer} = \dot{x} - {\textrm{i}} \dot{y}-\mathrm{i}\dot{\theta}[\bar{z}_b - ({x}- {\textrm{i}} {y})]. \end{equation}
Conservation of the fluid circulation implies that 
$\int _{l_b} \gamma (s,t) \, {\textrm {d}}s + \varGamma _{{w}} (t) = 0$.
The circulation parameter 
$\varGamma$ along the free vortex sheet 
$z_{{w}}(\varGamma ,t)$ is determined by the circulation shedding rates 
$\dot \varGamma _{{w}}$, according to the Kutta condition, which states that the fluid velocity at the trailing edge is finite and tangent to the flyer. The Kutta condition can be obtained from the Euler equations by enforcing that, at the trailing edge, the difference in pressure across the swimmer is zero. To this end, we integrate the balance of momentum equation for inviscid planar flow along a closed contour containing the vortex sheet and trailing edge,
\begin{equation} [p]_\mp(s) = p_-(s) - p_+(s) ={-} \frac{{\textrm{d}}\varGamma(s,t)}{{\textrm{d}}t} - \frac{1}{2}(u_-^{2} - u_+^{2}), \end{equation}
where 
$\varGamma (s,t) = \varGamma _{{w}}+\int _{-l}^{s} \gamma (s',t)\, {\textrm {d}}s'$, 
$-l \le s \le l$, is the circulation within the contour and 
$p_{\mp }(s,t)$ and 
$u_\mp (s,t)$ denote the limiting pressure and tangential slip velocities on both sides of the swimmer. Since the pressure difference across the free sheet is zero, it also vanishes at the trailing edge by continuity, which implies that
\begin{equation} \dot{\varGamma}_{{{w}}}={-}\tfrac{1}{2}(u_-^{2} - u_+^{2})|_{s={-} l}.\end{equation}
The values of 
$u_-$ and 
$u_+$ are obtained from the average tangential velocity component and from the velocity jump at the trailing edge, given by the sheet strength, evaluated at 
$s=- l$
\begin{equation} \bar{u}=\frac{u_+{+} u_-}{2}={\textrm{Im}}[(w-w_{swimmer})n], \quad u_-{-}u_+{=} \gamma .\end{equation}
Once shed, the vorticity in the free sheet moves with the flow. Thus the parameter 
$\varGamma$ assigned to each particle 
$z_{{w}}(\varGamma ,t)$ is the value of 
$\varGamma _{{w}}$ at the instant it is shed from the trailing edge. The evolution of the free vortex sheet 
$z_{{w}}$ is obtained by advecting it in time with the fluid velocity,
\begin{equation} \dot{\bar{z}}_{{{w}}} = w_{{w}}(z_{{{w}}},t) +w_{{b}}(z_{{{w}}},t).\end{equation}Appendix B. Forces and moments
The hydrodynamic force acting on the swimmer due to the pressure difference across the swimmer is given by,
\begin{equation} \int_{l_b} n[p]_\mp \, {\textrm{d}}s ={-}F \sin\theta + {\textrm{i}}F \cos\theta ,\end{equation}
where 
$F = \int _{l_b} [p]_\mp \, {\textrm {d}}s$. The hydrodynamic moment acting on the swimmer about its leading edge is given by
\begin{equation} M = {\textrm{Re}}\left[ \int_{l_b}\mathrm{i}\bar{n}(z_{le}-z_b)[p]_\mp \, {\textrm{d}}s \right], \end{equation}
where 
$z_{le}$ is position of the leading edge 
$s = \pm l$.
It is known that the strength of the bound vortex sheet exhibits an inverse square root singularity at the edges (Saffman Reference Saffman1992; Eldredge Reference Eldredge2019). The singularity at the trailing edge is regularized by enforcing the Kutta condition as discussed above. To regularize the singularity at the leading edge, we introduce a force parallel to the plate known as leading edge suction (Eldredge Reference Eldredge2019). Following the derivation provided in Eldredge (Reference Eldredge2019), we write the suction force, in dimensionless form as
\begin{equation} S = 2 {\rm \pi}\,\textrm{e}^{\textrm{i} \theta} \sigma^{2}, \end{equation}
where 
$\sigma$ is the suction parameter defined as
\begin{equation} \sigma = \frac{1}{2}(\dot{y} - l\dot{\theta} \cos \theta ) + \int_{l_b} \frac{\gamma (s, t)}{2 {\rm \pi}l} \textrm{Re} \left( \frac{ \tilde{z}(s, t) + l}{ \tilde{z}(s, t) - l} \right) ^{1/2} \, {\textrm{d}}s, \end{equation}
where 
$\tilde {z}(s,t) = z(s,t) - z \, \textrm {e}^{\textrm {i} \theta }$, is the complex position of any vortex sheet present in the fluid written in the plate's frame of reference. 
$\dot {y} - l\dot {\theta } \cos \theta$ is the velocity of the centre of the plate in the 
$y$-direction. Note that in (B3), the suction force is always positive (always a thrust force) and parallel to the plate.
Note that the majority of the suction force is due to the vertical motion of the leading edge relative to the surrounding fluid. For the pitching swimmer, since the leading edge has no vertical motion, the contribution of the leading edge suction force to the total thrust force of the swimmer is negligible. This is confirmed by our numerical experiments on a single pitching swimmer.
Last, we introduce a drag force 
$D$ that emulates the effect of skin friction due to fluid viscosity. This force is based on the Blasius laminar boundary layer theory as implemented by Fang (Reference Fang2016) in the context of the vortex sheet model. Blasius theory provides an empirical formula for skin friction on one side of a horizontal plate of length 
$2l$ placed in fluid of density 
$\rho _f$ and uniform velocity 
$U$. In dimensional form, Blasius formula is 
$D = \frac {1}{2} \rho _{{f}} (2l) (c_{{f}})U^{2}$, where the skin friction coefficient 
$C_{{f}} = {0.664}/{\sqrt {Re}}$ is given in terms of the Reynolds number 
$Re = {\rho _f U (2l)}/{\mu }$. Substituting back in the empirical formula leads to 
$D = C_d U^{3/2}$, where 
$C_d = 0.664 \sqrt { \rho _f \mu (2l)}$. Following Fang (Reference Fang2016), we write a modified expression of the drag force for a swimming plate
\begin{equation} D = C_d (\bar{U}_+ ^{3/2} + \bar{U}_- ^{3/2}), \end{equation}
where 
$\bar {U}_{\pm }$ are the spatially averaged tangential fluid velocities on the upper and lower side of the plate, respectively, relative to the swimming velocity 
$U$,
\begin{equation} \bar{U}_{{\pm}}(t) = \frac{1}{2l} \int_{{-}l} ^{l} u_{{\pm}} (s,t)\, {\textrm{d}}s - U. \end{equation}
We estimate 
$C_d$ to be approximately 
$0.02$ in the experiments of Ramananarivo et al. (Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016).
Appendix C. Numerical implementation
The bound vortex sheet is discretized by 
$2n+1$ point vortices at 
$z^{b}(t)$ with strength 
${\rm \Delta} \varGamma =\gamma {\rm \Delta} s$. These vortices are located at Chebyshev points that cluster at the two ends of the swimmer. Their strength is determined by enforcing no penetration at the midpoints between the vortices, together with conservation of circulation. The free vortex sheet is discretized by regularized point vortices at 
$z^{{{w}}}(t)$, that is released from the trailing edge at each timestep with circulation given by (A7). The free point vortices move with the discretized fluid velocity while the bound vortices move with the swimmer's velocity. The discretization of (2.2) and ((A7) and (A9)) yields a coupled system of ordinary differential evolution equations for the swimmer's position, the shed circulation and the free vorticity, that is integrated in time using the fourth-order Runge–Kutta scheme. The details of the shedding algorithm are given in Nitsche & Krasny (Reference Nitsche and Krasny1994). The numerical values of the timestep 
${\rm \Delta} t$, the number of bound vortices 
$n$ and the regularization parameter 
$\delta$ are chosen so that the solution changes little under further refinement.
Finally, to emulate the effect of viscosity, we allow the shed vortex sheets to decay gradually by dissipating each incremental point vortex after a finite time 
$T_{diss}$ from the time it is shed into the fluid. Larger 
$T_{diss}$ implies that the vortices stay in the fluid for longer times, mimicking the effect of lower fluid viscosity. For the results depicted in this study, we used 
$T_{diss} \in [1.5, 3.5]$ flapping period. We refer the reader to Huang et al. (Reference Huang, Ristroph, Luhar and Kanso2018) for a detailed analysis of the effect of dissipation time on the hydrodynamic forces on a stationary and moving plate in the vortex sheet model. Details of the numerical validation in comparison to Jones (Reference Jones2003) and Jones & Shelley (Reference Jones and Shelley2005) are provided in Huang et al. (Reference Huang, Nitsche and Kanso2016).
Appendix D. Swimming energetics
Heaving motions are produced by an active heaving force 
$F_{h}$ acting by the swimmer on the fluid in the 
$y$-direction. The value of 
$F_{h}$ is obtained from the balance of linear momentum on the swimmer in the 
$y$-direction,
\begin{equation} \textrm{Heaving:} \ m\ddot{y} = F_y + F_{h} . \end{equation}
Here, the hydrodynamic force 
$F_y$ acting on the swimmer in the 
$y$-direction is given by (B1).
Figure 7. (a) Schematic of the vortex sheet model for a two-dimensional flapping swimmer. (b) Depiction of the different hydrodynamic forces acting on the swimmer.
Figure 8. Hydrodynamic forces on the follower act as restoring forces. Snapshots of pairs of swimmers undergoing (a) heaving and (b) pitching motion during transient and steady-state formation. Green (thrust) and orange (drag) arrows represent period-averaged hydrodynamic forces acting on the follower. Right columns in (a,b) show the instantaneous thrust and skin drag (solid lines) and their period-averaged values (dashed lines) over one flapping period during transient and steady-state formation. The grey dashed lines denote the time instance of the snapshots shown to the left.
Pitching motions are produced by an active moment 
$M_{p}$ acting by the swimmer on the fluid about the leading edge. The value of 
$M_{p}$ is obtained from the balance of angular momentum about the swimmer's leading edge,
\begin{equation} \textrm{Pitching:} \ I \ddot{\theta} - {\textrm{Im}}[m(\dot{x}+{\textrm{i}}\dot{y})w_{l.e.}] = M + M_{p}, \end{equation}
Here, 
$I=m(2l)^{2}/3$ is the swimmer's moment of inertia about the leading edge, 
$w_{l.e.}$ is the swimmer's velocity at the leading edge, and 
$M$ is the hydrodynamic moment about the leading edge given in (B2).
The power input by the swimmer into the fluid due to heaving and pitching motions, respectively, is given by
$$\begin{gather} \textrm{Heaving:} \ P_{{h}} = F_{h} \dot{y}, \end{gather}$$
$$\begin{gather}\textrm{Pitching:} \ P_{{p}} = M_{p} \dot{\theta}. \end{gather}$$Note that, in both cases, the leading edge suction and skin drag forces do not contribute to the input power.
Appendix E. Hydrodynamic interaction forces
To further examine the hydrodynamic interactions between the two swimmers, we plot snapshots of the free vortex sheets and flow field for the pair of heaving and pitching swimmers in figure 8(a,b). We report two instances taken during the transient and steady-state motion, and from each regime, we report the hydrodynamic thrust and skin drag over one period of flapping and their time-period average. When the follower gets too close to the leader, the drag force dominates over thrust, causing the follower to decelerate and move further behind the leader. Conversely, when the distance between swimmers is larger than the steady-state spacing, the thrust force overcomes drag causing the follower to accelerate and the pair to move closer; see, e.g. top right of figures 8(a) and 8(b), respectively. The thrust and drag forces on the follower are balanced on average after steady state has been reached, effectively leading to zero acceleration and constant separation distance between the two swimmers.
