Joris Labarbe
email: joris.labarbe@northumbria.ac.uk

Oleg N Kirillov
email: oleg.kirillov@northumbria.ac.uk

Radiation-induced instability of a finite-chord Nemtsov membrane

Keywords: radiation-induced instabilities, dissipation through dispersion, radiation damping, anomalous Doppler effect, flow-structure interaction, nonlinear eigenvalue problem 2010 MSC: 70K50, 76E17, 35P30, 65N30 Phillips scattering theory [3, 4

 has shown that an infinite-chord membrane is susceptible to flutter instability due to excitation of long gravity waves on the free surface if the velocity of the flow exceeds the phase velocity of the waves and placed this phenomenon into the general physical context of the anomalous Doppler effect. In the present work we derive a full nonlinear eigenvalue problem for an integro-differential equation in the case of the finite-chord Nemtsov membrane in the finite-depth flow. In the shallow-and deep water limits we develop a perturbation theory in the small added mass ratio parameter acting as an effective dissipation parameter in the system, to find explicit analytical expressions for the frequencies and the growth rates of the membrane modes coupled to the surface waves. This result reveals a new intricate pattern of instability pockets in the parameter space and allows for its analytical description. The case of an arbitrary depth flow with free surface requires numerical solution of a new non-polynomial nonlinear eigenvalue problem. We propose an original approach combining methods of complex analysis and residue calculus, Galerkin discretization, Newton method and parallelization techniques implemented in MATLAB to produce high-accuracy stability diagrams within an unprecedentedly wide range of system's parameters. We believe that the Nemtsov membrane appears to play the same paradigmatic role for understanding radiation-induced instabilities as the famous Lamb oscillator coupled to a string has played for understanding radiation damping.

Introduction

Exactly 120 years ago [START_REF] Lamb | On a peculiarity of the wave-system due to the free vibrations of a nucleus in an extended medium[END_REF] had proposed a model of a one-dimensional harmonic oscillator without damping constrained to move in the vertical direction and coupled to a horizontally taut semiinfinite elastic string [START_REF] Lamb | On a peculiarity of the wave-system due to the free vibrations of a nucleus in an extended medium[END_REF]. Quite surprisingly, he had found that the emission of traveling waves in the continuum by the oscillating mass contributes an effective Rayleigh damping correction term to the oscillator equation yielding decay of its vertical motion [START_REF] Lamb | On a peculiarity of the wave-system due to the free vibrations of a nucleus in an extended medium[END_REF]. In the course of time the radiating Lamb oscillator became paradigmatic for understanding the radiation damping in open and damped subsystems of closed conservative systems and gave rise to a number of abstract models of dispersion of energy from a 'small', usually finitedimensional, subsystem to a 'large', infinite-dimensional wave field [START_REF] Love | Some illustrations of modes of decay of vibratory motions[END_REF][START_REF] Beck | On the physical interpretation of complex poles of the S-Matrix-I[END_REF][START_REF] Ursell | The decay of the free motion of a floating body[END_REF][START_REF]Causality and Dispersion Relations[END_REF][START_REF] Crighton | The 1988 Rayleigh medal lecture: Fluid loading-the interaction between sound and vibration?[END_REF][START_REF] Barbone | Vibrational modes of submerged elastic bodies[END_REF][START_REF] Soffer | Resonances, radiation damping and instability in Hamiltonian nonlinear wave equations[END_REF][START_REF] Soffer | Dissipation through dispersion[END_REF][START_REF] Fröhlich | Emission of Cherenkov radiation as a mechanism for Hamiltonian friction[END_REF][START_REF] Wayne | Dynamics of partial differential equations[END_REF][START_REF] Olver | Dispersive Lamb systems[END_REF][START_REF] Bertini | Dynamics and Lax-Phillips scattering for generalized Lamb models[END_REF].

Remarkably, deep understanding of the radiation damping (including Lamb's model) involves Lax-is a universal mechanism describing even asymptotic stability of solitary type solutions, when radiation going away from the solitons to infinity leaves them to move freely [START_REF] Soffer | Resonances, radiation damping and instability in Hamiltonian nonlinear wave equations[END_REF][START_REF] Soffer | Dissipation through dispersion[END_REF][START_REF] Wayne | Dynamics of partial differential equations[END_REF].

To illustrate radiating vibratory motions in dimensions higher than one [START_REF] Love | Some illustrations of modes of decay of vibratory motions[END_REF] presented several extensions of Lamb's model that included decay of electromagnetic oscillations of a perfectly conducting spherical antenna due to emission of electromagnetic waves and decay of mechanical vibrations of an elastic sphere emitting acoustic waves [START_REF] Love | Some illustrations of modes of decay of vibratory motions[END_REF][START_REF] Beck | On the physical interpretation of complex poles of the S-Matrix-I[END_REF]. In [START_REF] Ursell | The decay of the free motion of a floating body[END_REF][START_REF] Barbone | Vibrational modes of submerged elastic bodies[END_REF] in the frame of the resonant scattering theory an important question of interaction between the vibrational modes of a submerged solid and the scattering functions in the fluid for both light and heavy fluid loading has been discussed. It was established that modes of the fluid-solid system at zero fluid-loading can be identified with the solid whereas at infinite fluid-loading they correspond to scattering frequencies of the fluid alone with a soft boundary condition. The intermediate values of the fluid-loading parameter appear to be the most complicated as it is impossible to label the mode of the fluid-solid system as corresponding solely to a 'solid mode' or a 'fluid mode' [START_REF] Barbone | Vibrational modes of submerged elastic bodies[END_REF].

In [START_REF] Hagerty | Radiation induced instability[END_REF][START_REF] Bloch | Gyroscopically stabilized oscillators and heat baths[END_REF] gyroscopic versions of the Lamb model were proposed, rather artificial however, such as the spherical pendulum and a rigid body with internal rotors, coupled either to the classical non-dispersive wave equation or to a dispersive equation of Klein-Gordon form. In such systems, the gyroscopic Lamb oscillator is susceptible to instabilities induced by wave emission (the radiation-induced instabilities), to which such physically important effects belong as the famous Chandrasekhar-Friedman-Schutz (CFS) instability of rotating stars caused by emission of gravitational waves [START_REF] Chandrasekhar | Solutions of two problems in the theory of gravitational radiation[END_REF][START_REF] Lindblom | On the secular instabilities of the Maclaurin spheroids[END_REF][START_REF] Schutz | Perturbations and stability of rotating stars -III. Perturbation theory for eigenvalues[END_REF][START_REF] Chandrasekhar | On stars, their evolution and their stability[END_REF][START_REF] Andersson | Gravitational waves from instabilities in relativistic stars[END_REF], acoustic version of CFS instability [START_REF] Belyaev | Incompressible modes excited by supersonic shear in boundary layers: Acoustic CFS instability[END_REF], and the instability of vortices in a stratified rotating fluid due to emission of internal gravity waves [START_REF] Le Dizès | Radiative instability in stratified vortices[END_REF] as it happens, e.g., in the events of coalescence of lenticular vortices observed in recent experiments [START_REF] Orozco Estrada | Coalescence of lenticular anticyclones in a linearly stratified rotating fluid[END_REF].

In a recent work [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF] attention was paid to an overlooked classical model that appears to be a perfect candidate for the role of the Lamb oscillator in the field of radiation-induced instabilities. This is the model, proposed by Nemtsov in 1985, of an elastic membrane resting on the bottom of a uniformly flowing fluid layer of finite depth and loosing its stability due to emission of surface gravity waves [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF]. Nemtsov's membrane having infinitely many modes of free vibrations plays the role of a 'small' subsystem, the fluid with the free surface is the 'large' continuum supporting propagation of dispersive surface gravity waves, and the motion of the flow can contribute to a gyroscopic coupling [START_REF] Paidoussis | Fluid-Structure Interactions: Slender Structures and Axial Flow: 1 2nd edition[END_REF].

We remark that scattering of surface gravity waves even by rigid horizontal submerged plates already has numerous applications in marine and coastal engineering such as submerged breakwaters or underwater wave lenses that allow exchange of water and hardly disturb horizontal currents [START_REF] Porter | Linearised water wave problems involving submerged horizontal plates[END_REF][START_REF] Islam | Scattering and radiation of water waves by a submerged rigid disc in a two-layer fluid[END_REF][START_REF] Williams | The Wiener-Hopf and residue calculus solutions for a submerged semi-infinite elastic plate[END_REF]. The need for light, inexpensive and rapidly deployable wave barriers requires taking into consideration submerged horizontal flexible plates and membranes [START_REF] Williams | The Wiener-Hopf and residue calculus solutions for a submerged semi-infinite elastic plate[END_REF][START_REF] Cho | Interactions of a horizontal flexible membrane with oblique incident waves[END_REF]. Recent applications in energy harvesting exploit fluidstructure interaction, leading to the excitation (flutter) of an elastic plate or membrane, usually referred to as a flag [START_REF] Alben | The flapping-flag instability as a nonlinear eigenvalue problem[END_REF][START_REF] Shelley | Flapping and bending bodies interacting with fluid flows[END_REF], due to radiation of surface gravity waves when immersed in a moving flow with a free surface [START_REF] Mougel | Flutter and resonances of a flag near a free surface[END_REF]. Another relevant setting comes from the problem of turbulent friction reduction in a boundary layer by using compliant coatings. In particular, it involves studying propagation of waves in a layer of a viscoelastic material of finite thickness when a layer of an ideal incompressible fluid is moving over it [START_REF] Vedeneev | Propagation of waves in a layer of a viscoelastic material underlying a layer of a moving fluid[END_REF].

The scattering theory formalism is efficient for analytical derivation of such important quantities as the reflection and transmission coefficients and displacement of the free surface of the flow [START_REF] Williams | The Wiener-Hopf and residue calculus solutions for a submerged semi-infinite elastic plate[END_REF]. However, investigation of stability of a radiating object requires different methods. In a related set of problems on the stability of oscillations of moving wave emitters (e.g. radiation of elastic waves in rails by high-speed trains [START_REF] Metrikin | Unstable lateral oscillations of an object moving uniformly along an elastic guide as a result of an anomalous Doppler effect[END_REF] and emission of internal or surface gravity waves by a spherical body on an elastic spring moving parallel to the interface of two liquids [START_REF] Gaponov-Grekhov | The anomalous Doppler effect and the radiation instability of oscillator motion in hydrodynamics[END_REF][START_REF] Abramovich | Instability in the oscillations of a moving oscillator while it radiates surface and internal waves[END_REF]), the theory of Cherenkov radiation for structureless particles [START_REF] Fröhlich | Emission of Cherenkov radiation as a mechanism for Hamiltonian friction[END_REF][START_REF] Ginzburg | Radiation by uniformly moving sources (Vavilov-Cherenkov effect, transition radiation, and other phenomena)[END_REF] and its extension by Ginzburg and Frank [START_REF] Ginzburg | About Doppler effect at superlight velocity[END_REF] to the particles having internal degrees of freedom, provides important clues both for derivation of necessary and sufficient criteria for instability and for better understanding radiation-induced instabilities in the general physical context [START_REF] Nezlin | Negative-energy waves and the anomalous Doppler effect[END_REF].

Originally, Cherenkov radiation has been the name for the phenomenon that a charged particle, moving relativistically through a dielectric non-dispersive medium at constant speed v higher than the phase velocity of light v p < c in the medium, becomes a source of electromagnetic radiation [START_REF] Fröhlich | Emission of Cherenkov radiation as a mechanism for Hamiltonian friction[END_REF][START_REF] Tamm | General characteristics of Vavilov-Cherenkov radiation[END_REF][START_REF] Bekenstein | The many faces of superradiance[END_REF][START_REF] Carusotto | The Cerenkov effect revisited: from swimming ducks to zero modes in gravitational analogues[END_REF][START_REF] Shi | Superlight inverse Doppler effect[END_REF]. If the source has a natural frequency ω 0 in its own static frame, then in the observer static frame one receives the far-field angular distribution of radiation with frequency ω = ω 0 γ -1 /(1 -vv -1 p cos θ), where γ is the Lorentz factor, that turns out to be concentrated in the forward direction on a conical surface making an angle θ with the velocity vector [START_REF] Carusotto | The Cerenkov effect revisited: from swimming ducks to zero modes in gravitational analogues[END_REF][START_REF] Shi | Superlight inverse Doppler effect[END_REF]. The angle θ ch = arccos(v p /v), which is possible only if ω 0 = 0, defines the Cherenkov cone with the angular aperture 2θ ch , i.e. a locus in the space of wavenumbers of resonant modes into which the Cherenkov radiation from the structureless particle occurs [START_REF] Tamm | General characteristics of Vavilov-Cherenkov radiation[END_REF][START_REF] Carusotto | The Cerenkov effect revisited: from swimming ducks to zero modes in gravitational analogues[END_REF][START_REF] Shi | Superlight inverse Doppler effect[END_REF]. In this case, the radiated electromagnetic field is spatially concentrated on a wave front forming a Mach cone with the angular aperture π -2θ ch behind the source, which is the apex of the cone [START_REF] Ginzburg | Radiation by uniformly moving sources (Vavilov-Cherenkov effect, transition radiation, and other phenomena)[END_REF][START_REF] Carusotto | The Cerenkov effect revisited: from swimming ducks to zero modes in gravitational analogues[END_REF].

A source oscillating with the natural frequency ω 0 > 0 and moving at a velocity v < v p or at a superluminal velocity v > v p but emitting outside the Cherenkov cone, i.e. under the condition cos θ < v p /v, experiences conventional Doppler effect with the increase in ω > 0 while approaching the observer [START_REF] Carusotto | The Cerenkov effect revisited: from swimming ducks to zero modes in gravitational analogues[END_REF][START_REF] Shi | Superlight inverse Doppler effect[END_REF]. If the superluminal oscillator emits inside the Cherenkov cone then the condition cos θ > v p /v implies ω < 0, which means that the source becomes excited by passing from a lower energy level to an upper one during the emission process [START_REF] Ginzburg | Radiation by uniformly moving sources (Vavilov-Cherenkov effect, transition radiation, and other phenomena)[END_REF][START_REF] Ginzburg | About Doppler effect at superlight velocity[END_REF]. That is, the kinetic energy of the source supplies both the energy of the emitted photon and the positive increase in the internal energy of the source [START_REF] Tamm | General characteristics of Vavilov-Cherenkov radiation[END_REF][START_REF] Shi | Superlight inverse Doppler effect[END_REF]. Radiation inside the Cherenkov cone is known as the anomalous Doppler effect [START_REF] Ginzburg | Radiation by uniformly moving sources (Vavilov-Cherenkov effect, transition radiation, and other phenomena)[END_REF][START_REF] Ginzburg | About Doppler effect at superlight velocity[END_REF][START_REF] Nezlin | Negative-energy waves and the anomalous Doppler effect[END_REF][START_REF] Bekenstein | The many faces of superradiance[END_REF]. Recent work [START_REF] Shi | Superlight inverse Doppler effect[END_REF] further distinguishes between the range v p /v < cos θ < 2v p /v and the range cos θ > 2v p /v inside the Cherenkov cone as resulting in the superlight (superluminal) normal and inverse Doppler shift, respectively. It turns out that the major part of the change in the kinetic energy of the source contributes to a positive increase in the internal energy of the source (photon) for the superlight inverse (normal) Doppler effect [START_REF] Shi | Superlight inverse Doppler effect[END_REF].

The same basic process of generalized Cherenkov emission is characteristic of all the emission processes that take place when a uniformly moving source is coupled to some excitation field, even in the presence of dispersion that provides an individual Cherenkov cone for each frequency [START_REF] Ginzburg | Radiation by uniformly moving sources (Vavilov-Cherenkov effect, transition radiation, and other phenomena)[END_REF]: as soon as the source velocity exceeds the phase velocity of some mode of the field, the latter becomes continuously excited [START_REF] Carusotto | The Cerenkov effect revisited: from swimming ducks to zero modes in gravitational analogues[END_REF]. For the sources having internal degrees of freedom, this serves also as a necessary condition for the presence in the space of parameters of a domain of instability of the source due to the anomalous Doppler effect [START_REF] Metrikin | Unstable lateral oscillations of an object moving uniformly along an elastic guide as a result of an anomalous Doppler effect[END_REF]. A sufficient condition for the radiation-induced instability is the prevalence of reaction of waves emitted inside the Cherenkov cone over those emitted outside it [START_REF] Metrikin | Unstable lateral oscillations of an object moving uniformly along an elastic guide as a result of an anomalous Doppler effect[END_REF]. Therefore, it is not surprising that already Tamm in his Nobel lecture foresaw application of the anomalous Doppler effect to describe "self-excitation of some particular modes of vibrations of a supersonic airplane" [START_REF] Ginzburg | Radiation by uniformly moving sources (Vavilov-Cherenkov effect, transition radiation, and other phenomena)[END_REF][START_REF] Tamm | General characteristics of Vavilov-Cherenkov radiation[END_REF].

In [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF] Nemtsov considered stability of the membrane under the surface of the uniform flow in the limits of (i) shallow water and (ii) vanishing added mass ratio [START_REF] Minami | Added mass of a membrane vibrating at finite amplitude[END_REF] that measures coupling between the membrane and the flow and serves as an effective damping parameter. Since the surface gravity waves are non-dispersive in the shallow water approximation, the formalism of Cherenkov radiation and anomalous Doppler effect applied to the infinite-chord-length Nemtsov membrane predicts its instability in the range where the velocity of the flow is exceeding the phase velocity of the surface gravity waves (as a consequence, they appear to be traveling backward in the frame moving with the flow) and exactly when the phase of the induced surface gravity wave is equal to the phase of the elastic wave propagating in the membrane [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF]. However, Nemtsov's shallow water result, being effective in uncovering fundamental physical reasons for the membrane destabilization, could not answer to a question of practical importance, namely, what is the domain of instability when the parameters of the system are allowed to take arbitrary values?

In [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF] we extended analysis of Nemtsov's membrane with infinite chord to the case of dispersive surface gravity waves by allowing the fluid layer to have arbitrary depth and the added mass ratio parameter to take arbitrary non-negative values. New complete dispersion relation has been derived and analyzed with the perturbation theory for multiple roots of polynomials [START_REF] Kirillov | Nonconservative Stability Problems of Modern Physics[END_REF] to obtain an explicit analytical approximation to the critical flutter velocity that is in excellent agreement with the numerical computation of the full stability map. Moreover, we have identified in [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF] a new instability domain arising from a conical singularity in the parameter space that could not be detected in the restrictive assumptions of [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF]. 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. Finally, an elegant and applicable explicit expression for the total averaged energy has been derived by means of the direct integration and its reduction to the Cairns form involving the derivatives of the dispersion relation with respect to the frequency of oscillations has been proven. It was demonstrated that the radiation-induced instability of the membrane is the result of collision of modes of positive and negative energy and can be interpreted in terms of wave emission in the domain of the anomalous Doppler effect [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF].

In the present work, we consider a model of Nemtsov's membrane in its entirety, with a flow of finite depth and a membrane of finite chord. We formulate the dimensionless boundary value problem for this system and by means of Fourier analysis recover an explicit expression for the velocity potential in the form of an improper integral. Determining the domain of dependence from the Cherenkov condition written for dispersive surface gravity waves, we extract an integro-differential equation for the membrane displacement in the presence of the flow, which is our main object of investigation.

Following [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF] we first present a rigorous general treatment of the integro-differential equation in the shallow-water limit by means of the Laplace transform and complex analysis. We obtain an integral eigenvalue relation and develop a systematic procedure for its analysis based on the perturbation theory with respect to the small added mass ratio parameter and residue calculus. As a result, we find explicit expressions for the frequencies and growth rates of membrane's modes coupled to the free surface as a series in the small parameter. Analysing the first-order approximation we can treat membrane's destabilization as a classical dissipation-induced instability [START_REF] Kirillov | Nonconservative Stability Problems of Modern Physics[END_REF][START_REF] Kirillov | Paradoxes of dissipation-induced destabilization or who opened Whitney's umbrella?[END_REF][START_REF] Mackay | Movement of eigenvalues of Hamiltonian equilibria under non-Hamiltonian perturbation[END_REF][START_REF] Haller | Gyroscopic stability and its loss in systems with two essential coordinates[END_REF][START_REF] Bloch | Dissipation induced instabilities[END_REF]. Plotting neutral stability curves in the plane of velocity of the flow versus the speed of propagation of elastic waves along the membrane we uncover a new and intricate pattern of self-intersecting instability pockets reminding a similar phenomenon characteristic of some forms of the Hill equation [START_REF] Broer | Geometrical aspects of stability theory for Hill's equations[END_REF] and describe it explicitly in an analytical form. We derive an integral expression for the free surface of the flow which allows us to find and explore the fluid dynamical analogue to the superlight normal and inverse Doppler effects [START_REF] Shi | Superlight inverse Doppler effect[END_REF].

The case of the finite depth of the fluid layer requires numerical solution of the boundary eigenvalue problem for the original integro-differential equation. This investigation involves a thorough treatment of the improper integral using complex analysis with the subsequent Galerkin decomposition of the solution to generate an algebraic nonlinear and non-polynomial in the eigenvalue parameter eigenvalue problem. The non-polynomial dependence on the eigenvalue parameter arises in the coupling term with the added mass ratio parameter as a factor. Setting the latter parameter to zero, we obtain a standard quadratic eigenvalue problem determining the free membrane modes. Once taken into consideration, the coupling term 'turns on' the radiative instability mechanism tending to excite the flutter of the membrane.

We notice that non-polynomial eigenvalue problems frequently occur in the studies of fluid-structure interactions, see e.g. [START_REF] Alben | The flapping-flag instability as a nonlinear eigenvalue problem[END_REF][START_REF] Vedeneev | Panel flutter at low supersonic speeds[END_REF], and are notoriously hard to solve even numerically. The methods for their solution are a hot topic in modern numerical mathematics and linear algebra communities, see e.g. [START_REF] Mehrmann | Nonlinear eigenvalue problems: A challenge for modern eigenvalue methods[END_REF][START_REF] Bindel | Localization theorems for nonlinear eigenvalue problems[END_REF][START_REF] Güttel | The nonlinear eigenvalue problem[END_REF][START_REF] Mensah | Perturbation theory of nonlinear, non-self-adjoint eigenvalue problems: Simple eigenvalues[END_REF], where a broad range of approaches is discussed. Most of the methods presented are either based on Newton-Raphson iterative process or on contour integration and we are restricting ourselves to the former.

In order to reach an acceptable convergence rate of the Newton method, we derive the Jacobian in an analytical form using residue calculus and complex analysis instead of approximating it numerically. Since the domain of integration in our problem is one-dimensional and we intend to keep high accuracy of our numerical scheme, we use the Legendre-Gauss-Lobatto quadrature rule to approximate the integrals in the Galerkin discretization. Nodes and weights of this spectral collocation method are recovered using the Golub-Welsch algorithm [START_REF] Golub | Calculation of Gauss quadrature rules[END_REF], which is based on the inversion of a linear system obtained from the threeterm recurrence relation for Legendre polynomials. The computed eigenvalues correspond to the quadrature points, while the eigenvectors are used to recover the weights. This spectral quadrature is able to reach computer accuracy with less than 20 nodes of discretization and is used all over our code that has been fully implemented and parallelized in MATLAB using the Parallel Computing Toolbox available from the software and run on the High Performance Cluster at Northumbria University.

To the best of our knowledge, the approach developed in our work is original and making use of it, we are able to recover the eigenfrequencies of the complete system and hence, to perform an exhaustive stability analysis of the finite-chord Nemtsov membrane. From the numerically found growth rates we recover stability maps for the finite-chord membrane in the finite depth layer and compare with the shallow water approximation. We show that in the limit of infinite chord length the neutral stability boundaries perfectly correspond to the shallow water boundaries found in [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF]. We establish that there is a critical chord length such that the shorter membranes cannot be destabilized. The most intriguing finding is, however, the chains of intertwining instability pockets, which our method is able to resolve, thus confirming its excellent convergence and accuracy. We believe that our procedure is applicable to a broad class of fluid-structure interaction problems that require solving nonlinear eigenvalue problems.

Mathematical formulation

Following [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF], in a Cartesian coordinate system OXY Z, we consider an inextensible elastic rectangular membrane strip of constant thickness h, density ρ m , and tension T along the membrane chord in the Xdirection. The membrane has infinite span in the Y -direction and is held at Z = 0 at the leading edge (X = 0) and at the trailing edge (X = L) by simple supports.

The membrane is initially still and flat, immersed in a layer of inviscid, incompressible fluid of constant density ρ, 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. Therefore, the system is solved under the potential theory as it is the case in the original paper by Nemtsov [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF]. The bottom of the fluid layer at Z = 0 is supposed to be rigid and flat for X ∈ (-∞, 0] ∪ [L, +∞).

In contrast to Nemtsov [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF] who assumed that vacuum exists below the membrane, we suppose that a motionless incompressible medium of the same density ρ is present below the membrane with a pressure that is the same as the unperturbed pressure of the fluid [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF]. The system is in a uniform gravity field acting in the negative Z-direction with g standing for the gravity acceleration. Let w(X, t), where t is time, be a small vertical displacement of the membrane, u(X, t) the free surface elevation, and ϕ(X, Z, t) the potential of the fluid.

Following [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF] we choose the height of the fluid layer, H, as a length scale, and ω -1 0 , where ω 0 = g/H, as a time scale to introduce the dimensionless time and coordinates

τ = tω 0 , x = X H , y = Y H , z = Z H , (1)
the dimensionless variables

ξ = w H , η = u H , φ = ω 0 gH ϕ, (2)
the dimensionless parameters of the added mass ratio [START_REF] Minami | Added mass of a membrane vibrating at finite amplitude[END_REF] and membrane chord length

α = ρH ρ m h , Γ = L H , (3)
and two dimensionless numbers

M w = c √ gH , M = v √ gH , (4)
where c 2 = T /(ρ m h) is the squared speed of propagation of elastic waves in the membrane and √ gH is the speed of propagation of long surface gravity waves in the shallow water approximation. The chosen scale for the velocity explains our choice of notation in (4) because the Froude number v √ gH can be treated as a Mach number in the non-dispersive shallow water limit which simplifies comparison of our results with that of the supersonic aerodynamics [START_REF] Vedeneev | Panel flutter at low supersonic speeds[END_REF].

Denoting the fluid domain by Ω, the free surface, membrane, and rigid wall borders by ∂Ω 0 , ∂Ω 1 , and ∂Ω 2 , respectively, and assuming the time dependence for the velocity potential φ and the membrane displacement ξ in the form of e -iωτ we arrive at the dimensionless boundary value problem [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF]

∇ 2 φ = 0, in Ω (5a) ∇φ • n + (-iω + M ∂ x ) 2 φ = 0, on ∂Ω 0 (5b) ∇φ • n + (-iω + M ∂ x )ξ = 0, on ∂Ω 1 (5c) ∇φ • n = 0, on ∂Ω 2 (5d) ω 2 ξ + M 2 w ∂ 2 x ξ + α (-iω + M ∂ x ) φ = 0, on ∂Ω 1 (5e) ξ(0) = ξ(Γ) = 0, on ∂Ω 1 (5f)
where n is the vector of the outward normal to a surface and, for simplicity, we retain the same notation for the membrane displacement and the fluid potential after the separation of time.

As one can notice, the Laplace equation (5a) is supplemented with a combination of dynamic and kinematic free surface conditions (5b) and the impermeability conditions for the membrane (5c) and the walls (5d). The nonhomogeneous wave equation (5e) describes the physics along the membrane and is solved with the rigid boundary conditions (5f). This set of equations represents a boundary eigenvalue problem for the complex eigenfrequency ω and will be used for the stability analysis of the finite-chord Nemtsov membrane.

Velocity potential via inverse Fourier transform

As in the previous study [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF], since the fluid layer is assumed to have an infinite extension in the xdirection, we can write, respectively, the Fourier transform of the velocity potential φ and its inverse

φ(κ, z, ω) = +∞ -∞ φ(x, z, ω)e -iκx dx, φ(x, z, ω) = 1 2π +∞ -∞ φ(κ, z, ω)e iκx dκ, (6)
where κ is the dimensionless wavenumber. With ( 6) taken into account, the Laplace problem (5a-5c) yields the boundary value problem for φ in the Fourier space

∂ 2 z φ -κ 2 φ = 0, in Ω (7a) ∂ z φ -(ω -κM ) 2 φ = 0, on ∂Ω 0 (7b) ∂ z φ + Γ 0 (iωξ(x ) -M ∂ x ξ(x )) e -iκx dx = 0, on ∂Ω 1 (7c)
where x is a curvilinear abscissa along the membrane. Expression (7c) is the Fourier transform of the impermeability condition (5c).

The general solution to equation ( 7a) is known to be

φ(κ, z, ω) = A(κ, ω)e κz + B(κ, ω)e -κz , (8)
where the functions A and B are to be determined from the boundary conditions (7b) and (7c). This yields an expression for φ at z = 0

φ(κ, 0, ω) = κ -(ω -κM ) 2 tanh κ κ[(ω -κM ) 2 -κ tanh κ] Γ 0 (-iω + M ∂ x ) ξ(x )e -iκx dx . (9)
Returning to the physical space by means of the inverse Fourier transform [START_REF] Crighton | The 1988 Rayleigh medal lecture: Fluid loading-the interaction between sound and vibration?[END_REF] we recover an explicit form for the potential disturbance

φ(x, 0, ω) = 1 2π Γ 0 (-iω + M ∂ x ) ξ(x ) +∞ -∞ κ -(ω -κM ) 2 tanh κ κ[(ω -κM ) 2 -κ tanh κ] e iκ(x-x ) dκ dx . (10)

Integro-differential equation for membrane's deflection

Surface gravity waves on a non-moving finite-depth layer are in general dispersive with the phase speed σ(κ) = √ κ tanh κ/κ, which is varying between σ DW = 0 when κ → +∞ (deep water) and σ SW = 1 when κ → 0 (shallow water) [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF]. Suppose we have a uniform flow in the positive x-direction with the supercritical speed M > σ(κ) as shown in Fig. 1. Perturbation with the wavenumber κ of the flow surface from a point (x, y, ξ) of the membrane will spread from that point along concentric circles in the (x, y) plane with the phase speed σ(κ), which are washed downstream as M > σ(κ) [START_REF] Drela | Flight Vehicle Aerodynamics[END_REF][START_REF] Schiffer | Analytical theory of subsonic and supersonic flows[END_REF]. The point obstacle therefore affects only the flow pattern in the conical domain of influence [START_REF] Olver | Introduction to Partial Differential Equations[END_REF]. The boundary of this downstream half of the Mach cone is the envelope of the moving and expanding circles centred at the points with the horizontal coordinates x ≥ x [START_REF] Drela | Flight Vehicle Aerodynamics[END_REF][START_REF] Schiffer | Analytical theory of subsonic and supersonic flows[END_REF]. The presence of the point obstacle at x does not make itself felt at the points upstream (x < x) and outside of the half-cone [START_REF] Drela | Flight Vehicle Aerodynamics[END_REF][START_REF] Schiffer | Analytical theory of subsonic and supersonic flows[END_REF].

x z M ∞ ∞ 0 Γ 1 (x, y, ξ) (x , y, ξ) (a) M M x y (x, y, z) 0 Γ ∞ ∞ (b)
Figure 1: (a) Sketch of the Nemtsov system in the (x, z)-plane. (b) View of the system from above in the (x, y)-plane. A point obstacle at a position (x, y, ξ) in a uniform flow moving in the positive x-direction with a velocity M > σ(κ) affects only the pattern of the surface gravity waves with the wavenumber κ in the (white) downstream Mach cone representing the domain of influence [START_REF] Drela | Flight Vehicle Aerodynamics[END_REF][START_REF] Schiffer | Analytical theory of subsonic and supersonic flows[END_REF][START_REF] Olver | Introduction to Partial Differential Equations[END_REF]. The (gray) upstream half of the Mach cone is the domain of dependence [START_REF] Drela | Flight Vehicle Aerodynamics[END_REF][START_REF] Schiffer | Analytical theory of subsonic and supersonic flows[END_REF][START_REF] Olver | Introduction to Partial Differential Equations[END_REF] at the point (x, y, z).

On the other hand, the flow at the point (x, y, z) depends only on the flow in the upstream half of the same Mach cone (domain of dependence [START_REF] Olver | Introduction to Partial Differential Equations[END_REF]), see Fig. 1(b) [START_REF] Drela | Flight Vehicle Aerodynamics[END_REF][START_REF] Schiffer | Analytical theory of subsonic and supersonic flows[END_REF]. Therefore, for M > σ(κ) the perturbed potential [START_REF] Fröhlich | Emission of Cherenkov radiation as a mechanism for Hamiltonian friction[END_REF] along the membrane with the finite chord length Γ at a given y is composed of all the single sources with coordinates (x , ξ) where x ∈ [0, x] as shown in Fig. 1(a). Hence, the interval of integration in the first integral in (10) can be truncated from x ∈ [0, Γ] to x ∈ [0, x] to match the domain of dependence. Inserting the modified in this manner expression [START_REF] Fröhlich | Emission of Cherenkov radiation as a mechanism for Hamiltonian friction[END_REF] into (5e) yields the following integro-differential equation for the membrane displacement ξ

ω 2 ξ + M 2 w ∂ 2 x ξ = α 2π (iω -M ∂ x ) x 0 (-iω + M ∂ x ) ξ(x ) +∞ -∞ κ -(ω -κM ) 2 tanh κ κ[(ω -κM ) 2 -κ tanh κ]
e iκ(x-x ) dκ dx .

(11) The equation [START_REF] Wayne | Dynamics of partial differential equations[END_REF] has to be solved with respect to the eigenfrequency ω as a nonlinear eigenvalue problem, which is generally a highly challenging task due to the presence of an improper integral. By this reason, in the next section we begin the analysis of (11) in the shallow water limit that will allow us to apply perturbation theory and derive explicit approximation of the flutter domain. With the guidance provided by the analytical solutions in the shallow-and deep-water limits we finally present a numerical method that eventually results in the solution of the full problem.

Shallow water analysis of the finite-chord Nemtsov membrane

The goal of this section is to extend the result of [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF] by presenting a rigorous and general treatment of expression [START_REF] Wayne | Dynamics of partial differential equations[END_REF] in the shallow water approximation and provide new physical interpretation of the instability mechanism for a membrane of the finite chord in the finite-depth flow.

Velocity potential in the shallow water limit

Introducing the phase speed σ = ω/κ and re-writing the factor at e iκ(x-x ) in the integrand of the improper integral in [START_REF] Wayne | Dynamics of partial differential equations[END_REF] as

1 κ 2 1 -κ(σ -M ) 2 tanh κ [(σ -M ) 2 -tanh κ κ ] = 1 κ 2 1 (σ -M ) 2 -1 + O(κ 2 ) = 1 (ω -κM ) 2 -κ 2 + O(1) (12)
we find that in the long-wavelength (shallow-water) limit, κ → 0, the velocity potential simplifies

φ SW (x) = 1 2π x 0 (-iω + M ∂ x ) ξ(x ) +∞ -∞ e iκ(x-x ) dκ (ω -κM ) 2 -κ 2 dx (13)
and, after factorizing the denominator in the integrand, it can be further expressed in the equivalent form

+∞ -∞ e iκ(x-x ) dκ (ω -κM ) 2 -κ 2 = 1 M 2 -1 +∞ -∞ e iκ(x-x ) (κ -p 1 )(κ -p 2 ) dκ, (14)
where the pole p 1 = ω/(M + 1) corresponds to the wave travelling forward along the membrane and p 2 = ω/(M -1) to the wave travelling backward.

Explicit form of the velocity potential by means of residue calculus

For M > 1 the denominator in the expressions for the poles p 1,2 (ω) always remains real and positive, only the (complex, in general) frequencies ω define the location of the poles in the complex κ-plane. In this configuration, we shall focus on the frequencies with positive imaginary parts (Im(ω) > 0) to define a contour in the upper-half plane and integrate the expression [START_REF] Adam | Rays, Waves, and Scattering: Topics in Classical Mathematical Physics[END_REF].

We define the contour of integration C as a semi-circular and positively oriented curve, of radius R, closed with a segment along the real axis as follows

C = [-R, R] ∪ ∆ R , ∆ R = {Re it , 0 ≤ t ≤ π}. ( 15
)
Applying the Cauchy residue theorem around the contour (15) yields [START_REF] Ablowitz | Complex Variables. Introduction and Applications[END_REF] C

F (z)dz = R -R F (z)dz + ∆ R F (z)dz = 2πi 2 j=1 res (F (z), p j ) , (16)
where

F (z) = f (z)e iz(x-x ) and f (z) = [(ω -zM ) 2 -z 2 ] -1 .
According to Jordan's lemma, since the function f (z) is continuous for any z ∈ C , except at the poles p 1,2 (ω), and that we have lim R→+∞ |f (Re iθ )| = 0 for θ ∈ [0, π], the integral over ∆ R vanishes as we enlarge the radius. Therefore, the improper integral [START_REF] Adam | Rays, Waves, and Scattering: Topics in Classical Mathematical Physics[END_REF] reduces to [START_REF] Ablowitz | Complex Variables. Introduction and Applications[END_REF]

lim R→+∞ R -R F (κ)dκ = +∞ -∞ F (κ)dκ = 2πi 2 j=1 res(F (κ), p j ). ( 17
)
In general, the residue res(g(z), p) of a meromorphic function g(z) having a simple pole p and a factorized denominator for this pole, can be found as res(g(z), p) = lim z→p (z -p)g(z). Hence, for p 1 = ω/(M + 1), p 2 = ω/(M -1) and F (κ) in the form of ( 14), we have [START_REF] Ablowitz | Complex Variables. Introduction and Applications[END_REF] res(F (κ), p 1 ) = lim

κ→p1 1 (M 2 -1) e iκ(x-x ) (κ -p 2 ) = 1 (M 2 -1) e ip1(x-x ) (p 1 -p 2 ) = -e ip1(x-x ) 2ω , res(F (κ), p 2 ) = lim κ→p2 1 (M 2 -1) e iκ(x-x ) (κ -p 1 ) = 1 (M 2 -1) e ip2(x-x ) (p 2 -p 1 ) = e ip2(x-x ) 2ω . (18)
With the residues (18) in the expression [START_REF] Bloch | Gyroscopically stabilized oscillators and heat baths[END_REF], the velocity potential ( 13) takes an explicit form

φ SW (x) = x 0 (ω + iM ∂ x ) ξ(x ) 2 j=1 res(F (κ), p j )dx , = 1 2ω x 0 (ω + iM ∂ x ) ξ(x ) e ip2(x-x ) -e ip1(x-x ) dx , = -i 2ω x 0 V (x ) e ip1(x-x ) -e ip2(x-x ) dx , (19)
which reproduces the result by Nemtsov, if we denote V (x) = (-iω + M ∂ x )ξ(x) [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF]. Note that the term e ip1(x-x ) corresponds to the normal Doppler effect due to emission of surface gravity waves of positive energy whereas the term e ip2(x-x ) to the anomalous Doppler effect due to emission of surface gravity waves of negative energy.

Explicit expression for the membrane displacement by means of Laplace transform

Substituting solution [START_REF] Lindblom | On the secular instabilities of the Maclaurin spheroids[END_REF] with u(x) = -iV (x) and v(x) = e ip2x -e ip1x into the problem (11) yields

ω 2 ξ(x) + M 2 w ∂ 2 x ξ(x) + α 2ω (-iω + M ∂ x ) x 0 u(x )v(x -x ) dx = 0, (20)
which can be written as

ω 2 ξ(x) + M 2 w ∂ 2 x ξ(x) + α 2ω (-iω + M ∂ x ) (u * v) (x) = 0, (21)
where the symbol * denotes the operator of convolution of two functions supported on the interval [0, ∞). Equation ( 21) supplemented with the boundary conditions (5e) for the displacement ξ(x) is suitable to solve by the Laplace method. We recall that in the general case, the Laplace transform L of an arbitrary function g(x) is given as follows

ḡ(s) ≡ L [g(x)] = +∞ 0 g(x)e -sx dx, s ∈ C. ( 22
)
Also recall that if h(x) is a convolution

h(x) = (u * v)(x) = x 0 u(x )v(x -x )dx , then L [h(x)] = ū(s)v(s). ( 23
)
Using Leibniz integral rule

d dx h(x) = d dx x 0 u(x )v(x -x )dx = v(0)u(x) + x 0 dv(x -x ) dx u(x )dx , (24)
we find the Laplace transform of the derivative of the convolution to be

L [dh/dx] = ū(s)v(0) + ū(s)(sv(s) -v(0)) = sū(s)v(s). (25)
Following the definition in [START_REF] Andersson | Gravitational waves from instabilities in relativistic stars[END_REF], applying the standard differentiation and integration properties of the Laplace transform to [START_REF] Chandrasekhar | On stars, their evolution and their stability[END_REF] and taking into account [START_REF] Belyaev | Incompressible modes excited by supersonic shear in boundary layers: Acoustic CFS instability[END_REF] and [START_REF] Orozco Estrada | Coalescence of lenticular anticyclones in a linearly stratified rotating fluid[END_REF], we find

ω 2 ξ(s) + M 2 w s 2 ξ(s) -sξ(0) -ξ (0) - iα 2ω ū(s)v(s)(ω + isM ) = 0, (26)
where

ū(s) = ω ξ(s) + iM s ξ(s) -ξ(0) , v(s) = 1 s -ip 2 - 1 s -ip 1 = i(p 2 -p 1 ) (s -ip 1 )(s -ip 2 ) = - 2iω (ω + isM ) 2 + s 2 . ( 27
)
Applying the boundary condition ξ(0) = 0 to the expressions ( 26) and ( 27) yields

ω 2 + s 2 M 2 w ξ(s) -α (ω + isM ) 2 (ω + isM ) 2 + s 2 ξ(s) = M 2 w ξ (0). ( 28
)
Since the equation ( 28) is linear in ξ(s), we can isolate this term and invert the whole expression by means of Mellin's inverse formula to finally obtain

ξ(x) ≡ L -1 ξ(p) = 1 2πi lim T →∞ T -iν -T -iν M 2 w ξ (0)[(ω -pM ) 2 -p 2 ]e ipx (ω 2 -p 2 M 2 w )[(ω -pM ) 2 -p 2 ] -α(ω -pM ) 2 dp, (29)
where p = -is and ν is a real number greater than the imaginary part of all the poles. The bounds in integral (29) define a line in the complex plane that is usually closed with a portion of a circle, thus delimiting a closed contour C B (commonly known as the Bromwich contour [START_REF] Ablowitz | Complex Variables. Introduction and Applications[END_REF]). Using the same argument of Jordan's lemma as for [START_REF] Hagerty | Radiation induced instability[END_REF], we prove that contribution of the circular integral is negligible in the limit of infinite radius. Then, another application of the Cauchy residue theorem to the contour integral allows an explicit computation of the inverse.

Integral eigenfrequency relation

Requiring ξ(x) in the expression (29) to vanish at x = Γ in accordance with (5e), we obtain the following eigenfrequency relation

ξ(Γ) = M 2 w ξ (0) 2πi C B [(ω -pM ) 2 -p 2 ]e ipΓ (ω 2 -p 2 M 2 w )[(ω -pM ) 2 -p 2 ] -α(ω -pM ) 2 dp = 0, (30)
which can be written as follows

ξ(Γ) = M 2 w ξ (0) 2πi D(ω, α) = 0, (31)
where

D(ω, α) = C B e ipΓ D(ω, α, p) dp, (32)
and

D(ω, α, p) = ω 2 -p 2 M 2 w -α (ω -pM ) 2 (ω -pM ) 2 -p 2 , ( 33
)
is nothing else but the shallow water dispersion relation of the membrane of the infinite chord length in the case of a medium with constant pressure and the same density as that of the fluid being present below the membrane [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF].

An expression similar to [START_REF] Williams | The Wiener-Hopf and residue calculus solutions for a submerged semi-infinite elastic plate[END_REF], derived earlier by Nemtsov for the membrane of the finite chord length with vacuum below the membrane [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF], naturally has the corresponding shallow water dispersion relation in the denominator of the integrand, which slightly differs from the ours.

Perturbation of eigenfrequencies

In the case of α = 0, the eigenvalue relation [START_REF] Williams | The Wiener-Hopf and residue calculus solutions for a submerged semi-infinite elastic plate[END_REF] reduces to

D(ω, 0) = C B e ipΓ D(ω, 0, p) dp = 1 M 2 w C B -e ipΓ p 2 -p 2 0 dp = 0, ( 34
)
where p 0 = ω/M w . Applying the residue theorem to the last integral in [START_REF] Shelley | Flapping and bending bodies interacting with fluid flows[END_REF], we find

D(ω, 0) = 1 ωM w sin Γω M w = 0,
which yields frequencies of the free (decoupled from the flow) membrane

ω n = πn M w Γ = 0, n ∈ N. (35)
In contrast to Nemtsov, we develop a systematic approach based on the perturbation theory of simple eigenvalues ω n to find how they are affected by a weak coupling to the flow with free surface [START_REF] Mackay | Movement of eigenvalues of Hamiltonian equilibria under non-Hamiltonian perturbation[END_REF][START_REF] Haller | Gyroscopic stability and its loss in systems with two essential coordinates[END_REF][START_REF] Bloch | Dissipation induced instabilities[END_REF]. Then, simple roots ω(α) of the equation D(ω, α) = 0 can be represented as a series in α, 0 < α 1, as follows [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF][START_REF] Kirillov | Nonconservative Stability Problems of Modern Physics[END_REF][START_REF] Mensah | Perturbation theory of nonlinear, non-self-adjoint eigenvalue problems: Simple eigenvalues[END_REF][START_REF] Greenbaum | First-order perturbation theory for eigenvalues and eigenvectors[END_REF]]

ω = ω n -α ∂ α D ∂ ω D + o(α), ω(0) = ω n . ( 36
)
Computing the partial derivatives and evaluating them at α = 0 yields

∂ α D = 1 M 2 w (M 2 -1) C B e ipΓ (ω n -pM ) 2 p 2 -p 2 0,n 2 (p -p 1,n )(p -p 2,n ) dp, ∂ ω D = - 2ω n M 4 w C B e ipΓ (p 2 -p 2 0,n ) 2 dp, (37)
where

p 0,n = πn Γ , p 1,n = πn Γ M w M + 1 , p 2,n = πn Γ M w M -1 .
Applying the residue theorem to the integrals in [START_REF] Metrikin | Unstable lateral oscillations of an object moving uniformly along an elastic guide as a result of an anomalous Doppler effect[END_REF], we find

∂ ω D = -iΓ 2 nπM 3 w (-1) n , ∂ α D = (-1) n iΓ 3 2π 2 n 2 M 4 w (M 2 w -M 2 ) 2 -M 2 w -M 2 (M 2 w -M 2 -1) 2 -4M 2 + iΓ 3 2π 3 n 3 M 3 w    (M -1) 2 sin πnMw M -1 (M 2 w -(M -1) 2 ) 2 - (M + 1) 2 sin πnMw M +1 (M 2 w -(M + 1) 2 ) 2    + -Γ 3 2π 3 n 3 M 3 w    (M -1) 2 (-1) n -cos πnMw M -1 (M 2 w -(M -1) 2 ) 2 - (M + 1) 2 (-1) n -cos πnMw M +1 (M 2 w -(M + 1) 2 ) 2    . (38)
With the derivatives [START_REF] Gaponov-Grekhov | The anomalous Doppler effect and the radiation instability of oscillator motion in hydrodynamics[END_REF] the series expansion [START_REF] Vedeneev | Propagation of waves in a layer of a viscoelastic material underlying a layer of a moving fluid[END_REF] takes the form

ω = ω n + α 2ω n (M 2 w -M 2 ) 2 -M 2 w -M 2 (M 2 w -M 2 -1) 2 -4M 2 + (-1) n αΓ 2π 2 n 2    (M -1) 2 sin πnMw M -1 (M 2 w -(M -1) 2 ) 2 - (M + 1) 2 sin πnMw M +1 (M 2 w -(M + 1) 2 ) 2    + iαΓ 2π 2 n 2    (M -1) 2 1 -(-1) n cos πnMw M -1 (M 2 w -(M -1) 2 ) 2 - (M +1) 2 1 -(-1) n cos πnMw M +1 (M 2 w -(M + 1) 2 ) 2    +o(α). ( 39
)
In particular, from [START_REF] Abramovich | Instability in the oscillations of a moving oscillator while it radiates surface and internal waves[END_REF] we easily obtain the growth rate of the perturbed simple real eigenvalue ω n

Im (ω) = αΓ 2π 2 n 2    (M -1) 2 1 -(-1) n cos πn Mw M -1 (M 2 w -(M -1) 2 ) 2 - (M + 1) 2 1 -(-1) n cos πn Mw M +1 (M 2 w -(M + 1) 2 ) 2    , (40)
which, if re-written as follows

Im (ω) = α 2ΓM 2 w p 2 2,n [1 -(-1) n cos (p 2,n Γ)] p 2 2,n -p 2 0,n 2 - p 2 1,n [1 -(-1) n cos (p 1,n Γ)] p 2 1,n -p 2 0,n 2 , (41)
exactly reproduces the growth rate derived earlier by Nemtsov [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF].

Stability diagrams in the shallow water limit

Setting to zero the linear in α approximation to the growth rate (41) of the n-th mode of the membrane, we can find an approximation to the neutral stability curve for this mode that subdivides the plane of parameters M and M w into the domains of stability and flutter instability, Fig. 2.

First of all we observe a cluster of instability domains grouped in the region M > M w in Fig. 2. The threshold M = M w is visible in Fig. 2(c,d) as a black dotted line. For the Nemtsov membrane of infinite chord length the physical meaning of this threshold is the equality of the velocity of the flow to the phase speed of elastic waves propagating along the membrane, which is a consequence of the Cherenkov condition [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF][START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF]. However, in contrast to the infinite-chord membrane's stability map reported in [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF], there are infinitely many 'petals' of flutter instability for M > M w touching each other when n > 1, see Fig. 2(c,d). Quite surprisingly, the common points of the petals all belong to straight lines of the following form

M w = 2j n + 1 (M + 1) , j ∈ Z, ( 42
)
where j is a negative integer in the region M > M w . For instance, n = 4 and j = -1 yield M = 2M w -1, which is a line passing through the point with M w = M = 1 in Fig. 2(c,d). The growth rate along this line presented in Fig. 3(b) demonstrates vanishing to zero exactly at the common points of the instability regions. These points are located exactly at the intersections of the straight lines [START_REF] Nezlin | Negative-energy waves and the anomalous Doppler effect[END_REF] with the curves

M w = k n M 2 -1 , k ∈ N (43)
that are shown as red and dashed in Fig. 2. Solving equations ( 42) and ( 43) we obtain the coordinates of the crossing points

M = 2j + k + n k , M w = (2j + n)(2j + 2k + n) kn . ( 44
)
For instance, for j = -1, n = 4 this yields M = 1 + 2/k and M w = 1 + 1/k and for k = 2 results in M = 2 and M w = 3/2, see Fig. 2(c,d) and Fig. 3(b). Notice absence of instability domains in a wide gap centered at M w = M + 1 and corresponding to j = 0 in [START_REF] Nezlin | Negative-energy waves and the anomalous Doppler effect[END_REF], which is clearly visible in Fig. 2(a,c). In the case of an infinite-chord membrane with vanishing coupling parameter (α = 0) the relation M w = M + 1 corresponds to a crossing of dispersion curves of surface gravity waves and elastic waves in the membrane that unfolds into an avoided crossing (stability) for α > 0 [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF]. We can conclude therefore that this very property of an infinite-chord membrane manifests itself as a stability gap at M w = M + 1 for the finite-chord membrane. In a strike contrast to the infinite-chord membrane, stability diagrams of Fig. 2 display a regular pattern of intertwined instability tongues, each centered along a line (42) with j ≥ 1. All the tongues commence at M = 1 at the values of M w that are specified by [START_REF] Nezlin | Negative-energy waves and the anomalous Doppler effect[END_REF]. For instance, if n = 1, then the tongues grow from M w = 6, 10, 14, . . . 2(2j + 1), . . ., see Fig. 2(a). The growth rate along the line [START_REF] Nezlin | Negative-energy waves and the anomalous Doppler effect[END_REF] corresponding to the tongue with n = 1 and j = 1 is shown in Fig. 3(a). The growth rate vanishes at M = 1 + 3/k, k ≥ 1, i.e. exactly at the crossing points (44) that subdivide the instability tongue into a collection of infinitely many instability pockets, see Fig. 2. Note that similar intertwined resonance tongues with instability pockets are known for the Hill equation with some specific forms of periodic excitation [START_REF] Broer | Geometrical aspects of stability theory for Hill's equations[END_REF].

In order to highlight the periodic pattern of the instability pockets, finally we plot the stability map in the new coordinates in Fig. 4 by projecting the growth rate onto the (µ, M )-plane, where µ = M w /(M -1). Then, the lines (42) transform into the curves µ = (2j/n + 1)(M + 1)/(M -1) and the curves [START_REF] Tamm | General characteristics of Vavilov-Cherenkov radiation[END_REF] to the lines µ = k/n(M + 1), see Fig. 4.

Exploring fluid dynamics analogy to superlight normal and inverse Doppler effects

According to a recent study [START_REF] Shi | Superlight inverse Doppler effect[END_REF], a source with an internal structure moving in a medium at a velocity that exceeds speed of light in the medium can be excited due to emission of electromagnetic waves (Ginzburg-Frank anomalous Doppler effect [START_REF] Ginzburg | About Doppler effect at superlight velocity[END_REF]) with the Doppler shift of the emitted waves remaining normal at small superluminal velocities and becoming inverted beyond some critical superluminal velocity. As it was emphasized in [START_REF] Shi | Superlight inverse Doppler effect[END_REF], virtually any wave system in nature, including classical wave systems such as acoustic waves and surface waves can exhibit the analogous phenomena.

Non-dispersive character of surface gravity waves in the shallow water limit implies M = 1 as a critical value for the flow velocity to exceed the speed of surface gravity waves and thus as a necessary condition for existence of the anomalous Doppler effect [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF][START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF]. Indeed, flutter instability tongues in Fig. 2 and Fig. 4 exist at M > 1. We need to verify that successive surface gravity waves radiated by the Nemtsov membrane carry wavelengths larger (smaller) than a characteristic value when M > M c > 1 (1 < M < M c ). Using the Bernoulli integral at the free surface (z = 1), retaining only linear in φ terms [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF], and assuming the time dependence exp(-iωτ ), we obtain

η(x) = -(-iω + M ∂ x )φ(x, z = 1), ( 45
)
where the velocity potential φ(x, z) is obtained as the inverse Fourier transform ( 6) of the solution of the boundary value problem [START_REF] Barbone | Vibrational modes of submerged elastic bodies[END_REF], then evaluated at z = 1. Indeed, from the general expression ( 8), we find

φ(κ, 1, ω) = 1 (ω -κM ) 2 cosh κ -κ sinh κ Γ 0 iω -M d dx ξ(x )e -iκx dx . (46)
Inserting solution [START_REF] Shi | Superlight inverse Doppler effect[END_REF] into expression ( 6) and considering the shallow water limit (κ → 0), we arrive at the fluid potential at the free surface (z = 1)

φ(x, z = 1) = - 1 2π Γ 0 -iω + M d dx ξ(x ) +∞ -∞ e iκ(x-x ) (ω -κM ) 2 -κ 2 dκdx , ( 47
)
where membrane's frequency ω and displacement ξ are provided by expressions ( 39) and ( 29), respectively. With the help of ( 47) equation ( 45) yields

η(x) = - i 2π Γ 0 -iω + M d dx ξ(x ) +∞ -∞ (ω -κM )e iκ(x-x ) (ω -κM ) 2 -κ 2 dκdx . ( 48
)
Following a procedure that we used previously to derive expression (13), we apply the Cauchy residue theorem to [START_REF] Kirillov | Nonconservative Stability Problems of Modern Physics[END_REF] with exactly the same poles as in ( 14) and recover the surface wave 'mode' (for a given ω) as

η(x) = - 1 2(M 2 -1) Γ 0 -iω + M d dx ξ(x ) (M + 1)e iκ -(x-x ) + (M -1)e iκ + (x-x ) dx . ( 49
)
Plotting [START_REF] Kirillov | Paradoxes of dissipation-induced destabilization or who opened Whitney's umbrella?[END_REF] in Fig. 5, we observe the different regimes 1 < M < M c , M = M c and M > M c for a certain 'mode' of free surface, corresponding to the wavelengths that are shorter, equal and longer than Γ, respectively. To estimate the critical value M = M c we use the idea of phase synchronisation between modes of the membrane and the surface of the flow [START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF][START_REF] Metrikin | Unstable lateral oscillations of an object moving uniformly along an elastic guide as a result of an anomalous Doppler effect[END_REF] that in the infinite-chord membrane case gives a sufficient condition for the presence of an instability domain in the parameter space [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF]. For the finite-chord membrane we consider the frequencies ω n of vibration modes of the free membrane defined by equation [START_REF] Mougel | Flutter and resonances of a flag near a free surface[END_REF] and surface waves frequencies ω f that for simplicity we take from the shallow water dispersion relation of the decoupled system

D SW (ω, κ) = (ω f -κM ) 2 -κ 2 = 0, ω ± f = κ(M ± 1). ( 50
)
The modes ω - f = κ(M -1) correspond to slow surface gravity waves carrying negative energy and thus exciting flutter of the infinite-chord membrane (anomalous Doppler effect) [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF][START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF]. Choosing κ = 2π/Γ in the equality ω n = ω - f allows us to estimate the critical velocity of the flow as

M c = 1 + nM w 2 , ( 51
)
which agrees with the threshold observed in Fig. 5 for the surface modes computed from expression (49).

Deep water limit of the finite-chord Nemtsov membrane

Let us find the leading term of the factor at e iκ(x-x ) in the integrand of the improper integral in (10) in the limit κ → +∞, corresponding to the deep water approximation

1 κ 2 1 -κ(σ -M ) 2 tanh κ [(σ -M ) 2 -tanh κ κ ] = - 1 κ + o(κ -1 ). ( 52
)
This yields a simplified form of the potential (10)

φ DW (x, ω) = - 1 2π x 0 (-iω + M ∂ x ) ξ(x ) +∞ -∞ e iκ(x-x ) dκ κ dx . (53)
For x > x , the Dirichlet integral in ( 53) can be calculated around the pole κ = 0 in the sense of the Cauchy Principal Value (CPV) as follows

PV +∞ -∞ e iκ(x-x ) κ dκ = iπ, (54)
where PV := lim →0 R\[-, ] and is the radius of the sphere enclosing the singularity [START_REF] Ablowitz | Complex Variables. Introduction and Applications[END_REF]. Inserting ( 54) into ( 53) we recover the velocity potential in the deep water approximation

φ DW (x, ω) = - 1 2 x 0 (ω + iM ∂ x ) ξ(x ) dx . ( 55
)
With the potential (55) the equation ( 11) takes the form

ω 2 ξ + M 2 w ∂ 2 ξ + i α 2 ω + iM d dx x 0 ω + iM d dx ξ(x )dx = 0. ( 56
)
Using the same methodology as before and taking into account the boundary condition ξ(0) = 0, we find the Laplace transform of ( 56)

ω 2 + s 2 M 2 w + iα 2s (ω + isM ) 2 ξ(s) = M 2 w ξ (0). ( 57
)
Following once again the procedure described in the previous section, we take s = ip in the equation ( 57) and then inverse the whole expression to finally obtain the displacement as the Bromwich integral

ξ(x) = M 2 w ξ (0) 2πi lim T →∞ T -iν -T -iν 2pe ipx 2p(ω 2 -p 2 M 2 w ) + α(ω -pM ) 2 dp. ( 58
)
Evaluating ( 58) at x = Γ and taking into account the boundary condition ξ(Γ) = 0, we recover the eigenfrequency equation

ξ(Γ) = M 2 w ξ (0) 2πi C B 2pe ipΓ 2p(ω 2 -p 2 M 2 w ) + α(ω -pM ) 2 dp = 0, (59)
which can be written as follows

ξ(Γ) = M 2 w ξ (0)D(ω, α) = 0, ( 60
)
where

D(ω, α) = 1 2πi C B e ipΓ D(ω, α, p) dp ( 61
)
and

D(ω, α, p) = ω 2 -p 2 M 2 w + α (ω -pM ) 2 2p . ( 62
)
In the case of α = 0 the eigenvalue relation ( 60) reduces to

D(ω, 0) = 1 2πi C B e ipΓ D(ω, 0, p) dp = 1 M 2 w 1 2πi C B -e ipΓ p 2 -p 2 0 dp = 0, ( 63
)
where p 0 = ω/M w . Applying the residue theorem to the last integral in [START_REF] Ablowitz | Complex Variables. Introduction and Applications[END_REF], we reproduce the eigenfrequencies ω n given by equation [START_REF] Mougel | Flutter and resonances of a flag near a free surface[END_REF].

Simple roots ω(α) of the equation D(ω, α) = 0 can be represented as a series in α, 0 < α 1, as follows [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF][START_REF] Kirillov | Nonconservative Stability Problems of Modern Physics[END_REF]

ω = ω n -α ∂ α D ∂ ω D - α 2 2 ∂ 2 ω D ∂ ω D ∂ α D ∂ ω D 2 -2 ∂ 2 ωα D ∂ ω D ∂ α D ∂ ω D + ∂ 2 α D ∂ ω D + o(α 2 ), ω(0) = ω n . (64)
Computing the partial derivatives and evaluating them at α = 0 yields

∂ α D = -1 2πi C B e ipΓ (M p -ω n ) 2 2M 4 w p(p 2 -p 2 0,n ) 2 dp, ∂ ω D = -1 2πi C B 2e ipΓ ω n M 4 w (p 2 -p 2 0,n ) 2 dp, ( 65
)
where p 0,n = πn Γ .

Applying the residue theorem to the integrals in [START_REF] Delves | A numerical method for locating the zeros of an analytic function[END_REF], we find

∂ α D = iM (-1) n Γ 2 2M 3 w πn + M w Γ 2 [(-1) n -1] 2M 3 w π 2 n 2 , ∂ ω D = -i(-1) n Γ 2 M 3 w πn . ( 66
)
With the derivatives (66) the series expansion (64) takes the form

ω(α) = ω n + α M 2 + iα M w 2πn [(-1) n -1] + o(α). ( 67
)
Taking into account terms of the second order in α yields the following expression for the growth rate For all even n > 1 the growth rate (68) is negative up to the terms of higher order than α 2 :

Im(ω) = α M w 2πn [(-1) n -1] + α 2 Γ 32 2M 2 w (4 -π 2 n 2 ) π 4 n 4 [(-1) n -1] - M 2 (π 2 n 2 + 10(-1) n -14) π 2 n 2 - M 4 M 2 w + o(α 2 ). ( 68
)
Im(ω) = -α 2 M 2 Γ 32 1 - 4 π 2 n 2 + M 2 M 2 w + o(α 2 ), even n > 1.
In contrast, for odd n the growth rate can take positive values for some combinations of parameters. Equating the growth rate (68) to zero, we find the critical length of the membrane at the onset of flutter for odd n in the explicit form

Γ = 32M 3 w π 3 n 3 α -1 4M 4 w (π 2 n 2 -4) -M 2 M 2 w π 2 n 2 (π 2 n 2 -24) -M 4 π 4 n 4 , odd n ≥ 1. (69)
For instance, the minimal length of the membrane for n = 1 below which there is no flutter, is

Γ min = 128π 3 αM w (π 4 -32π 2 + 512) (70)
at

M min = M w √ 2 2π 24 -π 2 , (71)
as is shown in Fig. 6. Notice that Γ min quickly increases as α is tending to zero, thus reducing chances for a finite-chord membrane to be unstable in the deep water limit, if the coupling between the membrane and the flow is vanishingly weak.

The finite-chord Nemtsov membrane in the finite-depth layer

The case of a fluid layer with finite depth requires further attention in the derivation of eigenfrequency equation than the previously considered limits of shallow and deep water. The main difficulty is the nonpolynomial character of the dispersion relation for the finite depth flow [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF], which does not allow analytical solution via Cauchy residue theorem of the improper integral in [START_REF] Wayne | Dynamics of partial differential equations[END_REF] that therefore has to be now exclusively treated numerically. This major difference forces us first to count and localize the poles of the integrand using a combination of complex analysis and iterative solvers to make possible numerical implementation of the residue theorem. This allows us to discretize the integro-differential equation into an algebraic nonlinear eigenvalue problem for ω by means of Galerkin decomposition and solve the resulting equation with an appropriate numerical method. In this section, we present both the derivation of this nonlinear eigenvalue problem and the methods for its solution. We argue that our approach can easily be extended to a wide class of fluid-structure systems composed of a fixed elastic structure that interacts with a moving flow. 76) over the radius R for arbitrary pair of complex eigenfrequencies and Mach numbers. We used N = 20 nodes for the trapezoidal quadrature.

