1. Introduction
Flutter of membranes has been a classical subject for at least seven decades. Membranes submerged in a compressible gas flow occupying a space or a semi-space and their flutter at supersonic speeds have been considered already in the works by Miles (Reference Miles1947, Reference Miles1956), Goland & Luke (Reference Goland and Luke1954), Benjamin (Reference Benjamin1963) and Bolotin (Reference Bolotin1963).
Bolotin (Reference Bolotin1963), Spriggs, Messiter & Anderson (Reference Spriggs, Messiter and Anderson1969), Dowell & Ventres (Reference Dowell and Ventres1970) and Kornecki, Dowell & O'Brien (Reference Kornecki, Dowell and O'Brien1976) addressed the problem of the so-called membrane flutter paradox regarding the relation of stability criteria for an elastic plate to those for a membrane. Gislason (Reference Gislason1971) demonstrated both theoretically and experimentally that a membrane or elastic plate with a finite chord develops not only flutter but also a divergence instability.
Dowell (Reference Dowell1966), when critically appraising the study by Miles (Reference Miles1956) of an infinitely long, infinitely wide panel in a compressible flow occupying the upper semi-space, pointed out that the critical wavelength predicted in this study was infinite and the flutter velocity was zero, which was not physically meaningful. This observation has led him to the conclusion that the finite dimension of a membrane or a plate in the flow or a span direction is critical to the physically meaningful prediction of the instability (Dowell Reference Dowell1966). A similar effect of an elastic foundation was shown both theoretically and experimentally by Dugundji, Dowell & Perkin (Reference Dugundji, Dowell and Perkin1963).
Absolute and convective hydroelastic instabilities of slender elastic structures submerged in a uniform flow were discussed by Triantafyllou (Reference Triantafyllou1992). A comprehensive monograph by Dowell (Reference Dowell2015) is a standard reference in the field.
Recent works on membrane flutter are motivated by such diverse applications as stability of membrane roofs in civil engineering (Sygulski Reference Sygulski2007), flutter of travelling paper webs (Banichuk et al. Reference Banichuk, Jeronen, Neittaanmäki and Tuovinen2010, Reference Banichuk, Barsuk, Jeronen, Tuovinen and Neittaanmäki2019), aerodynamics of sails and membrane wings of natural flyers (Newman & Paidoussis Reference Newman and Paidoussis1991; Tiomkin & Raveh Reference Tiomkin and Raveh2017), as well as the design of piezoaeroelastic systems for energy harvesting (Mavroyiakoumou & Alben Reference Mavroyiakoumou and Alben2020).
Surface gravity waves on a motionless fluid of finite depth are a classical subject as well, going back to the seminal studies of Russell and Kelvin (Carusotto & Rousseaux Reference Carusotto, Rousseaux, Faccio, Belgiorno, Cacciatori, Gorini, Liberati and Moschella2013). Numerous generalizations are known taking into account, for instance, a uniform or a shear flow and surface tension (Maissa, Rousseaux & Stepanyants Reference Maissa, Rousseaux and Stepanyants2016), submerged solids (Smorodin Reference Smorodin1972; Arzhannikov & Kotelnikov Reference Arzhannikov and Kotelnikov2016) and hydrofoils (Faltinsen & Semenov Reference Faltinsen and Semenov2008), a flexible bottom (Mohapatra & Sahoo Reference Mohapatra and Sahoo2011) or a flexible plate resting on a free surface (Greenhill Reference Greenhill1886; Schulkes, Hosking & Sneyd Reference Schulkes, Hosking and Sneyd1987; Bochkarev, Lekomtsev & Matveenko Reference Bochkarev, Lekomtsev and Matveenko2016; Das, Sahoo & Meylan Reference Das, Sahoo and Meylan2018a,Reference Das, Sahoo and Meylanb; Das et al. Reference Das, Kar, Sahoo and Meylan2018). The latter setting has a straightforward motivation in the dynamics of sea ice and a less obvious application in the analogue gravity experiments (Barcelo, Liberati & Visser Reference Barcelo, Liberati and Visser2011; Weinfurtner et al. Reference Weinfurtner, Tedford, Penrice, Unruh and Lawrence2011; Carusotto & Rousseaux Reference Carusotto, Rousseaux, Faccio, Belgiorno, Cacciatori, Gorini, Liberati and Moschella2013). Recent work (Robertson & Rousseaux Reference Robertson and Rousseaux2018) discusses the effects of viscous dissipation of surface gravity waves to the analogue gravity.
Remarkably, another phenomenon that is analysed from the analogue gravity perspective is super-radiance (Barcelo et al. Reference Barcelo, Liberati and Visser2011; Carusotto & Rousseaux Reference Carusotto, Rousseaux, Faccio, Belgiorno, Cacciatori, Gorini, Liberati and Moschella2013; Brito, Cardoso & Pani Reference Brito, Cardoso and Pani2015) and its particular form, discovered by Ginzburg & Frank (Reference Ginzburg and Frank1947) and Ginzburg (Reference Ginzburg1996), known as the anomalous Doppler effect (ADE) (Nezlin Reference Nezlin1976; Nemtsov & Eidman Reference Nemtsov and Eidman1987; Bekenstein & Schiffer Reference Bekenstein and Schiffer1998). In electrodynamics, the ADE manifests itself when an electrically neutral overall particle, endowed with an internal structure, becomes excited and emits a photon during its uniform but superluminal motion through a medium, even if it started the motion in its ground state; the energy source is the bulk motion of the particle (Bekenstein & Schiffer Reference Bekenstein and Schiffer1998).
The anomalous Doppler effect in hydrodynamics was demonstrated for a mechanical oscillator with one degree of freedom, moving parallel to the border between two incompressible fluids of different densities (Gaponov-Grekhov, Dolina & Ostrovskii Reference Gaponov-Grekhov, Dolina and Ostrovskii1983). It was shown that the oscillator becomes excited due to radiation of internal gravity waves if it moves sufficiently fast. In Abramovich, Mareev & Nemtsov (Reference Abramovich, Mareev and Nemtsov1986) the ADE for such an oscillator was demonstrated due to radiation of surface gravity waves in a layer of an incompressible fluid.
Nemtsov (Reference Nemtsov1985) was the first who considered flutter of an elastic membrane resting at the bottom of a uniform horizontal flow of an inviscid and incompressible fluid as an anomalous Doppler effect due to emission of long surface gravity waves. In the shallow water approximation, he investigated both the case of a membrane that spreads infinitely far in both horizontal directions and the case when the length of the membrane in the direction of the flow (or the chord length) is finite whereas the span in the perpendicular direction is infinite. Nevertheless, the case of flow of arbitrary depth has not been studied in Nemtsov (Reference Nemtsov1985), and no numerical computation supporting the asymptotic results has been performed. Another issue that has not been addressed in Nemtsov (Reference Nemtsov1985) is the relation of stability domains for the membrane of finite length to those for the membrane of infinite length.
Vedeneev (Reference Vedeneev2004) studied flutter of an elastic plate of finite and infinite length at the bottom of a uniform horizontal flow of a compressible gas occupying the upper semi-space. He performed an analysis of the relation of the stability conditions for the finite plate to those for the infinite plate using the method of global stability analysis of Kulikovskii (Doaré & de Langre Reference Doaré and de Langre2006; Vedeneev Reference Vedeneev2016). A single-mode high-frequency flutter due to a negative aerodynamic damping and a binary flutter have been identified in Vedeneev (Reference Vedeneev2016). However, no connection has been made to the ADE and the concept of negative energy waves.
In the present work we reconsider the setting of Nemtsov in order to address the finite depth of the fluid layer, find flutter domains in the parameter space, analyse them using perturbation of multiple roots of the dispersion relation and investigate the flutter onset for the membrane of infinite chord length. We will explain the radiative instabilities via the interaction of positive and negative energy waves using an explicit expression for the averaged total energy derived rigorously from physical considerations and relate them to the anomalous Doppler effect. We believe that the Nemtsov membrane is as important for understanding the phenomenon of radiation-induced instabilities (Hagerty, Bloch & Weinstein Reference Hagerty, Bloch and Weinstein2003) as the famous Lamb oscillator coupled to a semi-infinite string was for understanding the radiative damping (Lamb Reference Lamb1900; Barbone & Crighton Reference Barbone and Crighton1994).
2. Model of a membrane interacting with a free surface
2.1. Physical system
In a Cartesian coordinate system 
$OXYZ$, consider an inextensible elastic rectangular membrane strip of constant thickness 
$h$ and density 
$\rho _m$, of infinite span in the 
$Y$-direction, held at 
$Z=0$ at the leading edge 
$(X=0)$ and at the trailing edge 
$(X=L)$ by simple supports, figure 1.
Figure 1. An elastic membrane with chord of length 
$L$ attached to two rigid walls along its leading 
$(X=0)$ and trailing 
$(X=L)$ edges on the bottom of a fluid layer of depth 
$H$ moving with the velocity 
$v$. 
$\Omega$ is the fluid domain and 
$\partial \Omega _0$, 
$\partial \Omega _1$ and 
$\partial \Omega _2$ are respectively the free surface, membrane and rigid wall boundaries.
The membrane is initially still and flat, immersed in a layer of inviscid, incompressible fluid of constant density 
$\rho$, with free surface at the height 
$Z=H$. The two-dimensional flow in the layer is supposed to be irrotational and moving steadily with velocity 
$v$ in the positive 
$X$-direction. The bottom of the fluid layer at 
$Z=0$ is supposed to be rigid and flat for 
$X\in (-\infty ,0] \cup [L,+\infty )$.
Nemtsov (Reference Nemtsov1985) assumed that a vacuum exists below the membrane. In the present study we prefer to consider that a motionless incompressible medium of the same density 
$\rho$ is present below the membrane with a pressure that is the same as the unperturbed pressure of the fluid (Vedeneev Reference Vedeneev2004, Reference Vedeneev2016).
Assuming small vertical displacement of the membrane 
$w(X,t)$, where 
$t$ is time, a constant tension 
$T$ along the membrane profile and neglecting viscous forces, we write the dimensional membrane dynamic equation as (Tiomkin & Raveh Reference Tiomkin and Raveh2017)
\begin{equation} \rho_m h \partial^2_{t}w=T\partial^2_{X}w-{\rm \Delta} P,\quad X\in[0,L], \ Z=0, \end{equation}
where 
${\rm \Delta} P(X,t)$ is the pressure difference across the interface 
$Z=0$. The simply supported boundary conditions for the membrane are
\begin{equation} w(0)=w(L)=0 \quad \text{at } Z=0. \end{equation} In general, to recover the pressure 
$P(X,Z,t)$ of the fluid we write the Euler equation for the vorticity-free flow (Carusotto & Rousseaux Reference Carusotto, Rousseaux, Faccio, Belgiorno, Cacciatori, Gorini, Liberati and Moschella2013; Maissa et al. Reference Maissa, Rousseaux and Stepanyants2016)
\begin{equation} \partial_t \boldsymbol{v}+\boldsymbol{\nabla}\left(\frac{\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{v}}{2}+\frac{P}{\rho}+gZ \right)=0 \end{equation} with 
$\boldsymbol {v}=v \boldsymbol {e}_X+\boldsymbol {\nabla } \varphi$, where 
$\varphi (X,Z,t)$ is the potential of the fluid, 
$\boldsymbol {e}_X$ is the unit vector in the 
$X$-direction and 
$g$ stands for the gravitational acceleration. This yields the integral of Bernoulli
\begin{equation} \frac{P}{\rho} + (\partial_t + v\partial_X)\varphi + \frac{1}{2}\boldsymbol{\nabla}\varphi\boldsymbol{\cdot}\boldsymbol{\nabla}\varphi + gZ = \text{const.} \end{equation}The incompressibility condition takes the form
\begin{equation} \nabla^2 \varphi=0. \end{equation} From (2.4) it follows that, in the case when a motionless medium of density 
$\rho$ is present below the membrane with its pressure equal to the unperturbed pressure of the fluid above the membrane, the linear in 
$\varphi$ expression for the pressure difference, 
${\rm \Delta} P(X,t)$, is
\begin{equation} {\rm \Delta} P(X,t) = -\rho(\partial_t + v\partial_X)\varphi(X,0,t). \end{equation}For the sake of completeness, we present also the analogous expression for the pressure difference for the case when there is a vacuum below the membrane (Nemtsov Reference Nemtsov1985)
\begin{equation} {\rm \Delta} P(X,t) = -\rho(\partial_t + v\partial_X)\varphi(X,0,t) - \rho g w(X,t). \end{equation}Impermeability of the rigid bottom implies the condition
\begin{equation} \boldsymbol{\nabla}\varphi\boldsymbol{\cdot}\boldsymbol{n} = 0 \quad \text{at } Z=0, \ X\in (-\infty,0] \cup [L,+\infty). \end{equation}The prescription of the normal velocity at the boundaries of moving surfaces allows us to express the kinematic condition for the membrane
\begin{equation} \boldsymbol{\nabla}\varphi\boldsymbol{\cdot}\boldsymbol{n} = -(\partial_t + v\partial_X)w \quad \text{at } Z=0, \ X\in [0,L], \end{equation}and to specify the same condition at the free surface
\begin{equation} \boldsymbol{\nabla}{\varphi} \boldsymbol{\cdot} \boldsymbol{n} = (\partial_t + v\partial_X)u, \end{equation}
where 
$u(X,t)$ is the free surface elevation and 
$\boldsymbol {n}$ is the vector of the outward normal to a surface. This implies that the projection of the vector 
$\boldsymbol {\nabla }\varphi$ to the normal will coincide with the positive 
$z$-direction for the free surface and have the opposite direction for the membrane, see figure 1.
Using the Bernoulli integral (2.4) at the free surface where 
$P=0$ and retaining only linear in 
$\varphi$ terms, we find
\begin{equation} g u = -(\partial_t + v \partial_X)\varphi. \end{equation}
Taking 
$u$ from (2.11) and substituting it into (2.10) we obtain the boundary condition at the free surface of the liquid that reads
\begin{equation} \boldsymbol{\nabla}\varphi\boldsymbol{\cdot}\boldsymbol{n} = -\frac{1}{g}(\partial_t + v\partial_X)^2\varphi \quad \text{at } Z=H. \end{equation}2.2. Dimensionless mathematical model
Let us choose the height of the fluid layer, 
$H$, as a length scale, and 
$\omega _0^{-1}$, where 
$\omega _0=\sqrt {g/H}$, as a time scale. Then, we can introduce the dimensionless time and coordinates
\begin{equation} \tau=t\omega_0,\quad x=\frac{X}{H}, \quad y=\frac{Y}{H}, \quad z=\frac{Z}{H}, \end{equation}the dimensionless variables
\begin{equation} \xi=\frac{w}{H},\quad \eta=\frac{u}{H},\quad \phi=\frac{\omega_0}{gH}\varphi, \end{equation}the dimensionless parameters of the added mass ratio (Minami Reference Minami1998) and membrane chord length
\begin{equation} \alpha=\frac{\rho H}{\rho_m h},\quad \Gamma=\frac{L}{H}, \end{equation}and the two Mach numbers (Vedeneev Reference Vedeneev2004, Reference Vedeneev2016)
\begin{equation} M_w=\frac{c}{\sqrt{g H}}, \quad M=\frac{v}{\sqrt{g H}}, \end{equation}
where 
$c^2=T/(\rho _m h)$ is the squared speed of propagation of elastic waves in the membrane and 
$\sqrt {g H}$ is the speed of propagation of long surface gravity waves in the shallow water approximation. The added mass ratio 
$\alpha$ is the ratio of the fluid to solid mass contained in the volume delimited by the dashed lines in figure 1 and in the membrane (Minami Reference Minami1998). In figure 1, 
$\Omega$ denotes the fluid domain and 
$\partial \Omega _0$, 
$\partial \Omega _1$ and 
$\partial \Omega _2$ stand, respectively, for the free surface, membrane and solid wall borders.
The dimensionless wave equation (2.1) is therefore
\begin{equation} \partial^2_{\tau}\xi-M_w^2\partial^2_{x}\xi=-\alpha \frac{{\rm \Delta} P}{\rho g H}, \quad x\in [0,\Gamma], \ z=0. \end{equation}Supplementing it with the expression (2.6), which in the dimensionless time and coordinates has the form
\begin{equation} \frac{{\rm \Delta} P}{\rho}=-\left(\omega_0\partial_{\tau}+\frac{v}{H}\partial_x\right)\varphi, \end{equation}we find
\begin{align} \partial^2_{\tau}\xi-M_w^2\partial^2_{x}\xi&= \alpha \left(\partial_{\tau}+\frac{v}{\omega_0 H}\partial_x \right)\frac{\omega_0}{g H}\varphi\nonumber\\ &= \alpha (\partial_{\tau}+M\partial_x )\phi. \end{align}The dimensionless boundary condition (2.9) is
\begin{equation} \boldsymbol{\nabla}\phi\boldsymbol{\cdot}\boldsymbol{n} = - (\partial_{\tau} + M\partial_x)\xi \quad \text{at } z=0, \ x\in [0,\Gamma], \end{equation}whereas the boundary condition (2.12) at the free surface in dimensionless form becomes
\begin{equation} \boldsymbol{\nabla}\phi\boldsymbol{\cdot}\boldsymbol{n} = - (\partial_{\tau} + M\partial_x)^2\phi \quad \text{at } z=1. \end{equation} Collecting together (2.19)–(2.21) and the obvious dimensionless versions of (2.5) and (2.8) and assuming a time dependence in the form of 
$\phi ,\xi \sim \textrm {e}^{-\textrm {i} \omega \tau }$ results in the following dimensionless set of equations and the boundary conditions for the case when a motionless medium is present below the membrane:
\begin{gather} \nabla^2\phi = 0,\quad \text{in } \Omega , \end{gather}
\begin{gather} \boldsymbol{\nabla}\phi\boldsymbol{\cdot}\boldsymbol{n} = -( -\textrm{i}\omega + M\partial_x )^2\phi, \quad \text{on } \partial\Omega_0, \end{gather}
\begin{gather} \boldsymbol{\nabla}\phi\boldsymbol{\cdot}\boldsymbol{n} = V(x), \quad \text{on } \partial\Omega_1, \end{gather}
\begin{gather} \boldsymbol{\nabla}\phi\boldsymbol{\cdot}\boldsymbol{n} = 0, \quad \text{on } \partial\Omega_2, \end{gather}
\begin{gather} {[}\omega^2 + M_w^2\partial_x^2]\xi = -\alpha(-\textrm{i}\omega + M\partial_x )\phi, \quad \text{on } \partial\Omega_1, \end{gather}
\begin{gather} \xi(0) = \xi(\Gamma) = 0, \quad \text{on } \partial\Omega_1, \end{gather}
where 
$V(x)=( \textrm {i}\omega - M\partial _x )\xi (x), x\in [0,\Gamma ]$ is the impermeability condition for the membrane. For simplicity, we retain the same notation for the membrane displacement and the fluid potential after the separation of time.
Therefore, due to the irrotational, incompressible and inviscid character of the fluid, our mathematical model (2.22) consists of the Laplace equation for the fluid potential (2.22a), supplemented by the kinematic conditions for the free surface (2.22b) and the membrane (2.22c). The pressure at the surface of the fluid is also prescribed as a dynamic condition and therefore closes the system of equations for the fluid in this model: the motion of the membrane is described by a non-homogeneous wave equation (2.22e) with the pressure of the fluid (recovered through the Bernoulli principle) as a source term. The membrane is supposed to be simply supported at its extremities as in (2.22f).
3. Methods and results
3.1. Membrane of infinite chord length
Our ultimate goal is to understand the fundamentals of the phenomenon of radiation-induced instabilities in the model (2.22) that we see as a reasonable analytically treatable substitute for the famous Lamb system (Lamb Reference Lamb1900; Barbone & Crighton Reference Barbone and Crighton1994; Hagerty et al. Reference Hagerty, Bloch and Weinstein2003). In this paper, as a first natural step, we analyse the case when the chord of the membrane is infinite, i.e. when the membrane extends from 
$-\infty$ to 
$+\infty$ in the 
$x$-direction.
The extension of the Nemtsov model to the case where the fluid layer presents a finite depth is our main concern. In the following we will show that even in the limit of infinite chord length the model (2.22) demonstrates physically meaningful radiation-induced flutter that sets in at finite values of the dimensionless flow velocity 
$M>M_w>0$, no matter what the values of the wavenumber 
$\kappa$ and the added mass ratio 
$\alpha$ are, in contrast to other known models discussed, e.g. in Miles (Reference Miles1956) and Dowell (Reference Dowell1966).
3.1.1. Dispersion relation for the fluid layer of arbitrary depth
Since the motion of the fluid is two-dimensional in the 
$(x,z)$-plane and the horizontal extension of the fluid layer is infinite in the 
$x$-direction too, we can represent the potential of the fluid 
$\phi$ in the physical space by means of the inverse Fourier transform of the potential 
$\hat \phi$ in the wavenumber space as
\begin{equation} \phi(x,z,\omega) = \frac{1}{2{\rm \pi}}\int_{-\infty}^{+\infty} \hat{\phi}(\kappa,z,\omega) \,\textrm{e}^{\textrm{i}\kappa x} \, \textrm{d} \kappa, \end{equation}
where 
$\kappa$ is the wavenumber and
\begin{equation} \hat\phi(\kappa,z,\omega) = \int_{-\infty}^{+\infty} \phi(x,z,\omega)\, \textrm{e}^{-\textrm{i}\kappa x} \, \textrm{d} x, \end{equation}
under the standard assumption that both 
$\phi (x)$ and 
$\hat \phi (\kappa )$ are absolutely integrable functions, implying they vanish at infinity.
Assuming that 
$\partial _x \phi$ is also absolutely integrable, which allows us to use twice the property 
$\widehat {\partial _x \phi }=\textrm {i}\kappa \hat \phi$, we find the Fourier transform of the Laplace equation (2.22a)
\begin{equation} \partial_z^2\hat{\phi} - \kappa^2\hat{\phi} = 0. \end{equation}The general solution to (3.3) is
\begin{equation} \hat{\phi}(\kappa,z,\omega) = A(\kappa,\omega)\, \textrm{e}^{\kappa z} + B(\kappa,\omega)\, \textrm{e}^{-\kappa z}, \end{equation}
where 
$A(\kappa ,\omega )$ and 
$B(\kappa ,\omega )$ are yet to be determined from the Fourier-transformed boundary conditions.
The boundary condition (2.22c), expressing the impermeability of the membrane at 
$z=0$, takes the form
\begin{equation} -\partial_z \phi=V, \end{equation}
because the outward direction of the normal vector 
$\boldsymbol {n}$ to the surface of the membrane is opposite to the positive 
$z$-direction, see figure 1. The Fourier transform of (3.5) reads
\begin{equation} \partial_z\hat{\phi} = -\hat{V}(\kappa,\omega), \end{equation}where
\begin{align} \hat{V}(\kappa,\omega) &= \int_{-\infty}^{+\infty} (\textrm{i}\omega\xi(s) - M\partial_{s}\xi(s))\, \textrm{e}^{-\textrm{i}\kappa s} \, \textrm{d} s\nonumber\\ &= \textrm{i}( \omega - \kappa M)\hat \xi. \end{align}
Substituting (3.4) into (3.6) yields at 
$z=0$
\begin{equation} \kappa ( A - B) = -\textrm{i}( \omega - \kappa M)\hat\xi. \end{equation}Similarly transforming the boundary condition (2.22b) at the free surface we find
\begin{equation} \partial_z\hat \phi = ( \omega -\kappa M )^2\hat\phi. \end{equation}
Substituting (3.4) into (3.9) yields at 
$z=1$
\begin{equation} \kappa ( A\textrm{e}^{\kappa} - B\textrm{e}^{-\kappa}) = ( \omega - \kappa M)^2 (A\textrm{e}^{\kappa} + B\textrm{e}^{-\kappa}). \end{equation}
Solving (3.8) and (3.10) simultaneously with respect to 
$A$ and 
$B$, we obtain
\begin{equation} \left.\begin{array}{c@{}} A(\kappa,\omega) = \dfrac{-\textrm{i}\hat \xi[(\omega - \kappa M)^2 + \kappa ](\omega - \kappa M)}{\kappa[(\omega - \kappa M)^2 - \kappa]\,\textrm{e}^{2\kappa} + \kappa[(\omega - \kappa M)^2 + \kappa]}, \\ B(\kappa,\omega) = \dfrac{\textrm{i}\hat \xi[(\omega - \kappa M)^2 - \kappa](\omega - \kappa M)}{\kappa[(\omega - \kappa M)^2 - \kappa] + \kappa [(\omega - \kappa M)^2 + \kappa]\,\textrm{e}^{-2\kappa}}. \end{array}\right\} \end{equation} The Fourier transform of the non-homogeneous wave equation (2.22e) for the membrane displacement evaluated at 
$z=0$ reads
\begin{equation} ( \omega^2 - \kappa^2 M_w^2 )\hat \xi - \textrm{i}\alpha(\omega - \kappa M)\hat \phi(\kappa,0,\omega) = 0. \end{equation} Inserting expression (3.4) for 
$\hat \phi$ with the coefficients (3.11) into (3.12), discarding 
$\hat \xi$ in the result and introducing new parameters, namely the phase velocity
\begin{equation} \sigma=\frac{\omega}{\kappa}, \end{equation}and the coupling parameter
\begin{equation} \beta=\frac{\alpha}{\kappa^2}, \end{equation}we obtain the following dispersion equation in the case where a medium with constant pressure is present below the membrane
\begin{equation} \beta = \frac{(M_w^2 - \sigma^2)[\kappa (\sigma -M)^2 - \tanh{\kappa}]}{\kappa (\sigma - M)^2[\kappa (\sigma - M)^2\tanh{\kappa} - 1]}. \end{equation} It is instructive to show another way of deriving the dispersion equation (3.15). For this, we notice that (3.6) and (3.7) allow us to express 
$\hat \xi$ by means of 
$\partial _z \hat \phi$. Using the result in (3.12), we can obtain a boundary condition for 
$\hat \phi (z)$ at 
$z=0$. This new boundary condition together with boundary condition (3.9) and equation (3.3) produce a closed-form boundary value problem for the Laplace equation with the Robin boundary conditions
\begin{equation} \left.\begin{array}{c@{}} \partial_z^2 \hat \phi-\kappa^2 \hat \phi = 0,\\ \partial_z \hat \phi ( \omega^2 - \kappa^2 M_w^2 ) - \alpha (\omega - \kappa M)^2\hat \phi = 0, \quad z=0,\\ \partial_z \hat \phi-(\omega - \kappa M)^2 \hat \phi= 0, \quad z=1. \end{array}\right\}\end{equation}
Substituting the general solution (3.4) into the boundary conditions of the problem (3.16) results in the system of two linear equations with respect to 
$A$ and 
$B$,
\begin{equation} \left.\begin{array}{c@{}} \kappa(A-B)( \omega^2 - \kappa^2 M_w^2 ) - \alpha (\omega - \kappa M)^2(A+B) = 0,\\ \kappa(A\textrm{e}^{\kappa} -B\textrm{e}^{-\kappa} )-(\omega - \kappa M)^2 (A\textrm{e}^{\kappa}+B\textrm{e}^{-\kappa}) = 0. \end{array}\right\} \end{equation}This system can be written in matrix form as
\begin{equation} (\omega^2 M_1 + \omega M_2 + M_3)\boldsymbol{f} = 0,\quad \boldsymbol{f}:=\left( \begin{array}{c} A \\ B\end{array} \right), \end{equation}
where the 
$2 \times 2$ matrices involved are
\begin{equation} \left.\begin{array}{c@{}} M_1 = -\left( \begin{matrix}\alpha-\kappa & \alpha+\kappa \\ \textrm{e}^{\kappa} & \textrm{e}^{-\kappa}\end{matrix} \right),\\ M_2= 2\kappa M\left( \begin{matrix} \alpha & \alpha \\ \textrm{e}^{\kappa} & \textrm{e}^{-\kappa} \\ \end{matrix} \right),\\ M_3= -\left( \begin{matrix} \kappa^2(M^2\alpha+M_w^2\kappa) & \kappa^2(M^2\alpha-M_w^2\kappa) \\ \kappa \textrm{e}^{\kappa}(M^2\kappa-1) & \kappa \textrm{e}^{-\kappa}(M^2\kappa+1) \\ \end{matrix} \right). \end{array}\right\} \end{equation}Computing the determinant of the matrix polynomial we arrive at the dispersion equation
\begin{align} D(\omega,\kappa)&= \det(\omega^2 M_1 + \omega M_2 + M_3)\nonumber\\ &= -\alpha(M\kappa-\omega)^2[(M\kappa-\omega)^2\tanh \kappa-\kappa] + \kappa (M_w^2\kappa^2-\omega^2) [(M \kappa-\omega)^2 - \kappa \tanh \kappa]\nonumber\\ &= 0, \end{align}
which, with the notation 
$\sigma =\omega /\kappa$ and 
$\beta =\alpha /\kappa ^2$, transforms exactly to (3.15).
For the sake of completeness we present also the dispersion relation for the system with a vacuum below the membrane
\begin{equation} \beta = \frac{(M_w^2 - \sigma^2)\left[\kappa (\sigma -M)^2 - \tanh{\kappa}\right]}{\left[\kappa^2 (\sigma - M)^4 - 1\right]\tanh{\kappa}}. \end{equation} In the shallow water approximation corresponding to the limit 
$\kappa \rightarrow 0$, the expression (3.21) reduces to
\begin{equation} \beta = (\sigma^2-M_w^2)\left( (\sigma -M)^2 - 1\right), \end{equation}which is nothing else but the shallow water dispersion relation derived by Nemtsov (Reference Nemtsov1985).
In order to get the dispersion relation (3.21), one must take the pressure difference (2.7), make it non-dimensional and use in the expression (2.17) which then reads as
\begin{equation} \partial^2_{\tau}\xi-M_w^2\partial^2_{x}\xi- \alpha \xi=\alpha (\partial_{\tau} + M\partial_x)\phi. \end{equation}After separation of time it reduces to the analogue of boundary condition (2.22e),
\begin{equation} [-\omega^2 - M_w^2\partial_x^2 - \alpha ]\xi - \alpha(-\textrm{i}\omega + M\partial_x)\phi(x,0,t) = 0, \end{equation}which has the following Fourier transform
\begin{equation} \left[\omega^2 -\kappa^2 M_w^2 + \alpha \right]\hat\xi -\textrm{i} \alpha\left(\omega -\kappa M\right)\hat\phi(\kappa,0,\omega) = 0. \end{equation}
Inserting the expression (3.4) for 
$\hat \phi$ with the coefficients (3.11) into (3.25) results, after familiar algebraic manipulations, in the dispersion relation (3.21).
3.1.2. Analysis of the dispersion equation
In the absence of coupling between the free surface and the membrane, i.e. for 
$\beta =0$, both the dispersion relation (3.15) and the dispersion relation (3.21) reduce to
\begin{equation} (\sigma^2-M_w^2)[\kappa (\sigma -M)^2 - \tanh{\kappa}]=0, \end{equation}
which yields the dispersion relation of the elastic waves in the free membrane 
$\sigma ^2=M_w^2$ and that of the surface gravity waves on a uniform flow: 
$\kappa (\sigma -M)^2 = \tanh {\kappa }$. The latter acquires a more familiar traditional form (Maissa et al. Reference Maissa, Rousseaux and Stepanyants2016)
\begin{equation} (\omega-\kappa Fr)^2=\kappa \tanh{\kappa} \end{equation}
after taking into account that 
$\sigma =\omega /\kappa$ and that 
$M$, as defined in (2.16a,b), can also be interpreted as the Froude number, 
$Fr$.
The roots of the decoupled dispersion equation (3.26) are real
\begin{equation} \sigma_1^{\pm}=\pm M_w, \quad \sigma_2^{\pm}=M\pm \sqrt{\frac{\tanh{\kappa}}{\kappa}}. \end{equation} If we consider the roots (3.28a,b) as functions of the fluid Mach number, 
$M$, we find that 
$\sigma _{1}^{\pm }$ are two horizontal straight lines and 
$\sigma _{2}^{\pm }$ are two straight lines with the slope equal to 
$1$, see figure 2(a). One can see that at 
$\beta =0$ the root branches intersect at four points forming the double roots 
$\sigma _0=M_w$ at
\begin{equation} M_0^{\pm}= M_w\pm\sqrt{\frac{\tanh \kappa}{\kappa}} \end{equation}
and the double roots 
$-\sigma _0$ at 
$-M_0^{\pm }$. The relation 
$M_w=M_0^+-\sqrt {(\tanh \kappa )/\kappa }=\sigma _2^-=\sigma _1^+=\sigma _0$ following from (3.29) and (3.28a,b) is the condition of ‘phase synchronism’ for the case of an arbitrary height of the fluid layer that extends the corresponding result obtained in Nemtsov (Reference Nemtsov1985) in the shallow water limit, 
$\kappa \rightarrow 0$.
Figure 2. Real (red, a–d) and imaginary (blue, e–h) parts of the roots of the dispersion relation (3.15) over the Mach number 
$M$ for 
$M_w=1$, 
$\kappa =1$ and (a,e) 
$\beta =0$, (b,f) 
$\beta =0.01$, (c,g) 
$\beta =0.1$ and (d,h) 
$\beta =1$.
With the increase in 
$\beta$ the roots 
$\pm \sigma _0$ situated at 
$M=\pm M_0^-$ split into simple real ones and this split is accompanied by unfolding of the crossings into avoided crossings, figure2.
Quite in contrast, the roots 
$\pm \sigma _0$ situated at 
$M=\pm M_0^+$ split into complex-conjugate pairs that form bubbles of instability at moderate values of 
$\beta$ that open up with the increase in 
$\beta$ to develop disconnected complex branches, as is seen in figure 2(d,h).
Let us re-write the dispersion relation (3.15) as follows
\begin{align} D(\sigma,M,\beta)&:=\beta\kappa (\sigma - M)^2\left[\kappa (\sigma - M)^2\tanh{\kappa} - 1\right]\nonumber\\ &\quad- {(M_w^2 - \sigma^2)\left[\kappa (\sigma -M)^2 - \tanh{\kappa}\right]}=0. \end{align}Then, we can apply to it the perturbation theory derived in appendix A.
Consider the double root 
$\sigma _0$ at 
$M= M_0^+$ and 
$\beta =\beta _0=0$. Adapt the approximate equation (A 18) to our model
\begin{align} &{\rm \Delta}\sigma(\partial^2_{\sigma M}D {\rm \Delta} M+\partial^2_{\sigma \beta}D {\rm \Delta} \beta) + \frac{1}{2}[ \partial_M^2 D ({\rm \Delta} M)^2+2\partial^2_{M\beta} D {\rm \Delta} M {\rm \Delta} \beta+ \partial_{\beta}^2 D ({\rm \Delta} \beta)^2] \nonumber\\ &\quad+ \frac{1}{2}\partial_{\sigma}^2 D({\rm \Delta}\sigma)^2 + \partial_M D {\rm \Delta} M+\partial_{\beta} D {\rm \Delta} \beta = 0, \end{align}
where 
${\rm \Delta} \sigma =\sigma -\sigma _0$, 
${\rm \Delta} M=M-M_0^+$ and 
${\rm \Delta} \beta =\beta$. Calculating the partial derivatives at 
$\sigma =\sigma _0$, 
$M= M_0^+$, and 
$\beta =\beta _0=0$, we find
\begin{equation} \left.\begin{array}{c@{}} \partial_{\sigma}^2 D= -8M_w \kappa\sqrt{\dfrac{\tanh \kappa}{\kappa}}, \quad \partial_{\sigma M}^2 D= 4M_w\kappa\sqrt{\dfrac{\tanh \kappa}{\kappa}},\\ \partial_{M}^2 D= 0, \quad \partial_{M\beta}^2 D=-\partial_{\sigma\beta}^2 D= 2\kappa(2(\tanh{\kappa})^2-1)\sqrt{\dfrac{\tanh \kappa}{\kappa}}, \quad \partial_{\beta}^2 D= 0, \\ \partial_{M} D= 0, \quad \partial_{\beta} D= (\tanh \kappa)^3-\tanh \kappa. \end{array}\right\} \end{equation}With the derivatives (3.32) the approximation (3.31) to the dispersion equation (3.30) near the crossing takes the form
\begin{equation} (\sigma-M_w)\left[\sigma-M+\sqrt{\frac{\tanh \kappa}{\kappa}}\right]= \beta\sqrt{\frac{\tanh \kappa}{\kappa}}\frac{(\tanh \kappa)^2-1}{4M_w}. \end{equation}
For any 
$\beta >0$ the crossing of the real roots 
$\sigma$ at 
$M= M_0^+$ unfolds into two hyperbolic branches of the real roots
\begin{align} &\beta\sqrt{\frac{\tanh \kappa}{\kappa}}\frac{[1-(\tanh \kappa)^2]}{4M_w} = \frac{1}{4}\left(M-M_w-\sqrt{\frac{\tanh \kappa}{\kappa}}\right)^2 \nonumber\\ &\quad - \left( \textrm{Re}\, \sigma-\frac{M_w+M}{2}+\frac{1}{2}\sqrt{\frac{\tanh \kappa}{\kappa}}\right)^2, \quad \textrm{Im} \,\sigma = 0 , \end{align}
that are connected to the ‘bubble’ of complex eigenvalues with the real parts 
$\textrm {Re}\, \sigma = \frac {1}{2}(M+M_w-\sqrt {(\tanh {\kappa })/\kappa })$ and with the imaginary parts that form an ellipse in the 
$(M, \textrm {Im}\, \sigma )$-plane
\begin{equation} (\textrm{Im}\, \sigma)^2+\frac{1}{4}\left(M-M_w-\sqrt{\frac{\tanh \kappa}{\kappa}}\right)^2=\beta\sqrt{\frac{\tanh \kappa}{\kappa}}\frac{[1-(\tanh \kappa)^2]}{4M_w }, \end{equation}
see figure 3. Equating to zero the discriminant of the quadratic in 
$\sigma$ in equation (3.33), we arrive at the following quadratic approximation to the neutral stability curve at the crossing point 
$M= M_0^+$
\begin{equation} \beta = M_w\frac{\left(M-M_w-\sqrt{(\tanh \kappa)/\kappa}\right)^2}{(1-(\tanh \kappa)^2)\sqrt{(\tanh \kappa)/\kappa}}. \end{equation}The bubble of instability (3.35) corresponds to the inner points of the instability domain bounded by (3.36).
Figure 3. Real and imaginary parts of the roots of the dispersion relation for 
$M_w=1$, 
$\kappa =1$ and 
$\beta =0.1$: (red) (3.15) and (blue, dashed) their approximations by (3.33) and (3.37) near the crossing points that exist at 
$\beta =0$, 
$M=M_0^{\pm }$, 
$\sigma =\sigma _0$. Notice an avoided crossing above the line 
$\textrm {Re}(\sigma )=M$ and the bubble of instability below this line.
Using the same methodology to approximate the avoided crossing close to 
$M=M_0^-$, 
$\sigma = \sigma _0$ and 
$\beta =\beta _0$ by (3.31), we obtain
\begin{equation} (\sigma-M_w)\left[\sigma-M-\sqrt{\frac{\tanh \kappa}{\kappa}}\right]=-\beta\sqrt{\frac{\tanh \kappa}{\kappa}}\frac{(\tanh \kappa)^2-1}{4M_w}. \end{equation}
Separating real and imaginary parts of 
$\sigma$ in (3.37) similarly to how it has been done in the previous case, one can see that the bubble of instability does not originate for 
$\beta >0$ in the unfolding of the crossing at 
$M=M_0^-$, see figure 3.
In figure 4 we show that the exact neutral stability boundaries obtained from equating the discriminant of the fourth-order polynomial (3.30) in 
$\sigma$ to zero and their approximation (3.36) calculated at the crossing point at 
$M=M_0^+$ are in a very good agreement.
Figure 4. Stability maps of the dispersion equation (3.30) given by its discriminant for (a) 
$M_w=1$ and 
$\kappa =1$, (b) 
$M_w=1$ and 
$\beta =0.1$, (c) 
$M_w=1$ and 
$M=1.6$ and (d) 
$\beta =0.5$ and 
$\kappa =1$. The regions of real phase speed 
$\sigma$ are shown in white (stability) and those of the complex 
$\sigma$ (temporal instability) in blue. The red dotted curve is the approximation (3.36). Notice the absence of instabilities for 
$M_w>M$ in panel (d).
It is instructive to change the point of view and to look at the critical values of parameters as functions of the Mach number 
$M_w$ of elastic waves in the membrane. In figure 5 we present stability maps of the dispersion equation (3.30) given by its discriminant in the 
$(M_w,\beta )$-plane for the fixed value of 
$M=M_0=2$ and increasing values of 
$\kappa$. We see that for all 
$\kappa$ the instability is possible only in the interval 
$|M_w|<M_0=2$, which agrees with figure 4. For 
$\beta =0$, the instability domain touches the 
$M_w$-axis at the points 
$M_w=M_0-\sqrt {(\tanh \kappa )/\kappa }$ and 
$M_w=-M_0+\sqrt {(\tanh \kappa )/\kappa }$. In the limit 
$\kappa \rightarrow 0$, the touching occurs at 
$M_w=M_0-1=1$ and 
$M_w=-M_0+1=-1$.
Figure 5. Stability maps of the dispersion equation (3.30) given by its discriminant for 
$M=M_0=2$ and: (a) 
$\kappa =0.5$, (b) 
$\kappa =0.55$, (c) 
$\kappa =0.58$, (d) 
$\kappa =0.8$, (e) 
$\kappa =1.5$, (f) 
$\kappa =3$. The regions of real phase speed 
$\sigma$ are shown in white (stability) and those of the complex 
$\sigma$ (temporal instability) in blue. The black dotted curve is the approximation (3.36) and the solid red ellipse is the conical approximation (3.43). When 
$\kappa \rightarrow \infty$, the central part of the instability domain approximated by (3.43) dominates over the side parts of the domain. Notice the absence of instabilities for 
$M_w>M_0$.
A qualitative change happens when 
$\kappa \ge \kappa _0$ where 
$\kappa _0>0$ is uniquely determined by 
$M_0>0$ from the equation
\begin{equation} \kappa_0 \tanh \kappa_0=\frac{1}{M_0^2}. \end{equation}
For instance, 
$M_0=2$ yields 
$\kappa _0\approx 0.5218134478$. At 
$\kappa =\kappa _0$ a new, isolated, domain of instability originates that touches the 
$M_w$-axis at 
$\beta =0$ and grows when 
$\kappa$ is further increased, figure 5. At some value of 
$\kappa$ the two domains touch each other and then form a unified domain. At 
$\kappa \rightarrow \infty$ the central part of the unified domain dominates over its side parts corresponding to the instability found by Nemtsov in the shallow water approximation when 
$\kappa \rightarrow 0$ and the coupling 
$\beta$ is weak, figure 5.
To understand the origin of the new instability, we plot the real and imaginary values of 
$\sigma$ as functions of 
$M_w$ in figure 6 for a given 
$M=M_0=2$. The central panel of figure6 corresponding to 
$\beta =0$ and 
$\kappa =\kappa _0$ shows four straight lines intersecting at five points, including the origin. The upper horizontal line corresponds to the fast surface gravity wave with 
$\sigma =M_0+\sqrt {(\tanh \kappa _0)/\kappa _0}\approx 3$, whereas the lower horizontal line to the slow surface gravity wave (Nemtsov Reference Nemtsov1985) with 
$\sigma =M_0-\sqrt {(\tanh \kappa _0)/\kappa _0}\approx 1$. The two inclined lines correspond to the forward and backward elastic waves in the membrane with 
$\sigma =\pm M_w$. When 
$\beta$ and 
$\kappa$ depart respectively from zero and 
$\kappa _0$, all the five crossings unfold either into avoided crossings (as elastic and fast surface gravity waves) or into bubbles of instability (as elastic and slow surface gravity waves) resulting in the high-frequency flutter due to radiation of long surface gravity waves. For 
$\beta >0$ the crossing at the origin transforms into an avoided crossing at 
$\kappa <\kappa _0$ or into the bubble of instability at 
$\kappa >\kappa _0$, which yields low-frequency flutter at short wavelengths 
$\kappa$.
Figure 6. Real (upper panels) and imaginary (lower panels) parts of the roots of the dispersion equation (3.30) for 
$M=M_0=2$ and: (a,d) 
$\beta =0.05$ and 
$\kappa =\kappa _0-0.1$, (b,e) 
$\beta =0$ and 
$\kappa =\kappa _0\approx 0.5218134478$, (c,f) 
$\beta =0.05$ and 
$\kappa =\kappa _0+0.3$. Notice that the bubbles of instability develop only for 
$\textrm {Re}(\sigma ) < M_0=2$.
Figure 7(a) allows us to track the evolution of the flutter domains as 
$\kappa$ varies from zero to infinity at 
$M_0=2$ and 
$\beta =0.03$. Nemtsov's radiation-induced flutter domain is the widest in the shallow water limit and evolves along the curves (shown as black solid lines in figure 7b)
\begin{equation} (M_0\pm M_w)^2=\frac{\tanh \kappa}{\kappa}, \end{equation}
to which the Nemtsov domains degenerate at 
$\beta =0$. Note that the Nemtsov flutter domain is perfectly approximated by formula (3.36) obtained from the unfolding of the eigenvalue crossing corresponding to the slow surface gravity wave and the elastic wave (dotted lines in figure 7a).
Figure 7. (a) Stability map of the dispersion equation (3.30) given by its discriminant for 
$M=M_0=2$ and 
$\beta =0.03$. The regions of real phase speed 
$\sigma$ are shown in white (stability) and those of the complex 
$\sigma$ (temporal instability) in blue. The black dotted curves correspond to the approximation (3.36) and the solid red line is the conical approximation (3.43). When 
$\beta =0$, the blue instability domains degenerate (central) to the ray 
$\kappa \ge \kappa _0\approx 0.5218134478$ and (sides) to the curves (3.39) shown as solid black lines. (b) Stability boundary with the conical singularity at 
$\kappa =\kappa _0$, 
$\beta =0$ and 
$M_w=0$, according to (blue, internal surface) the discriminant of the dispersion equation (3.30) and (red, external surface) to the approximation of the cone (3.43).
To understand the central instability domain shown in figure 7(a) for a given 
$\beta$ we plot it in the 
$(M_w,\kappa ,\beta )$-space in figure 7(b), given 
$M=M_0$. One can see that the domain is symmetric with respect to the plane 
$M_w=0$ and has a pronounced conical singularity at 
$\kappa =\kappa _0$ determined by (3.38) when 
$\beta =0$ and 
$M_w=0$. Equation (3.38) follows from the discriminant of the dispersion equation (3.30) at 
$\beta =0$ and 
$M_w=0$. The conical singularity of the stability boundary therefore exactly corresponds to the crossing of the eigenvalue curves at the origin in figure 6(b). Usually, the conical singularity of the stability boundary is associated with a double semi-simple eigenvalue with two linearly independent eigenvectors (Kirillov & Seyranian Reference Kirillov and Seyranian2004; Guenther & Kirillov Reference Guenther and Kirillov2006; Kirillov Reference Kirillov2009, Reference Kirillov2010, Reference Kirillov2013; Kirillov, Guenther & Stefani Reference Kirillov, Guenther and Stefani2009).
For this reason, we apply the perturbation theory of double eigenvalues presented in appendix A to the double zero eigenvalue 
$\sigma =\sigma _0=0$ at the crossing shown in figure6(b) and corresponding to the values of parameters 
$\beta =\beta _0=0$, 
$\kappa =\kappa _0$, 
$M=M_0$, 
$M_w=M_{w,0}=0$. A natural extension of the approximation formula (A 18) to the case of four parameters 
$\beta$, 
$\kappa$, 
$M_w$ and 
$M$ yields
\begin{align} &({\rm \Delta}\sigma)^2\frac{1}{2}\partial_{\sigma}^2 D+{\rm \Delta}\sigma\left(\partial^2_{\sigma \beta}D {\rm \Delta} \beta+\partial^2_{\sigma \kappa}D {\rm \Delta} \kappa+\partial^2_{\sigma M}D {\rm \Delta} M+\partial^2_{\sigma M_w}D {\rm \Delta} M_w\right) +\partial^2_{\beta \kappa} D {\rm \Delta} \beta {\rm \Delta} \kappa \nonumber\\ &\quad+\frac{1}{2}\left[ \partial_{\beta}^2 D ({\rm \Delta} \beta)^2+ \partial_{\kappa}^2 D ({\rm \Delta} \kappa)^2+ \partial_{M}^2 D ({\rm \Delta} M)^2+ \partial_{M_w}^2 D ({\rm \Delta} M_w)^2\right] +\partial^2_{M \kappa} D {\rm \Delta} M {\rm \Delta} \kappa \nonumber\\ &\quad+\partial^2_{M_w \kappa} D {\rm \Delta} M_w {\rm \Delta} \kappa+\partial^2_{\beta M} D {\rm \Delta} \beta {\rm \Delta} M+\partial^2_{\beta M_w} D {\rm \Delta} \beta {\rm \Delta} M_w+\partial^2_{\beta M} D {\rm \Delta} M_w {\rm \Delta} M \nonumber\\ &\quad+\partial_{\beta} D {\rm \Delta} \beta+\partial_{\kappa} D {\rm \Delta} {\kappa}+\partial_{M} D {\rm \Delta} {M}+\partial_{M_w} D {\rm \Delta} {M_w}=0. \end{align} Computing the corresponding partial derivatives of the left part of the dispersion equation (3.30), and evaluating them at 
$\beta =\beta _0=0$, 
$\kappa =\kappa _0$, 
$M=M_0$, 
$M_w=M_{w,0}=0$, where 
$M_0$ and 
$\kappa _0$ are related by (3.38), we find that the only non-zero derivatives are
\begin{equation} \left.\begin{array}{c@{}} \partial^2_{\sigma} D=-\partial^2_{M_w} D= 2 \kappa_0 M_0^2-\dfrac{2}{\kappa_0 M_0^2}, \quad \partial^2_{M \beta} D=-\partial^2_{\sigma \beta} D= 2 \kappa_0 M_0,\\ \partial^2_{\kappa \beta} D = M_0^4\kappa_0^2+M_0^2-1. \end{array}\right\} \end{equation}Taking this into account in (3.40), we find a simple approximation describing the unfolding of the double zero eigenvalue
\begin{align} &(M_0^4\kappa_0^2-1)(\sigma^2-M_w^2)-2\kappa_0^2M_0^3\sigma\beta+ \kappa_0M_0^2(M_0^4\kappa_0^2+M_0^2-1)(\kappa-\kappa_0)\beta\nonumber\\ &\quad +2\kappa_0^2M_0^3(M-M_0)\beta=0. \end{align} Let us further assume that 
$M=M_0$ is fixed. Then the last term in (3.42) vanishes, and the discriminant of the resulting quadratic polynomial in 
$\sigma$ produces the equation of a cone with the apex at 
$\kappa =\kappa _0$, 
$M_w=0$ and 
$\beta =0$
\begin{equation} M_0^6 \kappa_0^4 \beta^2-M_0^2 \kappa_0 (M_0^4\kappa_0^2+M_0^2-1)(M_0^4\kappa_0^2-1)(\kappa- \kappa_0)\beta+M_w^2(M_0^4\kappa_0^2-1)^2=0. \end{equation}
The cone (3.43) is shown in red in figure 7(b). With 
$\beta =0.03$, 
$M_0=2$ and 
$\kappa _0$ computed by means of (3.38), the approximation (3.43) fits the boundary of the exact instability domain with remarkable precision, as is evident in figure 7(a).
It is easy to see that in the plane 
$M_w=0$ the cone (3.43) defines the two lines
\begin{equation} \beta = \frac{(M_0^4\kappa_0^2-1)(M_0^4\kappa_0^2+M_0^2-1)}{M_0^4\kappa_0^3}(\kappa-\kappa_0),\quad \beta=0, \end{equation}
that approximate the instability domain near 
$\kappa =\kappa _0$, see figure 8(b). As soon as 
$M_w$ deviates from zero, the cone (3.43) again provides a very good fit to the actual stability boundary, figure 8(c). In the plane
\begin{equation} \beta = \frac{(M_0^4\kappa_0^2-1)(M_0^4\kappa_0^2+M_0^2-1)}{2M_0^4\kappa_0^3}(\kappa-\kappa_0), \end{equation}the cross-section of the cone (3.43) is described by the two lines
\begin{equation} \kappa = \kappa_0\pm M_w\frac{2\kappa_0M_0}{M_0^4\kappa_0^2+M_0^2-1}, \end{equation}that constitute a linear approximation to the stability boundary shown in figure 8(a).
Figure 8. For 
$M=M_0=2$ (a) cross-section of the instability domain with the conical singularity shown in figure 7(b) in the plane (3.45). The regions of real phase speed 
$\sigma$ are shown in white (stability) and those of the complex 
$\sigma$ (temporal instability) in blue. The red lines crossing at the apex of the cone at 
$\kappa = \kappa _0\approx 0.5218134478$ are linear approximations given by (3.46). (b) Cross-section in the plane 
$M_w=0$ of the instability domain and (red line) its linear approximation (3.44) at the conical point 
$\kappa =\kappa _0$. The black dotted line is given by (3.45). (c) Similar cross-section in the plane 
$M_w=0.1$ where the red curve is the approximation (3.43).
3.1.3. Wave energy of the Nemtsov system for membrane of infinite chord length
Let us use physical considerations to derive the expression for the averaged over the wave period energy of the Nemtsov system with the membrane of infinite chord length, by combining the approaches of the works by Maissa et al. (Reference Maissa, Rousseaux and Stepanyants2016) and Schulkes et al. (Reference Schulkes, Hosking and Sneyd1987).
In the linear wave theory, the energy is a function of the squared wave amplitude (Maissa et al. Reference Maissa, Rousseaux and Stepanyants2016). Therefore, the total energy per surface area of the membrane resulting both from the wave velocity of the structure and the elastic energy due to its tension is
\begin{equation} \mathcal{E}_m = \mathcal{K}_m + \mathcal{P}_m = \frac{1}{2} (\text{Re}[\partial_{\tau}\xi(x,\tau)])^2 + \frac{1}{2} M_w^2 (\text{Re}[\partial_x\xi(x,\tau)])^2, \end{equation}
where 
$\text {Re}$ stands for the real part of the vibration amplitude that is complex valued because of the assumed plane wave solution
\begin{equation} [\phi(x,z,\tau),\eta(x,\tau),\xi(x,\tau)] \sim [\hat{\phi}(z),\hat{\eta},\hat{\xi}] \textrm{e}^{\textrm{i}(\kappa x - \omega \tau)}. \end{equation}
Recall that 
$\hat {\phi }(z)$ is determined by the expression (3.4) with the coefficients (3.11) and 
$\hat {\eta }$, 
$\hat {\xi }$ are, respectively, displacement amplitudes of the free surface and the membrane.
The energy of the fluid depends on whether we assume a vacuum below the membrane (Nemtsov Reference Nemtsov1985) or a quiescent medium of the same density as the fluid above the membrane and with a pressure equal to the pressure of the unperturbed fluid (Vedeneev Reference Vedeneev2004, Reference Vedeneev2016). The gravitational potential energy of the free surface is the only term contributing to the total potential energy of the fluid in the latter context. Therefore,
\begin{equation} \mathcal{P}_f = \frac{1}{2} \alpha (\text{Re}[\eta(x,\tau)])^2. \end{equation} The kinetic energy of the flow per unit area is determined by the velocity field 
$\boldsymbol {u} = \boldsymbol {\nabla }\phi + M\boldsymbol {e}_x$, where 
$\boldsymbol {u}={\boldsymbol {v}}/{\sqrt {gH}}$, that needs to be directly integrated within the limits given by the surface of the membrane and the free surface of the fluid,
\begin{align} \mathcal{K}_f = \frac{1}{2} \alpha \int_{\text{Re}\, \xi}^{\text{Re}\, \eta} \vert\vert\text{Re}(\boldsymbol{u})\vert\vert^2 \, \textrm{d} z &= \frac{1}{2} \alpha \int_{\text{Re}\, \xi}^{\text{Re}\,\eta} \left[ (\text{Re}[\boldsymbol{\nabla}\phi])^2 + 2M\text{Re}[\partial_x\phi] + M^2 \right] \, \textrm{d} z \nonumber\\ &= \frac{1}{2} \alpha \int_{\text{Re}\, \xi}^{\text{Re}\, \eta} \left[ (\text{Re}[\partial_x\phi])^2 + (\text{Re}[\partial_z\phi])^2 \right] \, \textrm{d} z \nonumber\\ &\quad+ \alpha M \int_{\text{Re}\, \xi}^{\text{Re}\, \eta} \left[ \text{Re}[\partial_x\phi] + \frac{M}{2} \right] \, \textrm{d} z. \end{align} From assumption (3.48) and the explicit form of the complex amplitude 
$\hat {\phi }(z)$ determined by (3.4) with the coefficients (3.11), it follows that
\begin{equation} \left.\begin{array}{c@{}} \text{Re}[\partial_x\phi] = \textrm{i}\kappa\hat{\phi}(z)\cos{(\kappa x - \omega \tau)} , \quad \text{Re}[\partial_z\phi] = \textrm{i}\partial_z\hat{\phi}(z)\sin{(\kappa x - \omega \tau)}, \\ \text{Re}\, \xi = \hat{\xi}\cos{(\kappa x - \omega \tau)} , \quad \text{Re}\, \eta = \hat{\eta}\cos{(\kappa x - \omega \tau)}. \end{array}\right\} \end{equation}Taking into account the expressions (3.51) in (3.50), we find
\begin{equation} \int_{\text{Re}\, \xi}^{\text{Re}\, \eta} (\text{Re}[\partial_x\phi])^2 \, \textrm{d} z = -\kappa^2 \cos^2{(\kappa x - \omega \tau)} \int_{\text{Re}\, \xi}^{\text{Re}\, \eta} \hat{\phi}(z)^2 \, \textrm{d} z. \end{equation}Similarly, with the help of integration by parts, the Laplace equation (3.3), and expressions (3.51), we obtain
\begin{align} \int_{\text{Re} \, \xi}^{\text{Re}\, \eta} (\text{Re}[\partial_z\phi])^2 \, \textrm{d} z &= - \sin^2{(\kappa x - \omega \tau)} \left\{ [ \hat{\phi}\partial_z\hat{\phi} ]_{\text{Re}\, \xi}^{\text{Re}\, \eta} - \int_{\text{Re}\, \xi}^{\text{Re}\, \eta} \hat{\phi}(\partial_z^2\hat{\phi}) \, \textrm{d} z \right\} \nonumber\\ &= - \sin^2{(\kappa x - \omega \tau)} \left\{ [ \hat{\phi}\partial_z\hat{\phi} ]_{\text{Re}\, \xi}^{\text{Re}\, \eta} - \kappa^2 \int_{\text{Re}\,\xi}^{\text{Re}\,\eta} \hat{\phi}(z)^2 \, \textrm{d} z \right\} . \end{align} Finally, following Maissa et al. (Reference Maissa, Rousseaux and Stepanyants2016), we evaluate the last integral term in (3.50) with the help of the Lagrange mean value theorem, which is justified by the assumption that 
$\eta$ and 
$\xi$ are infinitesimally small perturbations of the surface boundaries 
$\partial \Omega _0$ and 
$\partial \Omega _1$. Performing this procedure, and then taking into account expressions (3.51), we obtain
\begin{align} &\int_{\text{Re}\, \xi}^{\text{Re}\, \eta} \left[ \text{Re}[\partial_x\phi] + \frac{M}{2} \right] \, \textrm{d} z \nonumber\\ &\quad= \int_{\text{Re}\, \xi}^{0} \text{Re}[\partial_x\phi] \, \textrm{d} z + \int_{0}^{1} \text{Re}[\partial_x\phi] \, \textrm{d} z + \int_{1}^{\text{Re}\, \eta} \text{Re}[\partial_x\phi] \, \textrm{d} z + \frac{1}{2} \int_{\text{Re}\,\xi}^{\text{Re}\, \eta} M \, \textrm{d} z \nonumber\\ &\quad= \text{Re}\, \eta\text{Re}[\partial_x\phi]|_{z=1} - \text{Re}\, \xi\text{Re}[\partial_x\phi]|_{z=0} + \int_{0}^{1} \text{Re}[\partial_x\phi] \, \textrm{d} z + \frac{M}{2} \text{Re}( \eta - \xi ) \nonumber\\ &\quad= \textrm{i}\kappa [ \hat{\eta}\hat{\phi}(1) - \hat{\xi}\hat{\phi}(0) ] \cos^2{(\kappa x - \omega \tau)} + \left[ \textrm{i} \kappa\int_0^1 \hat \phi(z)\, \textrm{d} z + \frac{M}{2} ( \hat{\eta} - \hat{\xi} ) \right] \cos{(\kappa x - \omega \tau)}. \end{align}
Note that the right-hand sides in the expressions (3.52)–(3.54) are 
$T$-periodic functions of time, where 
$T=2{\rm \pi} /\omega$. Averaging these expressions over the wave period 
$T$ according to the rule
\begin{equation} \langle\, f(\tau) \rangle = \frac{1}{T} \int_{0}^{T} f(\tau) \, \textrm{d} \tau, \end{equation}we deduce the mean kinetic energy of the fluid
\begin{equation} \langle \mathcal{K}_f \rangle = \frac{1}{4} \alpha \{ -[ \hat{\phi}\partial_z\hat{\phi} ]_{\partial \Omega_1}^{\partial \Omega_0} + 2\textrm{i}\kappa M [ \hat{\eta}\hat{\phi}(1) - \hat{\xi}\hat{\phi}(0) ] \}. \end{equation}
The term 
$\hat {\phi }\partial _z\hat {\phi }$ in (3.56) is evaluated with the help of the Bernoulli principle (2.11) and the free surface kinematic condition (2.10) at 
$\partial \Omega _0$, and the wave equation (2.22e) with the impermeability condition (2.22c) at 
$\partial \Omega _1$. This yields, respectively,
\begin{equation} \left.\begin{array}{c@{}} \hat{\phi}|_{\partial\Omega_0} = \dfrac{\hat{\eta}}{\textrm{i}(\omega - \kappa M)},\quad \hat{\phi}|_{\partial\Omega_1} = \dfrac{\omega^2 - \kappa^2 M_w^2}{\textrm{i}\alpha(\omega - \kappa M)}\hat{\xi},\\ \partial_z\hat{\phi}|_{\partial\Omega_0} = -\textrm{i}(\omega - \kappa M)\hat{\eta},\quad \partial_z\hat{\phi}|_{\partial\Omega_1} = -\textrm{i}(\omega - \kappa M)\hat{\xi}. \end{array}\right\} \end{equation}Substituting expressions (3.57) into (3.56) we obtain the final expression for the mean kinetic energy of the fluid
\begin{equation} \langle \mathcal{K}_f \rangle = \frac{1}{4} \{ \alpha\hat{\eta}^2 - (\omega^2 - \kappa^2 M_w^2) \hat{\xi}^2 + 2\textrm{i}\alpha\kappa M [ \hat{\eta}\hat{\phi}(1) - \hat{\xi}\hat{\phi}(0) ] \}. \end{equation}The other energies of the system, after taking into account (3.51) and time averaging (3.55), become
\begin{equation} \langle\mathcal{P}_m\rangle = \frac{1}{4}\kappa^2M_w^2\hat{\xi}^2, \quad \langle\mathcal{K}_m\rangle = \frac{1}{4}\omega^2\hat{\xi}^2, \quad \langle\mathcal{P}_f\rangle = \frac{1}{4}\alpha\hat{\eta}^2. \end{equation}
Notice that in the absence of the background flow 
$(M=0)$ the system respects the equipartition of energy
\begin{equation} \langle\mathcal{P}_f\rangle+\langle\mathcal{P}_m\rangle= \langle\mathcal{K}_f\rangle+\langle\mathcal{K}_m\rangle, \end{equation}in accordance with the virial theorem (Landau & Lifschitz Reference Landau and Lifschitz1987), because the flow is irrotational and thus derived from a fluid potential (Schulkes et al. Reference Schulkes, Hosking and Sneyd1987).
After summing up all the different terms given by (3.58) and (3.59a–c) we obtain the total averaged energy
\begin{equation} \langle \mathcal{E} \rangle = \frac{1}{2} \left\{ \kappa^2 M_w^2 \hat{\xi}^2 + \alpha \hat{\eta}^2 + \textrm{i}\alpha\kappa M \left[ \hat{\eta}\hat{\phi}(1) - \hat{\xi}\hat{\phi}(0) \right] \right\}, \end{equation}
thus providing an extension to the case when the velocity field contains a background flow (
$M\ne 0$).
A more suitable expression for the mean total energy can be obtained by expressing the different amplitudes of the system in (3.61) in terms of a unique one, for instance, 
$\hat {\xi }$. From the kinematic condition (2.10) on the free surface with the plane wave solution (3.48) and the coefficients (3.11), it is straightforward to express the surface amplitude 
$\hat {\eta }$ as
\begin{equation} \hat{\eta} = \frac{\textrm{i}\kappa}{\omega - \kappa M}\left[ A\textrm{e}^{\kappa} - B\textrm{e}^{-\kappa} \right] = \frac{(\omega - \kappa M)^2\hat{\xi}}{(\omega-\kappa M)^2\cosh{\kappa} - \kappa\sinh{\kappa}}. \end{equation} Substituting (3.62) into (3.61) and using the complex amplitude 
$\hat {\phi }(z)$ recovered from the boundary value problem (3.16), we find
\begin{align} \langle \mathcal{E} \rangle &= \frac{\hat{\xi}^2}{2}\left\{ \kappa^2M_w^2 + \alpha\frac{(\omega - \kappa M)^4(1-(\tanh{\kappa})^2)}{\left[(\omega-\kappa M)^2 - \kappa\tanh{\kappa}\right]^2} \right. \nonumber\\ &\quad\left. + \alpha M (\omega - \kappa M) \frac{\left[(\omega -\kappa M)^4 + \kappa^2\right]\tanh{\kappa} - 2\kappa(\omega -\kappa M)^2(\tanh{\kappa})^2}{\left[(\omega-\kappa M)^2 - \kappa\tanh{\kappa})\right]^2} \right\} . \end{align} Next, expressing the term 
$\kappa ^2M_w^2$ by means of the dispersion relation (3.20) and substituting the result into (3.63) yields a more compact formula for the total energy
\begin{equation} \langle \mathcal{E} \rangle =\frac{1}{4}\omega\left\{ 2\omega+\frac{\alpha}{\kappa}\frac{2(\omega-\kappa M)\tanh\kappa[(\omega-\kappa M)^4+\kappa^2-2\kappa(\omega-\kappa M)^2\tanh\kappa]}{[(\omega-\kappa M)^2-\kappa\tanh\kappa]^2}\right\}\hat{\xi}^2. \end{equation} Notice that the term in the braces in (3.64) is nothing else but the partial derivative 
${\partial \mathcal {D}}/{\partial \omega }$ of the dispersion relation (3.20) written in the following equivalent form
\begin{equation} \mathcal{D}(\omega,\kappa) := \mathcal{D}_m(\omega,\kappa) + \frac{\alpha}{\kappa} \frac{(\omega - \kappa M)^2[(\omega - \kappa M)^2\tanh{\kappa} - \kappa]}{\mathcal{D}_f(\omega,\kappa)}=0, \end{equation}
where 
$\mathcal {D}_m = \omega ^2-\kappa ^2M_w^2$ and 
$\mathcal {D}_f = [(\omega - \kappa M)^2 - \kappa \tanh {\kappa }]$ stand for the dispersion relation of, respectively, the free membrane and the free surface flow with a rigid boundary at the bottom. This proves that our total energy per unit area, averaged over the wave period, possesses the following simple representation in terms of the dispersion relation
\begin{equation} \langle \mathcal{E} \rangle = \frac{1}{4} \omega \frac{\partial\mathcal{D}}{\partial\omega} \hat{\xi}^2. \end{equation} The representation (3.66) can be found, e.g. in Cairns (Reference Cairns1979), and can be derived in the frame of the general Lagrangian variational approach (Ostrovskii, Rybak & Tsimring Reference Ostrovskii, Rybak and Tsimring1986; Whitham Reference Whitham1999), see also the recent work by Fukumoto, Hirota & Mie (Reference Fukumoto, Hirota, Mie, Kirillov and Pelinovsky2014) for historical notes and application to stability of vortices. Notice that according to (3.66) the energy vanishes at the points where 
$\omega =0$ or 
${\partial \mathcal {D}}/{\partial \omega }=0$, the latter condition corresponding to the existence of multiple roots of the dispersion relation. Correspondingly, the ratio 
$\langle \mathcal {E} \rangle / \omega$, which is the averaged wave action 
$\langle \mathcal {A} \rangle$ (Whitham Reference Whitham1999; Zhang et al. Reference Zhang, Qin, Davidson, Liu and Xiao2016), vanishes only at the locations of the multiple eigenvalues, cf. figures2 and 9. In the latter figure as well as in figure 10 we show several computations of the averaged wave energy and wave action over the fluid Mach number 
$M$, and, respectively, the membrane Mach number 
$M_w$, where 
$\omega$ is calculated with the use of the dispersion relation (3.20).
Figure 9. The averaged wave energy (a–d) 
$\langle \mathcal {E}\rangle$ given by the expression (3.64) and the action (e–h) 
$\langle \mathcal {A}\rangle = \langle \mathcal {E}\rangle / \omega$ over the Mach number 
$M$ evaluated for 
$M_w=1$, 
$\kappa =1$, 
$\hat {\xi }=0.01$ and: (a,e) 
$\alpha =0.1$, (b,f) 
$\alpha =0.5$, (c,g) 
$\alpha =1$ and (d,h) 
$\alpha =5$. Positive (respectively negative) energy/action is represented in red (respectively green).
Figure 10. The averaged wave energy (upper panels) 
$\langle \mathcal {E}\rangle$ given by the expression (3.64) and the action (lower panels) 
$\langle \mathcal {A}\rangle = \langle \mathcal {E}\rangle / \omega$ over the Mach number 
$M_w$ for 
$\hat \xi =0.01$, 
$M=M_0=2$, and (a,d) 
$\beta =0.05$ and 
$\kappa =\kappa _0-0.1$, (b,e) 
$\beta =10^{-3}$ and 
$\kappa =\kappa _0\approx 0.5218134478$, (c,f) 
$\beta =0.05$ and 
$\kappa =\kappa _0+0.3$.
4. Discussion
Comparing the eigenvalue plots of figures 2 and 6 with the averaged wave energy and wave action of each branch that are shown in figures 9 and 10, respectively, we notice that flutter instability is necessarily accompanied with the interaction of waves of opposite sign of energy/action. In contrast to the action, the energy changes sign also at the points where the phase velocity 
$\sigma$ changes sign, quite in accordance with (3.66).
Looking now at the roots (3.28a,b) of the decoupled dispersion equation (3.26), we conclude that the elastic waves 
$\sigma _1^{\pm }=\pm M_w$ propagating in the membrane always have positive energy whereas among the surface gravity waves 
$\sigma _2^{\pm }=M\pm \sqrt {(\tanh {\kappa })/{\kappa }}$ it is the energy of the slow wave 
$\sigma _2^-$ that becomes negative for 
$M>0$ as soon as 
$M>\sqrt {(\tanh {\kappa })/{\kappa }}$. Therefore at the crossing (3.29) corresponding to 
$M_0^{+} = M_w+\sqrt {(\tanh \kappa )/{\kappa }}$ the positive energy/action elastic wave meets the slow surface gravity wave that carries negative energy/action (Nemtsov Reference Nemtsov1985).
With 
$\beta$ increasing from zero, the crossing unfolds causing the eigenvalue branches to merge on the interval bounded by the points where 
$\partial _{\omega }D=0$. At these exceptional points (Kirillov Reference Kirillov2013) both the energy and the action change sign, see figures 9 and 10. On the interval the roots are complex and form the bubble of instability, see figures 2 and 6.
Since the fast surface gravity wave carries positive energy, one needs to add energy to the flow in order to excite this wave on the flow. In contrast, in order for the slow surface gravity wave carrying negative energy to build up on the flow, the energy must be extracted from the flow (Nezlin Reference Nezlin1976) via some mechanism for dissipation of its energy. In the Nemtsov problem, such a mechanism is the energy transfer from the slow surface gravity wave to an elastic wave associated with the membrane, which is a stationary medium and therefore has positive energy (Nezlin Reference Nezlin1976). One can say that this transferred energy yields flutter of the membrane due to emission of the slow surface gravity wave carrying negative energy.
In figures 2 and 6 as well as in figures 9 and 10 we observe that the flutter instability of the membrane occurs only if the velocity of the flow is higher than the phase velocity of the oscillations on the surface of the flow, 
$\sigma < M$, i.e. the flow moves faster than the waves it can excite (Nezlin Reference Nezlin1976; Nemtsov & Eidman Reference Nemtsov and Eidman1987). The condition 
$\omega =M\kappa$ or 
$\sigma =M$ is known as the Cerenkov condition for emission of radiation by a moving source (Ginzburg & Frank Reference Ginzburg and Frank1947; Nezlin Reference Nezlin1976; Ginzburg Reference Ginzburg1996; Bekenstein & Schiffer Reference Bekenstein and Schiffer1998; Carusotto & Rousseaux Reference Carusotto, Rousseaux, Faccio, Belgiorno, Cacciatori, Gorini, Liberati and Moschella2013). Substituted into a dispersion relation, the Cerenkov condition transforms the former into an expression defining a surface in the space of wavenumbers that determines the wake pattern behind the source (Schulkes et al. Reference Schulkes, Hosking and Sneyd1987; Carusotto & Rousseaux Reference Carusotto, Rousseaux, Faccio, Belgiorno, Cacciatori, Gorini, Liberati and Moschella2013). For the supercritical velocities 
$M>\sigma$ the surface in the space of wavenumbers develops a conical singularity known as the Cerenkov cone (Nemtsov Reference Nemtsov1985; Carusotto & Rousseaux Reference Carusotto, Rousseaux, Faccio, Belgiorno, Cacciatori, Gorini, Liberati and Moschella2013) with the angular aperture
\begin{equation} 2\arccos\left(\frac{\sigma}{M}\right). \end{equation}The anomalous Doppler effect (ADE) is the change in the sign of the field frequency radiated into the Cerenkov cone as compared with the field radiated outside this cone (Nezlin Reference Nezlin1976; Gaponov-Grekhov et al. Reference Gaponov-Grekhov, Dolina and Ostrovskii1983; Nemtsov Reference Nemtsov1985; Abramovich et al. Reference Abramovich, Mareev and Nemtsov1986; Carusotto & Rousseaux Reference Carusotto, Rousseaux, Faccio, Belgiorno, Cacciatori, Gorini, Liberati and Moschella2013). It is exactly the slow surface gravity wave that satisfies this condition
\begin{equation} \sigma_2^- - M= -\sqrt{\frac{\tanh \kappa}{\kappa}}<0. \end{equation}Hence, for the one-dimensional or, more precisely, plane two-dimensional waves, both the negative energy waves and the ADE correspond simply to waves with phase velocity lower than the flow velocity and wave vector pointing in the same direction as the flow (Nemtsov Reference Nemtsov1985; Ostrovskii et al. Reference Ostrovskii, Rybak and Tsimring1986). In our case, the radiated slow gravity wave increases the energy of oscillations of the membrane at the expense of the energy of the flow that supports this wave.
Finally, we plot the dispersion curves 
$\omega (\kappa )$ in figure 11, which show that the slow surface gravity wave branch and the membrane branch interact under the line 
$\textrm {Re}(\omega )=\kappa M$, if 
$\kappa >0$. Substituting the Cerenkov condition in the dispersion relation (3.20) we reduce it to 
$(M_w^2-M^2)\tanh \kappa =0$, thus providing a rationale for the absence of instabilities for 
$M_w^2>M^2$ that is evident in all our stability diagrams.
Figure 11. Dispersion curves ((red) real and (blue) imaginary parts of the roots 
$\omega$ of the dispersion relation (3.20)) for 
$M=2$, 
$\alpha \approx 0.0036725648$ and (a,d) 
$M_w=0$, (b,e) 
$M_w=0.09$, (c,f) 
$M_w=0.0967$, (g,j) 
$M_w=0.1$, (h,k) 
$M_w=0.5$, (i,l) 
$M_w=1$. Vertical dashed lines in the panels (a,d) correspond to 
$\kappa =\kappa _0\approx 0.5218134478$ and mark the onset of instability corresponding to the central instability zone in figure 7(a) and the conical instability zone in figure 7(b).
5. Conclusion
Through the revival of a classical hydrodynamical model performed in this work, we have been able to extend the stability analysis and to enhance knowledge of the underlying physics by making connections with the fundamental concepts such as negative energy waves and the ADE, supported by advanced mathematical tools.
Our continuation of Nemtsov's investigation of the radiation-induced flutter of a membrane in a uniform flow with the new derivation of the dispersion relation for a fluid layer of arbitrary depth and membrane of infinite chord length has led to a significant improvement in the computation of stability diagrams without any limitation on the range of the system parameters.
An exhaustive stability analysis has been performed using the original perturbation theory of multiple roots of the dispersion relation to compare with the exact stability domains, and both computations are proven to be in excellent agreement. More precisely, the crossings and avoided crossings of the dispersion curves are very well approximated by the simplified expressions for the phase speed of the membrane and fluid modes derived with the perturbation approach.
After computing the discriminant of the full dispersion relation, we have identified a new instability domain arising from a conical singularity in the parameter space that was not present in the early study of Nemtsov. This new domain is associated with a low-frequency flutter for short wavelengths and corresponds to the case when the velocity of propagation of elastic waves in the membrane is much smaller than the velocity of the flow.
Moreover, following the procedures used in previous studies on simplified hydrodynamical systems to calculate the averaged wave energy and after developing the method further to take into account the coupling between the free surface of the flow and the elastic membrane on the bottom, we have obtained an elegant and applicable expression for the total averaged energy. We have verified that, in the absence of the background flow, the system respects the equipartition of energy in accordance with the virial theorem, thus confirming that the existence of the negative energy waves can only be possible when the fluid is in motion.
We have shown that the formula for the total averaged energy recovered in our work by means of the direct integration of different physical fields is expressed via the derivative of the dispersion relation with respect to the frequency of oscillations and reduces exactly to the form described by Cairns (Reference Cairns1979).
The ADE is a direct consequence of the relative motion of an oscillator in a medium and more precisely, it occurs when the internal energy of the system increases due to the emission of negative energy waves (NEW). In our context, while the system is composed of a fluid layer and a membrane, such a phenomenon has been proved by Nemtsov to exist only when the conditions of phase synchronism and NEW emission are satisfied. The criterion for the phase synchronism in the system is easily identified in the computations of the dispersion curves as the crossings of the different branches that lead to the onset of positive growth rate and therefore to temporal instability. The latter phenomenon is a natural consequence of the highly excited state of energy that the system transits to due to the dominance of NEW over the waves carrying positive energy. Indeed, NEW emission is known as a process that increases the total energy of an oscillatory system while radiating energy away from the oscillator and, only when this gain in internal energy exceeds the losses from the contribution of positive energy modes, the total energy of the system starts growing in amplitude. Hence, it requires to have waves carrying energy of opposite signs that interact for the instability to develop.
Our expressions for the action and energy of the Nemtsov system demonstrated as expected the collision of waves carrying positive and negative energy as the onset for the radiative instability and the flutter of the membrane. Such a phenomenon is well known in the community of dynamical systems, but in this context, it is associated with the emission of NEW in the region of ADE. Hence, in addition to improving the stability analysis of the Nemtsov system and computing the averaged wave energy, our study provides a further, more detailed, examination of the ADE in hydrodynamics. Despite our problem being restricted to the study of planar waves, with the latter being emitted only in the horizontal direction, it is still sufficient for exploring the connection between the ADE and flutter theory.
An extension of this work to the case of a membrane with a chord of finite size, as described by the system of (2.22), is a promising necessary next step requiring asymptotic methods for the global stability analysis and numerical computations that we leave for future work.
Acknowledgements
We thank an anonymous referee for bringing important early works to our attention. We are grateful to Professor T.J. Bridges for helpful discussions. J.L. was supported by a PhD Scholarship from Northumbria University.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Sensitivity analysis of dispersion equations
In contrast to other works on frequency coalescence, e.g. Triantafyllou & Triantafyllou (Reference Triantafyllou and Triantafyllou1991), we adapt a more systematic multiparameter sensitivity analysis that can be found, e.g. in Kirillov & Seyranian (Reference Kirillov and Seyranian2002), Kirillov & Seyranian (Reference Kirillov and Seyranian2004), Kirillov (Reference Kirillov2007a), Kirillov (Reference Kirillov2007b), Kirillov (Reference Kirillov2009), Kirillov et al. (Reference Kirillov, Guenther and Stefani2009), Kirillov (Reference Kirillov2010) and Kirillov (Reference Kirillov2013).
Let us consider the dispersion equation
\begin{equation} D(\omega,p,q)=0, \end{equation}
where 
$D$ is a smooth function of scalar arguments 
$\omega$, 
$p$ and 
$q$. Assume that 
$D(\omega )$ is a polynomial of degree 
$n$ in 
$\omega$.
A.1. Sensitivity of simple roots
Let at 
$p=p_0$ and 
$q=q_0$ (A 1) have a simple root 
$\omega _0$ such that
\begin{equation} D_0:=D(\omega_0,p_0,q_0)=0, \end{equation}
where we use the symbol 
$:=$ to indicate a definition.
Following Kirillov & Seyranian (Reference Kirillov and Seyranian2002), Kirillov & Seyranian (Reference Kirillov and Seyranian2004), Kirillov (Reference Kirillov2007a), Kirillov (Reference Kirillov2007b), Kirillov (Reference Kirillov2010) and Kirillov (Reference Kirillov2013), we assume that 
$p=p(\varepsilon )$ and 
$q=q(\varepsilon )$. For 
$|\varepsilon |$ sufficiently small we can represent these functions as Taylor series
\begin{equation} \left.\begin{array}{c@{}} p(\varepsilon)= p_0+\varepsilon \dfrac{\textrm{d} p}{\textrm{d} \varepsilon}+\dfrac{\varepsilon^2}{2} \dfrac{\textrm{d}^2 p}{\textrm{d} \varepsilon^2}+o(\varepsilon^2),\\ q(\varepsilon)= q_0+\varepsilon \dfrac{\textrm{d} q}{\textrm{d} \varepsilon}+\dfrac{\varepsilon^2}{2} \dfrac{\textrm{d}^2 q}{\textrm{d} \varepsilon^2}+o(\varepsilon^2), \end{array}\right\} \end{equation}
with the derivatives evaluated at 
$\varepsilon =0$, and 
$p_0:=p(0)$ and 
$q_0:=q(0)$. Then, 
$\omega =\omega (\varepsilon )$ is also a root of (A 1), i.e. it satisfies the equation
\begin{equation} D_{\varepsilon}:=D(\omega(\varepsilon),p(\varepsilon),q(\varepsilon))=0. \end{equation}Differentiating (A 4), we find
\begin{equation} \frac{\textrm{d}}{\textrm{d} \varepsilon}D_{\varepsilon}=\partial_{\omega}D \frac{\textrm{d} \omega}{\textrm{d} \varepsilon}+\partial_{p}D \frac{\textrm{d} p}{\textrm{d} \varepsilon}+\partial_{q}D \frac{\textrm{d} q}{\textrm{d} \varepsilon}=0, \end{equation}
where the partial derivatives are evaluated at 
$\omega =\omega _0$, 
$q=q_0$, 
$p=p_0$.
Denoting 
${\rm \Delta} \omega =\varepsilon ({\textrm {d} \omega }/{\textrm {d} \varepsilon })\approx \omega -\omega _0$, 
${\rm \Delta} q=\varepsilon ({\textrm {d} q}/{\textrm {d} \varepsilon })\approx q-q_0$ and 
${\rm \Delta} p=\varepsilon ({\textrm {d} p}/{\textrm {d} \varepsilon })\approx p-p_0$, we find the expression for the increment of the simple root 
$\omega _0$ of (A 1) when the parameters depart from the values 
$q_0$ and 
$p_0$
\begin{equation} {\rm \Delta} \omega=- \frac{\partial_p D}{\partial_{\omega} D}{\rm \Delta} p- \frac{\partial_q D}{\partial_{\omega} D}{\rm \Delta} q+ o(|{\rm \Delta} p|, |{\rm \Delta} q|). \end{equation}A.2. Double root of the dispersion relation: generic case
Let at 
$p=p_0$ and 
$q=q_0$ the dispersion equation (A 1) have a double root 
$\omega _0$, which implies
\begin{equation} \left.\begin{array}{c@{}} D_0= 0,\\ \partial_{\omega}D_0:= \partial_{\omega}D(\omega_0,p_0,q_0)= 0. \end{array}\right\} \end{equation}
Assume that the perturbation of the parameters (A 3) causes splitting of the double root 
$\omega _0$ which generically is described by the Newton–Puiseux series (Kirillov & Seyranian Reference Kirillov and Seyranian2002, Reference Kirillov and Seyranian2004; Kirillov Reference Kirillov2007a,Reference Kirillovb, Reference Kirillov2010, Reference Kirillov2013)
\begin{equation} \omega(\varepsilon)=\omega_0+\omega_1\varepsilon^{1/2}+\omega_2\varepsilon+ \omega_3\varepsilon^{3/2}+\omega_4\varepsilon^{2}+o(\varepsilon^2). \end{equation} Expanding 
$D_{\varepsilon }$ as
\begin{equation} D_{\varepsilon}=\sum_{s=0}^n\frac{(\omega(\varepsilon)-\omega_0)^s}{s!} (\partial_{\omega}^s D+\varepsilon \partial_{\omega}^s D_1+\varepsilon^2\partial_{\omega}^s D_2 +o(\varepsilon^2) ), \end{equation}where
\begin{equation} \left.\begin{array}{c@{}} D_1:=\partial_pD \dfrac{\textrm{d} p}{\textrm{d} \varepsilon}+\partial_qD \dfrac{\textrm{d} q}{\textrm{d} \varepsilon},\\ D_2:=\dfrac{1}{2}\partial_p D \dfrac{\textrm{d}^2 p}{\textrm{d} \varepsilon^2}+\dfrac{1}{2}\partial_q D \dfrac{\textrm{d}^2 q}{\textrm{d} \varepsilon^2}+\dfrac{1}{2}\left( \partial_p^2 D \dfrac{\textrm{d}^2 p}{\textrm{d} \varepsilon^2}+2\partial^2_{pq} D \dfrac{\textrm{d} p}{\textrm{d} \varepsilon}\dfrac{\textrm{d} q}{\textrm{d} \varepsilon}+ \partial_q^2 D \dfrac{\textrm{d}^2 q}{\textrm{d} \varepsilon^2}\right), \end{array}\right\} \end{equation}
substituting expansion (A 8) into (A 9) and collecting the coefficients at the same powers of 
$\varepsilon$, we find
\begin{equation} \left.\begin{array}{c@{}} D_0= 0,\\ \omega_1 \partial_{\omega}D_0= 0,\\ D_1+\omega_1^2\dfrac{1}{2}\partial^2_{\omega}D+\omega_2 \partial_{\omega}D_0= 0. \end{array}\right\} \end{equation}
Looking for the coefficient 
$\omega _1\ne 0$, we see that the first two equations of (A 11) are satisfied in view of the fact that 
$\omega _0$ is a double root of the dispersion equation (A 1). Taking this into account, the last of (A 11) yields the expression for the coefficient 
$\omega _1$ in the expansion (A 8)
\begin{equation} \omega_1^2=-D_1 \left(\frac{1}{2}\partial_{\omega}^2 D \right)^{-1}, \end{equation}
where all the partial derivatives are calculated at 
$p = p_0$, 
$q = q_0$, 
$\omega = \omega _0$.
Therefore, if 
$D_1\ne 0$, the double root 
$\omega _0$ splits under variation of parameters (A 3) according to the formula
\begin{equation} \omega=\omega_0 \pm \sqrt{-\varepsilon D_1\left(\frac{1}{2}\partial_{\omega}^2D\right)^{-1}}+o(|\varepsilon|^{1/2}). \end{equation}In terms of the increments of the parameters, we can re-write (A 13) as
\begin{equation} {\rm \Delta} \omega= \pm \sqrt{-(\partial_pD {\rm \Delta} p+\partial_qD {\rm \Delta} q)\left(\frac{1}{2}\partial_{\omega}^2D\right)^{-1}}+ o(|{\rm \Delta} p|^{1/2}, |{\rm \Delta} q|^{1/2}). \end{equation}A.3. Double root of the dispersion relation: degenerate case
The case 
$D_1 = 0$ is degenerate, because the leading term in (A 13) of order 
$\varepsilon ^{1/2}$ vanishes and the expansion (A 8) is no longer valid, see e.g. Kirillov & Seyranian (Reference Kirillov and Seyranian2004). Substituting expansion (A 13) with 
$\omega _1=0$ into (A 9) and collecting coefficients of the same powers of 
$\varepsilon$, we obtain
\begin{equation} \left.\begin{array}{c@{}} D_1+\omega_2 \partial_{\omega} D_0= 0,\\ \omega_3 \partial_{\omega} D_0= 0,\\ D_2+\omega_2^2\dfrac{1}{2}\partial_{\omega}^2 D+\omega_2\partial_{\omega} D_1 +\omega_4\partial_{\omega} D_0= 0. \end{array}\right\} \end{equation}
Taking into account that 
$\partial _{\omega }D_0=0$ since 
$\omega _0$ is the double root and that 
$D_1 = 0$ due to our assumption, we conclude that the first two of (A 15) hold automatically. The third one simplifies as follows:
\begin{equation} \omega_2^2\frac{1}{2}\partial_{\omega}^2 D+\omega_2\partial_{\omega}D_1+D_2=0, \end{equation}
where all the derivatives are calculated at 
$\omega = \omega _0$, 
$p=p_0$ and 
$q = q_0$.
Therefore, the degeneracy, 
$D_1 = 0$, implies that the double root 
$\omega _0$ splits according to the formula
\begin{equation} \omega = \omega_0+\omega_2\varepsilon+o(\varepsilon), \end{equation}
where the coefficient 
$\omega _2$ is a root of the polynomial (A 16).
In combination with (A 10) and (A 17) the polynomial (A 16) transforms into
\begin{align} &({\rm \Delta} \omega)^2\frac{1}{2}\partial_{\omega}^2 D+{\rm \Delta} \omega(\partial^2_{\omega p}D {\rm \Delta} p+\partial^2_{\omega q}D {\rm \Delta} q)+\frac{1}{2}\left[ \partial_p^2 D ({\rm \Delta} p)^2+2\partial^2_{pq} D {\rm \Delta} p {\rm \Delta} q+ \partial_q^2 D ({\rm \Delta} q)^2\right] \nonumber\\ &\quad+ \partial_p D {\rm \Delta} p+\partial_q D {\rm \Delta} q=0. \end{align}Extension to the case of more than two parameters is straightforward, see e.g. Kirillov & Seyranian (Reference Kirillov and Seyranian2002), Kirillov & Seyranian (Reference Kirillov and Seyranian2004), Kirillov (Reference Kirillov2007a), Kirillov (Reference Kirillov2007b), Kirillov (Reference Kirillov2009), Kirillov et al. (Reference Kirillov, Guenther and Stefani2009), Kirillov (Reference Kirillov2010) and Kirillov (Reference Kirillov2013).
Open access
