2.1. Derivation of the Benney–Luke equation

As our study is focused on the weakly nonlinear evolution of small-amplitude waves, we proceed to simplify (2.1) in the case $0 < \epsilon \ll 1$ and $h = O(1)$. We begin by expanding the dynamic and kinematic boundary conditions (2.1b)–(2.1c) about $z = 0$ in powers of $\epsilon$, which, upon eliminating $\eta$, gives rise to the equation (Benney Reference Benney1962; Milewski & Keller Reference Milewski and Keller1996)

(2.2)\begin{align} \phi_{tt} + \phi_z = \epsilon(\partial_t (\phi_{tz}\phi_t) + \phi_{zz}\phi_t - \partial_t(|\boldsymbol{\nabla} \phi|^2) - \phi_z\phi_{zt}) + O(\epsilon^2) \quad\mathrm{for}\ \boldsymbol{x} \in \mathcal{D}, \quad z = 0. \end{align}

To reduce the fluid evolution to the dynamics arising on the linearised free surface, $z = 0$, we define the Dirichlet-to-Neumann operator, $\mathscr {L}$, so that $\mathscr {L} \phi |_{z = 0} = \phi _z|_{z = 0}$. Here, $\phi$ satisfies Laplace's equation (2.1a) over the linearised domain $-h < z < 0$, with $\partial _z \phi = 0$ on $z = -h$ (see (2.1e)) and $\boldsymbol {n}\boldsymbol {\cdot } \boldsymbol {\nabla } \phi = 0$ for $\boldsymbol {x} \in \partial \mathcal {D}$ (see (2.1d)). Notably, the Dirichlet-to-Neumann operator may be defined in terms of its spectral representation, as detailed in § 2.2. By denoting $u(\boldsymbol {x},t) = \phi (\boldsymbol {x},0,t)$, we finally obtain the finite-depth Benney–Luke equation (Benney Reference Benney1962; Benney & Luke Reference Benney and Luke1964; Milewski & Keller Reference Milewski and Keller1996)

(2.3)\begin{equation} u_{tt} + \mathscr{L} u + \epsilon\left( u_t(\mathscr{L}^2 + \varDelta)u + \frac{\partial {}}{\partial {t}} [(\mathscr{L} u)^2 +|\boldsymbol{\nabla} u|^2]\right) = O(\epsilon^2) \quad\mathrm{for}\ \boldsymbol{x} \in \mathcal{D}, \end{equation}

where we have simplified the nonlinear terms in (2.2) using $\phi _{zz} = -\Delta \phi$ and $u_{tt} = -\mathscr {L} u + O(\epsilon )$.

The remainder of our investigation will be focused on the evolution of resonant triads governed by the Benney–Luke equation (2.3). As resonant triads arising in confined geometries are governed primarily by quadratic nonlinearities, it is sufficient to neglect terms of size $O(\epsilon ^2)$ in (2.3); however, higher-order corrections to the Benney–Luke equation may be derived by following a similar expansion procedure (Benney Reference Benney1962; Milewski & Keller Reference Milewski and Keller1996; Berger & Milewski Reference Berger and Milewski2003). Although our investigation is mainly focused on the evolution of the velocity potential, $u$, one may recover the leading-order free-surface elevation from the dynamic boundary condition (2.1b), namely $\eta = -u_t + O(\epsilon )$.

2.2. Spectral representation of the Dirichlet-to-Neumann operator

The Dirichlet-to-Neumann operator, $\mathscr {L}$, may be understood in terms of the discrete set of orthogonal eigenfunctions of the horizontal Laplacian operator (Kreyszig Reference Kreyszig1989). Specifically, we consider the set of real-valued eigenfunctions, $\varPhi _n(\boldsymbol {x})$, satisfying

(2.4)\begin{equation} -\Delta \varPhi_n = k_n^2 \varPhi_n \quad \mathrm{for}\ n = 0, 1, \ldots, \end{equation}

where the corresponding eigenvalues, $k_n^2$, are ordered so that $0 = k_0 < k_1\leq k_2 \leq \ldots$. Moreover, each eigenfunction satisfies the boundary condition $\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {\nabla } \varPhi _n = 0$ on $\partial \mathcal {D}$, as motivated by the no-flux condition (2.1d). Finally, the orthogonal eigenfunctions are normalised so that $\langle \varPhi _m, \varPhi _n \rangle = \delta _{mn}$, where

(2.5)\begin{equation} \langle \,f, g\rangle = \frac{1}{S}\iint_\mathcal{D} f g \,\mathrm{d}A \end{equation}

defines an inner product for real functions $f$ and $g$, $S$ is the area of $\mathcal {D}$ and $\delta _{mn}$ is the Kronecker delta. Notably, $\varPhi _0(\boldsymbol {x}) = 1$ is the constant eigenfunction, with corresponding eigenvalue $k_0 = 0$.

To determine the Dirichlet-to-Neumann operator for sufficiently smooth $\phi$, we first substitute the series expansion $\phi (\boldsymbol {x},z) = \sum _{n = 0}^\infty \phi _n(z)\varPhi _n(\boldsymbol {x})$ into Laplace's equation (2.1a), where we have temporally omitted the time dependence. We then solve the resulting equation for $\phi _n(z)$ over the linearised domain $-h < z < 0$, in conjunction with the no-flux condition on $z = -h$ (see (2.1e)). It follows that $\partial _z\phi _n(0) = \hat {\mathscr {L}}_n \phi _n(0)$, where

(2.6)\begin{equation} \hat{\mathscr{L}}_n = k_n \tanh(k_n h) \end{equation}

is the spectral multiplier of the Dirichlet-to-Neumann operator, $\mathscr {L}$. By expressing the time-dependent free-surface velocity potential, $u = \phi |_{z = 0}$, in terms of the basis expansion $u(\boldsymbol {x}, t ) = \sum _{n = 0}^\infty u_n(t) \varPhi _n(\boldsymbol {x})$, it follows that the Dirichlet-to-Neumann map has the spectral representation $\mathscr {L} u = \sum _{n = 0}^\infty \hat {\mathscr {L}}_n u_n \varPhi _n$.

3.1. The existence of a critical depth

We proceed to determine necessary and sufficient conditions on the wavenumbers, $K_j$, for there to exist a depth, $h_c$, at which a resonant triad forms, where such a critical depth is unique. We summarise our results in terms of the following theorem:

Theorem 1 There exists a positive and finite value of $h$ such that $\varOmega _1 + \varOmega _2 = \varOmega _3$ if and only if

(3.6)\begin{equation} K_1 + K_2 < K_3 < (\sqrt{K_1} + \sqrt{K_2})^2. \end{equation}

When this pair of inequalities is satisfied, the corresponding value of $h$ is unique.

We briefly sketch the proof of Theorem 1, with full details presented in Appendix A. We first demonstrate that no solutions to $\varOmega _1 + \varOmega _2 = \varOmega _3$ are possible when the bounds in (3.6) are violated, i.e. when $K_1 + K_2 \geq K_3$ or when $\sqrt {K_1} + \sqrt {K_2} \leq \sqrt {K_3}$. We then consider the case where the inequalities (3.6) are satisfied and determine the existence of positive roots to the function $F(h) = (\varOmega _1(h) +\varOmega _2(h))/\varOmega _3(h) - 1$. In this case, we demonstrate that $\lim _{h \rightarrow 0} F(h) < 0$ and $\lim _{h\rightarrow \infty } F(h) > 0$, from which we conclude that $F(h)$ has at least one root (by continuity of $F$). Finally, we deduce that this root is unique by proving that $F(h)$ is a strictly monotonically increasing function of $h$ when the inequalities (3.6) are satisfied.

Two important conclusions may be deduced from Theorem 1. First, it follows from (3.6) that the wavenumber, $K_3$, corresponding to the largest angular frequency, $\varOmega _3$, is larger than both the other two wavenumbers ($K_1$ and $K_2$), but it cannot be arbitrarily large (as supplied by the upper bound). For a given pair of eigenmodes (say $\varPsi _1$ and $\varPsi _2$), we conclude that there are likely to be only finitely many eigenmodes that can resonate with this pair (indeed, that number might fairly small, or even zero). Second, when modes 1 and 2 coincide (a 1:2 resonance), one deduces that $\varOmega _1 = \varOmega _2$ and $K_1 = K_2$; as such, the existence bounds (3.6) simplify to $2K_1 < K_3 < 4K_1$, or $2 < K_3/K_1 < 4$ (Mack Reference Mack1962; Miles Reference Miles1984b).

3.2. Determining the critical depth

Although Theorem 1 determines necessary and sufficient conditions on the wavenumbers, $K_j$, for there to be a critical depth, $h_c$, at which a resonant triad exists, the critical depth remains to be determined. In general, the critical depth must be computed numerically (being the unique root of the nonlinear function $F(h)$); however, we demonstrate that useful quantitative and qualitative information may be obtained via asymptotic analysis. For the remainder of this section, we consider the rescaled wavenumbers, $\xi _1 = K_1/K_3$ and $\xi _2 = K_2/K_3$, and the rescaled depth, $\zeta = K_3 h$; it remains to determine the root, $\zeta _c$, of

(3.7)\begin{equation} F(\zeta) = \sqrt{\frac{ \xi_1 \tanh(\xi_1 \zeta )} {\tanh(\zeta)} } + \sqrt{\frac{ \xi_2 \tanh(\xi_2 \zeta )} {\tanh(\zeta)} } - 1, \end{equation}