Counting and localizing the poles of a non-polynomial integrand in the integro-differential equation

For the sake of clarity in this section, we reintroduce the wave equation of the finite depth problem (11) in terms of convolution, along with its boundary conditions

W (ξ; x, ω) = ω 2 + M 2 w ∂ 2 ∂x 2 ξ(x) + α -iω + M ∂ ∂x (u * v) (x) = 0, ξ(0) = ξ(Γ) = 0, (72)
where we have

u(x, ω) = -iω + M ∂ ∂x ξ(x), (73)
v(x, ω) = 1 2π +∞ -∞ [κ -(ω -κM ) 2 tanh κ]e iκx κ[(ω -κM ) 2 -κ tanh κ] dκ = 1 2π +∞ -∞ F (κ)dκ. ( 74
)
As stated before, the main issue in solving (72) is the presence of the improper integral (74). Indeed, in contrast with the shallow-and deep water cases, the difficulty in direct application of the Cauchy residue theorem to this integral is that it requires knowledge of the poles of the meromorphic function F (κ), or equivalently, the zeros of the analytic function

f (κ) = (ω -κM ) 2 -κ tanh κ, (75)
since the trivial pole κ = 0 of F (κ) is already known. As one can notice, expression (75) is nothing else but the dispersion relation of the surface gravity waves travelling along a fluid layer with finite depth and infinite extension, moving uniformly along a rigid bottom [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF].

