Benjamin-Feir instability of
Stokes waves in finite depth

A classical problem in fluid dynamics, pioneered by the famous work of Stokes [38] in 1847, concerns the spectral stability/instability of periodic traveling waves –called Stokes waves– of the gravity water waves equations in any depth.

Benjamin and Feir [3], Lighthill [32] and Zakharov [43, 45] discovered in the sixties, through experiments and formal arguments, that Stokes waves in deep water are unstable, proposing an heuristic mechanism which leads to the disintegration of wave trains. More precisely, these works predicted unstable eigenvalues of the linearized equations at the Stokes wave, near the origin of the complex plane, corresponding to small Floquet exponents $\mu$ or, equivalently, to long-wave perturbations. The same phenomenon was later predicted by Whitham [41] and Benjamin [2] for Stokes waves of wavelength $2\pi\kappa$, in finite depth $\mathtt{h}$, provided that $\kappa\mathtt{h}>1.363$ approximately. This phenomenon is nowadays called “Benjamin-Feir” –or modulational– instability, and it is supported by an enormous amount of physical observations and numerical simulations, see e.g. [15, 33]. We refer to [46] for an historical survey.

A serious difficulty for a rigorous mathematical proof of the Benjamin-Feir instability is that the perturbed eigenvalues bifurcate from the eigenvalue zero, which is defective, with multiplicity four. The first rigorous proof of a local branch of unstable eigenvalues close to zero for $\kappa\mathtt{h}$ larger than the Whitham-Benjamin threshold $1.363\ldots$ was obtained by Bridges-Mielke [9] in finite depth (see also the preprint by Hur-Yang [23]). Their method, based on a spatial dynamics and a center manifold reduction, breaks down in deep water. For dealing with this case Nguyen-Strauss [35] have recently developed a new approach, based on a Lyapunov-Schmidt decomposition. The novel spectral approach developed in Berti-Maspero-Ventura [6] allowed to fully describe, in deep water, the splitting of the four eigenvalues close to zero, as the Floquet exponent is turned on, proving in particular the conjecture that a pair of non-purely imaginary eigenvalues depicts a closed figure “8”, parameterized by the Floquet exponent.

The goal of this paper is to describe the full Benjamin-Feir instability phenomenon at any finite value of the depth $\mathtt{h}>0$. This analysis has fundamental physical importance since real-life experiments are performed in water tanks (for example the original Benjamin and Feir experiments, in Feltham’s National Physical Laboratory, had Stokes waves of wavelength 2.2 m and bottom’s depth of 7.62 m, see [2]). We also remark that the Benjamin-Feir instability mechanism is a possible responsible of the emergence of rogue waves in the ocean, we refer to [28, 29] and references therein for a vast physical literature. A first mathematically rigorous treatment of large waves is given in [18], via a probabilistic analysis, in the case of NLS.

Along this paper, with no loss of generality, we consider $2\pi$-periodic Stokes waves, i.e. with wave number $\kappa=1$. In Theorems 2.5 and 1.1 we prove the existence of a unique depth $\mathtt{h}_{\scriptscriptstyle{\textsc{WB}}}$, in perfect agreement with the Benjamin-Feir critical value 1.363…, such that:

• •

  Shallow water case: for any $0<\mathtt{h}<\mathtt{h}_{\scriptscriptstyle{\textsc{WB}}}$ the eigenvalues close to zero remain purely imaginary for Stokes waves of sufficiently small amplitude, see Figure 2(a)-left;

• •

  Sufficiently deep water case: for any $\mathtt{h}_{\scriptscriptstyle{\textsc{WB}}}<\mathtt{h}<\infty$, there exists a pair of non-purely imaginary eigenvalues which traces a complete closed figure “8” (as shown in Figure 2(a)-right) parameterized by the Floquet exponent $\mu$. By further increasing $\mu$, the eigenvalues recollide far from the origin on the imaginary axis where then they keep moving. As ${\mathtt{h}}\to\mathtt{h}_{\scriptscriptstyle{\textsc{WB}}}^{\,+}$ the set of unstable Floquet exponents shrinks to zero and the Benjamin-Feir unstable eigenvalues collapse to the origin, see Figure 3. This figure ‘8” was first numerically discovered by Deconink-Oliveras in [15].