when $\xi _1, \xi _2 > 0$ satisfy

(3.8)\begin{equation} \xi_1 + \xi_2 < 1 < \sqrt{\xi_1} + \sqrt{\xi_2}. \end{equation}

In figure 2(a), we present contours of the critical rescaled depth, $\zeta _c$, in the $(\xi _1, \xi _2)$-plane, restricted to the region demarcated by (3.8). Consistent with the limits $\lim _{\zeta \rightarrow 0}F(\zeta ) = \xi _1 + \xi _2 - 1$ and $\lim _{\zeta \rightarrow \infty }F(\zeta ) = \sqrt {\xi _1} + \sqrt {\xi _2} - 1$, we observe that the root, $\zeta _c$, tends to zero at the line $\xi _1 + \xi _2 = 1$, and approaches infinity at the curve $\sqrt {\xi _1} + \sqrt {\xi _2} = 1$. Furthermore, the uniqueness of the root of $F$ for given $(\xi _1, \xi _2)$ is reflected in the observation that the contours of $\zeta _c$ do not cross. Finally, we note that the contours are symmetric about the line $\xi _1 = \xi _2$, which is a direct consequence of the invariance of $F(\zeta )$ under the mapping $\xi _1 \leftrightarrow \xi _2$ (see (3.7)).

Figure 2. Contours of the rescaled critical depth, $\zeta _c = h_cK_3$, as a function of the rescaled wavenumbers, $\xi _1 = K_1/K_3$ and $\xi _2 = K_2/K_3$. (a) The contours computed numerically from (3.7). The black lines indicate the limiting cases of $\zeta _c \rightarrow 0$ (at $\xi _1 + \xi _2 = 1$) and $\zeta _c \rightarrow \infty$ (at $\sqrt {\xi _1} + \sqrt {\xi _2} = 1$). (b) The contours are overlaid by the leading-order approximation ((3.15); circles) and the higher-order correction ((3.16); diamonds) for $\zeta _c$ equal to 0.5, 1, 1.5 and 2.

Although we are primarily interested in the physically relevant case for which the cylinder's depth-to-width ratio, $h$, is of size $O(1)$, an informative analytic result may be obtained by considering $F(\zeta )$ in the limit $\zeta \ll 1$ (or $K_3 h \ll 1$). By utilising the Taylor expansion

(3.9)\begin{equation} \sqrt{\tanh(x)} \sim \sqrt{x}\left(1 - \frac{x^2}{6} + \frac{19}{360}x^4 + O(x^6)\right), \end{equation}

we obtain

(3.10)\begin{align} & \sqrt{\xi_1\tanh(\xi_1\zeta)} + \sqrt{\xi_2\tanh(\xi_2\zeta)} -\sqrt{\tanh(\zeta)} \nonumber\\ &\quad \sim \sqrt{\zeta}\left[\left(\xi_1 + \xi_2 - 1\right) - \frac{\zeta^2}{6}(\xi_1^3 + \xi_2^3 - 1) + O(\zeta^4)\right] , \end{align}

for $0 < \zeta \ll 1$. Whilst deriving (3.10), we have utilised the bound $\xi _1, \xi _2 < 1$ (see (3.8)), which additionally ensures that $0 < \xi _j\zeta \ll 1$ for $j = 1,2$. We note that the left-hand side of (3.10) is equal to $F(\zeta )\tanh (\zeta )$, so $\zeta _c$ satisfies

(3.11)\begin{equation} \xi_1 + \xi_2 - 1 - \frac{\zeta_c^2}{6}(\xi_1^3 + \xi_2^3 - 1) = O(\zeta_c^4), \end{equation}

provided that $0 < \zeta _c \ll 1$. By neglecting terms of size $O(\zeta _c^4)$ in (3.11), one may then easily solve for $\zeta _c$ in terms of $\xi _1$ and $\xi _2$.

Alternatively, a more succinct expression for $\zeta _c$ may be found by first noting that

(3.12)\begin{equation} \xi_1^3 + \xi_2^3 = (\xi_1 + \xi_2)^3 - 3\xi_1\xi_2(\xi_1 + \xi_2) = 1 - 3\xi_1\xi_2 + O(\zeta_c^2), \end{equation}

where we have utilised the leading-order approximation $\xi _1 + \xi _2 = 1 + O(\zeta _c^2)$ (see (3.11)) to determine the second equality. Upon substituting (3.12) into (3.11), we find that $\xi _1$, $\xi _2$ and $\zeta _c$ are now related by the notably simpler expression

(3.13)\begin{equation} \xi_1 + \xi_2 - 1 + \frac{\zeta_c^2}{2}\xi_1\xi_2 = O(\zeta_c^4). \end{equation}

By neglecting terms of $O(\zeta _c^4)$, the leading-order approximation for the rescaled critical depth, $\zeta _c$, is given by

(3.14)\begin{equation} \zeta_c \sim \sqrt{\frac{2(1 - \xi_1 - \xi_2)}{\xi_1\xi_2}}, \end{equation}

an expression valid when $0 < \zeta _c \ll 1$ and $\xi _1 + \xi _2 < 1$ (see (3.8)). Alternatively, one may deduce from (3.13) that the contours of $\zeta _c$ satisfy the approximate form

(3.15)\begin{equation} \xi_2 \sim \frac{1 - \xi_1}{1 + \dfrac{1}{2}\zeta_c^2\xi_1}, \end{equation}

where the term in the denominator is responsible for the increased ‘bending’ of the contours as $\zeta _c$ becomes progressively larger (see figure 2). We note that the additional simplification afforded by (3.12) allows for a far more tractable representation of the contours relative to solving (3.11) directly for $\xi _2$ given $\xi _1$ and $\zeta _c$.

Despite being derived under the assumption $0 < \zeta _c \ll 1$, we see in figure 2(b) that the contours given by (3.15) agree favourably with the numerical solution even up to $\zeta _c \approx 1$. However, it is readily verified from (3.14) that the asymptotic approximation of each contour crosses the boundary curve $\sqrt {\xi _1} + \sqrt {\xi _2} = 1$ at $\zeta _c = 4$ (for which $\xi _1 = \xi _2 = \frac {1}{4}$), thereby demonstrating that the reduced asymptotic form has limited applicability (even in a qualitative sense) for slightly larger values of $\zeta _c$. One may further improve the quantitative (and, to an extent, qualitative) agreement between the asymptotic analysis and numerical computation by including terms of size $O(\zeta ^4)$ in (3.10); indeed, an analogous calculation gives rise to the following higher-order correction to (3.13):

(3.16)\begin{equation} \xi_1 + \xi_2 - 1 + \frac{\zeta_c^2}{2}\xi_1\xi_2 + \frac{\zeta_c^4}{72} \xi_1\xi_2\left(\xi_1\xi_2 - 1\right) = O(\zeta_c^6). \end{equation}

Although one may then solve for $\zeta _c$ given $\xi _1$ and $\xi _2$ (or, alternatively, determine the contours of $\zeta _c$) by truncating terms of $O(\zeta _c^6)$ in (3.16), the resulting algebraic expressions yield little qualitative information. However, one may, in principle, use this reduced form as a reasonable initial guess for a numerical root-finding algorithm for determining the root of $F(\zeta )$, provided that $\zeta _c$ is not too large.

3.3. Example cavities

Our investigation into the emergence of resonant triads has been focused, thus far, on finite-depth cylinders with arbitrary horizontal cross-section. However, it is convenient to understand how the results of Theorem 1 influence the formation (or not) of resonant triads for some specific cross-sections, namely rectangular, circular and annular cylinders.

3.3.1. Rectangular cylinder

It is well known that resonant triads are impossible for plane gravity waves evolving across an unbounded horizontal domain of finite depth (Phillips Reference Phillips1960; Hasselmann Reference Hasselmann1961) (weak interactions are possible, however, in the shallow-water limit, $K_jh\rightarrow 0$, for which $\tanh (K_jh)$ in the dispersion relation (3.1) is replaced by its leading-order approximation, $K_j h$ (Phillips Reference Phillips1960; Bryant Reference Bryant1973; Miles Reference Miles1976)). We now utilise Theorem 1 to demonstrate a similar result: resonant triads are impossible for gravity waves evolving within a rectangular cylinder of finite depth. Our result generalises the special case of a 1:2 resonance, for which the impossibility of internal resonance in a rectangular cylinder was demonstrated by Miles (Reference Miles1976).

To proceed, we consider a rectangular cylinder with side lengths $L_x$ and $L_y$. By orientating the Cartesian coordinate system, $\boldsymbol {x} = (x,y)$, so that the cylinder cross-section is defined by the region $0 < x < L_x$ and $0 < y < L_y$, the eigenmodes are of the form

(3.17)\begin{equation} \varPhi_{mn}(x,y) = \frac{1}{\mathscr{N}_{mn}}\cos(p_m x)\cos(q_n y), \end{equation}

