3.1. Linear solutions

By solving the linearised water-wave equations, one obtains the linear standing-wave solution

(3.1) \begin{align} \begin{cases} \widetilde {\eta }_{k,l}(x,y,t) = \cos (k x)\cos (l y)\cos (\omega t), \\[5pt] \widetilde {\varphi }_{k,l}(x,y,t) = -\dfrac {\omega }{\kappa } \cos (kx)\cos (l y)\sin (\omega t), \end{cases} \end{align}

where $\kappa = \sqrt {k^2+l^2}$ , and $\omega$ satisfies the dispersion relation for deep-water gravity-capillary waves,

(3.2) \begin{align} \omega ^2 = \kappa + \kappa ^3. \end{align}

We use the following two assumptions when computing standing waves:

1. (i) $\varphi = 0$ at $t = 0,T/2,T,3T/2,\ldots$ ,

2. (ii) $\eta$ and $\varphi$ are even functions of $x$ and $y$ ,

where $T$ is the smallest temporal period. The first assumption implies that the kinetic energy vanishes every half-period. Although this also seems true for potential energy because $\widetilde {\eta }_{k,l}$ vanishes periodically in time, the potential energy of nonlinear standing waves never vanishes, as will be shown later. To construct standing waves in rectangular or square basins, we use the following three types of linear solutions.

1. (i) Case I (rectangular basin): $\eta = \epsilon \widetilde {\eta }_{k,l},\ \varphi = \epsilon \widetilde {\varphi }_{k,l}$ , $k,l\neq 0,\ k\neq l$ .

2. (ii) Case II (square basin): $\eta = \epsilon \widetilde {\eta }_{k,l} + \epsilon \widetilde {\eta }_{l,k},\ \varphi = \epsilon \widetilde {\varphi }_{k,l} + \epsilon \widetilde {\varphi }_{l,k}$ , $l = rk,\ r = 1,3,5,\ldots$ .

3. (iii) Case III (square basin): $\eta = \epsilon \widetilde {\eta }_{k,l} + \epsilon \widetilde {\eta }_{l,k},\ \varphi = \epsilon \widetilde {\varphi }_{k,l} + \epsilon \widetilde {\varphi }_{l,k}$ , $l = rk,\ r = 0,2,4,\ldots$ .

Here, $\epsilon$ is a small parameter. To support standing waves in square basins, it is required that $k$ and $l$ be rationally related, but we only study the case that $l$ is an integer multiple of $k$ . In general, one can also consider the following form of linear solution:

(3.3) \begin{align} \eta = \epsilon \widetilde {\eta }_{k,rk} - \epsilon \widetilde {\eta }_{rk,k},\quad\varphi = \epsilon \widetilde {\varphi }_{k,rk} - \epsilon \widetilde {\varphi }_{rk,k}. \end{align}

For $r = 0,2,4,\ldots$ and $y = Y+\pi /k$ , we have

(3.4) \begin{align} \cos (kx)\cos (rk y)-\cos (rkx)\cos (k y) = \cos (kx)\cos (rk Y)+\cos (rkx)\cos (k Y), \end{align}

thus (3.3) is equivalent to case III without specifying additional symmetries. On the other hand, for $r=1,3,5,\ldots$ , we find no convergent solution using (3.3). Note that

(3.5) \begin{align} \widetilde {\eta }_{k,l}(x,y,t) = \frac {1}{2}\widetilde {\eta }_{\kappa ,0}(\xi ,\zeta ,t) + \frac {1}{2}\widetilde {\eta }_{0,\kappa }(\xi ,\zeta ,t), \end{align}

where $\xi = (kx+ly)/\kappa$ and $\zeta = (kx-ly)/\kappa$ . This means that a three-dimensional standing wave can be decomposed into obliquely oriented two-dimensional standing waves, at least linearly. In particular, when $k=l$ , a case II standing wave can be transformed to a case III standing wave by rotating the $(x,y)$ -plane by $\pi /4$ .

3.2. Spatio-temporal symmetries

We consider triply periodic solutions in a three-dimensional cube

(3.6) \begin{align} \{(x,y,t)\mid -L_1\le x\le L_1,\, -L_2\le y\le L_2,\, 0\le t \le T\}, \end{align}

where $2L_1$ and $2L_2$ are the smallest periods in the $x$ - and $y$ -directions. Using the time-reversal argument, solutions must have a temporal symmetry about $t = T/2$ :

(3.7) \begin{align} \eta (x,y,T/2-t') &= \eta (x,y,T/2+t'), \\[-12pt]\nonumber \end{align}

(3.8) \begin{align} \varphi (x,y,T/2-t') &= -\varphi (x,y,T/2+t'), \end{align}

where $t'\in \mathbb R$ . Since $\varphi$ vanishes at $t = T/2$ , $\eta$ must returns to its initial state with a spatial shift (we do not consider other possibilities in this paper). Otherwise, $T/2$ becomes the smallest temporal period. To ensure consistency with linear solutions, there exist two types of shift.

1. (i) For cases I and II,

   (3.9) \begin{align} \eta (x,y,0) = \eta (x\pm L_1,y,T/2) = \eta (x,y\pm L_2,T/2). \end{align}

2. (ii) For case III,

   (3.10) \begin{align} \eta (x,y,0) = \eta (x\pm L_1,y\pm L_2,T/2). \end{align}

Using the time-reversal argument again, solutions must be symmetric about $t = T/4$ .

