2. Mathematical formulation for SCWs and their superharmonic normal forms

We consider surface gravity waves of an inviscid and incompressible fluid on deep water within a nonrotating medium. The motion is assumed to be irrotational, resulting in a velocity potential $\phi (x,y,z,t)$ , where $\eta (x,y,t)$ is the sea surface elevation with respect to the water level at rest. To formulate a generalized flow, we take $1/k$ and $(gk)^{-1/2}$ as our reference length and time, where k is the wave number of the incident wave and g is the gravitational acceleration [Reference Hsu, Tsuchiya and Silvester4, Reference Roberts10]. A new frame of reference ( $x^*,y^*,z^*,t^*,\phi ^*)$ moving with the wave and defined in [Reference Hsu, Tsuchiya and Silvester4, Reference Roberts10] is used to enhance the resolution: $x^*=x-ct$ , $y^*=y$ , $z^*=z$ , $t^*=t$ and $\phi ^*=\phi -cx^*$ , where c is the wave velocity equal to $\omega /m$ , $\omega $ is the angular frequency and $m=\sin \theta $ is the wave number in the seawall in the x-direction, $n=\cos \theta $ being the wave number in the transverse y-direction. Once the asterisks are omitted for the sake of simplicity, the flow equations are

(2.1) $$ \begin{align} \phi_{xx}+\phi_{yy}+\phi_{zz} &= 0 \quad\mbox{for } z\leqslant\eta(x,y), \end{align} $$

(2.2) $$ \begin{align} \phi_z &= 0 \quad\mbox{for } z\rightarrow-\infty, \end{align} $$

(2.3) $$ \begin{align} \eta_t+\phi_x\eta_x+\phi_y\eta_y-\phi_z&= 0\quad \mbox{on } z = \eta(x,y), \end{align} $$

(2.4) $$ \begin{align} \phi_t+\eta+\tfrac{1}{2}(\phi_x^2+\phi_y^2+\phi_z^2-c^2) &= 0\quad \mbox{on } z = \eta(x,y), \end{align} $$

where (2.1) is the linear Laplace continuity equation, (2.2) is the linear nonporous bottom condition, and (2.3) and (2.4) are the nonlinear kinematic and Bernoulli dynamical conditions at the free surface. To address both the SCW and the stability issues, we require the general solutions [Reference Ioualalen and Kharif5, Reference Ioualalen and Kharif6], with $t=t^*$ for simplicity, to be

(2.5) $$ \begin{align} \eta(x,y,t) &= \bar{\eta}(x,y)+\eta'(x,y,t), \end{align} $$

(2.6) $$ \begin{align} \phi(x,y,z,t) &= \bar{\phi}(x,y,z)+\phi'(x,y,z,t), \end{align} $$

where $\bar {\eta }(x,y)$ and $\bar {\phi }(x,y,z)$ , the SCW solutions of permanent form, are superimposed on the infinitesimal unsteady perturbations $\eta '(x,y,t)$ and $\phi '(x,y,z,t)$ , with $\eta '\ll \bar {\eta }$ and $\phi '\ll \bar {\phi }$ . After substituting forms (2.5) and (2.6) into the system of equations (2.1)–(2.4), then retaining only linear terms when perturbations and their derivatives are present, and finally ranging all components in the first two powers of the perturbation terms (powers 0 and 1), we obtain two distinct systems of equations: at the zeroth order of the perturbation elements (unperturbed wave field), we obtain the nonlinear SCWs equations for solutions ( $\bar {\eta },\bar {\phi }$ ); and at their first order, we obtain a linear system of equations solving the linear stability problem for solutions ( $\eta ', \phi '$ ), that is, the search of the normal forms [Reference Ioualalen and Kharif5, Reference Ioualalen and Kharif6].

At the zeroth order, the system of equations (2.1)–(2.4) is