where $\mathscr {N}_{mn} > 0$ is a normalisation constant. Notably, the wavenumbers $p_m = m{\rm \pi} /L_x$ and $q_n = n{\rm \pi} /L_y$ are chosen so that the no-flux condition is satisfied (see (2.1d)). For a triad determined by the non-negative integers $m_j$ and $n_j$ (for $j = 1, 2, 3$), the corresponding wavenumbers, $P_j = p_{m_j}$ and $Q_j = q_{n_j}$, must satisfy $P_1 + P_2 = P_3$ and $Q_1 + Q_2 = Q_3$ (under suitable reordering of the subscripts) in order for the eigenmode correlation condition (3.5) to be satisfied. By defining the wave vector $\boldsymbol {k}_j = (P_j, Q_j)$, the conditions on $P_j$ and $Q_j$ simplify to the single requirement $\boldsymbol {k}_1 + \boldsymbol {k}_2 = \boldsymbol {k}_3$, where the triangle inequality supplies that $|\boldsymbol {k}_3| \leq |\boldsymbol {k}_1| + |\boldsymbol {k}_2|$. As the eigenvalues, $K_j^2$, of the negative Laplacian operator are related to the wave vectors via $K_j = |\boldsymbol {k}_j|$, we deduce that $K_3 \leq K_1 + K_2$. Owing to the violation of the left-hand bound in (3.6), we conclude that resonant triads cannot exist in a rectangular cylinder of finite depth.

3.3.2. Circular cylinder

We consider a circular cylinder of unit radius in dimensionless variables (i.e. the dimensional radius is equal to $a$; see § 2). For polar coordinates $\boldsymbol {x} = (r,\theta )$, it is well known that the corresponding (complex-valued) eigenmodes may be expressed in the form

(3.18)\begin{equation} \varPhi_{mn}(r,\theta) = \frac{1}{\mathscr{N}_{mn}}\mathrm{J}_m(k_{mn}r)\mathrm{e}^{\mathrm{i} m \theta},\quad\mathrm{where}\ \mathscr{N}_{mn} = |\mathrm{J}_m(k_{mn})|\sqrt{1 - \frac{m^2}{k_{mn}^2}} \end{equation}

is the normalisation factor and $m$ is the azimuthal wavenumber (an integer). Furthermore, the no-flux condition (2.1d) determines that the radial wavenumbers, denoted $k_{mn}$, satisfy $\mathrm {J}_m'(k_{mn}) = 0$, where $0 < k_{m1} < k_{m2} < \cdots$ (we exclude $k_{00} = 0$ from consideration; see § 3). Notably, the eigenvalues of the negative Laplacian operator are precisely the squared wavenumbers, $k_{mn}^2$; consequently, the antinodes of each Bessel function play a pivotal role in determining the existence of resonant triads.

Akin to the rectangular cylinder, we find that the eigenmode correlation condition imparts an important restriction on the combination of eigenmodes that may resonate. For given $m_j$ and $n_j$ (for $j = 1,2,3$), we denote $K_j = k_{m_jn_j}$, $\varPsi _j = \varPhi _{m_j n_j}$ and $N_j = \mathscr {N}_{m_j n_j}$. Although the correlation condition given in (3.5) is defined for real eigenmodes, a similar condition holds for complex-valued eigenmodes, namely $\iint _{\mathcal {D}} \varPsi _1 \varPsi _2 \varPsi _3^*\,\mathrm {d}A \neq 0$. By considering the quantity

(3.19)\begin{align} \iint_{\mathcal{D}} \varPsi_1 \varPsi_2 \varPsi_3^* \,\mathrm{d}A &= \frac{1}{N_1 N_2 N_3}\left(\int_0^1 r\mathrm{J}_{m_1}(K_1 r) \mathrm{J}_{m_2}(K_2 r) \mathrm{J}_{m_3}(K_3 r)\,\mathrm{d}r\right) \nonumber\\ &\quad \times \left(\int_0^{2{\rm \pi}} \exp({\mathrm{i}(m_1 + m_2 - m_3)\theta})\,\mathrm{d}\theta\right), \end{align}

we deduce from the azimuthal integral that a necessary condition for the correlation integral to be non-zero is $m_1 + m_2 = m_3$ (Michel Reference Michel2019). This condition thus restricts the permissible combinations of azimuthal wavenumbers in a manner similar to the restriction on the permissible planar wavenumbers for the case of a rectangular cylinder. Unlike rectangular cylinders, however, we demonstrate that resonant triads are possible in a circular cylinder.

Despite the apparent restriction of the Bessel antinodes, $K_j$, and summation condition on the azimuthal wavenumbers, $m_j$, Theorem 1 determines that a vast array of resonant triads may be excited for judicious choices of the fluid depth. In table 1, we list a small number of resonant triads and each corresponding critical depth, $h_c$, subject to the restrictions $|m_j| \leq 3$ and $n_j \leq 3$; for larger values of $|m_j|$ and $n_j$, the corresponding wave field becomes increasingly oscillatory, to the extent that the effects of surface tension and dissipation might become appreciable. Moreover, even marginally relaxing the upper bounds on $|m_j|$ and $n_j$ vastly increases the number of resonant triads; indeed, the restriction $|m_j| \leq 4$ and $n_j \leq 4$ introduces 70 additional resonant triads relative to table 1. As the upper bounds for $|m_j|$ and $n_j$ are further increased, the typical difference between the various critical depths decreases and an increasingly large number of triads form at small values of the critical depth. Triads forming in shallow fluids (e.g. triads 21 to 25 in table 1) have physical relevance only at larger length scales (e.g. lakes) as dissipation could become a dominant factor at smaller scales.

Table 1. Combinations of the azimuthal wavenumbers, $m_j$, and radial mode indices, $n_j$, that form a resonant triad ($m_1 + m_2 = m_3$ and $\varOmega _1 + \varOmega _2 = \varOmega _3$) at critical depth, $h_c$, in a circular cylinder of unit radius. For each triad, the corresponding wavenumbers, $K_j = k_{m_j n_j}$, satisfy (3.6), and the correlation condition, $\iint _{\mathcal {D}} \varPsi _1\varPsi _2\varPsi _3^*\,\mathrm {d}A \neq 0$, is met. The list is restricted to resonant triads arising for $|m_j|, n_j \leq 3$, and we consider $m_1 \leq m_2$ and $m_3 \geq 0$ without loss of generality. We have omitted resonances that give rise to the same critical depth, but with the roles of modes 1 and 2 swapped. The triad numbers (left column) and shaded rows are referenced in the text.

Although the list of triads in table 1 is restricted to the lowest radial and azimuthal modes, we observe some general trends. In particular, we observe that the correlation integral, $\iint _{\mathcal {D}} \varPsi _1\varPsi _2\varPsi _3^*\,\mathrm {d}A$, generally decreases in magnitude as the fluid depth increases. Although the correlation integral remains non-zero (as is necessary to satisfy the correlation condition), its small value in some cases (e.g. triad 10) potentially corresponds to an elongation of the triad evolution time scale (see § 4.1). Moreover, we observe that the average wavenumber involved in the triad, $\bar {K} = \frac {1}{3}(K_1 + K_2 + K_3)$, is appreciably larger when the critical depth is very small. This correlation is consistent with the form of the corresponding angular frequency, $\varOmega _j = \sqrt {K_j \tanh (K_j h_c)}$, for which a small critical depth, $h_c$, is necessary for finite-depth effects to be appreciable when the typical wavenumber is large. To enumerate the myriad resonant triads arising in a circular cylinder when the upper bounds on $|m_j|$ and $n_j$ are relaxed, we provide MATLAB code in the supplementary material available at https://doi.org/10.1017/jfm.2023.441.

At this juncture, it is informative to assess how the triads listed in table 1 relate to the resonances explored in prior investigations. First, triad 8 in table 1 (dark grey row) was explored by Michel (Reference Michel2019) for a circular cylinder of radius 9.45 cm and an approximate fluid depth of 3 cm; it follows that the depth-to-radius ratio in Michel's experiment was approximately 0.317, close to the value of 0.30197 reported in table 1. Furthermore, table 1 (grey rows) incorporates two well-known examples of a 1:2 resonance, for which modes 1 and 2 coincide: (i) the critical depth $h_c = 0.83138$ (triad 2) corresponds to the second-harmonic resonance with the fundamental mode (Miles Reference Miles1976, Reference Miles1984a; Bryant Reference Bryant1989; Yang et al. Reference Yang, Dias, Liu and Liao2021); (ii) the critical depth $h_c = 0.19814$ (triad 14) corresponds to a standing wave composed of two resonant axisymmetric modes (Mack Reference Mack1962; Yang et al. Reference Yang, Dias, Liu and Liao2021). Finally, triads 1, 6, 16, 19 and 20 (table 1, light grey rows) form an interesting class of resonant triad, for which an axisymmetric mode ($m_3 = 0$) interacts with two identical counter-propagating non-axisymmetric modes ($m_1= -m_2 \neq 0$ and $n_1 = n_2$). In fact, our investigation in § 5.1 demonstrates that the axisymmetric mode is the so-called pump mode, and may thus excite the non-axisymmetric modes, even when the initial energy in each non-axisymmetric mode is negligible. We draw an analogy between this novel class of resonant triad and the excitation of beach edge waves (Guza & Davis Reference Guza and Davis1974) in § 6.