Since the characteristic equation ( 75) is not polynomial in κ we cannot say a priori how many zeroes it has in the complex κ-plane. We determine this number numerically using the standard expression [START_REF] Ablowitz | Complex Variables. Introduction and Applications[END_REF][START_REF] Delves | A numerical method for locating the zeros of an analytic function[END_REF]

ν = 1 2πi C f (κ) f (κ) dκ = R 2π 2π 0 f Re iθ f (Re iθ ) e iθ dθ (76)
that relates the integer number ν of zeros of f (κ) in the complex κ-plane inside an arbitrary closed contour C which, for numerical purposes, we choose to be a circle of radius R that is centered at the origin.

Solving (76) numerically using a standard trapezoidal quadrature with ω ∈ C and M ∈ R for increasing values of R, demonstrates convergence to ν = 2 with a good accuracy that is evident in Fig. 7. This result is supported by the fact that both in the shallow water limit (when at κ → 0 we have tanh(κ) ≈ κ) and

ω i = -0.1.
in the deep water limit (when tanh(κ) = 1 as κ → +∞), dispersion relation (75) reduces to a quadratic polynomial in κ which therefore has only 2 distinct roots. Now, when the number of zeros in expression ( 75) is established, we can use an iterative algorithm such as the standard Newton-Raphson method to locate precisely where these zeros lie in the complex κ-plane. In order for the algorithm to initiate, we need to provide an initial guess that is close enough to the exact value of the desired root. A natural choice is to use the roots of the dispersion relation (75) either in the shallow-or in the deep water limit