(2.7) $$ \begin{align} \bar{\phi}_{xx}+\bar{\phi}_{yy}+\bar{\phi}_{zz} &= 0 \quad\mbox{for } z\leqslant\bar{\eta}(x,y),\nonumber\\\bar{\phi}_z &= 0 \quad\mbox{for } z\rightarrow-\infty,\nonumber\\\bar{\phi}_x\bar{\eta}_x+\bar{\phi}_y\bar{\eta}_y-\bar{\phi}_z&= 0\quad \mbox{on } z = \bar{\eta}(x,y), \end{align} $$

(2.8) $$ \begin{align} \bar{\eta}+\tfrac{1}{2}(\bar{\phi}_x^2+\bar{\phi}_y^2+\bar{\phi}_z^2-c^2) &= 0\quad \mbox{on } z = \bar{\eta}(x,y). \end{align} $$

The wave steepness is defined as the half of the nondimensional peak-to-trough height:

$$ \begin{align*} h=\frac{\bar{\eta}(0,0)-\bar{\eta}(\pi/m,0)}{2}. \end{align*} $$

Two methods have been considered for the computation of the SCW solutions of permanent forms. The first is the asymptotic method, which was proposed by Roberts [Reference Roberts10], and was rebuilt and used by Ioualalen and Kharif [Reference Ioualalen and Kharif5]. The derived solutions are unique and, following Roberts [Reference Roberts10], the radius of convergence of their power series is increased beyond singularities (poles) by using the Padé approximation through Shanks’ transform [Reference Wynn12]. It is fair to say that the solution is no longer valid around the singularity, at the pole location. Unfortunately, this is precisely the location where Ioualalen and Kharif [Reference Ioualalen and Kharif5] found their sporadic instability domain. The second method is a direct numerical procedure proposed by Okamura [Reference Okamura9] who used a Newton procedure. The method has no technical limitations, so it reveals the entire set of SCW solutions.

At first order, the system of equations (2.1)–(2.4) writes

(2.9) $$ \begin{align} \phi^{\prime}_{xx}+\phi^{\prime}_{yy}+\phi^{\prime}_{zz} &= 0 \quad\mbox{for } z\leqslant\bar{\eta}(x,y), \end{align} $$

(2.10) $$ \begin{align} \phi^{\prime}_z &= 0 \quad\mbox{for } z\rightarrow-\infty, \end{align} $$

(2.11) $$ \begin{align} \eta'(\bar{\phi}_{zz}-\bar{\eta}_x\bar{\phi}_{xz}-\bar{\eta}_y\bar{\phi}_{yz})-\bar{\eta}_x\phi_x^{\prime}-\bar{\phi}_x\eta_x^{\prime} -\bar{\eta}_y\phi_y^{\prime}-\bar{\phi}_y\eta_y^{\prime}+\phi_z^{\prime}-\eta_t^{\prime}= 0\quad \mbox{on } z = \bar{\eta}(x,y), \end{align} $$

(2.12) $$ \begin{align} \phi_t^{\prime}+\bar{\phi}_x\phi_x^{\prime}+\bar{\phi}_y\phi_y^{\prime}+\bar{\phi}_z\phi_z^{\prime} +\eta'(1+\bar{\phi}_x\bar{\phi}_{xz}+\bar{\phi}_y\bar{\phi}_{yz}+\bar{\phi}_z\bar{\phi}_{zz}) =0\quad \mbox{on } z = \bar{\eta}(x,y),\end{align} $$

for which we look for nontrivial superharmonic solutions of the form,

(2.13) $$ \begin{align} \eta'(x,y,t) &= e^{-i\sigma t}\sum_{K=-\infty}^{\infty}\sum_{L=-\infty}^{\infty}a_{KL}e^{i(Kmx+Lny)}, \end{align} $$

(2.14) $$ \begin{align} \phi'(x,y,z,t) &= e^{-i\sigma t}\sum_{K=-\infty}^{\infty}\sum_{L=-\infty}^{\infty}b_{KL}e^{i(Kmx+Lny)} e^{\gamma_{KL}z}, \end{align} $$