We conclude our exploration of resonant triads arising in a circular cylinder by remarking that the fluid depth may, in some cases, be judiciously chosen so as to excite multiple triads. In general, the condition on the angular frequencies, $\varOmega _1 + \varOmega _2 = \varOmega _3$ (see (3.4)), cannot be satisfied for two distinct triads at the same fluid depth; however, nonlinear resonance may persist for both triads provided that each condition on the angular frequencies is approximately satisfied (Bretherton Reference Bretherton1964; McGoldrick Reference McGoldrick1972; Craik Reference Craik1986), at the cost of weak detuning (see § 4.3.1 for further details). Specifically, if triads 1 and 2 have critical depths $h_{c,1}$ and $h_{c,2}$, respectively, then there is potential excitement of both triads when the fluid depth, $h$, satisfies $|h - h_{c,j}| = O(\epsilon )$ for $j = 1,2$ (where $0 < \epsilon \ll 1$ is the typical wave slope; see § 2), giving rise to the approximation $\varOmega _1 + \varOmega _2 - \varOmega _3 = O(\epsilon )$ for each triad. For example, if $0 < h_{c,2} - h_{c,1} \ll 1$, then it may be sufficient to excite both triads at an intermediate depth, $h_{c,1} \leq h \leq h_{c,2}$. We note, however, that the excitation of multiple triads at a single fluid depth is not possible when the depth discrepancy, $|h - h_{c,j}|$, becomes too large (relative to the typical wave slope) for any of the triads under consideration.

To demonstrate the potential for the simultaneous excitation of two triads within a circular cylinder of finite depth, we consider two scenarios: (i) the excitation of two triads that share a common wave mode; and (ii) the excitation of two triads that do not share any common wave modes. Heuristically, case (ii) is more common than case (i) owing to the number of similar fluid depths in table 1; however, case (i) will likely generate a far richer set of dynamics owing to the nonlinear interaction between the two triads (McEwan, Mander & Smith Reference McEwan, Mander and Smith1972; Craik Reference Craik1986; Chow, Henderson & Segur Reference Chow, Henderson and Segur1996; Choi, Chabane & Taklo Reference Choi, Chabane and Taklo2021). As an example of case (i), we consider triads 11 and 12 in table 1, with nearby critical depths $h_{c,1} = 0.23678$ and $h_{c,2} = 0.21395$, respectively. As mode $(m_3,n_3) = (3,3)$ is common to both triads, inter-triad resonance may arise at an intermediate depth, e.g. $h = 0.225$. Furthermore, an example of case (ii) arises for triads 13 and 14 in table 1, with nearby critical depths $h_{c,1} = 0.19839$ and $h_{c,2} = 0.19814$. Neither of these triads share a common wave mode, so one would not expect the inter-triad energy exchange discussed in case (i). Nevertheless, one might anticipate a signature of these two triads to be visible in the surface evolution for an intermediate depth, e.g. $h = 0.19825$. The theoretical and numerical exploration of coupled triads in a circular cylinder will be the focus of future investigation.

3.3.3. Annular cylinder

A natural variation upon a circular cylinder is an annulus of inner radius $r_0 \in (0,1)$ and outer radius 1. By varying $r_0$, the annulus approaches a circular cylinder as $r_0 \rightarrow 0^+$, and a quasi-one-dimensional periodic ring as $r_0 \rightarrow 1^-$. Notably, resonant triads are impossible for a one-dimensional periodic ring, as can be shown by modifying the arguments presented for the case of a rectangular cylinder (see § 3.3.1). Thus, one might anticipate that the existence of triads in an annular cylinder depends critically on the inner radius, $r_0$. Rather than enumerating some possible triads for given values of $r_0$, we instead track the corresponding critical depth, $h_c$, for the triads identified for a circular cylinder (see table 1) as $r_0$ is progressively increased from zero. Of particular interest is determining whether a given triad exists for all $r_0 < 1$, or whether there is some critical inner radius, $r_c$, beyond which the triad ceases to exist, with either $h_c \rightarrow 0$ or $h_c \rightarrow \infty$ as $r_0 \rightarrow r_c^-$.

The (complex-valued) eigenmodes in an annular domain are cylinder functions of the form

(3.20)\begin{equation} \varPhi_{mn}(r,\theta) = \frac{1}{\mathscr{N}_{mn}}[\mathrm{J}_m(k_{mn}r) \cos(\gamma_{mn}{\rm \pi}) + \mathrm{Y}_m(k_{mn}r) \sin(\gamma_{mn}{\rm \pi})] \mathrm{e}^{\mathrm{i} m \theta}, \end{equation}

where $\mathscr {N}_{mn} > 0$ is a normalisation constant, $\mathrm {Y}_m$ is the Bessel function of the second kind with order $m$ (an integer) and $\gamma _{mn} \in [0,1]$ determines the weighting between the two Bessel functions. As shown in Appendix B, the no-flux condition (see (2.1d)) on the inner and outer walls determines that the wavenumbers, $k_{mn}(r_0)$, satisfy the equation

(3.21)\begin{equation} \mathrm{J}_m'(k_{mn} r_0) \mathrm{Y}_m'(k_{mn}) - \mathrm{J}_m'(k_{mn}) \mathrm{Y}_m'(k_{mn} r_0) = 0. \end{equation}

A formula for the corresponding value of $\gamma _{mn}$ is determined in Appendix B. Once again, the wavenumbers, $k_{mn}$, are ordered so that $0 < k_{m1} < k_{m2} < \cdots$ (excluding $k_{00} = 0$) and satisfy $-\Delta \varPhi _{mn} = k_{mn}^2\varPhi _{mn}$. Three correlated wave modes may form a resonant triad (for a judicious choice of the fluid depth) provided that the corresponding wavenumbers, $K_j$, which depend on the channel width, $1 - r_0$, satisfy the bounds given in Theorem 1.

Bifurcating from the limiting case of a circular cylinder, we track the critical depth (when such a depth exists) of different triads as $r_0$ is progressively increased. The predominant behaviour is characterised by the example presented in figure 3, for which we consider the triad whose critical depth is $h_c = 0.17266$ as $r_0 \rightarrow 0^+$ (see triad 15 in table 1). Given that $h_c$ is fairly small in this limit, one might anticipate that the triad ceases to exist with $h_c \rightarrow 0$; somewhat surprisingly, however, the opposite scenario arises, with $h_c \rightarrow \infty$ as $r_0 \rightarrow r_c^-$ ($r_c \approx 0.57$ in this example). It follows, therefore, that the triad may persist for narrow channels only when the fluid is sufficiently deep. We note, however, that there exist (at least) two relatively rare transitions for increasing $r_0$, which we briefly describe as follows: (i) the triad ceases to exist when $h_c \rightarrow 0$ as $r_0 \rightarrow r_c^-$, which may arise when bifurcating from a sufficiently shallow circular cylinder (e.g. triad 25 in table 1); and (ii) the triad continues to exist for all $r_0 < 1$, with $h_c \rightarrow 0$ and $K_j \rightarrow \infty$ as $r_0 \rightarrow 1$, yet the normalised depth, $h_c K_3$, remains finite, and the normalised wavenumbers, $K_1/K_3$ and $K_2/K_3$, remain within the triad existence region (e.g. triad 14 in table 1). Owing to the appreciable influence of viscous effects for relatively shallow fluids, the physical relevance of these latter two scenarios is somewhat nebulous, however.

Figure 3. The existence and predominant characteristics of a triad in an annular cylinder with inner radius $r_0$ and outer radius 1. The triad bifurcates from the critical depth $h_c = 0.17266$ as $r_0 \rightarrow 0$ (the limiting case of a circular cylinder), with corresponding wavenumbers presented in table 1 (see triad 15). (a) The critical depth, $h_c$ (blue curve), with $h_c \rightarrow \infty$ as $r_0 \rightarrow r_c^-$, where $r_c \approx 0.57$ (black line). (b) The corresponding wavenumbers, $K_j$, all of which remain finite for $r_0 < r_c$ (black line). (c) The normalised wavenumbers, $K_1/K_3$ and $K_2/K_3$, parametrised by increasing $r_0$ (blue arrow), with the limiting case $r_0 \rightarrow 0$ denoted by the white dot. The wavenumbers leave the triad existence region (see Theorem 1) via the left-hand boundary (black curve) as $r_0 \rightarrow r_c^-$.

4.1. Multiple-scale analysis

In order to determine the evolution of each of the three dominant wave modes involved in an exact resonance, we utilise the method of multiple scales (Kevorkian & Cole Reference Kevorkian and Cole1996; Strogatz Reference Strogatz2015). Specifically, we seek a perturbation solution to the Benney–Luke equation (2.3) of the form $u \sim u_0 + \epsilon u_1 + O(\epsilon ^2)$. The leading-order terms in (2.3) determine that $u_0$ satisfies $\partial _{tt}u_0 + \mathscr {L} u_0 = 0$; we choose to consider a leading-order solution comprised only of the three triad modes (all other modes are assumed to be smaller in magnitude and appear at higher order), giving rise to the leading-order form

(4.3)\begin{equation} u_0(\boldsymbol{x} ,t, \tau) = \sum_{j = 1}^3 [ A_j(\tau) \varPsi_j(\boldsymbol{x}) \exp({-\mathrm{i} \varOmega_j t}) + \mathrm{c.c.}]. \end{equation}

In (4.3), we have introduced the slow time scale $\tau = \epsilon t$, which governs the evolution of each complex amplitude, $A_j$. As $\epsilon$ and $t$ are both independent variables, we treat $\tau$ and $t$ as independent time scales, giving rise to the transformation of derivatives $\partial _t \mapsto \partial _t + \epsilon \partial _\tau$. Finally, $\mathrm {c.c.}$ denotes the complex conjugate of the preceding term, a contribution necessary for real $u_0$.