κ SW ± = ω M ± 1 , κ DW ± = 2ωM + 1 ± √ 4ωM + 1 2M 2 , (77)
which depends on the length Γ of the membrane. Indeed, since ω = ω n + iω i , where ω n = nπM w /Γ is the free membrane frequency [START_REF] Mougel | Flutter and resonances of a flag near a free surface[END_REF] and ω i ∈ R is fixed, we have a clear numerical evidence that the zeros of (75) tend either to κ SW ± as Γ → ∞ or to κ DW ± as Γ → 0, see Fig. 8. Since the location of the poles of the integrand in (74) is now determined, we continue our investigation by applying the residue theorem to this improper integral in the same way as we did in the shallow water limit. This time, due to the presence of a pole at the origin, the contour we consider is decomposed as

C = [-R, -] ∪ ∆ ∪ [ , R] ∪ ∆ R , ∆ = { e it , -π ≤ t ≤ 0}, ∆ R = {Re it , 0 ≤ t ≤ π}, (78)
where 0 < 1 and R 1. Using the oriented curve (78) in the contour integral (74) yields

C F (κ)dκ = - -R F (κ)dκ + ∆ F (κ)dκ + R F (κ)dκ + ∆ R F (κ)dκ. ( 79
)
Adopting the same argument as in the shallow water case and taking into account that F (z) is a continuous function for z ∈ ∆ R and that lim R→+∞ |g(Re iθ )| = 0, where g(κ) = F (κ)e -iκx and θ ∈ [0, π],