where $\gamma _{KL}=[(Km)^2+(Ln)^2]^{1/2}$ . Complex eigenvalues ( $\sigma $ ) and their respective eigenvector components ( $a_{KL}, b_{KL}$ ) are computed using the (spectral) Galerkin method of Ioualalen and Kharif [Reference Ioualalen and Kharif5, Reference Ioualalen and Kharif6]. For a given frequency $\Im \textrm {m}(-\sigma _{K,L})$ , $\Re e(-\sigma _{K,L})=0$ corresponds to a neutrally stable $(K, L)$ perturbation mode, $\Re e(-\sigma _{K,L})<0$ to a stable one and $\Re e(-\sigma _{K,L})>0$ applies for an instability.

3. SCW solutions. A brief review on harmonic resonances (HRs)

The semi-numerical asymptotic method of Roberts [Reference Roberts10] and the fully numerical Newton method of Okamura [Reference Okamura9] both look for N-truncated doubly periodic permanent solutions of the form,

$$ \begin{align*} \bar{\eta}(x,y)&=\sum_{I=0}^N\sum_{J=2-I.{\mathrm{mod}}(2)}^N\eta_{I,J}\cos(Imx)\cos(Jny),\\\bar{\phi}(x,y,z)&=\sum_{I=0}^N\sum_{J=2-I.{\mathrm{mod}}(2)}^N\phi_{I,J}\sin(Imx)\cos(Jny)e^{\gamma_{IJ}z}, \end{align*} $$

with $\gamma _{JK}=[(Jm)^2+(Kn)^2]^{1/2}$ . Harmonical coefficients $\eta _{I,J}$ and $\phi _{I,J}$ are obtained directly with the numerical method of Okamura [Reference Okamura9], which is also detailed by Ioualalen and Okamura [Reference Ioualalen and Okamura7]. Those of the asymptotic method of Roberts [Reference Roberts10] are obtained through successive powers of the small parameter h, previously defined as the wave steepness,

(3.1) $$ \begin{align} \bar{\eta}(x,y)=\sum_{r=1}^N\eta_rh^r , \quad \bar{\phi}(x,y,z)=\sum_{r=1}^N\phi_rh^r , \quad \omega=\sum_{r=0}^N\omega_rh^r. \end{align} $$

The symmetry of the wave propagation at the vertical wall imposes the first-order linear solution to be of the form,

$$ \begin{align*}\omega_0=1, \quad \eta_1=\cos(x)\cos(y) \quad \text{and} \quad \phi_1=\sin(x)\cos(y)e^z.\end{align*} $$

The flow symmetry also forces the general solutions to be of the following form, where harmonics I and J must be of same parity [Reference Roberts10]:

$$ \begin{align*} \eta_r&=\sum_{I,J}^r\eta_{r,I,J}\cos(Imx)\cos(Jny),\\\phi_r&=\sum_{I,J}^r\phi_{r,I,J}\sin(Imx)\cos(Jny)e^{\gamma_{IJ}z}. \end{align*} $$

The fully numerical method of Okamura [Reference Okamura9] is based on a Newton procedure. The system of equations (2.7) and (2.8) is solved at any collocation point of a computational domain covering one wavelength in both horizontal directions X and Y. The procedure starts with a first-guess solution (FGS) and Newton’s procedure operates several iterations until convergence is reached. Then, the next wanted solution, at a different h, is obtained by taking the previous converged solution as an FGS, that is, the choice of the FGS to solve the problem prescribes the future solution. This iterative procedure allows to calculate a segment/branch of the solution. However, in some circumstances, convergence is not reached, indicating the presence of a singular point. An alternative FGS is then necessary, permitting a jump from one branch solution to another around the singularity. Then, the iterative procedure continues on the second branch. This efficient iterative method may then reveal the multiple branches of the SCW solution.