So as to determine coupled evolution equations for each complex amplitude, $A_j$, we consider terms of $O(\epsilon )$ in the Benney–Luke equation (2.3). By substituting the leading-order solution, $u_0$, into the nonlinear terms and applying the triad condition for the angular frequencies (4.2), we obtain the following problem for $u_1$:

(4.4)\begin{equation} \partial_{tt}u_1 + \mathscr{L} u_1 ={-}\left[\sum_{j = 1}^3 f_j(\boldsymbol{x},\tau)\exp({-\mathrm{i} \varOmega_j t}) + \mathrm{c.c.}\right] + \mbox{non-resonant terms}. \end{equation}

As we will see below, each of the functions $f_j(\boldsymbol {x},\tau )$ appearing on the right-hand side of (4.4) will play a fundamental role when determining the amplitude equations; specifically,

(4.5) \begin{align} f_1 &={-}2\mathrm{i} \varOmega_1 \frac{\mathrm{d} {A_1}}{\mathrm{d} {\tau}} \varPsi_1 + \mathrm{i} A_2^* A_3^*([\varOmega_2 (L_3^2 - K_3^2) + \varOmega_3(L_2^2 - K_2^2) - 2\varOmega_1 L_2L_3]\varPsi_2\varPsi_3\nonumber\\ &\quad - 2\varOmega_1 \boldsymbol{\nabla}\varPsi_2 \boldsymbol{\cdot}\boldsymbol{\nabla} \varPsi_3), \end{align}

where $f_2$ and $f_3$ follow upon cyclic permutation of the subscripts $(1,2,3)$. Finally, we note that the ‘non-resonant terms’ in (4.4) are of the general form $p(\boldsymbol {x},\tau )\mathrm {e}^{\mathrm {i} \varsigma t}$, where we assume that the angular frequency, $\varsigma$, is not equal (or close) to any of the angular frequencies, $\pm \omega _n$, associated with linear wave modes (see § 3).

We proceed by projecting (4.4) onto each of the three wave modes, giving rise to differential equations of the form (for $j = 1,2,3$)

(4.6)\begin{equation} \partial_{tt}\hat{u}_{1,j} + L_j \hat{u}_{1,j} ={-}[\langle \varPsi_j, f_j \rangle \exp({-\mathrm{i} \varOmega_j t}) + \mathrm{c.c.}] + \mbox{non-resonant terms}, \end{equation}

where $\hat {u}_{1,j} = \langle \varPsi _j, u_1 \rangle$ is the projection of $u_1$ onto the mode $\varPsi _j$. By recalling that $L_j = \varOmega _j^2$, we immediately see that the term in square brackets in (4.6) is itself a solution to the linear operator $\partial _{tt} + \varOmega _j^2$. It follows that the solution of (4.6) comprises of particular solutions that have temporal dependence $t\exp ({\pm \mathrm {i} \varOmega _j t})$, leading to an ill-posed asymptotic expansion when $\epsilon t = O(1)$. The resolution to this problem is achieved via the solubility condition $\langle \varPsi _j, f_j\rangle = 0,$ which suppresses the secular growth.

By applying the solubility condition $\langle \varPsi _j, f_j\rangle = 0$ for $j = 1,2,3$, we conclude that the complex amplitude, $A_j(\tau )$, of each wave mode, $\varPsi _j(\boldsymbol {x}) \exp ({-\mathrm {i}\varOmega _j t})$, evolves according to the triad system of canonical form (Bretherton Reference Bretherton1964; Craik Reference Craik1986)

(4.7a–c)\begin{equation} \frac{\mathrm{d} {A_1}}{\mathrm{d} {\tau}} = \alpha_1 A_2^* A_3^*, \quad \frac{\mathrm{d} {A_2}}{\mathrm{d} {\tau}} = \alpha_2 A_1^* A_3^*, \quad \frac{\mathrm{d} {A_3}}{\mathrm{d} {\tau}} = \alpha_3 A_1^* A_2^*, \end{equation}

where

(4.8) \begin{equation} \alpha_1 = \frac{1}{2\varOmega_1}( [\varOmega_2 (L_3^2 - K_3^2) + \varOmega_3(L_2^2 - K_2^2) - 2\varOmega_1 L_2L_3] \mathscr{C} - 2\varOmega_1 \langle \varPsi_1, \boldsymbol{\nabla}\varPsi_2 \boldsymbol{\cdot}\boldsymbol{\nabla} \varPsi_3 \rangle), \end{equation}

while $\alpha _2$ and $\alpha _3$ follow by cyclic coefficient of the subscripts $(1,2,3)$. Furthermore, the correlation integral, $\mathscr {C}$, is defined

(4.9)\begin{equation} \mathscr{C} = \frac{1}{S}\iint_\mathcal{D} \varPsi_1\varPsi_2\varPsi_3\,\mathrm{d}A, \end{equation}

where we recall that $S$ is the area of the cylinder cross-section (see § 2.2). As the triad equations (4.7a–c) are valid for $\tau = O(1)$ (or $t = O(1/\epsilon )$), their dynamics yields an informative view of the long-time evolution of the resonant triad.

In order to assess the influence of the triad coefficients on the triad evolution (see § 4.2), we first simplify the algebraic form given in (4.8). As shown by Miles (Reference Miles1976), one may simplify the inner product $\langle \varPsi _1, \boldsymbol {\nabla } \varPsi _2 \boldsymbol {\cdot }\boldsymbol {\nabla } \varPsi _3 \rangle$ by repeated application of the divergence theorem and utilisation of the relationship $-\Delta \varPsi _j = K_j^2\varPsi _j$; it follows that

(4.10)\begin{equation} \langle \varPsi_1, \boldsymbol{\nabla}\varPsi_2 \boldsymbol{\cdot}\boldsymbol{\nabla} \varPsi_3 \rangle = \tfrac{1}{2}(K_2^2 + K_3^2 - K_1^2) \mathscr{C}, \end{equation}

where $\mathscr {C}$ is the correlation integral defined in (4.9). We then substitute (4.10) into (4.8) and simplify using the relation $\varOmega _1 + \varOmega _2 + \varOmega _3 = 0$. After some algebra, we derive the reduced expression

(4.11)\begin{equation} \alpha_1 = \frac{\mathscr{C}}{2\varOmega_1}\left( \varOmega_2 L_3^2 + \varOmega_3L_2^2 - 2\varOmega_1 L_2 L_3 + \sum_{l = 1}^3 \varOmega_l K_l^2 \right), \end{equation}

where $\alpha _2$ and $\alpha _3$ follow similarly. Finally, we demonstrate in Appendix C that the algebraic form of the triad coefficients may be further reduced to

(4.12)\begin{equation} \alpha_j = \frac{\mathscr{C}\beta}{2\varOmega_j} \quad\mathrm{for}\ j = 1,2,3, \end{equation}

where

(4.13)\begin{equation} \beta = \sum_{l= 1}^3 \varOmega_lK_l^2 - \frac{1}{2}\varOmega_1\varOmega_2\varOmega_3 (\varOmega_1^2 + \varOmega_2^2 + \varOmega_3^2). \end{equation}

Equations (4.7a–c), (4.9), (4.12) and (4.13) constitute the triad equations for resonant gravity waves confined to a cylinder of finite depth. Although the triad equations are of canonical form (Bretherton Reference Bretherton1964), the novelty of our investigation is the computation of the coefficients, $\alpha _j$, whose algebraic form is specific to our system.

4.2. Properties of the triad coefficients

The simplified form of the coefficients, $\alpha _j$ (4.12), allows for some important theoretical observations that were obfuscated by the more complicated expressions for $\alpha _j$ given in (4.8) and (4.11). In particular, as exactly two of the angular frequencies, $\varOmega _j$, have the same sign, we deduce from (4.12) that the two corresponding coefficients, $\alpha _j$, also have the same sign, with the third coefficient having the opposite sign. By utilising well-known results pertaining to the canonical triad equations, we conclude that all solutions to the triad equations (4.7a–c) are periodic in time, with solutions expressible in terms of elliptic functions (Ball Reference Ball1964; Bretherton Reference Bretherton1964; Simmons Reference Simmons1969; Craik Reference Craik1986). Typically, these solutions result in an exchange of energy between the comprising modes, although there is a class of periodic solution that, perhaps counter-intuitively, results in zero energy exchange for all time (Case & Chiu Reference Case and Chiu1977; Chabane & Choi Reference Chabane and Choi2019). Moreover, it is readily verified that the leading-order energy density, $\mathscr {E}_1 + \mathscr {E}_2 + \mathscr {E}_3$, is conserved, where $\mathscr {E}_j = \varOmega _j^2 |A_j|^2$, consistent with the Hamiltonian structure of the Euler equations (Bretherton Reference Bretherton1964; Craik Reference Craik1986). The reader is directed to the work of Craik (Reference Craik1986) for a more detailed account of the various properties of the canonical triad equations.