we find that the contribution of the arc integral over ∆ R vanishes as we enlarge the radius, according to Jordan's lemma [START_REF] Ablowitz | Complex Variables. Introduction and Applications[END_REF].

The contribution of the arc integral over ∆ in (79) also vanishes as → 0 because lim

→0 ∆ F (κ)dκ = lim →0 i 0 -π e iθ F ( e iθ )dθ = 0. ( 80
)
where the matrices in (87) can be explicitly recovered through the following expressions [START_REF] Vedeneev | Panel flutter at low supersonic speeds[END_REF] I ij = δ ij ,

K ij = jπM w Γ 2 δ ij , P ij = 2 Γ Γ 0 P (ξ (0) j ; x, ω)ξ (0) i dx, (88)
with δ ij standing for the Kronecker delta. The function P follows from (72) after application of the Leibniz rule [START_REF] Le Dizès | Radiative instability in stratified vortices[END_REF]:

P (ξ(x); x, ω) = iω (u * v) -M u * ∂v ∂x + u(x, ω)v(0, ω) , (89)
where v(x) is given by (82) and

∂v ∂x = i 2π +∞ -∞ [κ -(ω -κM ) 2 tanh κ]e iκx (ω -κM ) 2 -κ tanh κ dκ ( 90
)
should be computed separately with the similar approach. Indeed, integral (90) is nothing else but a slightly modified version of expression (74), without the pole located at the origin. Therefore, a similar analysis applied to (90) reduces it to