Roberts [Reference Roberts10] was at the beginning of the SCW studies, and revealed most of their characteristics and behaviours, thanks to his semi-numerical method. Harmonics ( $I,J$ ) are ranged and their respective coefficients are computed once for all for a given incident angle $\theta $ . Then, the wave expressions (3.1) are rapidly derived for any chosen value of the expansion parameter h, which is defined as the wave steepness. To any given value of h corresponds a unique solution of the wave height $\eta $ , velocity potential $\phi $ and frequency $\omega $ . It is inherent to a perturbation method [Reference Roberts10]. The asymptotic method used by Ioualalen and Kharif [Reference Ioualalen and Kharif6] is totally inspired by that of Roberts [Reference Roberts10]. He prescribed a nonsecular condition to describe the harmonic coefficients as power series of the wave height h instead of the arbitrary expansion parameter $\epsilon $ which has no physical meaning, while Ioualalen and Kharif [Reference Ioualalen and Kharif6] operate an inversion of the power series of the wave height h, that is, a function of $\epsilon $ , and then replace the parameter $\epsilon $ inside the coefficient expansions, to obtain a power series of h as well. The obtained solutions are then rigorously identical.

The combination of (2.7) and (2.8) leads to the well-known equation of a natural oscillator, the free surface. Roberts [Reference Roberts10] showed that its right-hand side may contain, through nonlinear interactions, a term proportional to $\sin (IX)\cos (IJ)e^{\gamma _{IJ}z}$ that is a solution of the homogeneous equation for some critical angles of resonance $\theta _{HR}$ , when $\cos ^2\theta _{HR}=(I^4-I^2)/(J^2-I^2)$ on deep water (Table 1, up to the order of truncation $N=11$ ).

Table 1

Resonance angles $\theta _{HR}(I,J)$ on deep water (in $^\circ $ ), up to order 11 of truncation [Reference Roberts10]. The flow symmetry forbids mixed parity for harmonics I and J [Reference Roberts10]. Thus, values $J=5, 7, 9, 11, \ldots $ , are not allowed for $I=2$ ; similarly for values $J=10$ and $I=3$ . In deep water, harmonic resonances (HRs) appear at $J=4$ and $J=9$ for $I=2$ and $I=3$ , respectively, at angle $0^\circ $ [Reference Roberts10].

To accurately compute HRs, Roberts [Reference Roberts10] analysed carefully the convergence of the power series (3.1). One must ensure that the expansion parameter h is smaller than the radius of convergence $R_c$ of the series and eventually increase $R_c$ by some techniques, like the Shanks transform [Reference Wynn12], if $h>R_c$ . To illustrate, let us consider angles $\theta =52.3^\circ $ and $\theta =53^\circ $ , which are close to the resonance angle $\theta _{HR}=52.2388^\circ $ (Table 1). At the angle $\theta _{HR}$ , a resonance occurs for a flat surface ( $h=0$ ). At angles $52.3^\circ $ and $53^\circ $ , it occurs at a nonzero wave steepness. Let us proceed with the computation of the d’Alembert criterion [Reference d’Alembert2] for the frequency $\omega (h)$ , that is, the ratio between two successive coefficients of its power expansion. Since $\omega $ retains only even orders of the expansion (3.1) for symmetry reasons [Reference Roberts10], it is developed in powers of $H = h^2$ . Thus, for simplicity, we report the ratio $(\omega _{n+2}/\omega _{n})^{1/2}$ in Figure 2.

Figure 2

Computed ratio $(\omega _{n+2}/\omega _{n})^{1/2}$ as a function of the SCW order of truncation n up to order $n = 73$ for angles $\theta =52.3^\circ $ ( $\circ $ ) and $\theta =53^\circ $ $(\vartriangle )$ ; $\omega _{n}$ being the component of the frequency at order $r=n$ (see (3.1)). The asymptotic method of Roberts [Reference Roberts10] is used.

At angle $\theta =52.3^\circ $ , where there is no other nearby resonance angles than $\theta _{HR}=52.2388^\circ $ up to order $n=73$ , the ratio converges into the single value $h\approx 0.054728$ (Figure 2). The process begins after order $n=6$ where HR (2,6) appears. According to the criterion, it can be considered as the “local” radius of convergence. The radius of convergence of the power series may also be approximated by using the Padé approximation [Reference Gilewicz3]. The pole of the (17,17) approximant is $h=0.054726$ which is in full agreement.