Of particular relevance to the evolution of the triad is the quantity $\beta$ (see (4.13)), which, together with $\mathscr {C}$, determines the time scale over which energy exchange arises. In particular, we present the form of $\beta$ in figure 4 for the case $\varOmega _1,\varOmega _2 > 0$ and $\varOmega _3 < 0$. As we will demonstrate below, $\beta < 0$ in this case; in general, the sign of $\beta$ is the same as the sign of the largest (in magnitude) angular frequency, $\varOmega _j$. Notably, $|\beta |$ decreases sharply towards zero as $K_1 + K_2 \rightarrow K_3$, corresponding to the limit $h_c \rightarrow 0$. Similarly, $|\beta |$ approaches zero in the limiting cases $K_1 \ll K_3$ or $K_2 \ll K_3$, corresponding to one low-oscillatory wave mode interacting with two highly oscillatory wave modes. Away from these limiting cases, however, $|\beta |$ depends only weakly on the wavenumbers, $K_j$, suggesting that the correlation integral, $\mathscr {C}$, predominantly controls the time scale of the triad evolution. Finally, we observe that $\beta$ is symmetric about the line $K_1 = K_2$, consistent with the invariance of (4.13) under the mapping $K_1\leftrightarrow K_2$ (and hence, $\varOmega _1 \leftrightarrow \varOmega _2$).

Figure 4. Contours of $\beta/K_3^{5/2}$ (see (4.13)) for the case $\varOmega _1, \varOmega _2 > 0$ and $\varOmega _3 < 0$ (with $\varOmega _1 + \varOmega _2 + \varOmega _3 = 0$), for which $\beta < 0$ (see (4.16)).

We conclude this section by proving that $\beta < 0$ in the case $\varOmega _1, \varOmega _2 > 0$ and $\varOmega _3 < 0$. By comparing the forms of (4.13) and (4.11), and then permuting the subscripts $(1,2,3) \mapsto (3,1,2)$, we first note that $\beta$ may be equivalently expressed as

(4.14)\begin{equation} \beta = \varOmega_1 L_2^2 + \varOmega_2L_1^2 - 2\varOmega_3L_1L_2 + \sum_{l = 1}^3\varOmega_l K_l^2, \end{equation}

or

(4.15)\begin{equation} \beta = \varOmega_1 (L_2^2 + K_1^2) + \varOmega_2(L_1^2 + K_2^2) + |\varOmega_3| ( 2L_1 L_2 - K_3^2). \end{equation}

By bounding $L_j = K_j \tanh (K_j h_c) < K_j$ for $0 < h_c < \infty$ and utilising the relation $\varOmega _1 + \varOmega _2 = |\varOmega _3|$, we obtain

(4.16) \begin{equation} \beta < |\varOmega_3|(K_1^2 + K_2^2 + 2K_1 K_2 - K_3^2) = |\varOmega_3| ((K_1 + K_2)^2 - K_3^2). \end{equation}

As resonant triads exist only when $K_1 + K_2 < K_3$ (see Theorem 1), we conclude that $\beta < 0$ in this case.

4.3. Summary

To summarise our theoretical developments, the velocity potential, $u$, at the fluid rest level ($z = 0$) evolves according to

(4.17)\begin{equation} u(\boldsymbol{x},t) \sim \sum_{j = 1}^3 [A_j(\tau) \varPsi_j(\boldsymbol{x}) \exp({-\mathrm{i} \varOmega_j t}) + \mathrm{c.c.}] + O(\epsilon), \end{equation}

where the complex amplitudes, $A_j(\tau )$, evolve over the slow time scale, $\tau = \epsilon t$, according to the triad equations (4.7a–c). In particular, the triad coefficients, $\alpha _j$ (see (4.12)), are defined in terms of the correlation integral, $\mathscr {C}$ (4.9), and the coefficient $\beta$ (4.13). Notably, we assume that $\mathscr {C}$ is non-zero; if this condition were violated then all three of the triad coefficients, $\alpha _j$, would be equal to zero, giving rise to non-interacting wave modes at leading order (contradicting the notion of a triad). Indeed, the condition $\mathscr {C} \neq 0$ is identical to the correlation condition detailed in (3.5), the origins of which we have now justified. Finally, the evolution of the free surface, $\eta$, may be recovered by recalling that $\eta = -u_t + O(\epsilon )$: we conclude that $\eta (\boldsymbol {x},t)$ has a similar leading-order form to $u(\boldsymbol {x},t)$, but each complex amplitude, $A_j(\tau )$, in (4.17) is replaced by $\mathrm {i} \varOmega _j A_j(\tau )$ (see (5.3) below).

