1. Introduction
A common locomotion strategy adopted by aquatic or flying animals (Gray Reference Gray1933; Shyy et al. Reference Shyy, Aono, Chimakurthi, Trizila, Kang, Cesnik and Liu2010; Wu Reference Wu2010), and more recently employed in the conception of large- and small-scale artificial swimmers (Barrett Reference Barrett1996; Ramananarivo, Godoy-Diana & Thiria Reference Ramananarivo, Godoy-Diana and Thiria2013; Williams et al. Reference Williams, Anand, Rajagopalan and Saif2014), is the flapping motion of appendages such as wings, tails and fins. A fundamental question that impacts the design of micro-swimmers and aerial vehicles (Williams et al. Reference Williams, Anand, Rajagopalan and Saif2014; Faux et al. Reference Faux, Thomas, Cattan and Grondel2018) is the critical size above which flapping-based propulsion remains efficient and applicable. Indeed, micro-organisms of very small scales, as cells or sperm, are known to exploit other locomotion strategies (Lauga Reference Lauga2011), ciliar or flagelar propulsion, respectively. As first stated by Purcell (Reference Purcell1977) in the so-called scallop theorem, a reciprocal motion of appendages, for which the paths during the two half-strokes are identical but time-reversal, does not allow a net thrust to be generated at those very small scales. This is due to the linearity and timeless nature of the surrounding flows which are entirely dominated by the viscous effects. An emblematic observation of the transition from ciliar to flapping propulsion has been achieved for the mollusc Clione Antartica (Childress & Dudley Reference Childress and Dudley2004) that disposes both of cilia and wings attached to its body. Whereas its cilia are always employed, the wings remain retracted to its body, being flapped exclusively after a critical velocity is reached. This switch of locomotion strategy was related to the evolution of the dynamical response of the surrounding flow as the Reynolds number, based on the swimming velocity, increases. The present paper aims at better understanding the emergence of flapping locomotion based on reciprocal motion.
To that aim, Vandenberghe, Zhang & Childress (Reference Vandenberghe, Zhang and Childress2004); Vandenberghe, Childress & Zhang (Reference Vandenberghe, Childress and Zhang2006) designed an experiment where flapping propulsion emerges exclusively from the flow and not from the motion asymmetry, which is generally explored in the animal world as to achieve a more efficient propulsion (Weis-Fogh Reference Weis-Fogh1973; Spagnolie et al. Reference Spagnolie, Moret, Shelley and Zhang2010). The experiment consists of a horizontal flat rectangular foil immersed in a cylindrical tank filled with still water and attached in its mid-span to a shaft. This shaft is vertically flapped with a sinusoidal motion and the foil is allowed to rotate, together with the shaft, around the vertical axis in the horizontal direction. Note that the foil is only heaving, not simultaneously heaving and pitching like in the experiments of Spagnolie et al. (Reference Spagnolie, Moret, Shelley and Zhang2010) where the foil besides the imposed heaving was allowed to passively pitch around its leading edge. For a small enough frequency, the flow induced by a heaving motion of fixed amplitude is left–right symmetric. Thus, no hydrodynamic force is generated over the foil in the horizontal direction for every instants, and the foil does not rotate. However, once a critical Stokes number 
$\beta =f^*(c^*)^2/\nu$ is attained (a non-dimensional number similar to the Reynolds number that uses the dimensional flapping frequency 
$f^*$ and the foil chord 
$c^*$ as characteristic time and length scales, as well as the fluid kinematic viscosity 
$\nu$), the surrounding flow breaks its initial symmetry and generates horizontal forces. The foil then achieves locomotion and eventually reaches a permanent forward regime in equilibrium with the fluid. Subsequent numerical studies were dedicated to understand how the transient dynamics and the self-propelled regimes of this model problem evolve with respect to its control parameters. These studies simplified Vandenberghe et al. (Reference Vandenberghe, Zhang and Childress2004) configuration, working with a two-dimensional cross-section of the experiment (imposed heaving and horizontally self-propelled foils), thus neglecting its rotational flow effects. Investigating the self-propulsion of elliptical foils in a two-dimensional incompressible flow under a fixed non-dimensional chord-based flapping amplitude 
$A=0.5$, Alben & Shelley (Reference Alben and Shelley2005) revealed that as the flapping frequency (equivalently the Stokes number) is increased, the foil motion transition between different self-propelled regimes that are a unidirectional propulsion (as in the experiments of Vandenberghe et al. Reference Vandenberghe, Zhang and Childress2004), a quasi-periodic back and forth motion around a fixed point in space (with a frequency remarkably lower than the flapping one) and even a chaotic motion. These authors have also shown that these self-propelled regimes are greatly impacted by the thickness-to-chord aspect ratio 
$h$ and the solid-to-fluid density ratio 
$\rho$. On one hand, thinner ellipses of aspect ratio 
$h=0.1$ are able to break symmetry at lower flapping frequencies and present thus a greater exponential growth of their horizontal velocity than thicker foils. On the other hand, for foils of greater density ratio (
$\rho >10$), the existence of non-coherent and chaotic regimes is greatly reduced or even suppressed when compared with their lighter equivalents. Lu & Liao (Reference Lu and Liao2006), for instance, have shown for a fixed flapping amplitude and frequency that a non-coherent state of motion can be suppressed thanks to the increase of the density ratio. A similar observation was made by Zhang et al. (Reference Zhang, Ni, Wang and He2009) while decreasing the aspect ratio of elliptical and rectangular foils, where in both cases smaller aspect ratios were found to be more prone to unidirectional locomotion than thicker ones. Lu & Liao (Reference Lu and Liao2006) have equally mapped in the plan flapping amplitude/frequency the transition between symmetric and non-symmetric flows (thus propelled ones) revealing that the transition occurs for smaller frequencies at higher flapping amplitudes. Later on, Deng & Caulfield (Reference Deng and Caulfield2016) established the same frontier for different aspect ratios, revealing that thinner foils break symmetry earlier in flapping amplitude and frequency. The authors also compared the frontier between symmetric and non-symmetric flows for propelled and non-propelled foils for two-dimensional ellipses (Deng & Caulfield Reference Deng and Caulfield2016) or three-dimensional oblate spheroids (Deng et al. Reference Deng, Teng, Caulfield and Mao2016, Reference Deng, Xue, Mao and Caulfield2017; Deng & Caulfield Reference Deng and Caulfield2018), indicating that the frontier of flow symmetry breaking is obtained in both two- and three-dimensional cases for smaller frequency and amplitude for self-propelled foils. Using two-dimensional numerical simulations, we revisit in this work the nonlinear regimes of locomotion for a thin rectangular foil (
$h^*=0.05 c^*)$ of density 
$\rho _{s}= 100 \rho _{f}$ and flapping with a fixed maximal amplitude 
$A^{*}=0.5 c^{*}$. By varying the Stokes numbers in the range 
$1 \leq \beta \leq 20$, we first aim at carefully characterizing and identifying the transition between various self-propelled regimes of the foil motion: the non-propulsive, the unidirectional propulsive and the back and forth motion. This parametric investigation is useful to accurately determine critical values of the Stokes numbers for which transition between these nonlinear regimes are observed. They will be used as the basis of comparison for the second objective of this work, i.e. predicting the emergence of these regimes using linear stability analysis of the periodic flows generated by the flapping foil.
Floquet analysis allows us to investigate the linear stability of periodic solutions (Floquet Reference Floquet1883). In hydrodynamics Barkley & Henderson (Reference Barkley and Henderson1996) first performed this analysis on a two-dimensional time-periodic wake flow to explain the onset of three-dimensional structures in the wake of a fixed circular cylinder. The Floquet analysis of time-periodic flows generated by flapping bodies has then been considered by Elston, Sheridan & Blackburn (Reference Elston, Sheridan and Blackburn2004); Elston, Blackburn & Sheridan (Reference Elston, Blackburn and Sheridan2006) for two-dimensional oscillating cylinder flows. They successfully explained the emergence of two- and three-dimensional flow asymmetries observed in experiments and simulations. More recently, Jallas, Marquet & Fabre (Reference Jallas, Marquet and Fabre2017) performed a Floquet analysis of a time-periodic propulsive wake generated by a pitching wing. They identify an unstable synchronous mode that successfully explains the lateral deviation of the propulsive vortex street observed when increasing the flapping frequency. To investigate the emergence of the self-propelled regimes described above, Deng & Caulfield (Reference Deng and Caulfield2016, Reference Deng and Caulfield2018) and Deng et al. (Reference Deng, Teng, Caulfield and Mao2016, Reference Deng, Xue, Mao and Caulfield2017) first proposed to consider a Floquet analysis for various flapping foil configurations. Based on the observation that flow symmetry breaking occurs prior to the self-propulsion of the foil in temporal simulations, Deng & Caulfield (Reference Deng and Caulfield2016, Reference Deng and Caulfield2018) applied a purely hydrodynamic analysis, that does not consider a perturbation of the foil speed in the propulsion direction. For a certain range of control parameters, results of this purely hydrodynamics analysis are in good agreement with nonlinear results. In particular, they identify unstable asynchronous modes at flapping frequencies where the foil exhibits a slow quasi-periodic back and forth motion. However, some disagreements between results of the purely hydrodynamic analysis and the nonlinear results of self-propelled simulations were also reported, as in the case of low aspect ratio 
$h=0.1$ ellipses, where linear analysis fails to predict the onset of unidirectional forward locomotion (Deng & Caulfield Reference Deng and Caulfield2016). In the present work, we introduce the so-called fluid–solid Floquet analysis that considers the foil speed as a perturbation variable and takes into account the inherent coupling between the flow perturbation and the rigid motion of the foil at the perturbation level. We will demonstrate, by comparison with nonlinear results, that this fluid–solid coupling is essential to correctly represent and predict the emergence of self-propelled regimes.
The importance of the fluid–solid coupling in linear stability analysis has a long history in aeroelasticity (see the review by Dowell et al. Reference Dowell, Curtiss, Scanlan and Sisto1989) that investigates the infinitesimal motion of structures immersed in high-Reynolds-number flows. Fluid–solid stability analyses for lower Reynolds number flows are more recent. To our knowledge, Cossu & Morino (Reference Cossu and Morino2000) first performed the fluid-stability analysis of the steady wake cylinder flow to explain the sub-critical vortex-induced vibration of the cylinder when mounted on a spring. The path of bodies freely rising or falling in fluids under the effect of gravity (see Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012 for a review) is another example where fluid–solid linear stability analyses successfully explained the emergence of various trajectories. Tchoufag, Fabre & Magnaudet (Reference Tchoufag, Fabre and Magnaudet2014a) first elucidate the path instability of buoyancy-driven disks/thin cylinders and then of freely rising spheroidal bubble (Tchoufag, Magnaudet & Fabre Reference Tchoufag, Magnaudet and Fabre2014b). Recently, Negi, Hanifi & Henningson (Reference Negi, Hanifi and Henningson2019) proposed a simplified formulation to handle the linearized fluid–structure interaction for rigid bodies. Fluid–solid stability analysis has also been extended to deformable (elastic) structures, to explain the dynamics of inverted flags in uniform flows (Goza, Colonius & Sader Reference Goza, Colonius and Sader2018) and of flexible splitter plates clamped to the rear of a cylinder (Pfister & Marquet Reference Pfister and Marquet2020). Note also that Tammisola, Lundell & Söderberg (Reference Tammisola, Lundell and Söderberg2012) investigated the global instability of planar jets and wakes in two immiscible fluids, focusing on the effect of surface tension. In all of these studies, the temporal evolution of perturbations over steady base flow solutions was considered. To our knowledge, the fluid–solid stability analysis of time-periodic flow solutions has never been addressed in the context of fluid–solid interaction. In the present study, we introduce the mathematical formalism of such analysis and apply it to explain the emergence of self-propelled flapping states. Additionally, a time-averaged analysis is proposed to highlight the role of the fluid–solid coupling in the destabilization of the Floquet modes. Such connections between linear modes and thrust efficiency have been, for instance, highlighted in the literature as key factors for an optimal frequency selection in flapping wings (Triantafyllou, Triantafyllou & Grosenbaugh Reference Triantafyllou, Triantafyllou and Grosenbaugh1993; Moored et al. Reference Moored, Dewey, Smits and Haj-Hariri2012).
This article is organized in two parts. In § 2 we investigate numerically the nonlinear regimes of locomotion for a self-propelled heaving foil. The configuration and non-dimensional parameters are introduced before describing the governing nonlinear equations and numerical methods. The self-propelled solutions obtained for a fixed flapping amplitude and density ratio are then carefully described for three values of the Stokes number. The transition between regimes of non-propulsive, unidirectional propulsive and back and forth motions are finally identified by varying the Stokes number in the range 
$2 \le \beta \le 20$. In § 3 we introduce first the fluid–solid Floquet stability analysis of self-propelled foils and then the time-averaged analysis that allows us to establish instability criteria based on the velocity and force of the Floquet mode. Results of this fluid–solid analysis, performed at 
$\rho =100$ for symmetric non-propulsive solutions, are first described by analysing the synchronous and asynchronous modes found unstable at different Stokes numbers. Those results are then compared with those obtained with the purely hydrodynamic Floquet analysis and with the nonlinear temporal simulations. The effect of the density ratio on the two unstable Floquet modes is finally described.
2. Problem formulation and self-propelled nonlinear solutions
We investigate the horizontally constrained locomotion of a vertically heaving foil of density 
$\rho _s$ immersed in an initially quiescent fluid of density 
$\rho _f$ and viscosity 
$\nu$. The foil, shown in figure 1, is similar to the one used in the experimental studies by Vandenberghe et al. (Reference Vandenberghe, Zhang and Childress2004, Reference Vandenberghe, Childress and Zhang2006). Its rectangular shape is characterized by the thickness 
$h^*$ and chord 
$c^*$ with rounded corners of diameter equal to the foil thickness. The periodic displacement imposed along the vertical axis 
$e_y$ is
\begin{equation} y_g^{*}(t)={-} A^{*} \cos(2{\rm \pi} f^{*} t^{*}), \end{equation}
where the superscript 
$*$ is used to indicate dimensional variables. Here 
$A^{*}$ is the maximal vertical amplitude and 
$f^{*}$ is the flapping frequency, and 
$T^{*}=1/f^{*}$ is the flapping period. The foil is free to move along the horizontal axis 
$e_x$ as a result of hydrodynamic forces acting on the solid–fluid interface 
$\varGamma _w$. This rigid-body fluid–structure interaction is characterized by four non-dimensional parameters, namely the frequency-based Stokes number 
$\beta$, the non-dimensional amplitude 
$A$, the solid–fluid density ratio 
$\rho$ and the non-dimensional thickness 
$h$, defined respectively as
\begin{equation} \beta = \frac{f^{*} (c^*)^2}{\nu}, \quad A=\frac{A^{*}}{c^*},\quad \rho=\frac{\rho_s}{\rho_f} \quad \text{and}\quad h=\frac{h^*}{c^*}. \end{equation}
These parameters are obtained by choosing the chord 
$c$ as characteristic length scale, the fluid density 
$\rho _{f}$ as characteristic mass scale and the flapping period 
$1/f^*$ as characteristic time. In the following, all variables are thus made non-dimensional using these scales. Note that the non-dimensional flapping period 
$T$ is thus equal to 
$1$ whatever the values of the Stokes number 
$\beta$, which is the only parameter or variable containing the dependency to the dimensional frequency 
$f^{*}$. Other choices of characteristic scale are also possible and made in the literature. For instance, Alben & Shelley (Reference Alben and Shelley2005) used the flapping velocity 
$A^{*} f^{*}$ as characteristic velocity, thus introducing the flapping amplitude based Stokes number 
$\beta _A = A^{*} f^{*} c ^*/\nu = A \beta$.
Figure 1. Sketch of the foil configuration under sinusoidal vertical motion and horizontal translation with velocity 
$\boldsymbol {u}_g$. The foil chord (trailing to leading edge) and its thickness, the frames of reference and the solid/fluid interface 
$\varGamma _w$ are indicated.
In the present study, the foil geometry and the flapping amplitude are fixed to 
$h=1/20$ and 
$A=0.5$, respectively. This aspect ratio is close to the experimental devices of Vandenberghe et al. (Reference Vandenberghe, Zhang and Childress2004), and this flapping amplitude equals the one adopted by Alben & Shelley (Reference Alben and Shelley2005). A discussion of the influence of these two parameters can be found in Zhang et al. (Reference Zhang, Ni, Wang and He2009) and Deng & Caulfield (Reference Deng and Caulfield2016). In this section we will investigate numerically the nonlinear dynamics of the foil for the fixed density ratio 
$\rho =100$ in the range of Stokes number 
$1 \le \beta \le 20$.
2.1. Governing nonlinear equations
The dynamics of the foil interacting with the surrounding fluid is described by the non-dimensional variable 
$\boldsymbol {q}=(\boldsymbol {u},p,u_g)^{T}$, where 
$\boldsymbol {u}=(u,v)$ is the two-dimensional fluid velocity field, 
$p$ is the pressure field and 
$u_g$ is the foil horizontal velocity. The fluid–solid variable is governed by the evolution equation
\begin{equation} \boldsymbol{\mathsf{B}}\frac{\partial\boldsymbol{q}}{\partial t} = \boldsymbol{\mathsf{R}}(\boldsymbol{q},v_g), \end{equation}
where 
$v_g(t)=2 {\rm \pi}A \sin (2 {\rm \pi}t)$ is the non-dimensional foil vertical velocity and the operators 
$\boldsymbol{\mathsf{B}}$ and 
$\boldsymbol{\mathsf{R}}$ are defined as
\begin{equation} \boldsymbol{\mathsf{B}}=\begin{bmatrix} \mathcal{I} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} , \quad \boldsymbol{\mathsf{R}}(\boldsymbol{q},v_g)=\begin{bmatrix} - \left( [\boldsymbol{u}-\boldsymbol{u}_g]\boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{u} - \boldsymbol{\nabla} p + \beta^{{-}1} \Delta \boldsymbol{u} \\ -\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} \\ \displaystyle (\rho S)^{{-}1} \displaystyle F_x(\boldsymbol{u},p) \end{bmatrix} . \end{equation}
The first and second lines are the incompressible Navier–Stokes equations written in a non-inertial frame of reference, denoted 
$(G,\boldsymbol {e}_x,\boldsymbol {e}_y)$ in figure 1, that translates at the foil speed 
$\boldsymbol {u}_g=(u_g,v_g)$ in the laboratory frame of reference 
$(O,\boldsymbol {e}_X,\boldsymbol {e}_Y)$. Note that both solid and fluid velocities are absolute velocities (Mougin & Magnaudet Reference Mougin and Magnaudet2003; Jenny & Dušek Reference Jenny and Dušek2004), the relative flow velocity 
$(\boldsymbol {u}-\boldsymbol {u}_g)$ appearing in the nonlinear term of the momentum equations. While the solid vertical velocity 
$v_g$ is imposed, the temporal evolution of the foil horizontal velocity 
$u_g$ is governed by Newton's second law, as stated by the third line in (2.3) and (2.4a,b). The horizontal acceleration is equal to the horizontal hydrodynamic force 
$F_x(\boldsymbol {u},p)$ weighted by the non-dimensional mass of the foil 
$\rho S$ (
$S=h(1-h)+{\rm \pi} h^2/4$ being its non-dimensional surface). This hydrodynamic force depends on the fluid velocity and pressure as
\begin{equation} F_x = \int_{\varGamma_w}([{-}p \boldsymbol{\mathcal{I}} + \beta^{{-}1} (\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^T)] \boldsymbol{\cdot} \boldsymbol{n}) \boldsymbol{\cdot} \boldsymbol{e_x} \, \mathrm{d} \varGamma_w , \end{equation}
where 
$\varGamma _w$ denotes the fluid–solid boundary. An additional coupling between the fluid and solid variables is due to the equality of velocities at the fluid–solid interface, i.e.
\begin{equation} \boldsymbol{u}\left(\varGamma_w,t\right)=\boldsymbol{u}_g(t)= [u_g(t), 2{\rm \pi} A \sin(2{\rm \pi} t)]^{{T}}. \end{equation}The fluid is at rest sufficiently far away from the foil.
2.2. Numerical methods
The system of (2.3)–(2.6) is discretized in time using the following 
$r$-order semi-explicit scheme:
\begin{equation} \left. \begin{gathered} \frac{\alpha_0}{\Delta t} \boldsymbol{u}^{n+1} + \boldsymbol{\nabla}{p}^{n+1} -\beta^{{-}1} \Delta \boldsymbol{u}^{n+1} = \boldsymbol{f}^{n+1},\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}^{n+1} = 0,\\ (u^{n+1},v^{n+1})(\varGamma_{w}) = \big(u_g^{n+1}, 2 {\rm \pi}\,A \, \sin(2 {\rm \pi}t^{n+1})\big),\\ \frac{\alpha_0}{\Delta t} u_g^{n+1} = (\rho S)^{{-}1}F_x(\boldsymbol{u}^{n+1},p^{n+1}) - \sum_{k=1}^r \frac{\alpha_k}{\Delta t} u_{g}^{n+1-k} ,\\ (u^{n+1},v^{n+1})(\varGamma_{e}) = (0,0). \end{gathered} \right\} \end{equation}
Here the right-hand side forcing term 
$\boldsymbol {f}^{n+1}$ in the momentum equation is defined as
\begin{equation} \boldsymbol{f}^{n+1} ={-} \sum_{k=1}^r \gamma_{k} (\boldsymbol{u}^{n+1-k}- \boldsymbol{u_g}^{n+1-k})\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^{n+1-k} - \sum_{k=1}^r \frac{\alpha_k}{\Delta t} \boldsymbol{u}^{n+1-k}, \end{equation}
with 
$\Delta t$ the time step and 
$(\boldsymbol {u}^{n+1},p^{n+1})$ the velocity and pressure at time 
$t_{n+1}=(n+1) \Delta t$. A quiescent fluid condition is applied in the external boundary 
$\varGamma _e$ of the computational domain, typically far away from the foil. The time derivatives are approximated by 
$r$-order backward differential formulae. The linear diffusion and pressure gradient terms are implicit, while the nonlinear convection terms are extrapolated with 
$r$-order formulae. A first-order scheme (
$r=1$, 
$\alpha _0=1$, 
$\alpha _1=-1$ and 
$\gamma _1=1$) is used for the first two temporal iterations (
$n \le 1$), before switching in the subsequent iterations (
$n>1$) to a second-order scheme (
$r=2$, 
$\alpha _0=3/2$, 
$\alpha _1=-2$, 
$\alpha _2=1/2$, 
$\gamma _1=2$ and 
$\gamma _2=-1$). To avoid severe time-step restrictions for small values of density ratio induced by an explicit coupling (Causin, Gerbeau & Nobile Reference Causin, Gerbeau and Nobile2005), the equality of fluid and solid velocity is treated implicitly. To allow the use of an existing fast implementation to solve the flow equations Jallas et al. (Reference Jallas, Marquet and Fabre2017), we use a segregated approach, proposed by Jenny & Dušek (Reference Jenny and Dušek2004) and detailed in appendix A, to solve the coupled fluid–solid problem. Typically, the time step is set to 
$\Delta t=10^{-2}$ for small values of the Stokes number (
$\beta =2$) and is decreased to 
$\Delta t=5\times 10^{-4}$ for larger values (
$\beta =19$), so as to ensure the numerical stability of this semi-explicit temporal scheme (Kress & Lötstedt Reference Kress and Lötstedt2006).
The linear equations (2.7) are discretized in space using a classical finite-element method. The flow velocity is discretized with quadratic elements (
$P2$) while the pressure is discretized with linear element (
$P1$). The implementation is based on the FreeFEM software (Hecht Reference Hecht2012). The computational domain, displayed in figure 2(a), is a circle of (non-dimensional) diameter 
$60$ centred at the foil centre of mass, the external boundary of this circular domain being 
$\varGamma _{e}$. A Delaunay triangulation of the computational domain results in mesh with typically 
$1.2 \times 10^{4}$ triangles. As spatially symmetric solutions are expected, particular attention was given to create a symmetric mesh and not artificially insert asymmetries in the flow. To create a mesh that is symmetric with respect to the 
$x$- and 
$y$-axis and refined in flow regions exhibiting large velocity gradients (see figure 2b), we have proceed as follows. Once a first solution has been computed, we adapt a mesh of a quarter domain to several instants of the periodic flow, using the hessian-based mesh adaptation implemented in FreeFEM. We refer to Fabre et al. (Reference Fabre, Citro, Ferreira Sabino, Bonnefis, Sierra, Giannetti and Pigou2018) for a practical review. After duplicating and appropriately rotating this quarter mesh, the full mesh can finally be assembled. The triangle size is typically of order 
${O}(10^{-2})$ close to the foil, and 
$1$ in the external part of the domain. Mesh refinement and domain size were chosen based on the convergence of the foil horizontal velocity and the vertical hydrodynamic force. Greater domains or mesh refinement have exhibited little influence over these values. The validation of this numerical method is detailed in appendix B by comparison with results of Spagnolie et al. (Reference Spagnolie, Moret, Shelley and Zhang2010).
Figure 2. Computational domain and mesh. (a) Full and (b) close-up views of a typical mesh adapted to the flow solution.
2.3. Results
Unsteady nonlinear simulations are performed for values of the Stokes numbers in the range 
$1 \le \beta \le 20$. The amplitude 
$A=0.5$ and the foil aspect ratio 
$h=1/20$ are kept fixed throughout this study. When increasing the Stokes number, three different types of solution are successively observed, called hereinafter the (i) symmetric non-propulsive, (ii) unidirectional propulsive and (iii) back and forth solutions. In § 2.3.1 we first describe these solutions for three representative values of Stokes number and for a fixed density ratio 
$\rho =100$, concluding the section by a summary of the Stokes numbers range for which these solutions are obtained. In § 2.3.2 these results are compared with those of a smaller density ratio closer to aquatic swimming (
$\rho =1$). These different type of solutions have already been experimentally or numerically observed in previous studies. The transition from non-propulsive to unidirectional propulsive solutions was investigated in the works of Vandenberghe et al. (Reference Vandenberghe, Zhang and Childress2004, Reference Vandenberghe, Childress and Zhang2006), while back and forth solutions have been computed numerically in Alben & Shelley (Reference Alben and Shelley2005), Lu & Liao (Reference Lu and Liao2006) and Deng & Caulfield (Reference Deng and Caulfield2018). Self-propelled regimes presented in this section are thus not new but aim to establish the transition route for comparison with the linear Floquet stability analysis performed in the next section.
2.3.1. Self-propelled regimes for 
$\rho =100$
A typical solution obtained for small values of the Stokes number is displayed for 
$\beta =2$ in figure 3. The flow induced by the flapping foil inherits the spatial symmetry of the foil and the temporal symmetry of the imposed vertical motion (Elston et al. Reference Elston, Sheridan and Blackburn2004). It satisfies the 
$x$-reflection spatial symmetry in the non-inertial frame of reference, i.e.
\begin{equation} (u,v,p,\omega_{z})(x,y,t)=({-}u,v,p,-\omega_{z})({-}x,y,t), \end{equation}and the spatio-temporal symmetry
\begin{equation} (u,v,p,\omega_{z})(x,y,t)=(u,-v,p,-\omega_{z})(x,-y,t+T/2), \end{equation}
which is the combination of the 
$y$-reflection symmetry and the 
$T/2$ time-reciprocal translation. The vorticity 
$\omega _{z}$, used to display the solution, in figures 3(a–d) and 3(e–h) respectively at the inertial and non-inertial frames of reference for four equally spaced instants of the period 
$T$ is clearly an odd function of the 
$x$ variable for every time instant. Physically, the spatial symmetry is seen by the vortices of equal shape but different sign shed one each side of the foil during its vertical motion. The spatio-temporal flow symmetry is observed by the inversion of the vorticity sign in opposite foil strokes. A direct consequence of the spatial flow symmetry is the absence of instantaneous hydrodynamic forces acting in the horizontal direction, i.e. 
$F_x(t)=0$. Consequently, the foil is not accelerated in that direction and its velocity remains equal to zero; hence, the name of symmetric non-propulsive solution.
Figure 3. Symmetric solution for 
$\beta =2$. The vorticity field is depicted in the laboratory frame of reference 
$(X,Y)$ (a–d) and in the non-inertial frame of reference 
$(x,y)$ (e–h) for four equally spaced instants of the unitary period. The initial time 
$t_0$ corresponds to the lowest vertical position of the foil. Time instants: (a) 
$t_0$; (b) 
$t_0+1/4$; (c) 
$t_0+1/2$; (d) 
$t_0+3/4$; (e) 
$t_0$; (f) 
$t_0+1/4$; (g) 
$t_0+1/2$; (h) 
$t_0+3/4$.
As the Stokes number (dimensional frequency) is increased, the flow breaks the spatial symmetry (2.9), as seen in figure 4(a–d) for 
$\beta =6$. Vortices shed on each side of the foil are of slightly different shape and intensity for all time instants. This asymmetric flow then induces an instantaneous horizontal force accelerating the foil. Figure 4(e) shows that an initial small perturbation of the horizontal velocity 
$u_g$ grows exponentially in time, before saturating for 
$t>200$ towards a periodic state, as shown in the close-up view displayed in figure 4(f). The amplitude of oscillation of the horizontal velocity is very weak compared with its time-averaged value, denoted hereinafter as 
$\langle u_g\rangle$. Being positive, the flapping foil self-propels in the positive 
$x$-direction. Solutions self-propelling in the negative 
$x$-direction can also be found by modifying the initial horizontal velocity. The effect of the initial condition on the symmetry breaking direction was investigated in Jallas et al. (Reference Jallas, Marquet and Fabre2017). The Fourier spectrum of the foil horizontal velocity, displayed in figure 5(a), shows that it oscillates at the (non-dimensional) frequency 
$f=2$, i.e. twice the (non-dimensional) flapping frequency equal to 
$f=1$ independently from the Stokes number. This frequency doubling of the horizontal velocity is related to the spatio-temporal flow symmetry (2.10). Over one flapping period, the horizontal force acting on the foil is identical during upward and downward strokes.
Figure 4. Unidirectional propulsive solution for 
$\beta =6$. (a–d) Vorticity flow field along a flapping period. Time evolution of the horizontal velocity 
$u_g$ (initially equal to zero) (e) over the whole simulated time and (f) restricted to time window indicated by the rectangle in (e). The instant 
$t_0$ is depicted with a vertical line in (e,f). The period of the horizontal velocity is 
$0.5$, half the vertical flapping period 
$1$. Time instants: (a) 
$t_0=324$; (b) 
$t_0+1/4$; (c) 
$t_0+1/2$; (d) 
$t_0+3/4$.
Figure 5. Frequency spectra of the horizontal foil speed for (a) the periodic unidirectional propulsive solution (
$\beta =6$) and (b) the quasi-periodic back and forth propulsive solution (
$\beta =13$). The fundamental frequency of the horizontal speed (vertical black lines) is twice the vertical flapping frequency. The low frequency of the quasi-periodic solution (figure 6e) is identified by the blue line.
For higher values of the Stokes number, the spatial flow symmetry is still broken but the propulsion is no longer unidirectional. The foil periodically reverses its propulsive direction, as if a restoring force was at play. A typical solution obtained for 
$\beta = 13$ is displayed in figure 6. The time evolution of the horizontal velocity shown in 6(e) clearly indicates that, after an initial exponential growth in time, the foil velocity slowly oscillates between positive and negative values, over a period approximately 
$50$ times larger than the flapping period. This solution is no longer periodic, but quasi-periodic, as clearly shown by the Fourier spectrum of the foil horizontal velocity displayed in figure 5(b). Two fundamental frequencies are obtained, one at 
$f=2$ corresponding to twice the flapping frequency, and one around 
$f=0.018$ corresponding to the slow period. The multiple peaks observed around each fundamental frequency are induced by nonlinear interactions. Coming back to the horizontal velocity of the foil, its time-average over the long period is zero. Thus, this solution is not a coherent (unidirectional) propulsive state (Alben & Shelley Reference Alben and Shelley2005). The foil oscillates back and forth around a fixed point in space. Nevertheless, the horizontal velocity, time-averaged along the (short) flapping period, is either positive or negative, as seen in figure 6(f). A propulsive effect is thus obtained at this time scale. The instantaneous vorticity fields displayed in figures 6(a–d) correspond to a flapping period (marked with the green dot in figure 6e) where the velocity of the foil is positive, while those shown in figure 6(g–j) correspond to a negative foil velocity (blue dot). In both cases, the leading-edge vortex is of smaller size and closer to the foil than the trailing-edge vortex. Interestingly, such vortex pattern is not observed for all phases of the long period, and in particular in the acceleration phases, marked with red dots in figure 6(e). The corresponding instantaneous vorticity fields are depicted in figure 7. In between the instants corresponding to figures 7(a) and 7(b), the foil accelerates and self-propels in the right direction but the leading-edge vortex (right) is now of larger size and further away from the foil, compared with the trailing-edge vortex (left). This suggests that the foil motion is induced by a suction of the leading-edge vortex. As the foil further accelerates, the leading- and trailing-edge vortices are progressively convected downstream until the more classical propulsive pattern is recovered when the foil reaches its maximal velocity (see figure 7d).
Figure 6. Back and forth solution for 
$\beta =13$. (a–d) and (g–j) Vorticity contours for a flapping period – starting from different time instants in (a,g). (e,f) Horizontal speed 
$u_g$ time evolution with dotted rectangle close-up in (f). Instants (a,g) are indicated by filled green and blue dots in (e). Long and small periods of 
$u_g(t)$ are respectively indicated in (e,f). Time instants: (a) 
$t_0=178$; (b) 
$t_0+1/4$; (c) 
$t_0+1/2$; (d) 
$t_0+3/4$; (g) 
$t_1=266$; (h) 
$t_1+1/4$; (i) 
$t_1+1/2$; (j) 
$t_1+3/4$.
Figure 7. Back and forth solution for 
$\beta =13$. (a–d) Vorticity contours of the time instants represented in figure 6(e) by filled red dots. An orange arrow indicates the horizontal velocity of these time instants. These time instants are equally spaced of three flapping periods.
The existence and characterization of the three propulsive solutions is displayed in figure 8 as a function of the Stokes number in the range 
$2\le \beta \le 20$. The regimes of non-propulsive, propulsive and back and forth solutions are identified with white, red and grey background colours, respectively. The horizontal velocity is depicted in figures 8(a) and 8(c), with black dots for the time-averaged value 
$\langle u_g\rangle$ and vertical bars for the fluctuation amplitude. For the back and forth solution, the long period is used for time-averaging. The frequencies 
$f$ of the foil velocities are shown in figures 8(b) and 8(d), the open circles denoting the vertical flapping frequency, while the filled circles correspond to the frequency of the horizontal velocity. Symmetric non-propulsive solutions exist for small Stokes numbers 
$\beta <4$ (region I) and for intermediate values in the range 
$9.53 \leq \beta \leq 11.25$ (region III). For these Stokes numbers, no locomotion is achieved and the foil remains in its position, producing a spatially symmetric periodic flow. Propulsive solutions appear for 
$\beta =4$. They are characterized by non-zero time-averaged horizontal velocities – both negative and positive depending on initial conditions – with very small amplitudes of fluctuations and a frequency of oscillation equal to 
$f=2$. As the Stokes number is increased, this frequency remains constant while the (absolute value of) time-averaged propulsive velocity continuously increases until 
$\beta \sim 8.5$. The mean velocity then decreases and abruptly (discontinuously) falls to zero for 
$\beta =9.58$. By decreasing again the Stokes number, we have identified a small range of Stokes numbers (
$9.53 \le \beta \le 9.58$), visible in figure 8(c), where non-propulsive and propulsive solutions co-exist. Therefore, the bifurcation from a propulsive to non-propulsive solution is sub-critical around 
$\beta =9.5$, unlike the transition from a non-propulsive to propulsive solution at 
$\beta =4$, which is super-critical. Finally, back and forth solutions are observed when increasing the Stokes number 
$\beta \ge 11.25$ (region IV). They are characterized by a zero time-averaged horizontal speed with a large amplitude of fluctuations. These quasi-periodic solutions are characterized by two fundamental frequencies, the high frequency (black dots) and the low frequency (blue dots). As seen in figure 8(d), the low frequency decreases towards zero when increasing the Stokes number.
Figure 8. Non-propulsive (white), unidirectional propulsive (red) and back and forth (grey) regimes as a function of the Stokes number 
$\beta$. (a,c) Time-averaged (circles) and oscillation amplitude (error bars) of the foil horizontal velocity. (b,d) Forcing frequency of the foil vertical velocity (open circles) and frequency of the horizontal foil velocity (filled circles). Panel (c) is a close-up view of (a) highlighting the transition between unidirectional propulsive, non-propulsive and back and forth solutions. Panel (d) is a close-up view of (b) showing the evolution of the low frequency in the back and forth solution as a function of the Stokes number. Parameters: 
$A = 0.5$ and 
$\rho = 100$.
2.3.2. Self-propelled regimes for $\rho =1$
Considering a density ratio closer to swimming organisms (
$\rho =1$), figure 9, the same self-propelled regimes and transitions are identified. In figure 9(a), as the previously presented case of 
$\rho =100$, symmetric non-propulsive solutions become unidirectional propulsive ones for the critical Stokes number 
$\beta =4$. The time-averaged velocity thus increases with an average value slightly slower than 
$\rho =100$ up to 
$\beta \sim 8.5$, decreasing beyond this point and abruptly falling to zero for 
$\beta =9.5$, a small change in the restabilization 
$\beta$ of the previous density ratio. By decreasing the Stokes number, co-existing unidirectional propulsive and non-propulsive solutions were again obtained, this time for 
$9.42 \leq \beta \leq 9.5$, visible in figure 9(b). Back and forth solutions are finally observed for 
$\beta \geq 10.6$. These solutions are significantly encouraged for smaller density ratios. Beside the smaller critical Stokes number, they present now a velocity four times greater in amplitude than for 
$\rho =100$. As presented by the two considered values of 
$\rho =1$ and 
$\rho =100$, the self-propelled regimes transition route appears to be robust to variations of the density ratio.
Figure 9. Self-propelled regimes for a density ratio 
$\rho = 1$. Background colours are the same as the previous figure. (a) Evolution of the time-averaged (circles) and amplitude (error bars) of the foil horizontal velocity with the Stokes number. (b) Close-up view highlighting the transition between the self-propelled regimes. Other control parameters: 
$A = 0.5$ and 
$h = 0.05$.
Self-propelled regimes and the transition route, for this Stokes number range, are also invariant with respect to the foil geometry. Simulations conducted for an elliptical foil of aspect ratio 
$h=0.1$, reported in appendix C, have presented similar results and transition route as the ones described for the rectangular foil with rounded edges.
3. Fluid–solid stability analysis of non-propulsive periodic solutions
The emergence of the propulsive solutions identified in the previous section with nonlinear unsteady simulations is now investigated by analysing the stability of non-propulsive periodic solutions. The Floquet stability analysis is introduced in § 3.1 by considering a perturbation of the horizontal foil velocity, in addition to the flow perturbation. The numerical method is then explained in § 3.3. Results of such stability analysis, that couples the fluid and solid perturbations, are presented in § 3.4 for the density ratio 
$\rho =100$. First the two synchronous and asynchronous modes found unstable are carefully described. Then these modes are discussed in light of the nonlinear results previously described. Finally, the influence of the density ratio on the linear results is discussed in § 3.5.
3.1. Fluid–solid Floquet stability analysis
The flow variables and the foil horizontal velocity are decomposed as
\begin{equation} (\boldsymbol{u},p,u_{g})=(\boldsymbol{u}_{b},p_{b},0) + \epsilon (\boldsymbol{u}',p',u_{g}'), \end{equation}
where 
$(\boldsymbol {u_b},p_b)$ denote the periodic base flow fields. Since it satisfies the spatial symmetry (2.9) at every instant of the flapping period, the foil horizontal velocity of the periodic base solution is equal to zero. Infinitesimal perturbations (
$\epsilon \ll 1$) are superimposed to the periodic base solution meaning that in addition to perturbing the base flow field 
$(\boldsymbol {u}',p')$ as in Deng & Caulfield (Reference Deng and Caulfield2016), the foil horizontal velocity 
$u_{g}'$ is perturbed. Note that no perturbation of the vertical velocity is considered since the flapping velocity 
$\boldsymbol {v}_{g}$ is imposed in the present analysis. By injecting the above decomposition into (2.3)–(2.6) and retaining the first-order term in 
$\epsilon$, we obtain the following system of equations governing the linear dynamic around the non-propulsive periodic solution:
\begin{align} \boldsymbol{\mathsf{B}} \, \frac{\partial}{\partial t} \underbrace{\begin{bmatrix} \boldsymbol{u}' \\ p' \\ u_g' \end{bmatrix}}_{\boldsymbol{q}'} = \underbrace{\begin{bmatrix} - [(\boldsymbol{u}_{b}(t)-v_g(t) \boldsymbol{e}_{y}) \boldsymbol{\cdot} \boldsymbol{\nabla}] - \boldsymbol{\nabla} \boldsymbol{u}_{b}(t) + \beta^{{-}1} \boldsymbol{\Delta} & -\boldsymbol{\nabla} & (\boldsymbol{\nabla} \boldsymbol{u}_{b}) \boldsymbol{\cdot} \boldsymbol{e}_{x} \\ -\boldsymbol{\nabla} & 0 & 0 \\ (\rho S)^{{-}1} \mathcal{F}_v & (\rho S)^{{-}1} \mathcal{F}_p & 0 \end{bmatrix}}_{\boldsymbol{\mathsf{L}}(\boldsymbol{u}_{b},p_{b})} \begin{bmatrix} \boldsymbol{u}' \\ p' \\ u_g' \end{bmatrix}. \end{align}
The first two rows are the linearized momentum and mass equations governing the fluid velocity and pressure perturbations. They are coupled to the foil velocity perturbation 
$u_g'$ via two terms. First, the bulk term 
$(\boldsymbol {\nabla } \boldsymbol {u}_{b}) \boldsymbol {\cdot } \boldsymbol {e}_{x}$ (block 
$(1,3)$ in the right-hand side matrix) that modifies the production of fluid perturbation in the momentum equation, and, second, the boundary conditions at the foil surface 
$\varGamma _{w}$, where the equality of fluid and solid perturbations hold:
\begin{equation} (u',v')(\boldsymbol{x}_w,t)=(u_g',0). \end{equation}
The third row indicates that the horizontal acceleration of the foil is equal to the horizontal force exerted by the flow perturbation, which is here separated into viscous 
$\mathcal {F}_v$ and pressure 
$\mathcal {F}_p$ contributions, respectively defined as
\begin{align} \mathcal{F}_v (\boldsymbol{u'}) = \beta^{{-}1} \displaystyle\int_{\varGamma_w} \left(\boldsymbol{\nabla}\boldsymbol{u'}+\boldsymbol{\nabla}\boldsymbol{u'}^{T}\right)\boldsymbol{\cdot} \boldsymbol{n}) \boldsymbol{\cdot} \boldsymbol{e_x} \, \mathrm{d}\varGamma_w,\quad \mathcal{F}_p (p')= \displaystyle\int_{\varGamma_w}(- p'\boldsymbol{n}) \boldsymbol{\cdot} \boldsymbol{e_x} \, \mathrm{d}\varGamma_w. \end{align}
We note that, in (3.2), the viscous and pressure forces are weighted by the invert of the foil mass 
$\rho S$ . Consequently, when the ratio of solid-to-fluid density increases, the effect of the flow on foil perturbations decreases. In the limit of infinite density ratio, i.e. 
$\rho \rightarrow \infty$, it even vanishes leading to the following one-way coupled fluid–solid system:
\begin{equation} \boldsymbol{\mathsf{B}} \, \frac{\partial}{\partial t} \begin{bmatrix} \boldsymbol{u}' \\ p' \\ u_g' \end{bmatrix} = \begin{bmatrix} - [(\boldsymbol{u}_{b}-v_g \boldsymbol{e}_{y}) \boldsymbol{\cdot} \boldsymbol{\nabla} ] - \boldsymbol{\nabla} \boldsymbol{u}_{b} + \beta^{{-}1} \boldsymbol{\Delta} & -\boldsymbol{\nabla} & (\boldsymbol{\nabla} \boldsymbol{u}_{b}) \boldsymbol{\cdot} \boldsymbol{e}_{x}\\ -\boldsymbol{\nabla} & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{u}' \\ p' \\ u_g' \end{bmatrix}. \end{equation}In that limit case, the horizontal acceleration of the foil is zero, but not its horizontal velocity. This velocity might still affect the flow perturbation via the coupling terms described above. This one-way coupling analysis is thus different from the hydrodynamic Floquet analysis, performed, for instance, by Elston et al. (Reference Elston, Sheridan and Blackburn2004), Elston et al. (Reference Elston, Blackburn and Sheridan2006) on a forced oscillating cylinder and more recently applied by Deng et al. (Reference Deng, Teng, Caulfield and Mao2016) and Deng & Caulfield (Reference Deng and Caulfield2016) respectively on the forced oscillation of an ellipsoid and on the self-propulsion of different aspect ratio oscillating ellipses.
In hydrodynamic Floquet analysis the horizontal perturbation velocity is assumed to be zero (
$u_g'=0$). The perturbation equations (3.2) then simplify to
\begin{equation} \boldsymbol{\mathsf{B}}_f \, \frac{\partial}{\partial t} \begin{bmatrix} \boldsymbol{u}_{f}' \\ p_{f}' \end{bmatrix} = \begin{bmatrix} - [(\boldsymbol{u}_{b}-v_g \boldsymbol{e}_{y}) \boldsymbol{\cdot} \boldsymbol{\nabla}] - \boldsymbol{\nabla} \boldsymbol{u}_{b} + \beta^{{-}1} \boldsymbol{\Delta} & -\boldsymbol{\nabla}\\ -\boldsymbol{\nabla} & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{u}_{f}' \\ p_{f}' \end{bmatrix}, \end{equation}
where 
$\boldsymbol{\mathsf{B}}_f$ is the portion of the operator 
$\boldsymbol{\mathsf{B}}$ related to the fluid variable, and the subscript 
$f$ is introduced to indicate that the perturbation is purely hydrodynamic. At the foil boundary, they satisfy the no-slip boundary condition
\begin{equation} (u_{f}',v_{f}')(\boldsymbol{x}_w,t)=(0,0). \end{equation}
All the above equations are closed by considering that the fluid perturbations vanish 
$\boldsymbol {u'}=0$ far away from the foil.
We would like to stress that the fluid–solid perturbation analysis encompasses the purely hydrodynamic perturbation analysis, since hydrodynamic perturbations should be retrieved in the fluid–solid analysis if the foil velocity perturbation is zero. This will be further discussed when presenting results in § 3.5 in the limit case 
$\rho \rightarrow \infty$.
Following Elston et al. (Reference Elston, Sheridan and Blackburn2004) or Jallas et al. (Reference Jallas, Marquet and Fabre2017), the perturbations are further decomposed in the form
\begin{equation} \boldsymbol{q'}(\boldsymbol{x},t)= \sum (\boldsymbol{\hat{q}}_j (\boldsymbol{x},t) \textrm{e}^{\lambda_j t} + \textrm{c.c.} ), \end{equation}
where 
$\boldsymbol {\hat {q}}_j$ are T-periodic functions, called the Floquet modes, associated to the complex numbers 
$\lambda _j$, called the Floquet exponents. The Floquet multiplier, defined as 
$\mu _j=\textrm {e}^{\lambda _j T}$, is rather used in the following. It represents the complex amplitude gain of the periodic Floquet mode over one period, i.e. 
$\boldsymbol {\hat {q}}_j(\boldsymbol {x},T)=\mu _j \boldsymbol {\hat {q}}_j(\boldsymbol {x},0)$. The polar decomposition of the Floquet multiplier is 
$\mu _j=|\mu _j| \textrm {e}^{\textrm {{i}} \phi _{\textrm {j}}}$, where the modulus 
$|\mu _j|$ quantifies the growth (or decay) of the corresponding Floquet mode over the period, and 
$\phi _{j}$ represents its phase shift over the same period. The stability of the periodic base solution is then addressed by considering the Floquet multiplier with the largest modulus. If its absolute value, denoted 
$|\mu _0|$, is greater than one, the corresponding Floquet mode will grow over one period and the periodic base solution is thus unstable. When the leading Floquet multiplier is real (
$\phi _{0}=0$), the Floquet mode is synchronous as the perturbation evolves in time with the period of the base flow. When the leading Floquet multiplier is complex (
$\phi _{0} \neq 0$), the Floquet mode is asynchronous and a frequency, denoted 
$f'$ in the following, related to the phase of the Floquet mode as 
$f' = \phi _{0} /(2 {\rm \pi})$ is introduced.
3.2. Time-averaged analysis of fluid–solid Floquet modes
To better understand how the periodic flow perturbation is related to the destabilisation of a fluid–solid Floquet mode, we examine the equation
\begin{equation} \lambda \hat{u}_{g} +\frac{\textrm{d}\hat{u}_{g}}{\textrm{d}t} = \frac{1}{\rho S} F_x(\hat{\boldsymbol{u}},\hat{p}) \end{equation}that expresses the instantaneous equilibrium between the horizontal force exerted by the flow component of the Floquet mode and the horizontal acceleration of the foil. The latter is composed of two terms, one related to the growth/decay of the mode, and one related to its instantaneous acceleration. Due to the periodicity of the Floquet mode, the latter disappears when time-averaging over a flapping period, yielding
\begin{equation} \lambda \langle \hat{u}_{g} \rangle = \frac{1}{\rho S} \langle \hat{F}_x \rangle, \end{equation}
where 
$\left \langle \cdot \right \rangle$ denotes the time-average over a flapping period.
For synchronous modes, the Floquet exponent and mode are real variables and the above expression gets
\begin{equation} \lambda_{r} = \frac{1}{\rho S} \frac{\langle \hat{F}_x\rangle}{\langle\hat{u}_g\rangle}. \end{equation}The growth rate of the Floquet mode is thus proportional to the ratio between the mean component of the force and the mean velocity of the Floquet mode. The Floquet mode is thus unstable (respectively stable) when the force and velocity are of same (different) sign. In the case of asynchronous modes, we introduce the polar decomposition of the time-averaged horizontal force and velocity in (3.10) to obtain the simple relation
\begin{equation} \lambda = \frac{1}{\rho S} \frac{|\langle \hat{F}_x\rangle|}{|\langle\hat{u}_g\rangle|} \textrm{e}^{\textrm{{i}} \psi} \end{equation}
showing that the growth rate (real part) is now also related to the phase difference 
$\psi$ between the time-averaged force and velocity, and not only to their ratio. The relations (3.11) and (3.12) will be used in § 3.4 for a physical discussion of the Floquet mode
3.3. Numerical method
The periodic non-propulsive solutions 
$(\boldsymbol {u_b},p_b)$ are computed by integrating in time the governing equations (2.3),(2.4a,b) using the same temporal and spatial discretization scheme as described in the previous section, but with the following boundary conditions:
\begin{gather} (u_b,v_b)(x_{w},y_{w},t) = (0,2 {\rm \pi}A \sin(2 {\rm \pi}t) ), \end{gather}
\begin{gather}(u_b,\partial_x v_b)(0,y,t) = (0,0). \end{gather}
The first set of boundary conditions, applied at any point 
$(x_w,y_w)$ of the foil surface, allows imposing the flapping motion of the foil in the vertical direction without any motion in the horizontal direction. The second set of boundary conditions, applied on the 
$y-$ axis, is used to enforce the spatial reflection symmetry of the flow characteristic of the non-propulsive solution. Typically, 
$50$ flapping periods are simulated to reach a periodic solution. Note that, for computational efficiency, the computational domain can be reduced to the left or right part of the full domain shown in figure 2, but this is not mandatory. In that case, the periodic base solution on the full domain is retrieved by using the spatial symmetry relation (2.9).
The strategy to compute Floquet modes is the one proposed by Barkley & Henderson (Reference Barkley and Henderson1996). The stability of an initial perturbation is assessed regarding the action of the propagator over one period 
$\boldsymbol {\varPhi }$, also known as the Monodromy matrix. The action of this Monodromy matrix over the perturbation at an arbitrary initial time 
$t_0$ is formally denoted by 
$\boldsymbol {q'}(\boldsymbol {x},t_0+T)=\boldsymbol {\varPhi } \, \boldsymbol {q'}(\boldsymbol {x},t_0)$. It is actually obtained by time integration along a period of the linearized equations (3.2) with boundary conditions (3.3), using the temporal and numerical discretization schemes previously described. An Arnoldi method with a modified Gram–Schmidt algorithm for the orthogonalization step (Saad Reference Saad2011) is implemented in the FreeFEM software (Hecht Reference Hecht2012) to approximate the Monodromy matrix in a low-dimensional space. The eigenvalues of this reduced matrix approximate the Floquet multiplier and its eigenvectors allow for reconstructing the Floquet modes at the initial time, i.e. 
$\hat {\boldsymbol {q}}(\boldsymbol {x},t_0)$. A minimal number of 
$30$ Arnoldi vectors is used in the following, this number being further increased in steps of 
$10$ when necessary in order for the dominant eigenvalue to converge to five significant digits. All computed modes are normalised by the kinetic energy of the coupled fluid–solid system. A validation of this method is detailed in appendix B.
Finally, as the Arnoldi method gives access to the Floquet mode at an initial time, the mode complete temporal evolution is obtained through time integration of the equation below over one flapping period,
\begin{equation} \boldsymbol{\mathsf{B}} \frac{\partial \hat{\boldsymbol{q}}}{\partial t} - \boldsymbol{\mathsf{L}}(\boldsymbol{u}_b,v_g) \hat{\boldsymbol{q}} ={-} \lambda \, \boldsymbol{\mathsf{B}} \, \hat{\boldsymbol{q}}, \end{equation}
using as an initial condition the Floquet mode obtained with the Arnoldi algorithm. The Floquet exponent 
$\lambda$ being known, the right-hand side term appropriately counteracts the temporal growth (respectively decay) of the unstable (respectively stable) Floquet mode.
3.4. Results of Floquet analyses for $\rho =100$
The stability analysis of non-propulsive periodic solutions has been performed for the flapping amplitude 
$A=0.5$ and Stokes numbers in the range 
$1 \le \beta \le 20$. The instantaneous vorticity fields of a typical base non-propulsive solution, obtained for 
$\beta =13$, are depicted in figure 10 at four instants in the flapping period. The spatial left–right symmetry (2.9) along the 
$y$-axis is clearly satisfied at every instant of the flapping period. By comparing the dipolar structure at time 
$t_0+1/2$ and 
$t_0$, this solution also satisfies the spatio-temporal symmetry (2.10).
Figure 10. Non-propulsive solution obtained for 
$A=0.5$ and 
$\beta =13$. The instantaneous vorticity field is depicted at four instants of the (unitary) flapping period. Time instants: (a) 
$t_0$; (b) 
$t_0+1/4$; (c) 
$t_0+1/2$; (d) 
$t_0+3/4$.
3.4.1. Floquet multipliers: fluid–solid versus hydrodynamic analysis
Results of the fluid–solid Floquet analysis performed for the density ratio 
$\rho =100$ are first depicted in figure 11, with the modulus and frequency of the leading Floquet multipliers as a function of 
$\beta$ in figures 11(a) and 11(b), respectively.We can clearly identify two ranges of Stokes number where the leading Floquet modes get unstable (
$|\mu _0|>1$) and they compare very well with regions II and IV, where unidirectional and back and forth self-propelled solutions were obtained in the nonlinear unsteady simulations.
Figure 11. Fluid–solid Floquet analysis for 
$\rho =100$. (a) Absolute value and (b) frequency of the dominant Floquet multiplier (black circles) as a function of the Stokes number 
$\beta$. The latin numbers indicate regimes of nonlinear solution identified in § 2, the red region corresponding to unidirectional self-propelled solutions, while the grey regions correspond to back and forth self-propelled solutions. In (b) the open circles indicate the low frequencies characteristic of the back and forth solutions.
In the range 
$4.00 \le \beta \le 9.53$ the unstable Floquet modes are synchronous (
$f'=0$). Thus, the perturbation does not break the periodicity of the base solution, in agreement with the unidirectional propulsive solution observed in region II. The frequency 
$f=2$ characterizing the horizontal speed of this propulsive solution (see figure 8b) is rather related to the spatio-temporal symmetry of the Floquet mode, as will be seen in § 3.4.2. A quantitative comparison of the thresholds is provided in table 1. The destabilisation of the synchronous mode at 
$\beta =4$ is in perfect agreement with the emergence of the unidirectional propulsive solution, i.e. the transition between region I and II. On the other hand, the threshold value 
$\beta =9.53$ corresponding to the stabilization of this mode is slightly different from the threshold 
$\beta =9.58$ above which the unidirectional propulsive solution disappears (transition from II to III). This is due to the sub-critical nature of the bifurcation at this threshold, clearly seen in figure 8(c). But the threshold value 
$\beta =9.53$ found with the linear stability corresponds perfectly to the disappearance of the symmetric solution when decreasing the Stokes number (transition from III to II).
Table 1. Critical thresholds obtained with the unsteady nonlinear simulations and the stability analyses. Here ‘
$\times$’ indicates that no unstable modes were obtained.
When further increasing the Stokes number 
$\beta \ge 11.25$, an asynchronous Floquet mode gets unstable, with a very low frequency (
$f'\sim 0.01$) compared with the flapping frequency (
$f=1$). The destabilization of this asynchronous mode occurs at the same value of the Stokes number for which the back and forth solution appears. The frequency of the asynchronous mode is compared with the frequency of this solution in figure 11(b). They compare very well at the threshold 
$\beta =11.25$, but when the Stokes number is increased, the agreement gets worse. Opposite trends are observed, with an increase of the Floquet mode frequency and a decrease of the back and forth frequency. The flow nonlinearities, i.e. which can be mean-flow distortion or higher-harmonics generation/interaction, clearly play a role in the frequency selection of this back and forth solution. Note that such discrepancy between the linear and nonlinear frequency has been observed for the unsteady wake of a fixed cylinder flow, and is predominantly due to mean-flow distortion in that case (Barkley Reference Barkley2006; Sipp & Lebedev Reference Sipp and Lebedev2007).
Before examining the unstable synchronous and asynchronous Floquet modes, we describe results of the purely hydrodynamic Floquet analysis, performed in the same range of Stokes number and displayed in figure 12 with black squares. First, we note that no unstable mode is found in the range of the Stokes number corresponding to region II. So, the purely hydrodynamic Floquet analysis cannot explain the emergence of the unidirectional propulsive solution. One unstable Floquet mode is found only for larger Stokes numbers 
$\beta \ge 12.20$. Just above the threshold, this mode is asynchronous, but when the Stokes number is increased to 
$\beta =13.3$, the pair of complex asynchronous modes collapses on the real axis becoming two real synchronous modes. One of these modes is further destabilized when increasing 
$\beta$, while the other one is stabilized for 
$\beta > 15$. The spatial structure of the unstable asynchronous modes found with the hydrodynamic stability analysis are very similar to the asynchronous mode obtained with the fluid–solid analysis, that are described in § 3.4.3, and thus will not be further described. Their frequency is much smaller, as displayed in figure 12(b), failing to predict the frequency of the back and forth solution even at the threshold. As indicated in table 1, this threshold is underpredicted by the purely hydrodynamic analysis. Unlike the fluid–solid Floquet analysis, the purely hydrodynamic Floquet analysis on one hand cannot explain the emergence of the unidirectional propulsive solutions, and on the other hand does not accurately predict the occurrence of the back and forth solutions. These two observations offer a possible explanation to observations made by Deng & Caulfield (Reference Deng and Caulfield2016) when comparing the results of unsteady nonlinear simulations and a purely hydrodynamic stability analysis. In their study, the purely hydrodynamic analysis did not estimate the enhancement (earlier transition) of the quasi-periodic nonlinear solutions for ellipses of aspect ratio 
$h=0.5$ and did not predict the unidirectional propulsion threshold for ellipses of aspect ratio 
$h=0.1$. Visibly, the onset of self-propulsion cannot be explained by the flow symmetry breaking instability alone and the self-propelled fluid–solid coupling is vital for its prediction.
Figure 12. Hydrodynamic Floquet analysis. Same legend as in previous figure.
3.4.2. Synchronous Floquet modes: emergence of unidirectional propulsion
Turning back to results of the fluid–solid Floquet analysis, the synchronous Floquet mode is depicted in figure 13 for 
$\beta =6$. The vorticity field of this periodic mode is displayed with coloured maps at four instants of the flapping period in figure 13(a–d), and the vorticity of the periodic base flow is superimposed using black (dashed) isolines for positive (negative) values. First we note that the synchronous Floquet mode breaks the left–right symmetry (2.9) of the base flow, since the perturbative vorticity is an even function of 
$x$ while the base vorticity is an odd function of 
$x$. But, as the base flow, it still satisfies the spatio-temporal symmetry (2.10), so that we can restrict our description of the mode to the upstroke phase 
$0\le t \le 1/2$. During the acceleration phase of this upstroke motion 
$(t<1/4)$, a patch of positive vorticity exists above the foil, in a region where the vorticity of the base flow is weak, since the latter is rather generated under the foil during the upstroke. This patch of vorticity corresponds to a shear region in the flow perturbation, that induces an increase in the horizontal forces exerted on the foil, as seen in figure 13(f). During the second half of the upstroke (
$1/4 < t < 1/2$) where the vertical velocity of the foil decreases, the patch of positive vorticity also decreases in size and amplitude. Meanwhile, a patch of negative vorticity appears under the foil, leading to a decrease of the horizontal force. This oscillation of the horizontal force results in an out-of-phase oscillation of the horizontal velocity, shown in figure 13(e). Interestingly, the average of the horizontal force and velocity over a flapping period (indicated with a dashed line in figures) are non-zero and positive here. Therefore, this synchronous Floquet mode is clearly at the origin of the propulsion of the foil in the horizontal direction. Note that the direction of propulsion is not determined by the Floquet mode, since its amplitude is arbitrarily positive (here) or negative, leading to a right (here) or left displacement of the foil.
Figure 13. Unstable synchronous Floquet mode for 
$\beta =6$. (a–d) Vorticity contours for a flapping period (base flow positive (respectively negative) values represented by solid (respectively dashed) lines). (e) Foil horizontal speed and (f) force. Time-averaged value along a flapping period of both values is represented by a dashed line in (e,f). In (e) the vertical speed 
$v_g$ of the base flow (light blue dashed line) is represented in the right axis. Time instants: (a) 
$t=0$; (b) 
$t=1/4$; (c) 
$t=1/2$; (d) 
$t=3/4$.
To stress again the role of the fluid–solid coupling in the destabilization of the mode, we consider the time-averaged analysis of the Floquet mode exposed in § 3.2. For synchronous modes, it was shown that their growth rate is given by
\begin{equation} \lambda_{r} = \frac{1}{\rho S} \frac{\langle \hat{F}_x\rangle}{\langle\hat{u}_g\rangle}, \end{equation}
i.e. the ratio between the time-averaged horizontal force and velocity, weighted by the foil mass. These quantities are plotted in figure 14 as a function of 
$\beta$. The time-averaged horizontal velocity, shown with black circles in figure 14(a), is positive for all values of the Stokes number. We note that its evolution is different from the time-averaged velocity of the foil computed with temporal simulation (open circles), indicating the importance of flow nonlinearities in the terminal foil velocity. Examining now the horizontal force in figure 14(b), its changes of sign correspond clearly to the destabilization and stabilization of the Floquet mode (figure 14c). When this force is positive (respectively negative), the growth rate is also positive (respectively negative), in agreement with the above relation (recalling that the horizontal velocity is always positive). Finally, we conclude that this synchronous Floquet mode is responsible for the emergence of the unidirectional propulsion solution obtained in region II, delimited in red in the figure. Retaining the fluid–solid coupling at the perturbation level is fundamental to explain the destabilization of this mode.
Figure 14. Time-averaged horizontal (a) velocity and (b) force for the synchronous Floquet mode as a function of 
$\beta$. (c) Real part of the Floquet exponent. In (a) the mean horizontal speed obtained with nonlinear simulations is represented with empty grey circles for comparison.
3.4.3. Asynchronous Floquet modes: emergence of back and forth solution
We now examine the asynchronous Floquet modes that become unstable for larger Stokes numbers. The complex mode, obtained at 
$\beta =13$ and displayed in figure 15, also breaks the left–right symmetry and satisfies the spatio-temporal symmetry (2.10). The instantaneous real (respectively imaginary) part of the vorticity is shown in figures 15(a–d) (respectively figure 15e–h) at four instants of the flapping period. The amplitude of the real part is noticeably larger than that of the imaginary part, and their spatial structures are quite different. The real part of the mode bears similarities with the synchronous Floquet mode found by Jallas et al. (Reference Jallas, Marquet and Fabre2017) to explain the deviation of propulsive wakes in flapping wings and the displacement modes of vortices (Fabre, Sipp & Jacquin Reference Fabre, Sipp and Jacquin2006; Brion, Sipp & Jacquin Reference Brion, Sipp and Jacquin2014). Let us consider the right solid dark line representing the base flow vortex of positive vorticity in figure 15(a). The perturbation has positive vorticity on the lower left and negative vorticity on the upper right part of the monopole. This superposition strengthens the lower left part of the monopole while weakening the upper right one, resulting in a net displacement of its centre to the lower left. The displacement of the dipolar vortex structure results in a horizontal force exerted on the foil whose temporal evolution is shown in figure 15(k). Due to the spatio-temporal flow symmetry, the frequency is twice the flapping frequency. Interestingly, the horizontal force strongly oscillates around a negative time-averaged value. Compared with the oscillation of the foil horizontal velocity displayed in figure 15(i), we first note that they are not in phase. As explained later, this phase difference is crucial to understand the destabilization of the asynchronous mode. Then, we also remark that the fluctuation of the horizontal velocity is much smaller and around a time-averaged value that is positive. Therefore, the real part of the asynchronous mode produces a mean resistive force over the flapping period, that decreases the foil horizontal velocity.
Figure 15. Unstable asynchronous Floquet mode for 
$\beta =13$. Vorticity contours of the (a–d) real and (e–h) imaginary parts of the mode at four instants of the flapping period. Positive and negative values of the base flow vorticity are depicted with solid and dashed contours. (i,j) Real and imaginary parts of the horizontal velocity 
$\hat {u}_{g}$. (k,l) Real and imaginary parts of the horizontal force 
$\hat {F}_{x}$. In (i–l) dashed lines represents the time-averaged value of the plotted quantity. Time instants: (a) 
$t=0$; (b) 
$t=1/4$; (c) 
$t=1/2$; (d) 
$t=3/4$; (e) 
$t=0$; (f) 
$t=1/4$; (g) 
$t=1/2$; (h) 
$t=3/4$.
Turning now to the imaginary part of the asynchronous mode depicted in figure 15(e–h), its spatial structure is of much smaller amplitude than for the real part. It looks like a combination between the synchronous mode (see figure 13), responsible for the unidirectional self-propulsion of the foil, and the real part of the asynchronous mode, which creates a mean resistive force during a flapping period. The temporal evolution of the foil velocity and horizontal force are displayed in figures 15(j) and 15(l), respectively. The fluctuation of the force (figure 15l) is now much smaller than for the real part (figure 15k). The real and imaginary horizontal forces are out-of-phase by 
$1/4$. During the upstroke of the foil (
$t<0.5$), the minimal and maximal values of the imaginary horizontal force are obtained at 
$t=1/2$ and 
$t=1/4$, respectively, while they are obtained at 
$t=1/8$ and 
$t=3/8$ for the real part. Interestingly, the time-averaged value of the imaginary part is now positive, as for the horizontal velocity (figure 15j). Therefore, the imaginary part of this asynchronous mode produces a mean propulsive force that increases the foil velocity.
To further understand the contrasting actions of the real and imaginary parts of the Floquet mode, we introduce the real quasi-periodic perturbation defined as 
$\boldsymbol {\tilde {q}} = \boldsymbol {q}' \textrm {e}^{-\lambda _r t}$. Compared with the real perturbation 
$\boldsymbol {q}'$, the exponential growth/decay given by the real part of the Floquet exponent is counteracted. The perturbation 
$\boldsymbol {\tilde {q}}$ is quasi-periodic as it retains the low frequency oscillation given by the imaginary part of the Floquet exponent, in addition to the high frequency flapping period. Using the Floquet decomposition (3.8), the temporal evolution of this quasi-periodic perturbation is simply given by
\begin{equation} \boldsymbol{\tilde{q}}(\boldsymbol{x},t,t_{\epsilon}) = \mathrm{Re}(\boldsymbol{\hat{q}})(\boldsymbol{x},t)\cos(2{\rm \pi} t_{\epsilon}) - \mathrm{Im}(\boldsymbol{\hat{q}})(\boldsymbol{x},t)\sin(2{\rm \pi} t_{\epsilon}), \end{equation}
where 
$t_{\epsilon } = f' t$ is a slow time scale, as the frequency of the Floquet mode is very small compared with the flapping frequency, i.e. 
$f' \ll 1$. For 
$t_{\epsilon } \sim 0$, the above expression shows that the quasi-periodic perturbation is dominated by the real part of the periodic Floquet mode 
$\mathrm {Re}(\boldsymbol {\hat {q}})(\boldsymbol {x},t)$, while for 
$t_{\epsilon } \sim 1/4$, it is dominated by its imaginary part 
$\mathrm {Im}(\boldsymbol {\hat {q}})(\boldsymbol {x},t)$. The quasi-periodic perturbation thus slowly evolves between the real and imaginary parts of the Floquet mode, on a time scale given by the low frequency of this asynchronous Floquet mode. This slow evolution is depicted in figure 16 for the unstable asynchronous mode at 
$\beta =13$. The contrasting actions of the real and imaginary parts of the Floquet mode are clearly visible in figure 16(a–d) which show the vorticity fields at four instants of the slow period. As expected from (3.17), the quasi-periodic perturbation is similar to the real part of the Floquet mode at time 
$t_{\epsilon }=0$ (compare figures 15a and 16a) or to its opposite at time 
$t_{\epsilon }=0.5$, while it is of much smaller amplitude and similar to the imaginary part of the Floquet mode at time 
$t_{\epsilon }=0.25$ and 
$t_{\epsilon }=0.75$. Let us now examine the quasi-periodic evolution of the horizontal force exerted on the foil displayed in figure 16(f), as well as the resulting foil speed and displacement shown in figures 16(e) and 16(g), respectively. Around the slow time 
$t_{\epsilon }=0$, the horizontal force time-averaged over the fast flapping period (thick line in figure 16f) is negative and opposite to the positive horizontal velocity (figure 16e). So the quasi-perturbation, shown in figure 16(a) and dominated by the real part of the Floquet mode, creates a resistive force. When it slowly evolves towards 
$t_{\epsilon }=0.25$, the horizontal force remains negative but it is then a propulsive force since the foil velocity is negative. For 
$t_{\epsilon }>0.32$, the horizontal force gets positive and is a resistive force to its horizontal motion until the sign of the foil velocity changes around 
$t_{\epsilon }=0.6$. This slow oscillation of the time-averaged force and velocity creates a back and forth displacement of the foil, depicted in figure 16(g), around a mean position value around 
$\tilde {x}_{g}=0.5$. Although this displacement is infinitesimally small, it is in agreement with the direction switching observed in the unsteady nonlinear solutions of regime IV. Note also the horizontal forces and speed time-averaged along the slow period 
$T_{\epsilon }=1/f'$ is zero, so that there is no unidirectional propulsion of the foil at that slow time scale.
Figure 16. Temporal evolution of the quasi-periodic perturbation 
$\boldsymbol {\tilde {q}}$ for 
$\beta =13$. (a–d) Vorticity of the perturbation (colours) and base flow (black lines) at four instants 
$t_{\epsilon }$ of the slow period 
$T_{\epsilon }=1/f'=33$. (e–g) Time evolution of the horizontal (e) velocity, (f) force and (g) position of the foil over the slow period, shown as a function of both time 
$t$ and 
$t_{\epsilon }$. The horizontal dashed lines are for the time-averaged value of the plotted quantity over the slow time period. In (f) the thick curve depicts the time-averaged value of the horizontal force over the flapping period. The vertical solid lines in (f,g) indicate the instants where the foil pass by its mean horizontal position. Time instants: (a) 
$t_{\epsilon }=0$; (b) 
$t_{\epsilon }=0.25$; (c) 
$t_{\epsilon }=0.5$; (d) 
$t_{\epsilon }=0.75$.
Finally, to understand the destabilisation of this asynchronous mode in light of the fluid–solid interaction, we consider again the time-averaged analysis. The asynchronous mode being complex, the polar decomposition of the time-averaged force and velocity is 
$\langle \hat {F}_x\rangle =|\langle \hat {F}_x\rangle | \textrm {e}^{\textrm {i} \phi _{\textrm {F}}}$ and 
$\langle \hat {u}_g\rangle =|\langle \hat {u}_g\rangle | \textrm {e}^{\textrm {i} \phi _\textrm {U}}$. As shown in § 3.2, the Floquet exponent then satisfies
\begin{equation} \lambda = \frac{1}{\rho S} \frac{|\langle \hat{F}_x\rangle|}{|\langle\hat{u}_g\rangle|} \textrm{e}^{\textrm{{i}} \psi}, \end{equation}
where 
$\psi =\phi _{F}-\phi _{U}$ is the phase difference between the force and velocity phases. The asynchronous mode is unstable (
$\mathrm {Re}(\lambda )>0$) when this phase difference lies in the interval 
$-{\rm \pi} /2 < \psi < {\rm \pi}/2$. The crucial role of the phase difference (rather than the force-to-velocity ratio) in the asynchronous mode destabilization is confirmed by examining figures 17(a) and 17(b) where both quantities are depicted as a function of the Stokes number. The mode gets unstable (grey area) precisely when the phase difference is 
$\phi <{\rm \pi} /2$. To better understand the physical meaning of this phase difference, the temporal evolution of the quasi-periodic perturbation of the time-averaged velocity 
$\langle \tilde {u}_g\rangle$ (solid) and force 
$\langle \tilde {F}_x\rangle$ (dashed) are plotted as a function of time in figure 17(c) for a stable mode and in figure 17(d) for an unstable mode. Note that the velocity and force are time-averaged over the (short) flapping period and their evolution is depicted over the long period (slow time scale 
$t_{\epsilon }$). In both figures, the green areas identify phases of motion where the hydrodynamic force is propulsive (since the force and velocity are of the same sign), while the white areas correspond to resistive phases of motion (force and velocity are of opposite signs). In the case of a stable mode (see figure 17c) the phase difference slightly above 
${\rm \pi} /2$ results in a mode with propulsive phases of motion that are shorter than resistive ones. On the other hand, in the case of an unstable mode (see figure 17d) the phase difference slightly under 
${\rm \pi} /2$ results in a mode with longer propulsive phases of motion than resistive ones. The phase difference between the horizontal force and velocity is thus related to the cumulative time of the propulsive phase over the resistive phase. When 
$-{\rm \pi} /2 < \psi < {\rm \pi}/2$, the propulsive phases last longer than resistive ones, and the asynchronous mode is unstable. A similar criterion was established by Navrose & Mittal (Reference Navrose and Mittal2016) for the instability threshold of a spring-mounted cylinder flow, based on the global stability analysis of the steady base flow solution. They showed that the phase difference between the vertical hydrodynamic force and displacement of the cylinder perturbation is related to the transfer of energy from the flow to the cylinder and drives the destabilization of the mode. The present criterion can be viewed as an extension to the instability of periodic fluid–solid interaction problems.
Figure 17. (a) Amplitude ratio and (b) phase difference of the (time-averaged) horizontal force and foil velocity of the asynchronous modes as a function of the Stokes number. The instability region is marked in grey. The dashed line corresponds to 
$\psi ={\rm \pi} /2$. (c,d) Temporal evolution of the quasi-periodic foil velocity 
$\langle \tilde {u}_g\rangle$ (solid curve) and force 
$\langle \tilde {F}_x \rangle$ (dashed curve) time-averaged along a flapping period for (c) a stable (
$\beta =10.5$) and (d) an unstable (
$\beta =13$) asynchronous mode, shown in (a,b) with black dots. The green area corresponds to propulsive phases (velocity and force have the same sign) while the white area corresponds to resistive phases (velocity and force are of opposite sign).
3.5. Effect of the fluid–solid density ratio
Before concluding, we investigate the effect of the fluid–solid density ratio 
$\rho$ on the results of the Floquet analysis. Two limit cases are considered in this section: high density ratios 
$\rho \gg 100$ that lead to a loosely coupled fluid–solid due to the vanishing action of the fluid over the solid problem (as presented in § 3.1), and the range of lower density ratios 
$\rho < 100$ that tend to the one of swimming organisms.
Let us first consider the high density ratio limit. The evolution of the Floquet exponent is shown in figures 18(a) and 18(b) for the synchronous and asynchronous modes, respectively. Starting from the value 
$\rho =100$ considered until now, and increasing the density ratio, the absolute value of the Floquet multipliers (black solid curve) decreases for both modes. However, their asymptotic behaviour in the limit 
$\rho \rightarrow \infty$, displayed with the dashed red curves in the two figures, is different. The synchronous mode evolves as 
$1 / \rho$ and, thus, only gets marginally stable. The asynchronous mode is stabilized for 
$\rho >10^3$ and its growth rate tends towards that of the purely hydrodynamic analysis.
Figure 18. Effect of the fluid–solid density ratio 
$\rho$ on (a,b) the Floquet exponent and (c,d) the perturbation of the foil horizontal speed for (a,c) the synchronous mode (
$\beta =6$) and (b,d) the asynchronous mode (
$\beta =12$). In (b,d) the real and imaginary parts of the Floquet exponent and the horizontal speed are respectively represented by solid and dash–dot lines. The solid curves correspond to results of the fluid–solid analysis. The red dashed curve corresponds to the asymptotic limits of the exponent: (a) 
$1/\rho$ curve and (b) the values of the purely hydrodynamic analysis 
$\mathrm {Re}(\lambda )=-0.01$ and 
$\mathrm {Im}(\lambda )=0.0043$. The dashed horizontal lines delimit in (a,b) 
$\mathrm {Re}(\lambda )=0$ and in (b) 
$\langle u_g\rangle =0$.
To further understand why the Floquet exponents of the synchronous and asynchronous mode exhibit different behaviours in that limit case, we propose to reconsider the time-averaged analysis. In the infinite density ratio limit, the relation (3.10), that links the growth rate to the mean value of the horizontal solid velocity and force, gets
\begin{equation} \lambda \langle \hat{u}_g\rangle = 0. \end{equation}
Thus, either the Floquet exponent or the mean horizontal velocity is zero. The synchronous mode corresponds clearly to the case 
$\lambda =0$ (figure 18a). The foil mean velocity does not necessarily vanish for high density ratios, as observed in figure 18(c) that displays the evolution of the time-averaged horizontal speed with the density ratio. The asynchronous modes correspond to the second case 
$\langle \hat {u}_g\rangle = 0$, as seen in figure 18(d). As a matter of fact, not only the mean horizontal velocity converges to zero, but so does the entire temporal evolution due to the negligible acceleration generated by high density ratios. In this case, the Floquet exponent does not tend to zero (see figure 18b), but rather to the value predicted by the purely hydrodynamic stability analysis (red dashed line). Indeed, in the limit 
$\rho \rightarrow \infty$, the fluid–solid linearized operator is block triangular as seen in (3.5) and the purely hydrodynamic Floquet multipliers are included in the fluid–solid Floquet spectrum. A similar asymptotic behaviour was observed by Fabre, Assemat & Magnaudet (Reference Fabre, Assemat and Magnaudet2011) when investigating the dynamics of free falling bodies in fluids using a fluid–solid stability analysis of the steady base solution (not periodic as in the present case). In the limit case, results of the purely hydrodynamic and fluid–solid stability analyses converged to the well know Von–Kármán vortex wake instability without an effect on the path of the falling body. The present result is somehow an extension to the fluid–solid stability analysis of a periodic solution.
Finally, figure 19 presents results of the fluid–solid stability analysis in the parameter space 
$(\beta ,\rho )$. The white area corresponds to stable regions, while the red and grey areas correspond to regions of unstable synchronous and asynchronous modes, respectively. We first consider the asymptotic limit of very large density ratios 
$\rho \gg 100$. As one can see in figure 19(a), the synchronous unstable mode (red region) thresholds are barely modified for high density ratios, always taking place for 
$\beta =4$ and 
$\beta =9.53$ for 
$\rho >100$. The instability threshold of asynchronous modes (grey region), on the other hand, is strongly modified with the increase of 
$\rho$ going from 
$\beta =11.25$ for 
$\rho =100$ to 
$\beta =12.2$ for 
$\rho =10 ^4$. For Stokes numbers between these two thresholds, the density ratio thus presents a stabilizing effect. Accordingly to our previous analysis, as 
$\rho \rightarrow \infty$, the asynchronous modes threshold tends to the purely hydrodynamic analysis one (displayed at the top of the figure). Synchronous modes, on the other hand, become marginally unstable and do not converge to the purely hydrodynamic analysis, their marginal instability explaining the minimal threshold changes for high density ratios.
Figure 19. (a) Asymptotic limit (
$1 \leq \rho \leq 10^9$) and (b) close-up on low density ratios (
$1 \leq \rho \leq 100$) of the instability regimes of the symmetric non-propulsive periodic solution identified with the fluid–solid Floquet analysis in the parameters plane (
$\beta ,\rho$). White regions correspond to stable solutions, while red and grey regions indicate unstable solutions to synchronous and asynchronous Floquet modes, respectively. Results of the purely hydrodynamic stability analysis are displayed at the top of (a) for comparison.
For small density ratios 
$1 \leq \rho \leq 100$ (figure 19b), the synchronous mode thresholds are again not strongly modified with the 
$\rho$ decrease. While the transition from stable to unstable synchronous modes remains constant at 
$\beta =4$, the stabilization threshold presents a significant variation going from 
$\beta =9.42$ for 
$\rho =1$ to 
$\beta =9.53$ for 
$\rho =100$. The instability threshold of the asynchronous mode is strongly modified, as for high density ratios. Its onset is now encouraged rather than delayed, varying from 
$\beta =11.25$ for 
$\rho =100$ to 
$\beta =10.6$ for 
$\rho =1$. These results indicate that swimming organisms (nearby a unity density ratio) are significantly more prone to non-coherent motions than flying organisms.
All instability thresholds and their trends are coherent with the unsteady nonlinear results previously reported. As shown for 
$\rho =1$ in table 2, the destabilisation of the synchronous mode and the onset of unidirectional propulsion both take place for 
$\beta =4$. The threshold value 
$\beta =9.42$ of the synchronous mode stabilization is again slightly different from the threshold 
$\beta =9.50$ of the transition between propulsive and non-propulsive solutions due to the bifurcation sub-critical nature (see figure 9b). The threshold value 
$\beta =9.42$ found with the linear stability corresponds, as in the case 
$\rho =100$, to the disappearance of the symmetric solution when decreasing the Stokes number (transition from III to II). The onset of back and forth solutions matches again the destabilization of the asynchronous Floquet mode for 
$\beta =10.6$.
Table 2. Critical thresholds obtained with the unsteady nonlinear simulations and the fluid–solid stability analysis for 
$\rho =1$.
4. Conclusions
The role of linear mechanisms in the transition to horizontal locomotion of a vertically flapping foil has been investigated. First, the occurrence of non-propulsive, unidirectional propulsive and back and forth solutions was established in the range of Stokes numbers 
$\beta \in [1,20]$ for a rectangular shaped foil with an aspect ratio 
$h=1/20$, a flapping amplitude 
$A=0.5$ and a solid–fluid density ratio 
$\rho =100$. Floquet stability analysis of the coupled fluid–solid system was thus performed over symmetry preserved non-propulsive solutions. Our study was finally concluded by analysing the effect of the solid–fluid density ratio on the stability analysis, this last study being compared with predictions of non-propelled foil stability analysis usually employed in the literature.
First, symmetric non-propulsive, unidirectional propulsive and back and forth solutions were obtained while raising 
$\beta$, as in Alben & Shelley (Reference Alben and Shelley2005) and Lu & Liao (Reference Lu and Liao2006). As expected for a low aspect ratio foil (Deng & Caulfield Reference Deng and Caulfield2016), non-propulsive solutions first transition to unidirectional propulsion. The results presented in this paper highlight the existence of a sub-critical transition between propulsive and again non-propulsive solutions (regimes II–III), with back and forth oscillations finally reached (regime IV) for higher 
$\beta$. The emergence of these nonlinear solutions was then investigated through a self-propelled stability analysis of non-propulsive solutions. This analysis revealed the existence of unstable synchronous and asynchronous Floquet modes in the region of unidirectional and back and forth solutions, respectively. We therefore studied the characteristics of these unstable modes, investigating their associated fluid flow, horizontal force and speed to finally relate these mechanisms to the locomotion regimes obtained in nonlinear simulations. The evolution of the modes mean horizontal force and speed with 
$\beta$ allowed us to establish a criterion of instability that link these quantities to the Floquet exponent.
In the case of synchronous modes, spatial symmetry breaking modes with non-zero force and speed are obtained, similar to unidirectional propulsion. Hydrodynamic forces that accelerate the horizontal speed lead to unstable Floquet exponents, and the transition to unstable modes is driven by the increase of the hydrodynamic force. The decay of this horizontal force leads subsequently to the mode stabilization and the re-emergence of symmetric non-propulsive solutions. Concerning asynchronous modes, the direction switching phenomenon observed in nonlinear solutions is explained by the competing action of the complex and real parts of these modes associated to the complex phase of their multiplier. This complex multiplier accurately predicts the low frequency of back and forth solutions at their onset. The temporal evolution of the quasi-periodic perturbation, resulting from the superposition of the real and imaginary parts of the Floquet mode, clearly shows how the horizontal force exerted by the flow perturbation is alternatively propulsive or resistive, i.e of the same or opposite sign to the foil velocity, acting as a restoring spring-like force on the foil position over the slow period. The destabilization of the asynchronous modes depends on the phase difference between the time-averaged force and velocity perturbation, as it measures how long this force is propulsive or resistive over the slow period. A generalization to three-dimensional foils and different flapping movements of this mode might offer an additional path to understand, for example, the snaking trajectory presented by Deng & Caulfield (Reference Deng and Caulfield2018) in the horizontal locomotion of oblate spheroids, first supposed to be connected to nonlinear effects, or the non-coherent motion of living organisms as the planktonic sea butterfly (Murphy et al. Reference Murphy, Adhikari, Webster and Yen2016).
The influence of the solid–fluid density ratio was finally investigated. For 
$\rho \in [1,10^9]$, we observed that the 
$\beta$ range of synchronous modes is not largely affected by the density ratio, whereas the transition to asynchronous modes is greatly impacted. To understand this behaviour, we have studied the coupled stability equations in the limit of high density ratios. We have shown that synchronous modes are asymptotically marginally unstable, whereas asynchronous modes converge to the uncoupled fluid system, since their horizontal speed converges to zero. These results explain the observed sensibility to the density ratio of the asynchronous modes threshold and the small variation of the transition to synchronous modes. Comparing these results with the non-propelled stability analysis, we have highlighted that the latter converges to the self-propelled analysis only for asynchronous modes in the limit of very large density ratios, since in this limit the fluid–solid coupling terms disappear.
We conclude thus that the studied fluid–solid Floquet stability analysis, that takes into account the inherently fluid–solid coupling of the studied self-propelled foil interacting with a viscous fluid, is best adapted to predict the onset of unidirectional and back and forth propulsion. Possible future applications of this coupled stability analysis are the return from back and forth to coherent unidirectional propulsion observed for higher stokes numbers (Alben & Shelley Reference Alben and Shelley2005), as well as bifurcations of self-propelled heaving foils passively pitching around their leading edge (Spagnolie et al. Reference Spagnolie, Moret, Shelley and Zhang2010) and of a self-propelled infinite array of flapping wings (Becker et al. Reference Becker, Masoud, Newbolt, Shelley and Ristroph2015).
Acknowledgements
This work was initiated during the PhD thesis of D. Jallas, co-supervised by D. Fabre and O.M., who gratefully acknowledges both of them. This project has received funding from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement No. 638307).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Segregated approach for solving the implicitly coupled fluid–solid problem
We describe here the segregated approach used to solve efficiently the temporally discretized (2.7). Manipulating the 
$r$-backward differential formula of the time derivative, we may split the unknown solid velocity 
$u_g^{n+1}$ as
\begin{equation} u_g^{n+1} = (\Delta t / \alpha_0) a_{g}^{n+1} + \hat{u}_{g} , \end{equation}
where the first term is an unknown component, proportional to the unknown acceleration 
$a_{g}^{n+1}=( \textrm {d} u_g /\textrm {d}t )^{n+1}$, and the second term is the known velocity component defined as
\begin{equation} \hat{u}_{g}={-}\sum_{k=1}^r (\alpha_k / \alpha_0) \, u_g^{n+1-k}. \end{equation}The above decomposition of the horizontal solid velocity is then used to split the fluid variables as
\begin{equation} (\boldsymbol{u}^{n+1},p^{n+1}) = \left( \frac{\Delta t}{\alpha_0} a_g^{n+1} \right) \left(\delta \boldsymbol{u}, \delta p\right) + (\hat{\boldsymbol{u}}^{n+1},\hat{p}^{n+1}), \end{equation}
where we have introduced the flow component 
$(\boldsymbol {\delta u},\delta p)$, proportional to the solid acceleration, and the flow component 
$(\boldsymbol {\hat {u}}^{n+1},\hat {p}^{n+1})$ that will depend on the solid velocity 
$\hat {u}_{g}$. Introducing the solid (A1) and fluid (A3) decomposition into the governing flow equations (2.7), we obtain two independent linear system of equations. The first one governs the flow component 
$(\boldsymbol {\hat {u}}^{n+1},\hat {p}^{n+1})$ as
\begin{equation} \left. \begin{gathered} \frac{\alpha_0}{\Delta t} \boldsymbol{\hat{u}}^{n+1} + \boldsymbol{\nabla}{\hat{p}}^{n+1} -\beta^{{-}1} \Delta \boldsymbol{\hat{u}}^{n+1} = \boldsymbol{f}^{n+1},\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\hat{u}}^{n+1} = 0, \\ (\hat{u}^{n+1},\hat{v}^{n+1})(\varGamma_w) = \big(\hat{u}_{g},2 {\rm \pi}A \sin(2 {\rm \pi}t^{n+1})\big).\end{gathered} \right\} \end{equation}
The boundary conditions at the fluid–solid interface 
$\varGamma _{w}$ is explicitly known. Therefore, any classical algorithm such as the Uzawa method or the projection splitting method can be used to obtain the solution of this discretized forced unsteady Stokes equation. Following Jallas et al. (Reference Jallas, Marquet and Fabre2017), we have used a preconditioned conjugate gradient algorithm (Cahouet & Chabard Reference Cahouet and Chabard1988) to impose the divergence-free condition. As it depends on the flow history through the right-hand side forcing terms 
$\boldsymbol {f}^{n+1}$ and the boundary conditions, it is solved at each temporal iteration.
The second problem governs the flow component 
$(\boldsymbol {\delta u},\delta p)$ that is proportional to the solid acceleration in the horizontal direction. It writes
\begin{equation} \left. \begin{gathered} \frac{\alpha_0}{\Delta t} \boldsymbol{\delta u} + \boldsymbol{\nabla}{\delta p}^{n+1} -\beta^{{-}1} \Delta \boldsymbol{\delta u}^{n+1} = 0,\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\delta u} = 0, \\ \left(\delta u,\delta v\right)\left(\varGamma_w\right) = (1,0). \end{gathered} \right\} \end{equation}
Again, the boundary conditions at the fluid–solid interface 
$\varGamma _{w}$ is explicitly known, but it is now independent from the temporal iteration. The solution can thus be obtained prior to the temporal iteration. It can be viewed as the short time response of a Stokes flow, initially at rest, to a unitary horizontal velocity.
The solution of the (above) two independent flow problems does not give access to 
$\left (\boldsymbol {u}^{n+1},p^{n+1}\right )$ in (A3) since the horizontal acceleration 
$a_{g}^{n+1}$ is still unknown. The final step of the algorithm is obtained by introducing this decomposition into the last equation of the governing equation (2.7), yielding
\begin{equation} \left(\rho S\right) \, a_{g}^{n+1} = \left( \frac{\Delta t}{\alpha_0} a_{g}^{n+1} \right) F_x\left(\boldsymbol{\delta u}, \delta p\right) + F_x\left(\hat{\boldsymbol{u}},\hat{p}\right). \end{equation}The horizontal acceleration is thus given by
\begin{equation} a_{g}^{n+1} = \left( \rho S - \frac{\Delta t}{\alpha_0} F_x(\boldsymbol{\delta u}, \delta p) \right)^{{-}1} F_x(\hat{\boldsymbol{u}},\hat{p}), \end{equation}and the velocity and pressure are obtained using (A3).
Appendix B. Validation of the nonlinear and linear fluid–solid solvers
The numerical method is primarily validated through the convergence of the time-averaged horizontal velocity and the maximal vertical force acting on the foil with the mesh number of elements and time step. In table 3 the influence of the time and spatial discretization is evaluated for a unidirectional propulsive solution (
$\beta =6$). A time step of 
$\Delta t = 0.005$ and a mesh of 17 116 triangles are sufficient to guarantee the convergence of the hydrodynamic force and horizontal velocity up to order 
${O}(10^{-3})$. This convergence was attained for all other Stokes numbers explored in this article. The mesh of 17 116 was used in all cases and the time step evolved from 
$\Delta t=10^{-2}$ for small values of the Stokes number (
$\beta =2$), being decreased to 
$\Delta t=5\times 10^{-4}$ for larger values (
$\beta =19$).
Table 3. Convergence of the maximal vertical force 
$\max (F_y)$ and the time-averaged horizontal velocity 
$\langle u_g\rangle$ with a time step and a mesh number of triangles for a unidirectional propulsive solution (
$\beta =6$). For the time-step convergence, a fixed mesh of 
$17116$ triangles was adopted, whereas for the mesh convergence, a fixed time step 
$\Delta t = 0.005$ was used.
The nonlinear solver employed in this study has then been validated by simulating the horizontal locomotion of a two-dimensional ellipsoid of aspect ratio 
$h=0.1$, flapping amplitude 
$A=0.25$, density ratio 
$\rho =10$ and Stokes number 
$\beta \in [2,15]$, as in the numerical study of Spagnolie et al. (Reference Spagnolie, Moret, Shelley and Zhang2010). In this work, the onset of locomotion is around 
$\beta >8$, and as seen in figure 20, both the emergence of propulsive solutions and their time-averaged horizontal velocities compare very well between the present numerical method and the values found on the reference.
Figure 20. Validation of the nonlinear solver. Time-averaged horizontal velocity of the foil as a function of the Stokes number 
$\beta$. Results of the numerical algorithm described in appendix A (filled dots) are compared with results (empty squares) extracted from figure 13 (clamped case) in Spagnolie et al. (Reference Spagnolie, Moret, Shelley and Zhang2010).
The linear solver is validated through the purely hydrodynamic Floquet stability analysis of the flow symmetry breaking around a heaving cylinder. Two distinct flapping amplitudes and Stokes numbers are considered (using the Keulegan–Carpenter number 
$KC=4{\rm \pi} A$): 
$(\beta ,KC)=(40,4.75)$ and 
$(\beta ,KC)=(100,3.65)$. We can see in the table 4 that the absolute value of the leading Floquet multiplier obtained in the two test cases is in good agreement with values of Elston et al. (Reference Elston, Blackburn and Sheridan2006).
Table 4. Linear solver validation: comparison of the leading Floquet multiplier 
$|\mu |$ with the values presented in table 3 of Elston et al. (Reference Elston, Blackburn and Sheridan2006).
Appendix C. Effect of the foil shape on self-propelled regimes and stability
To evaluate the influence of the foil shape and aspect ratio, the evolution with the Stokes number 
$\beta$ of the unsteady nonlinear dynamics and fluid–solid Floquet stability of an elliptical foil of aspect ratio 
$h=0.01$ is shown in figure 21. The flapping amplitude 
$A=0.5$ (identical to this work) and the density ratio 
$\rho =32$ have been fixed as to approach one of the configurations explored by Alben & Shelley (Reference Alben and Shelley2005).
Figure 21. Evolution of the (a) time-averaged horizontal velocity and (b) leading Floquet multiplier absolute value of an elliptical foil of minor/major axis aspect ratio 
$h=0.1$ with control parameters 
$A=0.5$ and 
$\rho =32$. In (a) – respectively (b) – white, red and grey background colours identify symmetric non-propulsive, unidirectional propulsive and back and forth regimes – respectively stable, synchronous unstable and asynchronous unstable multipliers.
The self-propelled regimes and their transition are similar to the rectangular foil with rounded edges. A remarkable difference is, nevertheless, the suppression of the intermediary non-propulsive regime III, between the unidirectional propulsive and the back and forth regimes. Apparently the increase of the aspect ratio favours, as the decrease of the density ratio (§ 2.3.2), the onset of non-coherent propulsion. The obtained onset of back and forth solutions (
$\beta >9.5$) closely matches the one of Alben & Shelley (Reference Alben and Shelley2005). The existence of a unidirectional propulsive regime prior to the back and forth one is, however, to our knowledge newly reported in the literature (apart from being briefly mentioned in the works of Deng & Caulfield Reference Deng and Caulfield2016). We suspect this regime has not yet been characterized due to two factors. On one hand, as illustrated in figure 4, this propulsive regime features a very long transitory regime of 
${\sim }200$ flapping periods. In addition, owing to the low Stokes number and its viscous nature, velocity perturbations added to the quiescent system initially decay. Following this initial decay, more than 
$100$ flapping periods are commonly needed for the horizontal velocity to grow above an initial perturbation 
$u_g=0.05$. This threshold is indeed beyond the one employed, for example, by Alben & Shelley (Reference Alben and Shelley2005) who differentiated non-propulsive and propulsive solutions by the growth of initial velocity perturbations after 
$80$ flapping periods.
Unstable synchronous and asynchronous modes, figure 21(b), are again obtained in the same 
$\beta$ range as the unidirectional propulsive and back and forth solutions. Despite the different geometry and the suppression of the intermediary non-propulsive regime, the fluid–solid Floquet stability analysis correctly predicts the onset of unidirectional propulsive and back and forth solutions.