∂v ∂x = iv(x) = - 2 j=1 res (F (κ), κ j ) , κ j ∈ C. (91)
With the expressions (82) and ( 91) we can recover an explicit form of all the matrices constituting the matrix pencil (87) of the nonlinear eigenvalue problem by direct numerical computation.

5.3. Jacobian of the nonlinear matrix pencil F (ω) Nonlinear eigenvalue problems constitute nowadays a challenging and ongoing research topic for a whole community of mathematicians. State-of-the-art reviews [START_REF] Mehrmann | Nonlinear eigenvalue problems: A challenge for modern eigenvalue methods[END_REF][START_REF] Bindel | Localization theorems for nonlinear eigenvalue problems[END_REF][START_REF] Güttel | The nonlinear eigenvalue problem[END_REF] identify and classify different classes of methods to solve them, depending on the nonlinearity. Most of the known methods are either based on the Newton-Raphson iterative process or use contour integration [START_REF] Bindel | Localization theorems for nonlinear eigenvalue problems[END_REF]. In the present paper we prefer the former.

The Newton-Raphson iterative process needs the derivative of the pencil (87) with respect to ω, or Jacobian of the system. In order to reach an acceptable convergence rate for the method, we compute the analytical form of the Jacobian instead of approximating it numerically.

From (87) we find