We briefly contrast our investigation of triad interaction with the early-time calculation of Michel (Reference Michel2019), who characterised the initial linear growth of a child mode induced by the nonlinear interaction of two parent modes (where all three modes comprise the triad). If modes 1 and 2 are the parent modes and mode 3 is the child mode, then the initial linear growth may be deduced directly from triad equations (4.7a–c) in the limit $|A_3| \ll |A_1| \sim |A_2|$. Specifically, the initial variation of $A_1$ and $A_2$ is slow relative to that of $A_3$, which has the approximate early-time form $A_3(\tau ) \approx \alpha _3 C_1^* C_2^* \tau + C_3$, where $C_j = A_j(0)$. Notably, the linear growth rate of the child mode depends on the corresponding triad coefficient, $\alpha _3$, and the product of the initial amplitudes of the two parent modes. However, our result for circular cylinders differs to that of Michel; we believe that the author neglected some important nonlinear contributions (compare Michel's equation (A2) with (2.4) and (2.4a) of Longuet-Higgins Reference Longuet-Higgins1962). As Michel's experiment verified the scaling of the interaction only up to a proportionality constant, this discrepancy was not captured.

4.3.1. The influence of weak detuning

As discussed earlier in § 4, the analysis in §§ 4.1 and 4.2 does not account for weak detuning of the angular frequencies, as might arise when the fluid depth, $h$, differs slightly from the critical depth, $h_c$. We now briefly consider the case of weak detuning, for which (4.2) is replaced by the condition $\varOmega _1 + \varOmega _2 + \varOmega _3 = \epsilon \sigma$ (see § 3.3.2); here, $\epsilon$ is the small parameter representative of the typical wave slope (see § 2) and $\sigma = O(1)$ determines the extent of the detuning (Bretherton Reference Bretherton1964; McGoldrick Reference McGoldrick1972). By following a very similar multiple-scale procedure to the case $\sigma = 0$, we obtain amplitude equations that are now augmented by a time-dependent modulation. Specifically, each complex amplitude now evolves according to

(4.18a–c)\begin{equation} \frac{\mathrm{d} {A_1}}{\mathrm{d} {\tau}} = \alpha_1 A_2^* A_3^* \mathrm{e}^{\mathrm{i} \sigma \tau}, \quad \frac{\mathrm{d} {A_2}}{\mathrm{d} {\tau}} = \alpha_2 A_1^* A_3^* \mathrm{e}^{\mathrm{i} \sigma \tau}, \quad \frac{\mathrm{d} {A_3}}{\mathrm{d} {\tau}} = \alpha_3 A_1^* A_2^* \mathrm{e}^{\mathrm{i} \sigma \tau}, \end{equation}

where each coefficient, $\alpha _j$, is defined in (4.12). Although detuning yields non-autonomous amplitude equations, autonomous equations may be derived by mapping $A_j(\tau ) \mapsto A_j(\tau ) \exp ({\mathrm {i} \sigma \tau /3})$ for all $j = 1,2,3$ (Craik Reference Craik1986). Finally, we note that the energy, $\mathscr {E}_1 + \mathscr {E}_2 + \mathscr {E}_3$, is not exactly conserved when considering the effects of detuning; instead, the energy slowly oscillates about a constant value (Craik Reference Craik1986).

4.3.2. The case of a 1:2 resonance

A 1:2 resonance is a resonant triad for which two modes comprising the triad coincide. For this case, we define two angular frequencies, $\varOmega _1$ and $\varOmega _2$, so that $\varOmega _2 = 2\varOmega _1$ (Miles Reference Miles1976), where the connection to resonant triads is clear when writing $\varOmega _1 + \varOmega _1 = \varOmega _2$. By following a very similar multiple-scale procedure to that outlined in § 4.1, we obtain

(4.19)\begin{equation} u(\boldsymbol{x},t) \sim \sum_{j = 1}^2 [A_j(\tau) \varPsi_j(\boldsymbol{x}) \exp({-\mathrm{i} \varOmega_j t}) + \mathrm{c.c.}] + O(\epsilon), \end{equation}

where

(4.20a,b)\begin{equation} \frac{\mathrm{d} {A_1}}{\mathrm{d} {\tau}} ={-}\gamma A_1^* A_2\quad \mathrm{and}\quad \frac{\mathrm{d} {A_2}}{\mathrm{d} {\tau}} = \frac{\gamma}{4} A_1^2. \end{equation}

In particular, the evolution of the amplitude equations (4.20a,b) depends on the coefficient $\gamma = \mathscr {C}(K_2^2 - K_1^2 - 3\varOmega _1^4)$, where $\mathscr {C} = ({1}/{S})\iint _{\mathcal {D}} \varPsi _1^2 \varPsi _2 \,\mathrm {d}A$ is the correlation integral. Indeed, the amplitude equations (4.20a,b) and coefficient, $\gamma$, are consistent with the results of Miles (Reference Miles1976) when expressing the evolution of each complex amplitude, $A_j$, in polar form (with appropriate rescaling). Finally, we note that a weak detuning (see § 4.3.1) may also be incorporated within the amplitude equations (4.20a,b), thereby accounting for a slight mismatch between the fluid depth, $h$, and the corresponding critical depth, $h_c$ (Miles Reference Miles1976).

Of particular interest is the evolution of weakly nonlinear waves steadily propagating around a circular cylinder of unit radius, focusing on the case where the fluid depth is precisely equal to the critical depth of a 1:2 resonance (Yang et al. Reference Yang, Dias, Liu and Liao2021). For the complex-valued eigenmodes defined in (3.18), the correlation condition, $\iint _{\mathcal {D}} \varPsi _1^2 \varPsi _2^* \,\mathrm {d}A \neq 0$, determines that the angular wavenumbers satisfy $m_2 = 2m_1$ (Chossat & Dias Reference Chossat and Dias1995; Yang et al. Reference Yang, Dias, Liu and Liao2021). By expressing the complex wave amplitudes in polar form, $A_j(\tau ) = a_j(\tau ) \exp ({\mathrm {i} \theta _j(\tau )})$ (for $j = 1,2$), (4.20a,b) may be recast as (Miles Reference Miles1976)

(4.21a–c)\begin{align} \frac{\mathrm{d} {a_1}}{\mathrm{d} {\tau}} ={-}\gamma a_1 a_2 \cos\varTheta, \quad \frac{\mathrm{d} {a_2}}{\mathrm{d} {\tau}} = \frac{\gamma}{4}a_1^2 \cos\varTheta, \quad \frac{\mathrm{d} {\varTheta}}{\mathrm{d} {\tau}} = 2\gamma a_2 \left[1 - \frac{a_1^2}{8a_2^2}\right] \sin\varTheta, \end{align}

where $\varTheta (\tau ) = \theta _2(\tau ) - 2\theta _1(\tau )$ is the time-dependent phase shift. Steadily propagating waves correspond to time-independent solutions for $a_1$, $a_2$ (both non-zero) and $\varTheta$, from which we deduce that $\cos \varTheta = 0$ and $a_1/a_2 = 2\sqrt {2}$. Indeed, it is remarkable that the amplitude ratio of the two dominant (normalised) wave modes is independent of the angular wavenumbers, $m_j$, the radial wavenumbers, $K_j$, and the corresponding angular frequencies, $\varOmega _j$ (see § 3.3.2 for details). Furthermore, one may readily determine the relationship between the angular velocity of the steady wave rotation and the corresponding wave amplitude, which may then be compared with the numerical solution of the full Euler equations (Yang et al. Reference Yang, Dias, Liu and Liao2021). This comparison, as well as a comparison with steadily propagating waves computed from various truncations of the Euler equations, will be the subject of future investigation.

5.1. Excitation via the triad pump mode

To first identify the triad pump mode and then characterise the resultant excitation, we consider the case for which $A_3$, say, is much larger in magnitude than the other two mode amplitudes, so $|A_1|, |A_2| \ll |A_3|$ (Davis & Acrivos Reference Davis and Acrivos1967; Hasselmann Reference Hasselmann1967; Simmons Reference Simmons1969). By linearising the triad equations (4.7a–c), we obtain

(5.1a–c)\begin{equation} \frac{\mathrm{d} {A_1}}{\mathrm{d} {\tau}} = \alpha_1 A_2^* A_3^*, \quad \frac{\mathrm{d} {A_2}}{\mathrm{d} {\tau}} = \alpha_2 A_1^* A_3^*, \quad \frac{\mathrm{d} {A_3}}{\mathrm{d} {\tau}} = 0, \end{equation}

from which we immediately conclude that $A_3$ is constant (whilst the linearisation assumption holds); we denote $A_3(\tau ) = C$ for some given complex number $C$. By considering second derivatives of $A_1$ and $A_2$, we deduce the linearised evolution equations (Craik Reference Craik1986)

(5.2a,b)\begin{equation} \frac{\mathrm{d} {{}^2A_1}}{\mathrm{d} {\tau^2}} = \alpha_1\alpha_2 |C|^2 A_1\quad\mathrm{and}\quad \frac{\mathrm{d} {{}^2A_2}}{\mathrm{d} {\tau^2}} = \alpha_1\alpha_2 |C|^2 A_2, \end{equation}

where $\alpha _1\alpha _2 = \mathscr {C}^2\beta ^2/(4\varOmega _1\varOmega _2)$ (see (4.12)). We conclude that $A_1(\tau )$ and $A_2(\tau )$ grow exponentially in time (whilst the linearisation approximation holds) when $\varOmega _1 \varOmega _2 > 0$, and exhibit sinusoidal oscillations when $\varOmega _1\varOmega _2 < 0$ (Davis & Acrivos Reference Davis and Acrivos1967; Hasselmann Reference Hasselmann1967; Craik Reference Craik1986). Thus, mode 3 may excite modes 1 and 2 when $\varOmega _1$ and $\varOmega _2$ have the same sign (and likewise for other mode permutations). As one angular frequency must have a different sign from the other two (so as to satisfy $\varOmega _1 + \varOmega _2 + \varOmega _3 = 0$), we conclude that the mode whose angular frequency is largest in magnitude (i.e. differs in sign) is the triad pump mode (Craik Reference Craik1986). Equivalently, the pump mode is the mode with largest wavenumber, $K_j$, providing a robust mechanism for an inverse energy cascade to lower wavenumbers (Annenkov & Shrira Reference Annenkov and Shrira2006).

To visualise the influence of the pump mode on the resultant free-surface pattern, we present the solution of the triad equations (4.7a–c) and the corresponding pump-mode approximation (5.1a–c) in figure 5. By recalling that the free surface satisfies $\eta = -u_t + O(\epsilon )$, we first deduce that

(5.3)\begin{equation} \eta(\boldsymbol{x},t) \sim \sum_{j = 1}^3 [\mathrm{i} \varOmega_j A_j(\tau) \varPsi_j(\boldsymbol{x}) \exp({-\mathrm{i} \varOmega_j t}) + \mathrm{c.c.}] + O(\epsilon). \end{equation}

For the case of a circular cylinder, we utilise the complex-valued eigenmodes defined in (3.18), corresponding to the superposition of steadily propagating waves for $m_j \neq 0$ (the rotation direction depends on the sign of $\varOmega _j /m_j$). Upon initialising the system so that the energy is primarily within the pump mode (mode 3), modes 1 and 2 are gradually excited due to nonlinear interaction, with exponential growth evident for $\tau \lesssim 10$. As time further increases, the dynamics departs from the pump-mode approximation: the energy in the pump mode appreciably decreases, whilst the energy in modes 1 and 2 saturates. The free surface varies qualitatively during this evolution, with an appreciable change in pattern structure visible by $\tau = 24$ (primarily a superposition of modes 1 and 2). Notably, the system evolution is periodic, which becomes apparent over longer time scales.

Figure 5. Excitation of a triad via its pump mode for the case of a circular cylinder. We consider triad 11 in table 1, but with $m_3 \mapsto -m_3$. We choose $\varOmega _1, \varOmega _2 > 0$ and $\varOmega _3 < 0$, so that mode 3 is the pump mode. (a) Evolution of the free surface, $\eta \sim -u_t$, over the slow time scale, $\tau = \epsilon t$, with $\epsilon = 10^{-3}$. (b) The evolution of the wave amplitudes, $|A_j|$, according to the triad equations ((4.7a–c), solid curves) and the pump-mode approximation ((5.1a–c), dashed-dotted curves). Insets: modes 1 (blue), 2 (red) and 3 (gold) at $\tau = 0$; all three modes rotate counter-clockwise. The simulations were initialised from $A_1(0) = 0.01$ and $A_2(0) = 0.01\mathrm {i}$, where $A_3(0)$ was chosen to be the positive real number satisfying $\mathscr {E}_1 + \mathscr {E}_2 + \mathscr {E}_3 = 1$, with $\mathscr {E}_j = \varOmega _j^2 |A_j|^2$ (see § 4.2).

5.2. Excitation via an applied pressure source

Based on the ideas of the previous section, we consider a methodology for exciting the pump mode of a triad, which will subsequently excite the remaining two modes (provided that the initial disturbance of each of the remaining modes is non-zero). Notably, several methods for exciting internal resonances have been considered in prior investigations, primarily focusing on imposed motion of the fluid vessel via horizontal (Miles Reference Miles1976, Reference Miles1984c) or vertical vibration (Miles Reference Miles1976, Reference Miles1984b; Miles & Henderson Reference Miles and Henderson1990; Henderson & Miles Reference Henderson and Miles1991). Furthermore, one may, in principle, utilise sinusoidal paddles or plungers to excite a particular triad's pump mode for a given geometry (similar wave makers are used in rectangular wave tanks McGoldrick Reference McGoldrick1970a; Henderson & Hammack Reference Henderson and Hammack1987). However, for large-scale fluid tanks, imposed motion of the vessel may be impractical (if the tank were set in a concrete base, for example), and it may be challenging to determine the correct paddle motion necessary to excite a chosen pump mode for geometrically complex cylinders. We choose, therefore, to consider a slightly different approach: we instead excite the pump mode via a pulsating pressure source located just above the free surface (e.g. an air blower).

In order to incorporate a pressure source within our mathematical framework, we first reformulate the dimensionless dynamic boundary condition (2.1b) as

(5.4)\begin{equation} \phi_t + \eta + \frac{\epsilon}{2}(|\boldsymbol{\nabla} \phi|^2 + \phi_z^2) + \epsilon P(\boldsymbol{x}, t) = 0 \quad \mbox{for}\ \boldsymbol{x} \in \mathcal{D}, \quad z = \epsilon \eta, \end{equation}

where the dimensional pressure is $\epsilon ^2 a \rho g P$ for fluid density $\rho$ ($P = 0$ corresponds to atmospheric pressure). The pressure source is chosen to be small in magnitude so that the resultant wave excitation arises over the slow time scale, $\tau = \epsilon t$, and may thus be saturated by weakly nonlinear effects. By modifying the developments outlined in § 2.1, we derive the forced Benney–Luke equation

(5.5)\begin{equation} u_{tt} + \mathscr{L} u + \epsilon\left( u_t(\mathscr{L}^2 + \varDelta)u + \frac{\partial {}}{\partial {t}}[(\mathscr{L} u)^2 + |\boldsymbol{\nabla} u|^2] + \partial_t P\right) = O(\epsilon^2) \quad\mathrm{for}\ \boldsymbol{x} \in \mathcal{D}, \end{equation}

which will be the starting point for the asymptotic analysis.

Before proceeding further, we first describe two forms of the pressure source relevant to our investigation. For a stationary pressure source oscillating periodically over the fast time scale, $t$, we express $P(\boldsymbol {x},t) = f(\tau ) s(\boldsymbol {x}) \exp ({-\mathrm {i} \varOmega _p t}) + \mathrm {c.c.}$, where $s(\boldsymbol {x})$ is a fixed spatial profile (generally spanning the cavity), $f(\tau )$ accounts for a slow modulation in the magnitude of the pressure, and $\varOmega _p$ is the pulsation angular frequency. We choose $\varOmega _p$ to be close to the angular frequency of the pump mode, which, without loss of generality, we assume to be mode 3 (i.e. $\varOmega _3$ has the opposite sign from $\varOmega _1$ and $\varOmega _2$). We denote, therefore, $\varOmega _p = \varOmega _3 + \epsilon \mu$, where $\mu = O(1)$ determines the extent of the frequency mismatch. For a pressure source orbiting the centre of a circular cylinder at a constant angular velocity, we instead posit that $P$ has the form $P(r,\theta, t) = f(\tau ) s(r, \theta - \varOmega _p t)$, where $\varOmega _p = (\varOmega _3 + \epsilon \mu )/m_3$ is the angular velocity of the pressure source (assuming that the pump mode is non-axisymmetric, i.e. $m_3 \neq 0$).

For both standing and orbiting pressure sources, we now follow a similar multiple-scale procedure to that outlined in § 4.1, starting from the forced Benney–Luke equation (5.5). So as to discount the possibility that the pressure source excites more than one mode in the triad, we assume that neither $|\varOmega _1|$ or $|\varOmega _2|$ are close to $|\varOmega _3|$. Furthermore, we incorporate a weak detuning in the triad angular frequencies, denoting $\varOmega _1 + \varOmega _2 + \varOmega _3 = \epsilon \sigma$ (see § 4.3.1). It follows that each complex amplitude, $A_j(\tau )$, evolves according to

(5.6a–c)\begin{equation} \left.\begin{gathered} \frac{\mathrm{d} {A_1}}{\mathrm{d} {\tau}} = \alpha_1 A_2^* A_3^* \mathrm{e}^{\mathrm{i} \sigma \tau}, \quad \frac{\mathrm{d} {A_2}}{\mathrm{d} {\tau}} = \alpha_2 A_1^* A_3^* \mathrm{e}^{\mathrm{i} \sigma \tau}, \\ \frac{\mathrm{d} {A_3}}{\mathrm{d} {\tau}} = \alpha_3 A_1^* A_2^* \mathrm{e}^{\mathrm{i} \sigma \tau} - \varOmega_3s_3f(\tau)\exp({-\mathrm{i} \mu \tau}), \end{gathered}\right\} \end{equation}