At angle $53^\circ $ , the ratio drops into successive thresholds due to the presence of nearby resonance angles. The first is the same $\theta _{HR}=52.2388^\circ $ , appearing at order 6 and triggering a first convergence at $h\approx 0.191135$ . The value h is larger than the value for angle $\theta =52.3^\circ $ because angle $53^\circ $ is farther from $\theta _{HR}=52.2388^\circ $ . The ratio then trends to a second value $h\approx 0.0749$ due to HR (4,26) appearing at order 26 ( $\theta _{HR}=52.9133^\circ $ , not shown in Table 1), but operating at higher order; and finally around $h\approx 0.0231$ due to HR (5,41) ( $\theta _{HR}=52.9918^\circ $ ). The successive drops may potentially result in a zero-ratio, which is the effective radius of convergence of all power series (3.1). Furthermore, we may expect that SCW solutions will not converge because of the occurrence of HRs for every angle $\theta $ . Following Roberts [Reference Roberts10], Ioualalen and Kharif [Reference Ioualalen and Kharif5] explored low-order HRs ( $I < 4$ ) by using the Shanks transform [Reference Wynn12] to increase their radius of convergence. They ignored high-order resonances ( $I> 3$ ), conjecturing that high-order HRs would be damped through viscosity in real ocean waves.

The coefficient $\phi _{2,6}$ of the resonant component (2,6) is plotted in Figure 3 versus the wave steepness h at the nontrivial angle of resonance $\theta =52.3^\circ $ , which is close to the resonance angle 52.2388 of HR (2,6) (Table 1). The fully numerical Newton method of Okamura [Reference Okamura9] and the semi-numerical asymptotic procedure of Roberts [Reference Roberts10] are used. Using Newton’s method, branch (3) is obtained starting from a flat surface and increasing the wave steepness h (upper panel of Figure 3, in disks). A “severe” positive jump is then obtained around the critical singularity $h_c\approx 0.0547$ , which is precisely the radius of convergence and the pole described above using the asymptotic method. Below this critical value of h, branch (3) is similar to the solution obtained with the asymptotic method (in blue and red lines in Figure 3): both are weak and positive. The two solutions differ around the pole value. As expected, the asymptotic solution, without any processing, diverges at the singularity (in dashed blue without $\varepsilon $ -algorithm). However, when the radius of convergence is extended using the Shanks transform [Reference Wynn12], here the $\varepsilon $ -algorithm ( $\varepsilon $ -alg.), the processed solution operates a slight jump into weak negative values (red line, with $\varepsilon $ -alg.). The passage from positive to negative values happens around the pole. The asymptotic solution cannot be reliable around this point.

Figure 3

(a) Coefficient $\phi _{2,6}$ of the harmonics (2,6) of the velocity potential versus the wave steepness h at the incidence angle $\theta =52.3^\circ $ using Newton’s method. The branches of the solution are numbered (1)–(3), branches (1) and (2) being connected through a turning point (TP). Note branch (2) is negative, although not visible at this scale. Branch (3), spanning from point ( $P_1$ ) to point ( $P_2$ ), is plotted in small disks when stable, large disks when unstable and plus signs when normal form solutions have not converged. Similarly, for branches (1) and (2), there are small and large triangles and plus signs extending from point ( $P_4$ ) to point ( $P_3$ ) through the TP. For reference, there are also plotted solutions obtained with an asymptotic method with (solid red) and without (dashed blue) use of the Shanks transform ( $\varepsilon $ -algorithm). (b) Enlargement around the TP for better visibility. In particular, the passage between stability regimes for branches (1) and (3) and around TP is clearly delineated (sizes of the disks and triangles for ease to read).

Figure 4

The coefficient $\phi _{1,1}$ of the harmonies (1,1) of (a) the velocity potential and (b) the wave frequency $\omega $ versus the wave steepness h at incidence angle $\theta =52.3^\circ $ , and their three solution branches (1), (2) and (3).

Using the Newton method past the pole, the remaining branches (2) and (3) are obtained; thanks to the choice of a judicious first-guess solution in the iterative process. Past the bifurcation at $h_c$ , the complete multiple solution is thus obtained. Branch (3) escapes from branches (2) and (1), while those latter are connected through a turning point (TP). It is interesting to note the similarity between branch (2) and the processed asymptotic solution: both are weak, negative and have the same trend. Around the TP, the complete and the processed asymptotic solutions differ significantly: both solutions operate a jump from positive to negative value, but to a larger extent for the former one.