We remark that our approach fully describes all the eigenvalues close to $0$, providing a necessary and sufficient condition for the existence of unstable eigenvalues, i.e. the positivity of the Benjamin-Feir discriminant function $\Delta_{\scriptscriptstyle{\textsc{BF}}}(\mathtt{h};\mu,\epsilon)$ defined in (1.6).

The results of Theorems 2.5 and 1.1 are complementary to those of [6]. In following the natural spectral approach developed in [6], we encounter a major difference in the proof, that we now anticipate. In the infinitely deep water ideal case it turns out that the “reduced” $4\times 4$ matrix obtained by the Kato reduction procedure is a small perturbation of a block-diagonal matrix which possesses yet the correct Benjamin-Feir unstable eigenvalues. This is not the case in finite depth. The correct eigenvalues of the “reduced” $4\times 4$ matrix emerge only after one non-perturbative step of block diagonalization. We shall explain in detail this point after the statement of Theorem 2.5. This is related with the fact that, in infinite deep water, among the four eigenvalues close to zero of the linearized operator at the Stokes wave, two are $\mathcal{O}(\mu)$, whereas the other two have much larger size $\mathcal{O}(\sqrt{\mu})$, whereas in finite depth all four eigenvalues have size $\mathcal{O}(\mu)$. In addition, along the whole proof, one needs to carefully track the explicit dependence with respect to $\mathtt{h}$ of the entries of the reduced $4\times 4$ matrix.

Let us now present rigorously our results.

In this section we present in detail the complete spectral Theorem 2.5. We first introduce the pure gravity water waves equations and the Stokes waves solutions.
The water waves equations. We consider the Euler equations for a 2-dimensional incompressible, irrotational fluid under the action of gravity. The fluid fills the region

$${\mathcal{D}}_{\eta}:=\left\{(x,y)\in\mathbb{T}\times\mathbb{R}\;:\;-\mathtt{h}\leq y<\eta(t,x)\right\}\,,\quad\mathbb{T}:=\mathbb{R}/2\pi\mathbb{Z}\,,$$

with finite depth and space periodic boundary conditions. The irrotational velocity field is the gradient of a harmonic scalar potential $\Phi=\Phi(t,x,y)$ determined by its trace $\psi(t,x)=\Phi(t,x,\eta(t,x))$ at the free surface $y=\eta(t,x)$. Actually $\Phi$ is the unique solution of the elliptic equation $\Delta\Phi=0$ in ${\mathcal{D}}_{\eta}$ with Dirichlet datum $\Phi(t,x,\eta(t,x))=\psi(t,x)$ and $\Phi_{y}(t,x,y)=0$ at $y=-\mathtt{h}$.

The time evolution of the fluid is determined by two boundary conditions at the free surface. The first is that the fluid particles remain, along the evolution, on the free surface (kinematic boundary condition), and the second one is that the pressure of the fluid is equal, at the free surface, to the constant atmospheric pressure (dynamic boundary condition). Then, as shown by Zakharov [44] and Craig-Sulem [13], the time evolution of the fluid is determined by the following equations for the unknowns $(\eta(t,x),\psi(t,x))$,

$$\eta_{t}=G(\eta)\psi\,,\quad\psi_{t}=-g\eta-\dfrac{\psi_{x}^{2}}{2}+\dfrac{1}{2(1+\eta_{x}^{2})}\big{(}G(\eta)\psi+\eta_{x}\psi_{x}\big{)}^{2}\,,$$

(2.1)

where $g>0$ is the gravity constant and $G(\eta):=G(\eta,\mathtt{h})$ denotes the Dirichlet-Neumann operator $[G(\eta)\psi](x):=\Phi_{y}(x,\eta(x))-\Phi_{x}(x,\eta(x))\eta_{x}(x)$. In the sequel, with no loss of generality, we set the gravity constant $g=1$, see Remark 2.4.

The equations (2.1) are the Hamiltonian system

$$\partial_{t}\begin{bmatrix}\eta\\ \psi\end{bmatrix}=\mathcal{J}\begin{bmatrix}\nabla_{\eta}\mathcal{H}\\ \nabla_{\psi}\mathcal{H}\end{bmatrix},\quad\quad\mathcal{J}:=\begin{bmatrix}0&\mathrm{Id}\\ -\mathrm{Id}&0\end{bmatrix},$$

(2.2)