where the coefficients, $\alpha _j$, are defined in (4.12). Notably, the pump mode may only be excited provided that the corresponding eigenmode is non-orthogonal to the pressure source, corresponding to a non-zero projection, i.e. $s_3\neq 0$, where $s_3 = \langle \varPsi _3, s\rangle$. Similar equations describing the evolution of forced resonant triads have been explored by McEwan et al. (Reference McEwan, Mander and Smith1972) (with the inclusion of linear damping) and Raupp & Silva Dias (Reference Raupp and Silva Dias2009).

In the special case of time-independent forcing ($\,f$ constant) and no frequency detuning ($\sigma = \mu = 0$), the dynamics of the forced triad equations has been analysed by Harris, Bustamante & Connaughton (Reference Harris, Bustamante and Connaughton2012), with both periodic and quasi-periodic dynamics reported. We also consider this case, leaving the effects of detuning and variable forcing for future investigation. In this setting, when $|A_1|$, $|A_2|$ and $|A_3|$ are initially small relative to the magnitude of the forcing, $|\varOmega _3 s_3 f|$, the initial growth in $A_3$ is approximately linear (see figure 6a). As mode 3 is the pump mode, the growth in $A_3$ excites $A_1$ and $A_2$, thus activating the triad. The conservation laws of the forced triad equations (Harris et al. Reference Harris, Bustamante and Connaughton2012) result in a temporary diminution of mode 3, which is later augmented by the external forcing; whence the process repeats. In some parameter regimes, the resulting evolution of the forced triad is periodic in time (see figure 6(b) and Raupp & Silva Dias Reference Raupp and Silva Dias2009); in contrast to the findings of Harris et al. (Reference Harris, Bustamante and Connaughton2012), however, we also identify initial conditions (with all other parameters unchanged) that result in hitherto unidentified chaotic dynamics (see figure 6c).

Figure 6. Evolution of the forced triad equations (5.6a–c) for $\sigma = \mu = 0$ and constant $f$. We consider the same triad as figure 5, with $s_3 f = 0.1$. In all four panels, $A_1(0) = 0.02\mathrm {i}$ and $A_2(0) = 0.01$. For $A_3(0) = 0.01$, we observe (a) the initial excitation of the triad and (b) the resultant periodic dynamics (the initial growth is highlighted within the grey box). (c) For $A_3(0) = 0.01\mathrm {i}$, the triad evolution appears chaotic. (d) The separation distance, $\delta (\tau )$, of a trajectory randomly perturbed at time $\tau _{{pert}}$, with $\delta (\tau _{{pert}}) = 10^{-10}$ and the same initialisation (at $\tau = 0$) as (c). The black curves bound $\delta (\tau )$ for 30 equally spaced values of $\tau _{{pert}}$ in the interval $100 \leq \tau _{{pert}} \leq 115$; the blue curve corresponds to $\tau _{{pert}} = 100$.

To verify the chaotic nature of this latter example, we consider the separation distance of two initially adjacent trajectories in phase space, for which we observe exponential divergence in time (see figure 6d). This exponential divergence is indicative of a positive maximal Lyapunov exponent (Strogatz Reference Strogatz2015), which characterises the sensitivity to initial conditions exhibited by chaotic systems. Specifically, for a solution, $A_j(\tau )$, and its perturbation, $A_{j,{pert}}(\tau )$, we consider the evolution of the separation distance, defined as

(5.7)\begin{equation} \delta(\tau) = \sqrt{\sum_{j = 1}^3 |A_j(\tau) - A_{j,{pert}}(\tau)|^2 }. \end{equation}

To compute $\delta (\tau )$, we first simulate the forced triad equations (5.6a–c) on the interval $0 \leq \tau \leq \tau _{{pert}}$, with $A_j(0) = A_{j,{pert}}(0)$ and $\tau _{{pert}}$ chosen to be sufficiently large so as to ensure that the chaotic attractor (should one exist) be approached. At $\tau = \tau _{{pert}}$, the real and imaginary parts of $A_{j,{pert}}$ are both randomly perturbed according to a uniform distribution on the interval $(-1, 1)$, with the resulting complex perturbation scaled so that $\delta (\tau _{{pert}}) = 10^{-10}$. We then evolve $A_j(\tau )$ and $A_{j,{pert}}(\tau )$ for $\tau \geq \tau _{{pert}}$, giving rise to initial exponential growth of $\delta (\tau )$, with saturation when the perturbation distance is comparable to the ‘diameter’ of the chaotic attractor (see figure 6d).

To confirm that the exponential growth was not specific to a perturbation about a particular point on the chaotic attractor (Strogatz Reference Strogatz2015), we considered 30 equally spaced values of $\tau _{{pert}}$ on the interval $100 \leq \tau _{{pert}} \leq 115$, roughly corresponding to the time taken for one ‘loop’ of the chaotic attractor to take place (see figure 6c). Each simulation was computed with a fourth-order Runge–Kutta method and a time step of 0.005. For each value of $\tau _{{pert}}$, we observed similar exponential divergence of trajectories (see figure 6d); moreover, the evolution of $\delta (t)$ during the growth phase remained unchanged when the time step was decreased to 0.001, with numerical errors only accumulating over longer time scales. Our results thus provide strong evidence that there is a positive maximal Lyapunov exponent in this particular portion of parameter space, indicative of a chaotic dynamics.