To summarize, using the fully numerical method, we obtain the effective set of solutions, that is, a single branch prior to the TP, branch (3), and a multi-branch solution subsequent to the TP, which is composed of the continuation of branch (3) and the two connected branches (1) and (2). The asymptotic method, with or without the use of the $\varepsilon $ -alg., first provides the accurate segment of branch (3) for h smaller than the radius of convergence of the coefficient expansion. Then, around the radius of convergence, that is, at the TP, when complemented with the Shanks transform, the transform connects continuously the lower segment of branch (3) (before the singularity) to an approximate of branch (1) after the singularity, resulting in an overall single solution [Reference Gilewicz3, Reference Wynn12]. This single solution can thus be considered as “hybrid”, because the transform operates a matching between two distinct branches. The use of the Shanks transform seems efficient around and past the singularity by increasing the radius of convergence of the power series. However, it still remains an artificial procedure because it uses one particular Padé approximant, bearing in mind that Padé approximants may differ drastically [Reference Gilewicz3]. Moreover, Gilewicz [Reference Gilewicz3] also proposed other alternatives such as his so-called best Padé approximant (BAP) to increase radii of convergence. Thus, the use of the “hybrid” solution is an interesting alternative, but should be avoided when an exact set of solutions is available, (like that of Okamura) [Reference Okamura9] which we use hereafter. This is precisely what was done by Roberts [Reference Roberts10], by also developing a fully numerical method after his asymptotic one [Reference Roberts and Schwartz11].

Once the resonant component (2,6) forces a multiple solution past the TP, the entire set of coefficients of both the wave heights, the velocity potential and the wave frequency behave similarly (Figure 4). The full resonance of the SCW field can thus be described as a “contamination” of the entire SCW spectrum by the resonant harmonics (2,6) through interactions between frequency components that are present in nonlinear terms of (2.7) and (2.8). Once the SCW field is obtained, its mode of occurrence and persistency is further analysed through the computation of its normal forms.

4. Stability of resonant SCWs: their normal forms

We consider the resonance (2,6) on deep water versus the wave steepness h. As described by Ioualalen and Kharif [Reference Ioualalen and Kharif5], a harmonic resonance $(\pm K,L)$ corresponds to a superharmonic instability that may occur only when the two eigenvalues $\sigma _{K,L}^s=\sigma _{-K,L}^{-s}$ of opposite signatures coalesce. The collision of the two eigenmodes $(\pm K,L)$ is then interpreted as Ioualalen and Kharif’s [Reference Ioualalen and Kharif6] class Ia $(K,L)$ instability. It corresponds to a resonance between the two eigenmodes $(\pm K,L)$ and the $2K$ modes $(1,\pm 1)$ of the basic SCW of permanent form,

$$ \begin{align*} \Omega_1 &= -\Omega_2+K\Omega_{01}+K\Omega_{02},\\ \boldsymbol{k}_1 &= \boldsymbol{k}_2+K\boldsymbol{k}_{01}+K\boldsymbol{k}_{02}, \end{align*} $$

where $\Omega _i=|{\boldsymbol {k}}_i|^{1/2}$ , $\Omega _{0i}= |{\boldsymbol {k}}_{0i}|^{1/2}$ for $i=1,2$ , ${\boldsymbol {k}}_1=(mK,nL)$ and ${\boldsymbol {k}}_2=(-mK,nL)$ being the wave numbers of the perturbation modes, and ${\boldsymbol {k}}_{01}=(m,n)$ and ${\boldsymbol {k}}_{02}=(m,-n)$ being those of the basic wave.