where $\nabla$ denote the $L^{2}$-gradient, and the Hamiltonian $\mathcal{H}(\eta,\psi):=\frac{1}{2}\int_{\mathbb{T}}\left(\psi\,G(\eta)\psi+\eta^{2}\right)\mathrm{d}x$ is the sum of the kinetic and potential energy of the fluid. In addition of being Hamiltonian, the water waves system (2.1) possesses other important symmetries. First of all it is time reversible with respect to the involution

$$\rho\begin{bmatrix}\eta(x)\\ \psi(x)\end{bmatrix}:=\begin{bmatrix}\eta(-x)\\ -\psi(-x)\end{bmatrix},\quad\text{i.e. }\mathcal{H}\circ\rho=\mathcal{H}\,.$$

(2.3)

Moreover, the equation (2.1) is space invariant, since, being the bottom flat,

$$\tau_{\theta}G(\eta)\psi=G(\tau_{\theta}\eta)[\tau_{\theta}\psi]\,,\quad\forall\theta\in\mathbb{R}\,,\quad\text{where}\quad\tau_{\theta}u(x):=u(x+\theta)\,.$$

In addition, the Dirichlet-Neumann operator satisfies $G(\eta+m,\mathtt{h})=G(\eta,\mathtt{h}+m)$, for any $m\in\mathbb{R}$.
Stokes waves. The Stokes waves are traveling solutions of (2.1) of the form $\eta(t,x)=\breve{\eta}(x-ct)$ and $\psi(t,x)=\breve{\psi}(x-ct)$ for some real $c$ and $2\pi$-periodic functions $(\breve{\eta}(x),\breve{\psi}(x))$. In a reference frame in translational motion with constant speed $c$, the water waves equations (2.1) become

$$\eta_{t}=c\eta_{x}+G(\eta)\psi\,,\quad\psi_{t}=c\psi_{x}-\eta-\dfrac{\psi_{x}^{2}}{2}+\dfrac{1}{2(1+\eta_{x}^{2})}\big{(}G(\eta)\psi+\eta_{x}\psi_{x}\big{)}^{2}$$

(2.4)

and the Stokes waves $(\breve{\eta},\breve{\psi})$ are equilibrium steady solutions of (2.4).

The bifurcation result of small amplitude of Stokes waves is due to Struik [39] in finite depth, and Levi-Civita [31], and Nekrasov [34] in infinite depth. We denote by $B(r):=\{x\in\mathbb{R}\colon\ |x|<r\}$ the real ball with center 0 and radius $r$.

The expansions (2.6)-(2.8) are derived in the Appendix B for completeness, although present in the literature (they coincide with [42, section 13, chapter 13] and [2, section 2]). Note that in the shallow water regime $\mathtt{h}\to 0^{+}$ the expansions (2.6)-(2.8) become singular. For the analiticity properties of the maps stated in Theorem 2.1 we refer to [8].

We also mention that more general time quasi-periodic traveling Stokes waves – which are nonlinear superpositions of multiple Stokes waves traveling with rationally independent speeds – have been recently proved for (2.1) in [5] in finite depth, in [16] in infinite depth, and in [4] for capillary-gravity water waves in any depth.
Linearization at the Stokes waves. In order to determine the stability/instability of the Stokes waves given by Theorem 2.1, we linearize the water waves equations (2.4) with $c=c_{\epsilon}$ at $(\eta_{\epsilon}(x),\psi_{\epsilon}(x))$. In the sequel we closely follow [6] pointing out the differences of the finite depth case.

By using the shape derivative formula for the differential $\mathrm{d}_{\eta}G(\eta)[\hat{\eta}]$ of the Dirichlet-Neumann operator one obtains the autonomous real linear system

$$\begin{bmatrix}\hat{\eta}_{t}\\ \hat{\psi}_{t}\end{bmatrix}=\begin{bmatrix}-G(\eta_{\epsilon})B-\partial_{x}\circ(V-c_{\epsilon})&G(\eta_{\epsilon})\\ -1+B(V-c_{\epsilon})\partial_{x}-B\partial_{x}\circ(V-c_{\epsilon})-BG(\eta_{\epsilon})\circ B&-(V-c_{\epsilon})\partial_{x}+BG(\eta_{\epsilon})\end{bmatrix}\begin{bmatrix}\hat{\eta}\\ \hat{\psi}\end{bmatrix}$$

(2.9)

where