dF dω = -2ωI + α dP dω , (92)
where the derivative dP dω involves the derivative of the function P (x, ω) as defined in (89). Applying the Leibniz rule [START_REF] Le Dizès | Radiative instability in stratified vortices[END_REF] once again, we obtain

∂P ∂ω = i (u * v) + iω ∂u ∂ω * v + u * ∂v ∂ω -M ∂u ∂ω * ∂v ∂x + u * ∂ 2 v ∂x∂ω + ∂u ∂ω (x, ω)v(0, ω) + u(x, ω) ∂v ∂ω (0, ω) , ( 93
)
where u is given by (73), v by (82), and ∂v ∂x by (91). As in the previous case, the higher-order derivatives of (74) involve slightly different improper integrals.

Their explicit expressions are recovered as

∂v ∂ω = 1 π +∞ -∞ (ω -κM )(tanh 2 κ -1)e iκx [(ω -κM ) 2 -κ tanh κ] 2 dκ, ∂ 2 v ∂x∂ω = i 2π +∞ -∞ κ(ω -κM )(tanh 2 κ -1)e iκx [(ω -κM ) 2 -κ tanh κ] 2 dκ. ( 94
)
Despite the numerator in the integrand of (94) is different from (74), the poles remain identical to (90), with the only difference that they are of second order. Therefore, we can apply the same procedure involving the residue theorem as we did in the derivation of (82) (with a particular attention to the pre-factors of the integrals). As the poles are no longer simple, we need to use the expression (84) for the residues when computing derivatives (94).

Numerical evaluation of the integral in (88)

In our approach, the integrals in the expression (88) and its derivative, as well as in all the convolutions, are integrated using numerical quadrature. Since the domain of integration is one-dimensional and we want to keep high accuracy in our numerical scheme, we prefer to use the Legendre-Gauss-Lobatto quadrature rule to approximate every integral. The nodes and weights of this spectral collocation method are recovered using the Golub-Welsch algorithm [START_REF] Golub | Calculation of Gauss quadrature rules[END_REF], which is based on the inversion of a linear system obtained from the threeterm recurrence relation for Legendre polynomials. The computed eigenvalues correspond to the quadrature points, while the eigenvectors are used to recover the weights. In this method, the nodes are defined on an interval [-1, 1] before being mapped, using a linear transformation, to the interval of integration [0, Γ]. This spectral quadrature is able to reach computer accuracy with less than 20 nodes of discretization of the interval [0, Γ] and is thus used all over our code.

Newton-like numerical method for solving the nonlinear eigenvalue problem

Now when every element of the matrix pencil (87) and its derivative (92) can be recovered by direct numerical computation, we shall introduce the method that we will be using throughout the stability analysis of the nonlinear eigenvalue problem. This numerical method is introduced in the recent review [START_REF] Mehrmann | Nonlinear eigenvalue problems: A challenge for modern eigenvalue methods[END_REF] and is designed to solve the characteristic equation det F (ω) = 0 from the inversion of successive linear problems. As always in Newton-like methods, it has to start with an initial guess that is close enough to the exact solution ω. In our case, due to strong nonlinearity in ω in the term (89), we shall restrict ourselves to reasonably low values of α. Indeed, this parameter acts as a linear factor at the nonlinear in ω operator and keeping it sensibly small prevents us from departing too far from the free membrane solution (corresponding to α = 0). This allows us to initiate our algorithm by choosing the eigenfrequency of the free membrane [START_REF] Mougel | Flutter and resonances of a flag near a free surface[END_REF] as an initial starting point. With this first guess ω (0) = ω n , the method of successive linear problems is an iterative routine, where the p-th iteration requires to solve a linear eigenvalue problem

F (ω (p) )u = θF (ω (p) )u. ( 95
)
After inversion of expression (95), we re-initiate the method as follows

ω (p+1) = ω (p) -θ, (96)
where θ is chosen to be the smallest eigenvalue of (95) in the absolute value. As expected from the method, we can easily reach quadratic convergence in ω. The computational method presented through this section has been fully implemented and parallelized in MATLAB using the Parallel Computing Toolbox available from the software. All stability maps in the parameter spaces that will be presented in the next section are recovered from the direct computation of growth rates Im(ω) of the nonlinear matrix pencil (87).

Im(ω)

Finite Depth Shallow Water

Stability maps for the finite-chord Nemtsov membrane in the finite-depth flow

Applying the computational method described above to the nonlinear eigenvalue problem defined by the matrix pencil (87), we find frequencies and growth rates of the finite-chord Nemtsov membrane coupled to the free surface flow of finite depth. First, we are benchmarking our method against the analytical solution (41) in the shallow-water limit, corresponding to Γ → ∞, and shown in Fig. 2. In Fig. 9 we show an analogue of Fig. 2 for the membrane with the chord length Γ = 10 and α = 10 -4 that demonstrates all the structural characteristics that are present in the shallow water stability map including the tongues of intertwining flutter regions and the wide stability gap. Nevertheless, one can observe that in the case of the finite-chord membrane in the finite-depth flow some of the tongues are either separated to individual instability islands or merged into continuous instability belts. These new effects are caused by the finite values of the chord length of the membrane and the finite depth of the fluid layer.

In a similar way, we produce stability maps in the (µ, M )-plane computed from the numerical solution of the algebraic nonlinear matrix pencil (87) for Γ = 10 and α = 10 -4 . Comparing the results shown in Fig. 10 with the analytical solutions in the shallow water limit that are visualized in Fig. 4, we notice good qualitative and quantitative agreement of the two approaches.

A drawback of our first-order in α analytical expression (41) for the growth rates in the shallow water approximation was that the size of the membrane Γ in it was only playing the role of a scaling factor and thus did not change the shape of the instability domains as it is the case for the finite depth solution. In contrast, the numerical solution of the nonlinear eigenvalue problem with the pencil (87) allows us to investigate zones of the radiation-induced flutter in a broad range of variation of the chord length Γ.

In Fig. 11(a,b) we present stability maps in the (Γ, M )-plane for the membrane modes with n = 1, 2 that show stability close to the critical values M = 1 and Γ = 0 with instability (blue domain) everywhere else. Notice that as Γ → ∞, the lower boundary of the flutter domain for the finite depth layer tends to the horizontal neutral stability curve (shown as a red dotted line) that follows from the first-order in α expansion of the growth rate (41) in the shallow water approximation. On the other hand, in the opposite limit of Γ → 0, corresponding to the deep water approximation, the finite-chord-length Nemtsov membrane is stable, in accordance with the perturbation analysis in section 4. These observations confirm that our numerical results for the finite depth problem are in a very good agreement with the analytical treatment presented in the previous sections and that our numerical procedure converges to correct eigenvalues.