The stability problem is solved using the procedure of Ioualalen and Kharif [Reference Ioualalen and Kharif5]. The perturbation forms (2.13) and (2.14) are truncated at order M and the SCWs of permanent form $\bar {\eta }(x,y,z)$ and $\bar {\phi }(x,y)$ are computed using the Newton method of Okamura [Reference Okamura9]. Once both expressions are substituted into the first-order system of equations (2.9)–(2.12), we obtain a generalized eigenvalue problem of the form ${\boldsymbol {A}}{\boldsymbol {u}}=i\sigma {\boldsymbol {B}}{\boldsymbol {u}}$ that we solve with the QZ algorithm [Reference Moler and Stewart8], where $\sigma $ are the eigenvalues and ${\boldsymbol {u}}=(a_{jk},b_{jk})^{T}$ their corresponding eigenvectors. Complex matrices ${\boldsymbol {A}}$ and ${\boldsymbol {B}}$ are functions of the basic flow. The Galerkin numerical procedure of Ioualalen and Kharif [Reference Ioualalen and Kharif5] is used to solve the eigenvalue problem. Equations (2.11) and (2.12) are numerically integrated over one wavelength in the two horizontal directions using fast Fourier transforms (FFTs) over $\nu \times \mu $ points of coordinates $\alpha x_u=2\pi u/\nu , u=0,\ldots ,\nu -1$ and $\beta y_v=2\pi v/\mu , v=0,\ldots ,\mu -1$ . A generalized eigenvalue problem is obtained at the free surface $z=\bar {\eta }(x,y)$ ,

$$ \begin{align*} &\sum_{K=-M}^{M} \sum_{L=-M}^{M} F_{K-r,L-s}\{E_{KL}^{(1)}\} a_{KL} + \sum_{K=-M}^{M} \sum_{L=-M}^{M} F_{K-r,L-s}\{G_{KL}^{(1)}\}b_{KL} =i\sigma a_{rs}, \\ &\sum_{K=-M}^{M} \sum_{L=-M}^{M} F_{K-r,L-s}\{E_{KL}^{(2)}\} a_{KL} + \sum_{K=-M}^{M} \sum_{L=-M}^{M} F_{K-rL-s}\{G_{KL}^{(2)}\} b_{KL} \\ &\quad=i\sigma \sum_{K=-M}^{M} \sum_{L=-M}^{M} F_{K-r,L-s}\{H_{KL}^{(2)}\} b_{KL}, \end{align*} $$

with

$$ \begin{align*} E_{KL}^{(1)} &=-\bar{\phi}_{zz}+imK\bar{\phi}_x+inL\bar{\phi}_y+ \bar{\eta}_x\bar{\phi}_{xz}+ \bar{\eta}_y\bar{\phi}_{yz}, \\ G_{KL}^{(1)} &=(imK\bar{\eta}_x +inL\bar{\eta}_y-\gamma_{KL}) e^{\gamma_{KL}\bar{\eta}}, \\ E_{KL}^{(2)} &=1+\bar{\phi}_x\bar{\phi}_{xz}+\bar{\phi}_y\bar{\phi}_{yz}+\bar{\phi}_z\bar{\phi}_{zz},\\ G_{KL}^{(2)} &=(imK\bar{\phi}_x +inL\bar{\phi}_y+\gamma_{KL}\bar{\phi}_z) e^{\gamma_{KL}\bar{\eta}}, \\ H_{KL} &=e^{\gamma_{KL}\bar{\eta}}. \end{align*} $$

The functions

$$ \begin{align*} F_{K-r,L-s}\{f_{KL}\}=\sum_{u=0}^{\nu-1}\sum_{v=0}^{\mu-1}f_{KL}e^{i\alpha(K-r)x_u}e^{i\alpha(L-s)y_v}, \end{align*} $$

where $r=-M,\ldots ,M$ and $s=-M,\ldots ,M$ , are computed using two-dimensional FFTs. Convergence of the eigenvalues $\sigma $ is obtained by increasing M. In the following, all computations are carried out up to order $M=49$ if necessary.