$$V:=V(x):=-B(\eta_{\epsilon})_{x}+(\psi_{\epsilon})_{x}\,,\ \ B:=B(x):=\frac{G(\eta_{\epsilon})\psi_{\epsilon}+(\psi_{\epsilon})_{x}(\eta_{\epsilon})_{x}}{1+(\eta_{\epsilon})_{x}^{2}}=\frac{(\psi_{\epsilon})_{x}-c_{\epsilon}}{1+(\eta_{\epsilon})_{x}^{2}}(\eta_{\epsilon})_{x}\,.$$

(2.10)

The functions $(V,B)$ are the horizontal and vertical components of the velocity field $(\Phi_{x},\Phi_{y})$ at the free surface. Moreover $\epsilon\mapsto(V,B)$ is analytic as a map $B(\epsilon_{0})\to H^{\sigma,s-1}(\mathbb{T})\times H^{\sigma,s-1}(\mathbb{T})$.

The real system (2.9) is Hamiltonian, i.e. of the form $\mathcal{J}\mathcal{A}$ for a symmetric operator $\mathcal{A}=\mathcal{A}^{\top}$, where $\mathcal{A}^{\top}$ is the transposed operator with respect the standard real scalar product of $L^{2}(\mathbb{T},\mathbb{R})\times L^{2}(\mathbb{T},\mathbb{R})$.

Moreover, since $\eta_{\epsilon}$ is even in $x$ and $\psi_{\epsilon}$ is odd in $x$, then the functions $(V,B)$ are respectively even and odd in $x$, and the linear operator in (2.9) is reversible, i.e. it anti-commutes with the involution $\rho$ in (2.3).

Under the time-independent “good unknown of Alinhac” linear transformation

$$\begin{bmatrix}\hat{\eta}\\ \hat{\psi}\end{bmatrix}:=Z\begin{bmatrix}u\\ v\end{bmatrix}\,,\qquad Z=\begin{bmatrix}1&0\\ B&1\end{bmatrix},\quad Z^{-1}=\begin{bmatrix}1&0\\ -B&1\end{bmatrix},$$

(2.11)

the system (2.9) assumes the simpler form

$$\begin{bmatrix}u_{t}\\ v_{t}\end{bmatrix}=\widetilde{\mathcal{L}}_{\epsilon}\begin{bmatrix}u\\ v\end{bmatrix},\qquad\widetilde{\mathcal{L}}_{\epsilon}:=\begin{bmatrix}-\partial_{x}\circ(V-c_{\epsilon})&G(\eta_{\epsilon})\\ -1-(V-c_{\epsilon})B_{x}&-(V-c_{\epsilon})\partial_{x}\end{bmatrix}\,.$$

(2.12)

Note that, since the transformation $Z$ is symplectic, i.e. $Z^{\top}\mathcal{J}Z=\mathcal{J}$, and reversibility preserving, i.e. $Z\circ\rho=\rho\circ Z$, the linear system (2.12) is Hamiltonian and reversible as (2.9).

Next we perform a conformal change of variables to flatten the water surface. Here the finite depth case induces a modification with respect to the deep water case. By [1, Appendix A], there exists a diffeomorphism of $\mathbb{T}$, $x\mapsto x+\mathfrak{p}(x)$, with a small $2\pi$-periodic function $\mathfrak{p}(x)$, and a small constant $\mathtt{f}_{\epsilon}$, such that, by defining the associated composition operator $(\mathfrak{P}u)(x):=u(x+\mathfrak{p}(x))$, the Dirichlet-Neumann operator writes as [1, Lemma A.5]

$$G(\eta_{\epsilon})=\partial_{x}\circ\mathfrak{P}^{-1}\circ{\mathcal{H}}\circ\tanh\big{(}(\mathtt{h}+\mathtt{f}_{\epsilon})|D|\big{)}\circ\mathfrak{P}\,,$$

(2.13)

where ${\mathcal{H}}$ is the Hilbert transform, i.e. the Fourier multiplier operator

$$\mathcal{H}(e^{\mathrm{i}\,jx}):=-\mathrm{i}\,\textup{sign}(j)e^{\mathrm{i}\,jx}\,,\quad\forall j\in\mathbb{Z}\setminus\{0\}\,,\quad\mathcal{H}(1):=0\,.$$

The function $\mathfrak{p}(x)$ and the constant $\mathtt{f}_{\epsilon}$ are determined as a fixed point of (see [1, formula (A.15)])