1. (i) For cases I and II,

   (3.11) \begin{align} \eta (x,y,T/4-t') &= \eta (x\pm L_1,y,T/4+t') = \eta (x,y\pm L_2,T/4+t'), \\[-12pt]\nonumber \end{align}
   (3.12) \begin{align} \varphi (x,y,T/4-t') &= -\varphi (x\pm L_1,y,T/4+t') = -\varphi (x,y\pm L_2,T/4+t'). \end{align}

2. (ii) For case III,

   (3.13) \begin{align} \eta (x,y,T/4-t') &= \eta (x\pm L_1,y\pm L_2,T/4+t'), \\[-12pt]\nonumber \end{align}
   (3.14) \begin{align} \varphi (x,y,T/4-t') &= -\varphi (x\pm L_1,y\pm L_2,T/4+t'). \end{align}

Consequently, we only need to consider a quarter temporal period, say $0\le t\le T/4$ , in computation.

Owing to the even symmetry with respect to $x$ and $y$ (termed fourfold symmetry hereafter), i.e.

(3.15) \begin{align} \eta (x,y,t) &= \eta (-x,y,t) = \eta (-x,-y,t) = \eta (x,-y,t), \\[-12pt]\nonumber \end{align}

(3.16) \begin{align} \varphi (x,y,t) &= \varphi (-x,y,t) = \varphi (-x,-y,t) = \varphi (x,-y,t), \end{align}

we only need to consider the values of $\eta$ and $\varphi$ within one quadrant, say $\{(x,y,t)\mid -L_1\le x\le 0,\ {-L_2}\le y\le 0,\ 0\le t\le T/4\}$ . For cases I and II, (3.11)–(3.12) imply

(3.17) \begin{align} \eta (x,y,t) &= \eta (x\pm L_1,y\pm L_2,t), \\[-12pt]\nonumber \end{align}

(3.18) \begin{align} \varphi (x,y,t) &= \varphi (x\pm L_1,y\pm L_2,t). \end{align}

Combining these with the fourfold symmetry, we have

(3.19) \begin{align} \eta (-x,-y,t) &= \eta (x\pm L_1,y\pm L_2,t), \\[-12pt]\nonumber \end{align}

(3.20) \begin{align} \varphi (-x,-y,t) &= \varphi (x\pm L_1,y\pm L_2,t). \end{align}

This property is used for case I standing waves to reduce the number of unknowns, as shown in figure 1.

Figure 1.

Typical contours of $\eta$ for a case I standing wave. The rectangles surrounded by the black dashed and solid lines represent a unit periodic cell and a quadrant on the $(x,y)$ -plane, respectively. The shaded region denotes the real computational domain by using (3.19) and (3.20).

For cases II and III, solutions have a stronger eightfold symmetry, i.e. fourfold symmetry plus even symmetry about the two diagonals $x=\pm y$ :

(3.21) \begin{align} \eta (x,y,t) &= \eta (y,x,t)= \eta (-y,-x,t), \\[-12pt]\nonumber \end{align}

(3.22) \begin{align} \varphi (x,y,t) &= \varphi (y,x,t)= \varphi (-y,-x,t). \end{align}

Combining these with (3.17) and (3.18), we have the following identities for case II:

(3.23) \begin{align} \eta (x,y,t) &= \eta (y\pm L,x\mp L,t) = \eta (-y\pm L,-x\mp L,t), \\[-12pt]\nonumber \end{align}

(3.24) \begin{align} \varphi (x,y,t) &= \varphi (y\pm L,x\mp L,t) = \varphi (-y\pm L,-x\mp L,t), \end{align}

where $L = L_1 = L_2$ . This means that case II solutions are also symmetric about the four straight lines $x+y = \pm L$ , $x-y = \pm L$ . We use these symmetries to reduce the computational domain from a quadrant to one of its quarter sub-triangles, as shown in figure 2. For case III, however, solutions are only symmetric about the two diagonals, as shown in figure 3. Moreover, there exist additional spatial symmetries for all three cases at $t=T/4$ . Setting $t' = 0$ in (3.11)–(3.14) and using the fourfold symmetry, we have the following.

Figure 2.

Typical contours of $\eta$ for a case II standing wave. The squares surrounded by the black dashed and solid lines represent a unit periodic cell and a quadrant on the $(x,y)$ -plane, respectively. The blue straight lines show the two diagonals $x = \pm y$ , and the four red straight lines denote $x+y = \pm L$ and $x-y = \pm L$ . The shaded region denotes the real computational domain by using (3.23) and (3.24).

Figure 3.

Typical contours of $\eta$ for a case III standing wave. The squares surrounded by the black dashed and solid lines represent a unit periodic cell and a quadrant on the $(x,y)$ -plane, respectively. The blue straight lines show the two diagonals $x = \pm y$ . The shaded region denotes the real computational domain by using (3.21) and (3.22).

Figure 3. Long description

Panel A: A contour plot representing a case III standing wave. The plot features blue and red contour lines indicating wave patterns. The black dashed and solid lines form a square representing a unit periodic cell and a quadrant on the x-y plane, respectively. The blue straight lines depict the two diagonals of the square. The axes are labeled x and y. Panel B: A zoomed-in view of the shaded region from Panel A, showing detailed contour lines within the real computational domain. The axes are also labeled x and y.

1. (i) For cases I and II,

   (3.25) \begin{align} \eta (-x,y,T/4) &= \eta (x\pm L_1,y,T/4) = \eta (x,y\pm L_2,T/4) = \eta (x,-y,T/4), \\[-12pt]\nonumber \end{align}
   (3.26) \begin{align} \varphi (-x,y,T/4) &= -\varphi (x\pm L_1,y,T/4) = -\varphi (x,y\pm L_2,T/4) = \varphi (x,-y,T/4). \end{align}

2. (ii) For case III,

   (3.27) \begin{align} \eta (-x,-y,T/4) &= \eta (x\pm L_1,y\pm L_2,T/4), \\[-12pt]\nonumber \end{align}
   (3.28) \begin{align} \varphi (-x,-y,T/4) &= -\varphi (x\pm L_1,y\pm L_2,T/4). \end{align}

This means that at $t = T/4$ , case I and case II solutions are symmetric about $x = \pm L_1/2$ and $y = \pm L_2/2$ within each quadrant, and case III solutions are symmetric about $x+y = \pm L$ and $x-y = \pm L$ , as shown in figure 4.

Figure 4.

Typical contours of $\eta (x,y,T/4)$ for (a) case I, (b) case II and (c) case III standing waves. The rectangle and squares surrounded by the solid black lines each represent a quadrant on the $(x,y)$ -plane. The blue straight lines are the two diagonals $x = \pm y$ . The red straight lines denote $x+y = \pm L$ and $x-y = \pm L$ . The green lines indicate $x = \pm L_1/2$ and $y = \pm L_2/2$ . The shaded regions denote the real computational domain.

Note that the reflection symmetry imposed for standing waves is an assumption. Although it is clear that there exist asymmetric travelling-standing waves (Wilkening Reference Wilkening2021), it remains an open question regarding the existence of asymmetric standing waves.

3.3. Numerical implementation

We introduce $I+1,\ J+1$ and $N+1$ uniform collocation points for $x,\ y$ and $t$ , i.e.

(3.29) \begin{align} x_i &= \frac{2(i-1)L_1}{I}-L_1,\quad i = 1,2,\ldots ,I+1, \\[-12pt]\nonumber \end{align}

(3.30) \begin{align} y_{\!j} &= \frac{2(j-1)L_2}{J}-L_2,\quad j = 1,2,\ldots ,J+1, \\[-12pt]\nonumber \end{align}

(3.31) \begin{align} t_n &= \frac{(n-1)T}{N},\quad n = 1,2,\ldots ,N+1, \end{align}

and our task is to solve the function values of $\eta$ and $\varphi$ at these points, along with the frequency $\omega$ , given wave amplitude $H$ (this can be the wave amplitude in physical space, the amplitude of a fundamental Fourier mode, or a measure of energy). Based on the spatio-temporal symmetries, we have the following unknowns to solve for the case I standing waves.

1. (i) At $t=0$ , $\eta (x_i,y_{\!j},0)$ for $i=1,2,\ldots,I/4+1$ and $j=1,2,\ldots,J/2+1$ .

2. (ii) At $t = t_n$ ( $n=2,3,\ldots, N/4$ ), $\eta (x_i,y_{\!j},t_n)$ and $\varphi (x_i,y_{\!j},t_n)$ for $i=1,2,\ldots ,I/4+1$ and $j=1,2,\ldots ,J/2+1$ .

3. (iii) At $t=T/4$ , $\eta (x_i,y_{\!j},T/4)$ for $i=1,2,\ldots ,I/4+1$ and $j=1,2,\ldots ,J/4+1$ , $\varphi (x_i,y_{\!j},T/4)$ for $i=1,2,\ldots ,I/4$ and $j=1,2,\ldots ,J/4$ .

Figure 5.

Schematic of the reduced computational domain by using spatio-temporal symmetries.

All together, these give rise to $(I/4+1)(J/2+1)N/2 - J/4 \sim IJN/16$ unknowns for $\eta$ and $\varphi$ . For cases II and III, it is convenient to let $I=J$ . The total numbers of unknowns for $\eta$ and $\varphi$ are reduced to $(I/4+1)^2N/2 \sim I^2N/32$ for case II, and $(I/4+1)(I/2+1)N/2 \sim I^2N/16$ for case III. A schematic is drawn in figure 5 to show the reduced computational domain for all three cases. Equations (2.6) and (2.7) give rise to the same number of equations after using the spatio-temporal symmetries. To close the system, we use the wave-amplitude equation

(3.32) \begin{align} \eta (0,0,0) = H, \end{align}

for given values of $H$ .

The system is solved via Newton’s method by setting the iterative tolerance to $10^{-10}$ . The actual residuals for convergent solutions are usually several orders smaller than this value. All derivatives and operations involving the DNOs are calculated via a pseudo-spectral method and fast Fourier transform. Once we get a convergent solution, we use the continuation method to search for other solutions along the solution branch by gradually increasing $H$ . In the presence of turning points, we alternatively use $H$ and $\omega$ as the bifurcation parameters.

4.1. Comparisons with existing results

We first compare our numerical solutions with some existing results. For two-dimensional gravity-capillary standing waves in finite water depth $h$ , Concus (Reference Concus1962) derived a third-order Stokes expansion and the nonlinear dependence of $\omega$ on $\epsilon$ , the amplitude of the fundamental Fourier mode. For $k=1$ , this relation reads (using our scalings)

(4.1) \begin{align} \omega \sim \omega _0+\frac {1}{2}\epsilon ^2\omega _2 + O(\epsilon ^3), \end{align}

where

(4.2) \begin{align} \omega _0 =\sqrt {2\tanh (h)},\quad \omega _2 = \frac {-2\omega _0^5-39\omega _0-66\omega _0^{-3}+2592\omega _0^{-7}}{320(1-6\omega _0^{-4})}, \end{align}

and $\epsilon$ represents the amplitude of the $\cos (x)\cos (\omega t)$ mode. Figure 6 compares Concus’s asymptotic result and our numerical results of boundary-integral formulation (see Appendix B), the cubic model and the quintic model, for $k = 1, h = \infty$ . Figure 6(a) exhibits the bifurcation curve plotted under $\omega$ and $\epsilon$ . All results show perfect consistency and are indistinguishable when $\epsilon \lt 0.2$ . As $\epsilon$ gradually increases towards $0.4$ , the cubic model agrees well with Concus’s analytic result, while the quintic model gives almost identical results to the boundary-integral method, indicated by both the frequency curves and the wave profiles in figure 6(b).

Figure 6.

Comparisons between Concus’s asymptotic result (4.1), and the numerical results using the boundary-integral method, the cubic model and the quintic model for $k=1, h=\infty$ . (a) Plot of $\omega$ versus $\epsilon$ . (b) Initial wave profiles with crest-to-trough amplitude $0.8$ ( $\epsilon = 0.394$ ).

Figure 6. Long description

Panel A: A line graph plots omega against epsilon. The x-axis represents epsilon ranging from 0 to 0.4, and the y-axis represents omega ranging from 1.35 to 1.42. The graph includes four data series: Boundary integrals, Concus, Cubic model, and Quintic model. All data series show a decreasing trend as epsilon increases. Panel B: A line graph plots eta(x, 0) against x. The x-axis represents x ranging from negative pi to pi, and the y-axis represents eta(x, 0) ranging from negative 0.5 to 0.5. The graph includes four data series: Boundary integrals, Concus, Cubic model, and Quintic model. All data series show a wave-like pattern with a peak at x equals 0 and troughs at x equals negative pi and pi. An inset zooms in on the region from negative 0.5 to 0.5 on the x-axis and from 0.31 to 0.34 on the y-axis, showing detailed differences among the models.

Another comparison is made for the case II three-dimensional gravity standing waves in finite-depth water. Verma & Keller (Reference Verma and Keller1962) derived the second-order Stokes-type expansion and the asymptotic expression for the nonlinear frequency, which has the same form as (4.1). For $k=l=1$ , the coefficients are

(4.3) \begin{align} \omega _0 = \sqrt {\sqrt {2}\tanh (\sqrt {2}h)},\quad \omega _2 = -\frac {(3\omega _0^4-2)^2}{16\omega _0(\tanh (2h)-2\omega _0^2)} + \frac {36\omega _0^{-7}-24\omega _0^{-3}+5\omega _0-23\omega _0^{5}}{64}, \end{align}

and $\epsilon$ represents the amplitude of $\cos (x)\cos (y)\cos (\omega t)$ . Figure 7(a) compares the numerical results obtained from the the cubic and quintic models, the asymptotic result of Verma & Keller (Reference Verma and Keller1962), and the numerical result of Rycroft & Wilkening (Reference Rycroft and Wilkening2013), for $k = l = 1,\ h = \pi$ . The cubic model slightly deviates from other results when $\epsilon \gt 0.2$ , but still exhibits a qualitative agreement. In contrast, the quintic model agrees remarkably with the existing asymptotic and numerical results. In figure 7(b), we exhibit the initial wave profile for $H = 0.46$ ( $\epsilon = 0.3675$ ) obtained from the quintic model. Interested readers may compare it with those numerical solutions shown in Rycroft & Wilkening (Reference Rycroft and Wilkening2013) to see the similarity.

In the following sections we will only exhibit the numerical results obtained from the quintic model owing to its excellent consistency with the original full potential-flow formulation.

Figure 7.

Comparisons between the asymptotic result of Verma & Keller (Reference Verma and Keller1962), numerical results of Rycroft & Wilkening (Reference Rycroft and Wilkening2013), the cubic model and the quintic model for case II gravity standing waves with $k=l=1, h=\pi$ . (a) Plot of $\omega$ versus $\epsilon$ . (b) An initial profile obtained from the quintic model with $H=0.46$ ( $\epsilon =0.3675$ ).

4.2. Standing waves and bifurcations

Figure 8.

Surface profiles of a case II standing wave for $k = l = 1,\ H=0.7$ at $t = 0,T/4,T/2$ .

Figure 8. Long description

The image contains three 3D surface plots labeled as Panel A, Panel B, and Panel C. Each plot shows the surface profile of a standing wave with different parameters. Panel A: This 3D surface plot shows the surface profile of a standing wave with the function η(x,y,0) plotted against x and y axes. The x and y axes range from -π to π, and the z-axis ranges from approximately -1.5 to 0.5. The surface exhibits a complex wave pattern with peaks and troughs, and contour lines are drawn on the base to indicate the wave’s amplitude. Panel B: This 3D surface plot shows the surface profile of a standing wave with the function η(x,y,T/4) plotted against x and y axes. The x and y axes range from -π to π, and the z-axis ranges from approximately -0.6 to 0.5. The surface exhibits a different wave pattern compared to Panel A, with distinct peaks and troughs, and contour lines are drawn on the base. Panel C: This 3D surface plot shows the surface profile of a standing wave with the function η(x,y,T/2) plotted against x and y axes. The x and y axes range from -π to π, and the z-axis ranges from approximately -1.5 to 0.5. The surface exhibits another distinct wave pattern with peaks and troughs, and contour lines are drawn on the base.

Figure 9.

Energy distribution over one temporal period for the standing wave shown in figure 8. Here, $\mathcal K, \mathcal G, \mathcal S, \widetilde{\mathcal{H}}$ represent kinetic, gravitational, surface-tension and total energies.

We first show representative case II standing waves with $k = l = 1$ . In contrast to the gravity standing waves featuring pyramidal-shaped crests, gravity-capillary standing waves exhibit rounded crests and troughs. Figure 8 displays a solution with $H=0.7, \omega = 1.867$ , calculated using $128\times 128\times 128$ grid points. Figures 8(a–c) correspond to the wave profiles at $t=0,T/4,T/2$ , respectively. At $t=T/4$ , the surface develops multiple crests and troughs, with the maximum wave amplitude decreasing to $0.16$ . Figure 9 exhibits the energy distribution among the kinetic, gravitational and surface-tension parts over one temporal period. The kinetic energy vanishes every half-period, while the potential energy reaches its minimum $0.497$ at $t = T/4,3T/4,\ldots$ . The total energy remains at a constant level $7.661$ perfectly, with maximum relative error $7\times 10^{-8}$ .

Figure 10.

(a) Initial cross-sections of case II standing waves for $H=0.1,0.2,0.4,0.6,0.7$ . The red curve corresponds to $0.7\cos (x)\cos (y)$ . (b) Initial diagonal cross-sections of the same waves as in (a).

Figure 10 shows the initial cross-sections of case II standing waves with successively increasing wave amplitudes $H=0.1,0.2,0.4,0.6,0.7$ . Figure 10(a) corresponds to the cross-sections on the plane $y=0$ , which are nearly sinusoidal. As shown by the comparison with $0.7\cos (x)\cos (y)$ , the crest of the nonlinear solution is slightly broader than the linear solution. This is also shown in figure 10(b), where we plot the cross-sections along the diagonal $x=y=\xi /\sqrt {2}$ . Note that $\xi =\pm \pi /\sqrt {2}\approx \pm 2.22$ are nodes of linear standing waves, but the nonlinear solutions exhibit negative surface elevations at these points.

Figure 11.

(a–e) Surface profiles of a resonant case II standing wave with $k = l = 1,\ H=0.33$ at $t = 0,T/16,3T/16,7T/32,T/4$ . ( f) Amplitudes of the first $100$ Fourier coefficients for $\eta (x,y,0)$ .

A common feature for two-dimensional standing waves is the existence of various nonlinear resonant solutions (Smith & Roberts Reference Smith and Roberts1999; Wilkening Reference Wilkening2011; Wilkening & Yu Reference Wilkening and Yu2012; Abassi & Wilkening Reference Abassi and Wilkening2025). In three dimensions, Rycroft & Wilkening (Reference Rycroft and Wilkening2013) first computed resonant gravity standing waves. Here, we observe some similar phenomena. Figure 11 shows a representative case II resonant standing wave with $k=l=1,\ H=0.33, \omega = 2.0192$ , calculated using $128\times 128\times 128$ grid points. In contrast to the smooth surface profiles shown in figure 8, the resonant standing wave exhibits complex surface wrinkles. A local hump and dip alternately emerge on top of the crest at a faster oscillating frequency relative to the fundamental Fourier mode, as shown in figures 11(a–c). The surface develops a volcano-shaped appearance with four saddles along the diagonals at $t=7T/32$ , as shown in figure 11(d). At $t=T/4$ , the wave amplitude drops to $0.0247$ approximately, and the surface profile shows a similar pattern to that in figure 8(b), but has much flatter crests and four sharper spikes at $(\pm \pi /2,\pm \pi /2)$ . Figure 11( f) shows the Fourier amplitudes associated with the $\cos (nx)\cos (my)$ mode for $\eta (x,y,0)$ . Note that all Fourier modes with odd-valued sum $n+m$ vanish, as a consequence of (3.11) and (3.12). Owing to the eightfold symmetry, the spectrum is invariant under the interchange of $n$ and $m$ . As can be clearly seen, the three dominant Fourier components are $\cos (x)\cos (y)$ , $\cos (3x)\cos (5y)$ and $\cos (5x)\cos (3y)$ . According to the dispersion relation (3.2), the two latter modes have the same linear frequency $\sqrt {35\sqrt {34}} \approx 7.075 \omega$ , thus yielding a resonance of the seventh temporal harmonic. The eighth temporal harmonic, dominated by the $\cos (2x)\cos (6y)$ and $\cos (6x)\cos (2y)$ modes, is also nearly resonant, as indicated by linear frequency $\sqrt {41\sqrt {40}} \approx 7.975 \omega$ .

Figure 12.

(a) Cross-sections of the standing wave shown in figure 11 at $t = 0,T/32,T/16,\ldots ,T/4$ from bottom to top. The profiles are vertically shifted for the clarity. (b) Plot of $\eta (0,0,t)$ over one temporal period. (c) Amplitudes of the temporal Fourier modes $\cos (p\omega t)$ from $p = 0$ to $p=20$ .

Figure 12. Long description

Panel A: A line graph displays cross-sections of a standing wave at different vertical shifts for clarity. The horizontal axis represents the spatial coordinate x ranging from negative pi to pi, and the vertical axis represents the wave profile η(x, 0, t) ranging from negative 0.4 to 1.0. The lines show the wave profiles at various heights, with each profile shifted vertically for clarity. Panel B: A line graph shows the temporal variation of the wave profile at a specific point over one temporal period. The horizontal axis represents time t ranging from 0 to T, and the vertical axis represents the wave profile η(0, 0, t) ranging from negative 0.4 to 0.2. The line shows the periodic nature of the wave profile over time. Panel C: A bar graph displays the amplitudes of the temporal Fourier modes. The horizontal axis represents the Fourier mode index p ranging from 0 to 20, and the vertical axis represents the amplitude of cos(pwt) on a logarithmic scale ranging from 10^-5 to 10^0. The bars show the amplitudes of different Fourier modes, indicating the contribution of each mode to the wave profile.

The same resonant wave is shown in figure 12. Figure 12(a) presents the cross-sections of the wave profile on the plane $y=0$ over a quarter temporal period, exhibiting five shorter waves clearly. Figure 12(b) shows the time history of $\eta (0,0,t)$ over a temporal period, which indicates seven minor oscillations. To measure the amplitude of each temporal harmonic, we take the Fourier transform of $\eta (x,y,t)$ with respect to $t$ , and display the infinity norm associated with $\cos (p\omega t)$ . As can be clearly seen, the five dominant components correspond to $p = 1, 7, 2, 0,8$ .

Figure 13.

Bifurcations of case II standing waves for $k = l = 1$ , and representative initial wave profiles.

Figure 13 shows the bifurcation of case II standing waves for $k = l = 1$ , which grows from the infinitesimal linear solution with frequency $\omega = 2.0598$ . Following the primary branch where no resonances occur, solutions are characterised by smooth surfaces with rounded crests and troughs. A large-amplitude solution, which is labelled by the triangle  $D$ , is identical to the standing wave shown in figure 8. The nonlinear frequency $\omega$ monotonically decreases with $H$ , until $\omega \approx 2.026$ when the bifurcation breaks into two disjoint branches. The local structure is shown in the small box nearby. Moving along each branch in a specific direction, the resonance between the fundamental and seventh temporal harmonics gradually becomes significant. Two representative resonant solutions are labelled $A$ and $B$ . The latter solution, which features a local hump on top of the crest, is the same solution as shown in figure 11. In contrast, solution $A$ exhibits a local dip on the crest, revealing that the two branches of resonant solutions are out of phase in the seventh temporal harmonic. When $\omega \approx 1.93$ , we observe a second break-up. A resonant solution with $H = 0.59$ is shown, labelled $C$ . This break-up is caused by resonances among multiple modes, including the seventh, eighth and ninth temporal harmonics. While there exist more break-ups following the bifurcation, we do not present them because of the similar structure.

Figure 14.

Surface profiles of case III standing waves for $k = 1,\ l = 0$ . (a–c) A non-resonant solution with $H=0.3$ at $t = 0,T/4,T/2$ . (d) A resonant solution with $H=0.29$ at $t = 0$ .

Figure 14. Long description

Panel A: A 3D surface plot showing a non-resonant solution with eta(x, y, 0). The x-axis ranges from -π to π, and the y-axis ranges from -π to π. The z-axis represents the surface elevation eta(x, y, 0) with values ranging approximately from -0.6 to 0.2. The plot features a smooth surface with varying elevations, indicated by color gradients from blue to red. Panel B: A 3D surface plot showing a non-resonant solution with eta(x, y, T/4). The x-axis ranges from -π to π, and the y-axis ranges from -π to π. The z-axis represents the surface elevation eta(x, y, T/4) with values ranging approximately from -0.05 to 0.05. The plot features a more complex surface with multiple peaks and troughs, indicated by color gradients from blue to red. Panel C: A 3D surface plot showing a non-resonant solution with eta(x, y, T/2). The x-axis ranges from -π to π, and the y-axis ranges from -π to π. The z-axis represents the surface elevation eta(x, y, T/2) with values ranging approximately from -0.6 to 0.2. The plot features a smooth surface with varying elevations, indicated by color gradients from blue to red. Panel D: A 3D surface plot showing a resonant solution with eta(x, y, 0). The x-axis ranges from -π to π, and the y-axis ranges from -π to π. The z-axis represents the surface elevation eta(x, y, 0) with values ranging approximately from -0.6 to 0.2. The plot features a smooth surface with varying elevations, indicated by color gradients from blue to red.

Figure 14 shows typical case III standing waves for $k = 1,\ l = 0$ . (This does not mean that solutions are independent of $y$ according to the explanation in § 3.1.) The surface profiles of a non-resonant solution with $H = 0.3,\ \omega = 1.3979$ are shown in figures 14(a–c), corresponding to $t = 0,T/4,T/2$ . Note that this solution is essentially the same as the case II standing wave with $k=l=\sqrt {2}/2,\ H=0.3$ after a rotation of the $(x,y)$ -plane by $\pi /4$ . This is also reflected by the fact that the shift connecting $\eta (x,y,0)$ and $\eta (x,y,T/2)$ is along the diagonal direction, which becomes the $x$ - or $y$ -direction after the rotation. Figure 14(d) exhibits a resonant standing wave with $H=0.29,\ \omega = 1.3994$ . The resonance is primarily between the fundamental and ninth temporal harmonics, dominated by the $\cos (2x)\cos (5y)+\cos (5x)\cos (2y)$ modes.

Figure 15.

Typical level sets of the left-hand side of (4.4) for (a) $\alpha = \beta = 1$ , and (b) $\alpha = -1,\ \beta = 1$ . The green dots denote particular values of $\boldsymbol{\kappa}$ supporting resonant standing waves.

In previous studies of gravity-capillary waves, it is known that the three-wave resonance accounts for the classical Wilton ripples, travelling waves where the fundamental and second-order Fourier modes move at the same speed. (Other types of Wilton ripples exist, but correspond to high-order resonances.) To further investigate resonant standing waves, we consider the same mechanism here. Given a triad of wavenumbers $\boldsymbol{\kappa }$ , $\widetilde {\boldsymbol{\kappa }}$ and $\breve {\boldsymbol{\kappa }} = \boldsymbol{\kappa }+\widetilde {\boldsymbol{\kappa }}$ , the condition of three-wave resonance reads

(4.4) \begin{align} \omega (\breve {\boldsymbol{\kappa }}) - \omega (\boldsymbol{\kappa }) - \omega (\widetilde {\boldsymbol{\kappa }}) = 0, \end{align}

where $\omega (\boldsymbol{\kappa })$ is evaluated from (3.2) with $\kappa =|\boldsymbol{\kappa }|$ . Here, we let $\boldsymbol{\kappa } = (k,l)$ and $\widetilde {\boldsymbol{\kappa }} = (\alpha k, \beta l)$ , with $\alpha$ and $\beta$ being rational numbers. Figure 15(a) shows the left-hand side of (4.4) for $\alpha = \beta = 1$ . The zero level set is identical to a circle of radius $\sqrt {2}/2$ , and each point on it corresponds to a collinear resonant triad. To construct case II standing waves, we take

(4.5) \begin{align} \begin{cases} \boldsymbol{\kappa } = \widetilde {\boldsymbol{\kappa }} = (0.5,0.5), \quad \breve {\boldsymbol{\kappa }} = (1,1),\\[7pt] \omega (\boldsymbol{\kappa }) = \omega (\widetilde {\boldsymbol{\kappa }}) = \dfrac {\sqrt {3\sqrt {2}}}{2}\approx 1.0299, \quad \omega (\breve {\boldsymbol{\kappa }}) = \sqrt {3\sqrt {2}} \approx 2.0598, \end{cases} \end{align}

along with their transverse counterparts. These give rise to resonant standing waves featuring square-shaped patterns. Other possible wavenumbers can be found by calculating the intersections of the resonant circle and the straight line $l = rk$ , where $r$ is a rational number. To find non-collinear resonant triads, we let $\alpha = -1,\ \beta = 1$ . The left-hand side of (4.4) is shown in figure 15(b). In particular, we find the following resonant triad numerically (accurate up to the sixth decimal place):

(4.6) \begin{align} \begin{cases} \boldsymbol{\kappa } = (1,1.531164),\quad \widetilde {\boldsymbol{\kappa }} = (-1,1.531164), \quad \breve {\boldsymbol{\kappa }} = (0, 3.062328),\\[3pt] \omega (\boldsymbol{\kappa }) =\, \omega (\widetilde {\boldsymbol{\kappa }}) = 2.818705, \quad \omega (\breve {\boldsymbol{\kappa }}) = 5.637410. \end{cases} \end{align}

Figure 16.

Typical level sets of the left-hand side of (4.4) for $\alpha = 2,\ \beta = 4$ . The dashed line represents the zero level set of the the left-hand side of (4.7) for $\gamma = 2$ . The green dot represents a particular value of $\boldsymbol{\kappa}$ supporting resonant standing waves.

In general, (4.4) alone is not sufficient to support resonant standing waves. The other necessary condition is

(4.7) \begin{align} \omega (\widetilde {\boldsymbol{\kappa }}) - \gamma \omega (\boldsymbol{\kappa }) =0, \quad \gamma \in \mathbb Q, \end{align}

which means that the frequencies of the resonant triads must be rationally related. The previous two cases correspond to $ \gamma = 1$ . Figure 16 shows the left-hand side of (4.4) for $\alpha = 2,\ \beta = 4$ . The dashed line, which represents the zero level set of the left hand-side of (4.7) with $\gamma = 2$ , intersects the resonant curve at $k = 0.459181,\ l = 0.145664$ (accurate up to the sixth decimal place). Thus a particular choice supporting resonant standing waves is

(4.8) \begin{align} \begin{cases} \boldsymbol{\kappa } = (0.459181,0.145664),\quad \widetilde {\boldsymbol{\kappa }} = (0.918362,0.582656), \quad \breve {\boldsymbol{\kappa }} = (1.377543, 0.728320),\\[3pt] \omega (\boldsymbol{\kappa }) = 0.770405,\quad \omega (\widetilde {\boldsymbol{\kappa }}) = 1.540811, \quad \omega (\breve {\boldsymbol{\kappa }}) = 2.311216. \end{cases} \end{align}

Figure 17.

Two solution branches of case II standing waves for $k = l = 0.5$ : (a) $\omega$ versus $H$ ; (b) amplitude ratio versus $\epsilon$ .

Figure 18.

A case II standing wave with $k = l = 0.5,\ H=0.73$ on branch $1$ . (a–c) Top views of the solutions at $t=0,T/8,3T/16$ . A periodic cell is surrounded by the dashed lines. (d–f) The corresponding wave profiles.

We first examine the standing waves with $k=l=0.5$ . Figure 17 shows two fundamental solution branches (labelled $1$ and $2$ ) bifurcating from infinitesimal linear solutions at the frequency $\omega =1.0299$ . For branch $2$ , $\omega$ monotonically decreases with $H$ . For branch $1$ , $\omega$ initially grows for $H \lt 0.35$ , and subsequently decreases with $H$ . At small amplitude, $\eta$ can be written as

(4.9) \begin{align} \eta \sim \epsilon \cos (0.5 x)\cos (0.5 y)\cos (\omega t) + a\epsilon \cos (x)\cos (y)\cos (2\omega t) + O(\epsilon ^2), \end{align}

where $a$ measures the relative importance of the two resonant modes, thus influencing the wave morphology. Figure 17(b) shows that $a$ takes opposite signs on the two branches, and tends to $\pm 0.5$ as $\epsilon \rightarrow 0$ . Figure 18 displays a representative standing wave with $H=0.73\ (\epsilon =0.4361),\ \omega = 1.02800$ on branch $1$ . Figures 18(a–c) show the top views of the solutions at $t = 0,T/8,3T/16$ , exhibiting waffle-shaped patterns clearly. The corresponding wave profiles within a periodic cell (the region surrounded by the dashed lines) are shown in figures 18(d–f). The solution is calculated using $128\times 128\times 128$ grid points. Similarly, figure 19 shows the top views and wave profiles of a branch $2$ standing wave with $H=0.135\ (\epsilon =0.7095),\ \omega = 0.98008$ at $t = 0,T/8,3T/16$ . Figures 20(a) and 20(b) exhibit the cross-sections of the two standing waves over a quarter temporal period, resembling the classical Wilton ripples when the fundamental and second-order Fourier components are resonant.

Figure 19.

A case II standing wave with $k = l = 0.5,\ H=0.19$ on branch $2$ . (a–c) Top views of the solutions at $t=0,T/8,3T/16$ . A periodic cell is surrounded by the dashed lines. (d–f) The corresponding wave profiles.

Figure 19. Long description

The image contains six panels labeled (a) through (f), each depicting different aspects of water wave analysis. Panel A: A contour plot with x and y axes ranging from -4π to 4π. The plot shows a repeating pattern of contours in red and blue. Panel B: Another contour plot similar to Panel A but with a slightly different pattern. Panel C: A contour plot with a more intricate and dense pattern compared to Panels A and B. Panel D: A 3D surface plot with x and y axes ranging from -2π to 2π, and a z-axis representing η(x, y, 0). The surface shows peaks and troughs in red and blue. Panel E: A 3D surface plot similar to Panel D but with a different time parameter T/8. The surface shows a more complex pattern. Panel F: Another 3D surface plot with a time parameter 3T/16, showing a different wave profile with peaks and troughs.

Figure 20.

Cross-sections on the plane $y=0$ : (a) the branch $1$ solution with $H=0.73$ ; (b) the branch $2$ solution with $H=0.19$ . From bottom to top, $t = 0, T/32,T/16,\ldots,T/4$ .

Figure 21.

Bifurcations of case I standing waves for $k = 1,\ l = 1.5312$ . (a) Plot of $\omega$ versus $H$ . The blue and red curves correspond to three-dimensional solution branches, and the black curve represents the two-dimensional solution branch for $k = 0,\ l = 3.0624$ . (b) Plot of amplitude ratio versus $\epsilon$ . The blue and red dashed lines represent two fitted analytic relations for small $\epsilon$ .

Figure 21. Long description

Panel A: A line graph plots omega against H. The blue and red curves represent three-dimensional solution branches, while the black curve shows the two-dimensional solution branch for a specific parameter. An inset zooms in on a small section of the graph, highlighting detailed differences between the branches. Panel B: Another line graph plots the amplitude ratio against epsilon. The blue and red dashed lines indicate fitted analytic relations for small epsilon. The graph shows how the amplitude ratio changes with varying epsilon, with specific annotations marking key values.

Next, we examine the standing waves generated by non-collinear resonant triads. Figure 21 shows the two solution branches (also labelled $1$ and $2$ ) for the case I standing waves with $k = 1,\ l = 1.5312$ . Branch $1$ bifurcates from a two-dimensional standing wave with $k=0,\ l=3.0624$ at $H = 7\times 10^{-5},\ \omega = 2.81880$ , while branch $2$ bifurcates from an infinitesimal three-dimensional standing wave at the linear frequency $\omega = 2.81876$ . The small frequency gap between the two solution branches at $H=0$ is a consequence of the slight deviation at the fourth decimal place of $l$ from the exact value supporting resonance. The two frequency curves are similar to those shown in figure 17(a): one branch exhibits a monotonic decrease of frequency with wave amplitude, while the other displays a slight frequency upshift followed by a subsequent downshift. Following (4.9), we write $\eta$ as

(4.10) \begin{align} \eta \sim \epsilon \cos (x)\cos (1.5312 y)\cos (\omega t) + a\epsilon \cos (3.0624 y)\cos (2\omega t) + O(\epsilon ^2). \end{align}

Figure 21(b) shows that $a$ takes opposite signs on the two solution branches, and the two amplitude-ratio curves are approximately symmetric in the log scale. For $\epsilon \gt 10^{-3}$ , $|a|$ stabilises to a nearly constant level $\pm 0.324$ . Because of the different bifurcation mechanisms of the two branches, their amplitude ratios exhibit totally different behaviours when $\epsilon \lt 10^{-4}$ : for branch $2$ , $a \sim 7\times 10^{-5}/\epsilon$ , while for branch $1$ , $a \sim -0.324^2 \epsilon /{}(7\times 10^{-5})\approx -1500\epsilon$ , which is implied by the symmetry of the two curves.

Figure 22.

Profiles of small-amplitude case I standing waves for $k = 1$ and $l = 1.5312$ on branch $1$ : (a) $H = 8\times 10^{-5}$ , (b) $H = 10^{-3}$ .

Figure 22 shows the initial profiles of two small-amplitude branch $1$ solutions, with $H=8\times 10^{-5}$ in figure 22(a), and $H=0.001$ in figure 22(b). The former is close to the bifurcation point, thus featuring a nearly two-dimensional shape, while the latter develops a three-dimensional form consisting of both resonant modes. At larger wave amplitude, branch $1$ solutions exhibit similar wave profiles owing to the nearly constant value of $a$ . Figure 23 shows the solution with $H=0.4\ (\epsilon =0.3480),\ \omega = 2.7182$ on branch $1$ . Figures 23(a–c) correspond to the top views of the solutions at $t = 0,3T/16,T/4$ , which clearly show the alternative appearance of hexagonal and quasi-two-dimensional patterns. The three-dimensional views of the wave profiles within a periodic cell are plotted in figures 23(d–f). In contrast, figure 24 exhibits the top views and wave profiles of a branch $2$ standing wave with $H=0.135\ (\epsilon = 0.3402),\ \omega = 2.7222$ . Both solutions are calculated using $96\times 96\times 96$ points. Figures 25(a) and 25(b) show the cross-sections on the planes $y=0$ and $x=0$ for the branch $1$ solution over a quarter temporal period. The former resembles the non-resonant two-dimensional standing wave profiles, while the latter is similar to the resonant profiles shown in figure 20.

Figure 23.

A standing wave with $k = 1,\ l = 1.5312,\ H=0.4$ on branch $1$ . (a–c) Top views at $t=0,{}3T/16,T/4$ . (d–f) The corresponding profiles within the periodic cell surrounded by the dashed lines.

Figure 23. Long description

Panel A: A contour plot shows the top view of a standing wave with x and y axes ranging from -2π to 2π and -4.10 to 4.10 respectively. The plot features alternating red and blue contours indicating wave patterns. Panel B: Another contour plot similar to Panel A but with a different wave pattern. Panel C: A third contour plot with yet another distinct wave pattern. Panel D: A 3D surface plot displays the profile of the standing wave within a periodic cell. The x and y axes range from -π/2 to π/2 and -2.05 to 2.05 respectively, with the z-axis representing the wave height. The plot shows a peak in the center. Panel E: A 3D surface plot similar to Panel D but with a different profile and time frame. Panel F: Another 3D surface plot with a distinct profile and time frame.

Figure 24.

A standing wave with $k = 1,\ l = 1.5312,\ H=0.135$ on branch $2$ . (a–c) Top views at $t=0,{}3T/16,T/4$ . (d–f) The corresponding three-dimensional profiles within a periodic cell (regions surrounded by the dashed lines).

Figure 25.

(a) Cross-sections on the plane $y=0$ for the branch $1$ solution with $H=0.4$ . (b) Cross-sections on the plane $x=0$ for the same solution. From bottom to top, $t = 0, T/32,T/16,\ldots,T/4$ .

Figure 26.

Four bifurcations for standing waves with $k=0.459181,\ l=0.145664$ : (a) $\omega$ versus $H$ ; (b–e) amplitude ratio versus $\epsilon$ .

Figure 27.

Standing waves with $k=0.459181,\ l=0.145664$ at $t=0,T/8,T/4$ . (a–c) A branch $1$ solution with $H = 0.2$ . (d–f) A branch $2$ solution with $H = 0.07$ . (g–i) A branch $3$ solution with $H = 0.195$ . (j–l) A branch $4$ solution with $H = 0.185$ . For the convenience of visualisation, the figures are stretches in the $x$ - and $z$ -directions.

Figure 27. Long description

Panel A: A 3D surface plot shows standing water waves at initial time with axes labeled x, y, and η(x, y, 0). The plot features a wavy surface with red and blue regions indicating peaks and troughs, respectively. Panel B: A 3D surface plot displays standing water waves at T/8 with axes labeled x, y, and η(x, y, T/8). The plot shows a similar wavy pattern with red and blue regions. Panel C: A 3D surface plot illustrates standing water waves at T/4 with axes labeled x, y, and η(x, y, T/4). The plot features a more complex wavy surface with red and blue regions. Panel D: A 3D surface plot shows standing water waves at initial time with axes labeled x, y, and η(x, y, 0). The plot features a wavy surface with red and blue regions indicating peaks and troughs, respectively. Panel E: A 3D surface plot displays standing water waves at T/8 with axes labeled x, y, and η(x, y, T/8). The plot shows a similar wavy pattern with red and blue regions. Panel F: A 3D surface plot illustrates standing water waves at T/4 with axes labeled x, y, and η(x, y, T/4). The plot features a more complex wavy surface with red and blue regions. Panel G: A 3D surface plot shows standing water waves at initial time with axes labeled x, y, and η(x, y, 0). The plot features a wavy surface with red and blue regions indicating peaks and troughs, respectively. Panel H: A 3D surface plot displays standing water waves at T/8 with axes labeled x, y, and η(x, y, T/8). The plot shows a similar wavy pattern with red and blue regions. Panel I: A 3D surface plot illustrates standing water waves at T/4 with axes labeled x, y, and η(x, y, T/4). The plot features a more complex wavy surface with red and blue regions. Panel J: A 3D surface plot shows standing water waves at initial time with axes labeled x, y, and η(x, y, 0). The plot features a wavy surface with red and blue regions indicating peaks and troughs, respectively. Panel K: A 3D surface plot displays standing water waves at T/8 with axes labeled x, y, and η(x, y, T/8). The plot shows a similar wavy pattern with red and blue regions. Panel L: A 3D surface plot illustrates standing water waves at T/4 with axes labeled x, y, and η(x, y, T/4). The plot features a more complex wavy surface with red and blue regions.

Finally, we show the standing waves with $k=0.459181,\ l=0.145664$ . In the leading order, $\eta$ can be written as a combination of three components:

(4.11) \begin{align} \eta \sim {}& \epsilon \cos (0.459181 x)\cos (0.145664 y)\cos (\omega t)\nonumber \\ &{}+ a_1\epsilon \cos (0.918362 x)\cos (0.582656y)\cos (2\omega t)\nonumber \\ &{}+ a_2\epsilon \cos (1.377543 x)\cos (0.728320y)\cos (3\omega t) + O(\epsilon ^2). \end{align}

It is expected that there are four branches of solutions, corresponding to the combinations of the signs of $a_1$ and $a_2$ : $(+,+),(+,-),(-,+),(-,-)$ . Figure 26(a) shows the frequency curves of the four bifurcations, all stemming from the infinitesimal linear solutions at the frequency $\omega = 0.770405$ . Branches 1, 2 and 4 exhibit frequency downshift, while branch $3$ shows a frequency upshift. On the other hand, only three branches (2, 3, 4) exhibit resonant standing waves. This can be viewed in figures 26(b–e), where we plot the amplitude ratios versus $\epsilon$ along the four bifurcations. For branch $1$ , it is clear that $a_1\propto -\epsilon ^2$ and $a_2\propto -\epsilon ^3$ for $\epsilon \ll 1$ . Consequently, the small-amplitude standing waves on this branch are non-resonant. In contrast, the amplitude ratios take much greater values on the other three branches. When $\epsilon \ll 1$ , $a_1\approx -0.75,\ a_2\approx 0.6$ on branch $2$ , $a_1\approx 0.75, a_2\approx 0.6$ on branch $3$ , and $a_1\approx 0.75,\ a_2\approx -0.6$ on branch $4$ . No solution is found for the case $a_1\approx -0.75,\ a_2\approx -0.6$ when $\epsilon \ll 1$ . Based on our numerical experiments using various initial guesses, it is believed that there are only four solution branches growing from the bifurcation point, excluding the trivial case that a single periodic wave repeats itself multiple times within one computational domain. The profiles at $t=0,T/4,T/2$ of a non-resonant solution and three flower-like resonant solutions are shown in figure 27, where the four rows from top to bottom correspond to the branches $1$ to $4$ , respectively.