In Figure 5, we have plotted the frequencies [ $-\Re \textrm {e}(\sigma _{\pm 2,6})$ ] of the eigenvalues $\sigma _{\pm 2,6}$ along with their growth rates [ $-\Im \textrm {m}(\sigma _{\pm 2,6})$ ] for branches (1) to (3) of the SCW field. Let us first consider branch (3). As expected, as h increases from point (P1) towards the TP (Figures 3 and 5), the frequencies of the harmonics (−2,6) and (2,6) converge to each other while remaining neutrally stable. Then, they coalesce when they both vanish and give rise to an instability beyond the TP. As described by Ioualalen and Kharif [Reference Ioualalen and Kharif5], a coalescence at zero-frequency is consistent with a resonance, since both harmonic components propagate at the same phase speed as the basic SCW field. They are phase-locked, which is precisely the definition of a resonance phenomenon.

Figure 5

Frequency [ $-\Re \textrm {e}(\sigma _{\pm 2,6})$ ] (dots) and growth rate [ $-\Im \textrm {m}(\sigma _{\pm 2,6})$ ] (+ signs) versus the wave steepness h for angle $\theta =52.3^\circ $ . Positive values of the frequencies and growth rates correspond to eigenvalues $\sigma _{2,6}$ and negative values correspond to $\sigma _{-2,6}$ . (a) Solution branch (3), and (b) branches (2) and (1). Points ( $P_1$ ) to ( $P_4$ ) and TP are the same as in Figure 3.

Beyond a certain value of the wave steepness h, Ioualalen and Kharif [Reference Ioualalen and Kharif5], using an asymptotic method, observed a re-stabilization. They described the h-range of instability as a “bubble”-like structure with a specific maximum growth rate. In contrast, with our complete SCW solution, and in particular for our branch (3), we do not observe any re-stabilization, at least as far as the normal form problem converges, which is not the case past a certain value of the wave steepness (Table 2). We even find an inflection point from $h\approx 0.10$ enhancing the growth rate instead of reaching a maximum (Figure 5).

For the other branches, starting from point (P3) towards the TP all along branch (2), the frequencies converge to each other and also give rise to an instability around the TP at the beginning of branch (1) (Figures 3 and 5). Within branch (1) from the TP to point (P4), an instability arises with no re-stabilization of the converged solutions (Table 2, Figure 5), in contrast to that observed by Ioualalen and Kharif [Reference Ioualalen and Kharif5]. Similarly to the unstable segment of branch (3), an inflection point is also observed from the same wave steepness $h\approx 0.10$ , with very similar growth rates.

We propose the following reason for the re-stabilization computed by Ioualalen and Kharif [Reference Ioualalen and Kharif5], that is, the existence of their “bubble”-like structure. Their hybrid SCW solution operates a jump from the h-left segment of our branch (3) to branch (2). As we show here, both parts are stable. Around the jump, the solution operates a jump from the beginning of the unstable branch (3) to the beginning of the unstable branch (1). The passage automatically requires a narrow h-range of instability, that is, their “bubble”. It is thus due to the hybrid but incomplete structure of their solution.

Ultimately, in our case study, the stability analysis indicates a possible passage from partially unstable segments of branches (1) and (3) onto the always-neutrally stable branch (2) for any given wave steepness h (Figures 3 and 6). Regardless of the wave steepness h, a subcritical bifurcation may always operate from an unstable branch onto a stable one (Figures 3 and 6). There exists always a pre-existing (neutrally) stable branch solution upon which a transition from an unstable segment may occur (Figure 6). Therefore, the physical appearance of HRs is still questionable, despite the possible large growth rate of their associated instability.

Table 2

Convergence of the growth rates $-\Im \textrm {m}(\sigma _{\pm 2,6})$ versus the wave steepness h at incidence angle $\theta =52.3^\circ $ . Three values of the parameter h are taken for each branch (3) ( $h_3$ ) and (1) ( $h_1$ ), $h=0.1386$ and $h=0.1129$ being respectively the highest values where convergence is reached.

Figure 6

Schematic representation of the SCW three-solution branches (1), (2) and (3) for coefficient $\phi _{2,6}$ of the velocity potential versus the wave steepness h for angle $\theta =52.3^\circ $ (Figure 3). Dashed lines apply for unstable segments while solid lines are for neutrally stable ones. Arrows represent possible jumps form unstable to stable branch segments.