$$\mathfrak{p}=\frac{\mathcal{H}}{\tanh\big{(}(\mathtt{h}+\mathtt{f}_{\epsilon})|D|\big{)}}[\eta_{\epsilon}(x+\mathfrak{p}(x))]\,,\qquad\mathtt{f}_{\epsilon}:=\frac{1}{2\pi}\int_{\mathbb{T}}\eta_{\epsilon}(x+\mathfrak{p}(x))\mathrm{d}x\,.$$

(2.14)

By the analyticity of the map $\epsilon\to\eta_{\epsilon}\in H^{\sigma,s}$, $\sigma>0$, $s>1/2$, the analytic implicit function theorem implies the existence of a solution $\epsilon\mapsto\mathfrak{p}(x):=\mathfrak{p}_{\epsilon}(x)$, $\epsilon\mapsto\mathtt{f}_{\epsilon}$, analytic as a map $B(\epsilon_{0})\to H^{s}(\mathbb{T})\times\mathbb{R}$. Moreover, since $\eta_{\epsilon}$ is even, the function $\mathfrak{p}(x)$ is odd. In Appendix B we prove the expansion

$$\mathfrak{p}(x)=\epsilon{\mathtt{c}}_{\mathtt{h}}^{-2}\sin(x)+\epsilon^{2}\frac{(1+{\mathtt{c}}_{\mathtt{h}}^{4})(3+{\mathtt{c}}_{\mathtt{h}}^{4})}{8{\mathtt{c}}_{\mathtt{h}}^{8}}\sin(2x)+\mathcal{O}(\epsilon^{3})\,,\quad\mathtt{f}_{\epsilon}=\epsilon^{2}\frac{{\mathtt{c}}_{\mathtt{h}}^{4}-3}{4{\mathtt{c}}_{\mathtt{h}}^{2}}+\mathcal{O}(\epsilon^{3})\,.$$

(2.15)

Under the symplectic and reversibility-preserving map

$$\mathcal{P}:=\begin{bmatrix}(1+\mathfrak{p}_{x})\mathfrak{P}&0\\ 0&\mathfrak{P}\end{bmatrix}\,,$$

(2.16)

the system (2.12) transforms, by (2.13), into the linear system $h_{t}=\mathcal{L}_{\epsilon}h$ where $\mathcal{L}_{\epsilon}$ is the Hamiltonian and reversible real operator

	$\displaystyle\mathcal{L}_{\epsilon}:=\mathcal{P}\,\widetilde{\mathcal{L}}_{\epsilon}\,\mathcal{P}^{-1}$	$\displaystyle=\begin{bmatrix}\partial_{x}\circ({\mathtt{c}}_{\mathtt{h}}+p_{\epsilon}(x))&|D|\tanh((\mathtt{h}+\mathtt{f}_{\epsilon})|D|)\\ -(1+a_{\epsilon}(x))&({\mathtt{c}}_{\mathtt{h}}+p_{\epsilon}(x))\partial_{x}\end{bmatrix}$		(2.17)
		$\displaystyle=\mathcal{J}\begin{bmatrix}1+a_{\epsilon}(x)&-({\mathtt{c}}_{\mathtt{h}}+p_{\epsilon}(x))\partial_{x}\\ \partial_{x}\circ({\mathtt{c}}_{\mathtt{h}}+p_{\epsilon}(x))&|D|\tanh((\mathtt{h}+\mathtt{f}_{\epsilon})|D|)\end{bmatrix}$

where

$${\mathtt{c}}_{\mathtt{h}}+p_{\epsilon}(x):=\displaystyle{\frac{c_{\epsilon}-V(x+\mathfrak{p}(x))}{1+\mathfrak{p}_{x}(x)}}\,,\quad 1+a_{\epsilon}(x):=\displaystyle{\frac{1+(V(x+\mathfrak{p}(x))-c_{\epsilon})B_{x}(x+\mathfrak{p}(x))}{1+\mathfrak{p}_{x}(x)}}\,.$$

(2.18)

By the analiticity results of the functions $V,B,\mathfrak{p}(x)$ given above, the functions $p_{\epsilon}$ and $a_{\epsilon}$ are analytic in $\epsilon$ as maps $B(\epsilon_{0})\to H^{s}(\mathbb{T})$. In the Appendix B we prove the following expansions.