Exploring the stability map in the (Γ, M )-plane further for the higher-order membrane modes with n = 3 unveils even more intriguing pattern shown in Fig. 11(c,d). First, Fig. 11(c) highlights a curious structure of two prominent subdomains in the parameter plane with the stability gap corresponding to a sector of stability also visible in the left part of Fig. 9(c). Second, the lower subdomain in Fig. 11(c) decomposes to a bunch of instability tongues spreading along the Γ-axis, see Fig. 9(d). Third, for every tongue, even for that commencing very close to the M -axis (Fig. 9(e)), there is a critical value of the chord length Γ such that the shorter membranes are stable, quite in accordance with the analytical results for the deep water reported in section 4, cf. Fig. 6. Fourth, and, probably the most rewarding, is the evidence that in the limit of Γ → ∞ the boundaries of the instability tongues of Fig. 9(d) converge to the shallow water solution shown by red dotted horizontal lines in Fig. 9, which is also confirmed by the convergence of the corresponding growth rates shown in Fig. 12. Therefore, the pattern of instability tongues that we discovered first in the shallow water approximation manifests itself also in the general case of the finite-chord membrane in the finite-depth fluid flow. With the increase in n, the instability tongues start to break and intertwine, making the pattern even more intriguing, see Fig. 13 for the stability maps corresponding to n = 4.

As is evident in Fig. 14 modification of the added mass ratio parameter α by an order of magnitude from α = 10 -4 to α = 10 -3 deforms the pattern of instability regions. Although all the qualitative features remain in place, the behaviour of the stability boundaries at large values of Γ does not demonstrate a perfect convergence to the shallow water solution as one can see in the closer views of the stability domains in Fig. 14(e,f). Notice, however, that the red dotted lines are obtained from the first order in α perturbation expansion [START_REF] Ginzburg | About Doppler effect at superlight velocity[END_REF] of the full dispersion relation [START_REF] Williams | The Wiener-Hopf and residue calculus solutions for a submerged semi-infinite elastic plate[END_REF]. This fact explains the discrepancy and suggests that sensitivity of stability of the Nemtsov membrane to α has been proven to be rather important, quite in accordance with the remark of Barbone and Crighton [START_REF] Barbone | Vibrational modes of submerged elastic bodies[END_REF] that the intermediate values of the fluid-loading parameter appear to be the most complicated as it is impossible to label the mode of the fluid-solid system as corresponding solely to a 'solid mode' or a 'fluid mode'. This complication manifests itself in the fact that α is a factor of the nonlinear and non-polynomial in ω operator in (87). By this reason, as soon as α departs from the origin, the choice of the free membrane mode [START_REF] Mougel | Flutter and resonances of a flag near a free surface[END_REF] is less suitable as a first guess to initiate the Newton-like iterative process (96), which therefore works less stably at larger values of α. Parameter continuation proposed, e.g. in [START_REF] Alben | The flapping-flag instability as a nonlinear eigenvalue problem[END_REF] in a different setting, only slightly improved performance of our method, reflecting the fact that nonlinear eigenvalue problems of fluid-structure interaction are notoriously hard. Nevertheless, our method allowed to obtain new results in the classical problem in the unprecedentedly broad range of all other important parameters, including, first of all the chord length, dimensionless velocity of the flow, and speed of propagation of elastic waves in the membrane.

Concluding Remarks

In this paper we studied conditions for the onset of a radiation-induced instability of the Nemtsov membrane in a uniform flow with free surface. In contrast to previous works [START_REF] Labarbe | Membrane flutter induced by radiation of surface gravity waves on a uniform flow[END_REF][START_REF] Nemtsov | Flutter effect and emission in the region of anomalous and normal Doppler effects[END_REF] that were limited either by the shallow water approximation or by the assumption that the membrane has infinite chord length, we consider the problem in its entirety and take into account both the finite chord of the membrane and the finite depth of the fluid layer.

First, we derive a new integro-differential equation for the deflection of the finite-chord membrane that is coupled to the finite-depth flow.

Then, we develop an analytical procedure allowing to find the eigenvalues of the membrane that is weakly coupled to the flow in the shallow-and deep water approximations. Our original contribution is a systematic procedure that combines Laplace transform, residue calculus, and perturbation of eigenvalues. The analytical solution allowed us to plot detailed stability maps and find a new pattern of intertwining instability tongues that to the best of our knowledge has not been previously reported in the literature. Furthermore, we were able to find analytically the geometrical structure that governs position, orientation and self-intersections of the instability tongues.

Next, we developed an original numerical method to treat the finite-chord membrane in the finite-depth uniform flow with the free surface. With a combination of complex analysis and Galerkin discretization we reduced the boundary eigenvalue problem for an integro-differential equation to an algebraic non-polynomial nonlinear eigenvalue problem and solved it with a Newton-like method. This approach allowed us to explore the onset of instability with respect to the chord length of the membrane, velocity of the flow, and speed of elastic waves propagating along the membrane at small but finite values of the added mass ratio parameter that plays a role of an effective damping due to radiation of surface gravity waves.

We believe we have made a convincing case that the Nemtsov membrane is able to play the same paradigmatic role for understanding radiation-induced instabilities as the famous Lamb oscillator coupled to a string has played for understanding radiation damping. We believe that our procedure is applicable to a broad class of fluid-structure interaction problems that require solving nonlinear eigenvalue problems. An extension of it allowing for a significant continuation with respect to the coupling parameter has proven to be a harder topic that we leave beyond the scope of this paper.

Figure 2 :

 2 Figure 2: (Filled blue) Instability regions for (a,b) the first (n = 1) and (c,d) the fourth (n = 4) modes of the finite-chord Nemtsov membrane weakly coupled to the flow in the shallow water limit. The black dotted line is equation M = Mw. The green straight lines are given by Mw = (2j/n + 1)(M + 1) for j ∈ Z and the red dashed curves by Mw = k/n(M 2 -1) for k ∈ N. Notice absence of instability domains in a wide gap centered at Mw = M + 1 and corresponding to j = 0.

Figure 3 :

 3 Figure 3: Growth rates of perturbed frequencies of the membrane along the line (42) with α = 10 -3 , Γ = 10 and (a) n = 1 and j = 1 with zeros at M = 1 + 3/k, k ≥ 1 and (b) n = 4 and j = -1 with zeros at M = 1 + 2/k, k > 1.

Figure 4 :

 4 Figure 4: Stability maps with logarithmic scale for the growth rate (41) over µ = Mw/(M -1) and M . The dashed white curves and red lines represent the lines (42) and curves (43), respectively, for Γ = 1, α = 10 -4 and n according to the caption.

Figure 5 :

 5 Figure5: Surface 'modes' as given by expression[START_REF] Kirillov | Paradoxes of dissipation-induced destabilization or who opened Whitney's umbrella?[END_REF] for a fixed ω recovered by first-order in α approximation[START_REF] Abramovich | Instability in the oscillations of a moving oscillator while it radiates surface and internal waves[END_REF], with parameters n = 1, Mw = 1, Γ = 10, α = 10 -4 and M according to the legend. Dashed vertical lines specify x-coordinates of membrane's edges.

Figure 6 :

 6 Figure 6: (a): (Blue) Instability domain computed from expression (69) for n = 1, Mw = 2, and α = 0.04, with (dotted lines) the minimum (70) of Γ at the minimizer given by (71). (b): Growth rate for the left figure at a fixed value of Γ = 250.

  (ω, M ) = (3π -4i, 4.5)

Figure 7 :

 7 Figure 7: Convergence of expression (76) over the radius R for arbitrary pair of complex eigenfrequencies and Mach numbers. We used N = 20 nodes for the trapezoidal quadrature.

Figure 8 :

 8 Figure 8: Convergence of zeros of (75) to the roots (77) over Γ, for ω = ωn + ω i with fixed values of imaginary parts (a): ω i = 0.001 and (b): ω i = -0.1.

(a) n = 1 (b) n = 2 (

 12 c) n = 3 (d) n = 4

Figure 9 :

 9 Figure 9: Stability maps in the (Mw, M )-plane from the numerical solution of the matrix pencil (87) for Γ = 10, α = 10 -4 and n according to the legend.

Figure 10 :

 10 Figure 10: Stability maps in the (µ, M )-plane from the numerical solution of the matrix pencil (87) for Γ = 10, α = 10 -4 and n according to the legend. We use a logarithmic scale for the growth rate.

Figure 11 :

 11 Figure 11: Instability domains (blue) of the finite depth system and finite-chord membrane for the eigenfrequencies ω recovered from the matrix pencil (87). Parameters used are Mw = 1, α = 10 -4 and n according to the legend. Neutral stability line obtained from the first-order expansion of the shallow water growth rate (41) are shown as dotted red.

Figure 12 :Figure 13 :

 1213 Figure 12: Finite depth growth rates converge to the shallow water ones as Γ → ∞ for Mw = 1, α = 10 -4 , n = 3, and Γ = 100.

(a) n = 1 ( 4 Figure 14 :

 1414 Figure 14: Instability domains (blue) of the finite depth system and finite-chord membrane for the eigenfrequencies ω recovered from the matrix pencil (87). Parameters used are Mw = 1, α = 10 -3 and n according to the legend. The dotted red lines represent the neutral stability curves obtained from the first-order expansion of the shallow water growth rate (41).

Aknowledgments

We are grateful to Prof. M. L. Overton from the Courant Institute for helpful discussions. We thank the London Mathematical Society for supporting Prof. Overton's visit to Northumbria through the Scheme 4 Research in Pairs grant No 41820. J. L. was supported by a Ph.D. Scholarship from Northumbria University that also provided him with an opportunity to run the code on the High Performance Cluster.

Finally, taking the two limits of R → +∞ and → 0 simultaneously and then applying the residue theorem to (79), we obtain lim

where the improper integral has to be taken in the sense of Cauchy Principal Value and where κ j are the zeros of (75) lying in the upper-half plane.

A similar argument works as well for an oriented contour C in the lower half of the complex κ-plane, to enclose the poles with the negative imaginary parts. Now we are prepared to recover the function v(x, ω) from the general expression (74) as

which constitutes the main result of this section. Note that owing to the fact that the poles are computed numerically and that F in (82) cannot be expanded in the Laurent series, we need to calculate the residues in (82) with an alternative, however equivalent, expression to that used for obtaining [START_REF] Chandrasekhar | Solutions of two problems in the theory of gravitational radiation[END_REF]. Namely, if we consider a meromorphic function M (z) = a(z)/b(z) with a simple pole z 1 that is a root of b(z), then the residue for M (z) at z 1 reads as [START_REF] Ablowitz | Complex Variables. Introduction and Applications[END_REF] res (M (z),

If the function M (z) has a double pole z 2 , we find similarly that [START_REF] Ablowitz | Complex Variables. Introduction and Applications[END_REF] res (M (z),

Expressions ( 83) and (84) will be utilized in the numerical treatment of equation ( 72), for instance, in the computation of the coresponding Jacobian.

Galerkin discretization and reduction to an algebraic nonlinear eigenvalue problem

The last step in solving the integro-differential equation (72) numerically is to introduce a modal form for the displacement ξ(x) that respects the boundary conditions ξ(0) = ξ(Γ) = 0. For this purpose, we introduce the following Galerkin decomposition, based on a superposition of modes of a free membrane vibrating in vacuum

Substituting (85) into (72) and using the orthogonality of the modes (85) while integrating over the membrane chord, we find

The expression (86) represents the j-th scalar equation of the algebraic nonlinear eigenvalue problem in ω, which matrix pencil can be written as