1. Introduction
Wave turbulence is a state of motion in systems characterised by nonlinear interactions among many waves at different scales. The long-term statistical properties of such systems are described theoretically in the framework of wave turbulence theory. In this framework, a kinetic equation can be derived to model the evolution of the wave spectrum under wave–wave interactions. Such interactions have to satisfy the resonance conditions (or approximations in the case of quasi-resonances)
\begin{align} \boldsymbol {k}_1\pm \boldsymbol {k}_2\pm \cdots \pm\boldsymbol {k}_N=0,\end{align}
\begin{align} \omega _1\pm \omega _2\pm \cdots \pm \omega _N=0,\end{align}
where
$N$
is the number of modes in the interaction (
$N=4$
for a quartet, and
$N=3$
for a triad), and
$\boldsymbol {k}_i$
is the wavenumber vector, which is related to frequency
$\omega _i$
by the dispersion relation.
For surface gravity waves in deep water, it is widely accepted that the relevant interactions in the kinetic equation are quartet resonant interactions with
$N=4$
. This is associated with the fact that triad resonant interactions are prohibited by the dispersion relation
$\omega \sim |\boldsymbol {k}|^{1/2}$
. As a result, a normal form transformation (see the detailed formulation in Krasitskii Reference Krasitskii1994) can be implemented on the dynamical equations to remove the quadraticnonlinearity terms (corresponding to non-resonant triad interactions). The resulted equation, named the Zakharov equation (Zakharov Reference Zakharov1968), contains only the cubic terms, based on which the kinetic equation with
$N=4$
can be further derived. This process implies that the effect of non-resonant triad interactions, instead of being null, can be merged into and understood at the level of quartet interactions (see other studies in the Fermi–Pasta–Ulam–Tsingou system, e.g. Ganapa (Reference Ganapa2023)). This argument about triad interactions, however, remains on a theoretical level without any direct numerical demonstration.
Another effect of triad interactions, which is discussed more often in the literature, is the generation of bound modes. While such an effect is absent from the Zakharov equation, it can indeed be expected from the original dynamical equation with quadratic nonlinearity. Given two free modes
$\boldsymbol {k}_1$
and
$\boldsymbol {k}_2$
, the triad non-resonant interaction is supposed to generate another mode
$\boldsymbol {k}_3$
with
$\omega _3$
not satisfying the dispersion relation, hence named a bound mode. Such bound modes can be observed in the wavenumber–frequency spectrum as branches away from the dispersion relation curves, which are reported and discussed in many experimental and numerical studies (Lvov, Nazarenko & Pokorni Reference Lvov, Nazarenko and Pokorni2006; Krogstad & Trulsen Reference Krogstad and Trulsen2010; Cobelli et al. Reference Cobelli, Przadka, Petitjeans, Lagubeau, Pagneux and Maurel2011; Taklo et al. Reference Taklo, Trulsen, Gramstad, Krogstad and Jensen2015; Aubourg et al. Reference Aubourg, Campagne, Peureux, Ardhuin, Sommeria, Viboud and Mordant2017; Campagne et al. Reference Campagne, Hassaini, Redor, Valran, Viboud, Sommeria and Mordant2019; Zhang & Pan Reference Zhang and Pan2022b
). While the bound mode has been proposed as a building block connecting two triads to form a resonant quartet (see figure 1 and the discussion in the previous paragraph), a single triad is considered as purely non-resonant, i.e. it is capable of generating only bound modes, thus not contributing to the free-wave portion of the spectrum.
Figure 1. A diagrammatic representation of a resonant quartet (
$\boldsymbol {k}_1, \boldsymbol {k}_2, \boldsymbol {k}_3, \boldsymbol {k}_4$
) formed by the connection of two non-resonant triads satisfying
$\boldsymbol {k}_1+\boldsymbol {k}_2=\boldsymbol {k}_c=\boldsymbol {k}_3+\boldsymbol {k}_4$
,
$\omega _1+\omega _2=\omega _c=\omega _3+\omega _4$
.
Despite the above-mentioned understanding on triad interactions, their effects on spectral evolution of gravity waves have not been completely understood, upon which questions often arise, say, in interpreting experimental data (Aubourg et al. Reference Aubourg, Campagne, Peureux, Ardhuin, Sommeria, Viboud and Mordant2017). In this work, we aim to provide a comprehensive study of this problem. The core of this study is a numerical algorithm that we developed to directly track the contributions of quadratic and cubic terms in the dynamical equation (hence triad and quartet interactions) to spectral evolution.
Through an analysis of the evolution of a typical gravity wave spectrum in a time interval, we find that the contributions from triad and quartet interactions share similar magnitude and trend for most wavenumbers, except that the former shows a much higher peak at small wavenumbers, where the initial energy level is low. We point out that the comparable contributions from triad and quartet interactions at most wavenumbers is a direct manifestation of the normal form transformation, i.e. the effect of triad interactions understood at the quartet level. We further study the rapid spectral evolution at low wavenumbers, which according to our results is due to non-resonant triad interactions. However, an analysis of the spectral content at low wavenumbers shows that there is significant energy in free waves in addition to bound waves, which seems to contradict the previous understanding of non-resonant interactions (generating only bound modes). We reconcile this contradiction by providing a new understanding: whenever a bound mode is generated, it is always accompanied by a free mode at the same wavenumber. In nature, this is similar to the scenario of a frictionless pendulum forced at non-natural frequency, which leads to an oscillation at both forced and natural frequencies (the latter being a homogeneous solution). Based on this understanding, we analytically derive the energy ratio between bound and free modes in non-resonant triad interactions through a perturbation analysis. The obtained formula is finally validated favourably through simulations of a range of different triad interactions.
The paper is organised as follows. In § 2, we first review the problem formulation and normal form transformation. Then we present our algorithm to decompose quadratic and cubic terms in the gravity-wave dynamical equation, and the method to track their contributions in spectral evolution. In § 3, we analyse the evolution from a (tail-damped) JONSWAP spectrum, and elucidate the effect of triad interactions in this process, regarding the connection to form quartets and non-resonant interactions to distribute energy in generated bound and free modes. In § 4, we derive the analytical formula for the distribution of energy in triad interactions, with validations from direct simulations. Discussions and conclusions are provided in § 5.
2. Formulation and methodology
We consider gravity waves on a two-dimensional (2-D) free surface of an incompressible and irrotational ideal fluid with infinite depth. The flow field can be described by a velocity potential
$\phi (\boldsymbol {x},z,t)$
satisfying Laplace’s equation, where
$\boldsymbol {x}=(x,y)$
denotes the horizontal coordinates,
$z$
is the vertical coordinate, and
$t$
is time. The evolution of the surface elevation
$\eta (\boldsymbol {x},t)$
and surface velocity potential
$\psi (\boldsymbol {x},t)=\phi (\boldsymbol {x},z,t)|_{z=\eta }$
satisfy the Euler equations in Zakharov form (Zakharov Reference Zakharov1968; see also Mei, Stiassnie & Yue Reference Mei, Stiassnie and Yue2005; Korotkevich Reference Korotkevich2023)
\begin{equation} \frac {\partial \eta }{\partial t}+{\boldsymbol \nabla}_{\boldsymbol {x}}\eta {\boldsymbol \cdot \nabla} _{\boldsymbol {x}}\psi -(1+{\boldsymbol \nabla} _{\boldsymbol {x}}\eta {\boldsymbol \cdot \nabla} _{\boldsymbol {x}}\eta )\phi _z=0, \end{equation}
\begin{equation} \frac {\partial \psi }{\partial t}+\eta +\frac {1}{2}{\boldsymbol \nabla} _{\boldsymbol {x}}\psi {\boldsymbol \cdot \nabla} _{\boldsymbol {x}}\psi -\frac {1}{2}(1+{\boldsymbol \nabla} _{\boldsymbol {x}}\eta {\boldsymbol \cdot \nabla} _{\boldsymbol {x}}\eta )\phi _z^2=0, \end{equation}
where
$\phi _z(\boldsymbol {x},t)=\partial \phi /\partial z|_{z=\eta }$
is the vertical velocity evaluated at the free surface, and
${\boldsymbol \nabla} _{\boldsymbol {x}}=(\partial /\partial x,\partial /\partial y)$
denotes the horizontal gradient. At each time instant,
$\phi _z$
can be determined as a solution to a Dirichlet-to-Neumann problem involving the boundary value
$\psi$
on
$z=\eta$
. With proper choice of mass and time units, the density and gravitational acceleration are both set to be unity so that they do not appear in (2.2).
2.1. Normal form transformation
Here we review the normal form transformation to remove the quadratic terms in (2.1) and (2.2). The canonical variable for gravity waves is defined as
\begin{equation} a_{\boldsymbol {k}}=\frac {1}{\sqrt {2}}\left(k^{-1/4} \hat {\eta }_{\boldsymbol {k}}+\mathrm {i}k^{1/4} \hat {\psi }_{\boldsymbol {k}}\right), \end{equation}
where a caret denotes the Fourier component in the wavenumber domain, and
$k=|\boldsymbol {k}|$
. Writing (2.1) and (2.2) in terms of
$a_{\boldsymbol {k}}$
, we have (truncated up to cubic nonlinearity)
\begin{equation} \begin{aligned} \mathrm {i}\frac {\partial a_{\boldsymbol {k}}}{\partial t}={}&\omega _{k}a_{\boldsymbol {k}} + \sum _{\boldsymbol {k}_1,\boldsymbol {k}_2} V^{[1]}_{012}a_1a_2\delta _{0-1-2} +\sum _{\boldsymbol {k}_1,\boldsymbol {k}_2} V^{[2]}_{012}a_1^*a_2\delta _{0+1-2} +\sum _{\boldsymbol {k}_1,\boldsymbol {k}_2} V^{[3]}_{012}a_1^*a_2^*\delta _{0+1+2} \\ &+\sum _{\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_3} W^{[1]}_{0123}a_1a_2a_3\delta _{0-1-2-3} +\sum _{\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_3} W^{[2]}_{0123}a_1^*a_2a_3\delta _{0+1-2-3} \\ &+\sum _{\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_3} W^{[3]}_{0123}a_1^*a_2^*a_3\delta _{0+1+2-3} +\sum _{\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_3} W^{[4]}_{0123}a_1^*a_2^*a_3^*\delta _{0+1+2+3}, \end{aligned} \end{equation}
where
$V^{[n]}$
,
$W^{[n]}$
are interaction coefficients, and
$\delta$
is the compact form of the Kronecker delta function satisfying
$\delta _{0\pm 1\pm 2}=\delta (\boldsymbol {k}\pm \boldsymbol {k}_1\pm \boldsymbol {k}_2)$
and
$\delta _{0\pm 1\pm 2\pm 3}=\delta (\boldsymbol {k}\pm \boldsymbol {k}_1\pm \boldsymbol {k}_2\pm \boldsymbol {k}_3)$
.
The goal is to introduce a near-identity transformation
$a_{\boldsymbol {k}}\rightarrow b_{\boldsymbol {k}}$
in the form
\begin{equation} \begin{aligned} a_{\boldsymbol {k}}={}&b_{\boldsymbol {k}}+\sum _{\boldsymbol {k}_1,\boldsymbol {k}_2} A^{[1]}_{012}b_1b_2\delta _{0-1-2} + \sum _{\boldsymbol {k}_1,\boldsymbol {k}_2} A^{[2]}_{012}b_1^*b_2\delta _{0+1-2} + \sum _{\boldsymbol {k}_1,\boldsymbol {k}_2} A^{[3]}_{012}b_1^*b_2^*\delta _{0+1+2} \\ &+ \sum _{\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_3} B^{[1]}_{0123}b_1b_2b_3\delta _{0-1-2-3} + \sum _{\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_3} B^{[2]}_{0123}b_1^*b_2b_3\delta _{0+1-2-3} \\ &+ \sum _{\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_3} B^{[3]}_{0123}b_1^*b_2^*b_3\delta _{0+1+2-3} + \sum _{\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_3} B^{[4]}_{0123}b_1^*b_2^*b_3^*\delta _{0+1+2+3}, \end{aligned} \end{equation}
so that when we write the equation in terms of
$b_{\boldsymbol {k}}$
, it yields the form
\begin{equation} \mathrm {i}\frac {\partial b_{\boldsymbol {k}}}{\partial t} = \omega _kb_{\boldsymbol {k}} + \sum _{\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_3} T_{0123}b_1^*b_2b_3\delta _{0+1-2-3}. \end{equation}
Equation (2.6) is known as the Zakharov equation, based on which the kinetic equation associated only with quartet interactions can be derived. In order to obtain the Zakharov equation, the coefficients
$A^{[n]}$
and
$B^{[n]}$
in (2.5) need to be chosen in a way such that the quadratic terms vanish under the transformation. The formulation and explicit expressions for these coefficients, together with those in (2.4) and (2.6), are well documented and can be found in e.g. Krasitskii (Reference Krasitskii1994) and Janssen & Onorato (Reference Janssen and Onorato2007). Here, we provide an example of
$A^{[n]}$
:
\begin{equation} A^{[2]}_{012}= -\frac {V^{[2]}_{012}}{\omega _k+\omega _1-\omega _2}\quad \text {for}\ \boldsymbol {k}_2=\boldsymbol {k}+\boldsymbol {k}_1, \end{equation}
where
\begin{equation} V^{[2]}_{012}=\frac {1}{2^{5/2}\pi }\left [\frac {k_2^{1/4}(\boldsymbol {k}\cdot \boldsymbol {k}_1+kk_1)}{(kk_1)^{1/4}} + \frac {k_1^{1/4}(\boldsymbol {k}\cdot \boldsymbol {k}_2-kk_2)}{(kk_2)^{1/4}} + \frac {k^{1/4}(\boldsymbol {k}_1\cdot \boldsymbol {k}_2-k_1k_2)}{(k_1k_2)^{1/4}} \right ]. \end{equation}
We remark that the denominator in (2.7) remains non-zero since there is no solution to (1.1) and (1.2) with
$N=3$
. This is the critical point allowing the normal form transformation to be implemented.
2.2. Algorithm to decompose contributions from quadratic and cubic terms
Our goal here is to compute numerically the contributions from quadratic and cubic terms in (2.4) to the evolution of modal energy
$e_{\boldsymbol {k}}=\omega _ka_{\boldsymbol {k}}a_{\boldsymbol {k}}^*$
. However, a direct computation based on (2.4) is very expensive with the convolutions, i.e.
$O(N^2)$
and
$O(N^3)$
computational complexity for quadratic and cubic terms, with
$N$
being the number of modes. To achieve an efficient
$O(N\log N)$
computation, we need to go back to (2.1) and (2.2) in the spectral domain where all terms can be computed by making use of the fast Fourier transform. For this purpose, we decompose (2.1) and (2.2) in the spectral domain as (the decomposed equations are equivalent to (2.4))
\begin{equation} \frac {\partial \hat {\eta }_{\boldsymbol {k}}}{\partial t}=\hat {L}_{\eta }+\hat {N}_{\eta }^{(2)}+\hat {N}_{\eta }^{(3)}+O(\epsilon ^4), \end{equation}
\begin{equation} \frac {\partial \hat {\psi }_{\boldsymbol {k}}}{\partial t}=\hat {L}_{\psi }+\hat {N}_{\psi }^{(2)}+\hat {N}_{\psi }^{(3)}+O(\epsilon ^4), \end{equation}
with linear terms
\begin{align} \hat {L}_{\eta }=\widehat {\phi _z^{(1)}},\end{align}
\begin{align} \hat {L}_{\psi }=-\hat {\eta }_{\boldsymbol {k}},\end{align}
and nonlinear terms (up to cubic nonlinearity)
\begin{equation} \hat {N}_{\eta }^{(2)}=\widehat {\phi _z^{(2)}}-\widehat {\big({\boldsymbol \nabla} _{\boldsymbol {x}}\eta {\boldsymbol \cdot \nabla} _{\boldsymbol {x}}\psi \big)}, \end{equation}
\begin{equation} \hat {N}_{\psi }^{(2)}=\frac {1}{2}\widehat {\big[\phi _z^{(1)} \big]^2}-\frac {1}{2}\widehat {\big({\boldsymbol \nabla} _{\boldsymbol {x}}\psi {\boldsymbol \cdot \nabla} _{\boldsymbol {x}}\psi \big)}, \end{equation}
\begin{equation} \hat {N}_{\eta }^{(3)}=\widehat {\phi _z^{(3)}}+\widehat {\left({\boldsymbol \nabla} _{\boldsymbol {x}}\eta {\boldsymbol \cdot \nabla} _{\boldsymbol {x}}\eta \phi _z^{(1)}\right)}, \end{equation}
\begin{equation} \hat {N}_{\psi }^{(3)}=\widehat {\phi _z^{(1)}\phi _z^{(2)}}. \end{equation}
In the above equations,
$\phi _z^{(m)}$
represents
$m$
th-order nonlinearity terms in
$\phi _z$
when it is expressed as power series in
$\eta$
and
$\psi$
, say, using the method in Dommermuth & Yue (Reference Dommermuth and Yue1987) and Pan, Liu & Yue (Reference Pan, Liu and Yue2018). With such expressions of
$\phi _z^{(m)}$
, all terms in (2.9) and (2.10) can be evaluated straightforwardly, employing the fast Fourier transform for the calculation of derivatives.
We next connect (2.9) and (2.10) to the spectral energy evolution rate
$r_{\boldsymbol {k}}=\partial e_{\boldsymbol {k}}/\partial t$
, the quantity of our interest in the study. In terms of
$\eta$
and
$\psi$
,
$e_{\boldsymbol {k}}$
can be written as
\begin{equation} e_{\boldsymbol {k}}=\frac {1}{2}(k\hat {\psi }_{\boldsymbol {k}}\hat {\psi }^*_{\boldsymbol {k}}+\hat {\eta }_{\boldsymbol {k}}\hat {\eta }^*_{\boldsymbol {k}}), \end{equation}
and
$r_{\boldsymbol {k}}$
can be written as
\begin{equation} r_{\boldsymbol {k}}=\frac {1}{2}\left [k\left (\frac {\partial \hat {\psi }_{\boldsymbol {k}}}{\partial t}\hat {\psi }_{\boldsymbol {k}}^*+\frac {\partial \hat {\psi }_{\boldsymbol {k}}^*}{\partial t}\hat {\psi }_{\boldsymbol {k}}\right ) + \left (\frac {\partial \hat {\eta }_{\boldsymbol {k}}}{\partial t}\hat {\eta }_{\boldsymbol {k}}^*+\frac {\partial \hat {\eta }_{\boldsymbol {k}}^*}{\partial t}\hat {\eta }_{\boldsymbol {k}}\right )\right ] \end{equation}
(see a similar implementation in the Majda–McLaughlin–Tabak model in Hrabski & Pan (Reference Hrabski and Pan2022)). By substituting (2.9) and (2.10) in (2.18), we obtain
\begin{equation} r^{(m)}_{\boldsymbol {k}}=\operatorname {Re}[k\hat {N}_{\psi }^{(m)}\hat {\psi }_{\boldsymbol {k}}^*+\hat {N}_{\eta }^{(m)}\hat {\eta }_{\boldsymbol {k}}^*],\quad m=2,3,\ldots, \end{equation}
where
$r^{(m)}_{\boldsymbol {k}}$
gives the
$k$
-mode energy evolution rate due to
$m$
th-order nonlinearity in (2.9) and (2.10). By definition,
$r^{(1)}_{\boldsymbol {k}}=0$
, while
$r^{(2)}_{\boldsymbol {k}}$
and
$r^{(3)}_{\boldsymbol {k}}$
correspond to the contributions from quadratic and cubic terms.
2.3. Numerical set-up
In this work, we simulate the evolution of a spectrum via (2.1) and (2.2), with (2.19) implemented to track the contributions of triad and quartet interactions to spectral evolution. The initial condition for the simulation is chosen as a directional tail-damped JONSWAP spectrum (Hasselmann et al. Reference Hasselmann1973). Specifically, the spectrum in the frequency–angle domain is written as
$S(\omega ,\theta )=G(\omega )\,D(\theta )$
, where
$\theta$
is the angle with respect to the positive
$x$
direction. We set
$G(\omega )=J(\omega )\,H(\omega)$
, where
$J(\omega )$
is a JONSWAP spectrum with the peak enhancement factor
$\gamma =6$
(peak period
$T_p$
and significant wave height
$H_s$
specified later), and
$H(\omega )$
is a tail-damped function in the form
\begin{equation} H(\omega )= \left \{ \begin{array}{lc} 1,& \quad\omega ^2\lt k_a,\\ {\rm e}^{\lambda (\omega ^2-k_a)},& \quad\omega ^2\geqslant k_a, \end{array} \right . \end{equation}
where
$\lambda$
and
$k_a$
are parameters controlling the damping rate and the effective range of damping, chosen as
$k_a=40$
and
$\lambda =-0.21$
in this study. The purpose of including
$H(\omega )$
is to allow more room for evolution at the tail of the spectrum (since the spectrum
$J(\omega )$
is relatively close to the stationary state at the tail). The angle spreading function
$D(\theta )$
is chosen as a cosine-squared function in the form
\begin{equation} D(\theta )= \left \{ \begin{array}{lc} \dfrac {2}{\pi }\cos ^2\theta ,& \quad |\theta |\leqslant \pi /2,\\ 0,& \quad |\theta |\gt \pi /2, \end{array} \right . \end{equation}
which satisfies
$\int_{-\pi }^{\pi }D(\theta )\,{\rm d}\theta =1$
. The initial phase of each mode is set as a random number uniformly distributed in
$[0,2\pi )$
, which rules out otherwise possible spectral adjustment due to phase randomisation in the initial stage of the simulation. We also note that while results in § 3 are based on the initial JONSWAP spectrum, the physical mechanisms of interest do not depend on the specific choice of the initial spectrum.
The simulation of (2.1) and (2.2) is conducted in a doubly periodic domain of size
$2\pi \times 2\pi$
(corresponding to a fundamental wavenumber
$k_0=1$
) with
$512\times 512$
modes. An order-consistent high-order spectral (HOS) method (West et al. Reference West, Brueckner, Janda, Milder and Milton1987) is used to compute the nonlinear terms up to cubic nonlinearity. For the time integration, we apply an integration factor scheme to solve for the linear terms analytically, and use a fourth-order Runge–Kutta method to integrate the nonlinear terms explicitly. (A detailed description of these time-marching schemes can be found in Pan (Reference Pan2020) in the context of capillary waves.)
To stabilise the simulation, we add dissipation terms
$D_{\eta }(k)=\gamma _k\eta$
and
$D_{\psi }(k)=\gamma _k\psi$
to the right-hand sides of (2.1) and (2.2), respectively, with the dissipation coefficients defined as
$\gamma _k=\gamma _0(k/k_d)^{\nu }$
. Such a dissipation scheme serves as a low-pass filter modelling the dissipation due to wave breaking (Xiao et al. Reference Xiao, Liu, Wu and Yue2013), which has been tested against wave tank experimental data and also applied in previous work (e.g. Zhang & Pan Reference Zhang and Pan2022b
). Here, we use
$\gamma _0=-50$
,
$k_d=115$
and
$\nu =30$
, and have tested that the resultant dissipation is negligible compared to the evolution due to nonlinear terms, so the evaluated
$r^{(2)}_{\boldsymbol {k}}$
and
$r^{(3)}_{\boldsymbol {k}}$
provide a correct view of the total evolution of the spectrum.
3. Results
We start by defining the angle-integrated modal energy and its evolution rate:
\begin{equation} E(k_r)=\sum _{|k-k_r|\lt \delta k}e_{\boldsymbol {k}}, \end{equation}
\begin{equation} R^{(m)}(k_r)=\sum _{|k-k_r|\lt \delta k}r_{\boldsymbol {k}}^{(m)}, \end{equation}
where
$\delta k$
, chosen as 1, is a parameter characterising the width of the (ring-shaped) region for the summation. Our interest is in the normalised and averaged energy evolution rate
$|\overline {R^{(m)}}|/|\overline {E}|$
obtained from a time interval of spectral evolution, with the overbar denoting the time average in this interval. The denominator
$|\overline {E}|$
accounts for the significant variation of energy level at different scales, which enables a fair comparison over all wavenumbers. We also remark that the quantity
$|\overline {R^{(m)}}|/|\overline {E}|$
has dimensions of inverse time
$[t]^{-1}$
and can be considered as an estimate of the reciprocal of the nonlinear time scale.
Figure 2. Spectra of surface elevation obtained at
$t=0$
(dashed orange line) and
$t=80T_p$
(solid black line) for (a)
$\epsilon =0.1259$
and (b)
$\epsilon =0.0629$
. Plots of
$|\overline {R^{(2)}}|/|\overline {E}|$
(red) and
$|\overline {R^{(3)}}|/|\overline {E}|$
(blue) as functions of
$k_r$
for (c)
$\epsilon =0.1259$
and (d)
$\epsilon =0.0629$
.
Starting from the tail-damped JONSWAP spectrum introduced in § 2.3, we perform simulations up to
$80T_p$
, where
$T_p$
is the peak wave period corresponding to the peak wavenumber
$k_p=20$
. This simulation time is found to be sufficient to capture the physics of interest in this work. The nonlinearity level
$\epsilon =H_sk_p/2$
is controlled by varying
$H_s$
in the initial spectrum. Figure 2 shows the evolution of the surface-elevation spectrum
$S_{\eta }(k)$
and the associated
$|\overline {R^{(m)}}|/|\overline {E}|$
(with
$m=2,3$
) computed over
$[0,80T_p]$
at two nonlinearity levels,
$\epsilon =0.1259$
and
$0.0629$
. From figures 2(a) and 2(b), we see significant spectral evolution occurring at small and large scales for both nonlinearity levels, with the evolution at higher nonlinearity being more pronounced (a phenomenon related to the finite size effect, which limits the spectral evolution at lower nonlinearity; see e.g. Pan & Yue Reference Pan and Yue2014; Zhang & Pan Reference Zhang and Pan2022b
). The quantities
$|\overline {R^{(m)}}|/|\overline {E}|$
shown in figures 2(c) and 2(d) reveal two observations that are true for both nonlinearities. (i) For most wavenumbers, especially closer to the tails of the spectrum,
$|\overline {R^{(2)}}|/|\overline {E}|$
and
$|\overline {R^{(3)}}|/|\overline {E}|$
show a comparable amplitude and similar trend. (ii) At small wavenumbers,
$|\overline {R^{(2)}}|/|\overline {E}|$
is dominant to drive the spectral evolution. These two observations actually hold for more nonlinearity levels that we have tested. To see this, we define another quantity
\begin{equation} Q=\frac {|\overline {R^{(2)}}|}{|\overline {R^{(2)}}|+|\overline {R^{(3)}}|}, \end{equation}
to measure the relative intensity of triad interactions. Figure 3 plots
$Q(k_r)$
as a function of
$\epsilon$
for
$k_r=3$
,
$61$
and
$99$
. It is clear that both facts can be observed, in terms of
$Q\approx 1$
for
$k_r=3$
for observation (i), and
$Q=0.3{-}0.6$
for
$k_r=61$
and
$99$
for observation (ii). We next discuss in detail the mechanisms associated with the two facts.
Figure 3. Plots of
$Q$
as functions of
$\epsilon$
at
$k_r=3$
(red),
$k_r=61$
(green) and
$k_r=99$
(blue).
3.1. Discussion for observation (i)
The observation that
$|\overline {R^{(2)}}|\sim |\overline {R^{(3)}}|$
for most wavenumbers is in fact a manifestation of the normal form transformation discussed in § 2.1. In particular, for any resonant quartet, there are two ways for energy transfer to occur. The first way is obviously from cubic terms in the governing (2.1) and (2.2). The second way is through quadratic terms, whose effect on spectral evolution is equivalent to the corresponding part in
$T_{0123}$
in (2.6) obtained from the normal form transformation. As a result, the effects of both quadratic and cubic terms can be seen on the quartet level, leading to comparable trend and magnitude between
$|\overline {R^{(2)}}|$
and
$|\overline {R^{(3)}}|$
in most wavenumbers.
The effect of the quadratic term on the quartet level can also be understood in a diagrammatic way. As shown in figure 1, for any resonant quartet consisting of modes
$\boldsymbol {k}_1$
,
$\boldsymbol {k}_2$
,
$\boldsymbol {k}_3$
and
$\boldsymbol {k}_4$
, one can always find a bridging mode
$\boldsymbol {k}_c=\boldsymbol {k}_1+\boldsymbol {k}_2=\boldsymbol {k}_3+\boldsymbol {k}_4$
, so that the resonant quartet is decomposed into two non-resonant triads
$(\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_c)$
and
$(\boldsymbol {k}_3,\boldsymbol {k}_4,\boldsymbol {k}_c)$
. Moreover, the bridging mode is associated with the same bound frequency in the two triads, simply because
$\omega _c=\omega _1+\omega _2=\omega _3+\omega _4$
. This ensures that the (frequency) oscillation of the bridging mode generated from the first triad can be perfectly transferred to the second triad, which is a critical condition for the transfer to be established. We note that this condition is satisfied only when
$\boldsymbol {k}_1$
,
$\boldsymbol {k}_2$
,
$\boldsymbol {k}_3$
and
$\boldsymbol {k}_4$
form a resonant quartet, which is the only situation (as opposed to non-resonant quartets) where such triad-induced transfer becomes effective.
Furthermore, as seen in figure 2, the relative contributions of
$|\overline {R^{(2)}}|$
and
$|\overline {R^{(3)}}|$
to the tail of the spectra remain similar at the two nonlinearity levels even though the finite size effect limits the spectral evolution at lower nonlinearity. This suggests that both cubic- and quadratic-term-induced transfers on the quartet level are subject to the interplay between nonlinear broadening and wavenumber discreteness. When the nonlinearity level is sufficiently low in a finite domain, we can expect that quasi-resonances corresponding to both cubic and quadratic terms are depleted, leading to a wave field where the nonlinear interactions are dominated by exact resonances (Kartashova, Nazarenko & Rudenko Reference Kartashova, Nazarenko and Rudenko2008; Zhang & Pan Reference Zhang and Pan2022a
).
3.2. Discussion for observation (ii)
We now turn to the spectral evolution at small wavenumbers induced by triad interactions dominantly. Although the evolutions seen in figures 2(a) and 2(b) look like an ‘inverse cascade’, the mechanism is completely different from the well-studied inverse cascade due to quartet resonant interactions, i.e. those resulting in the Kolmogorov-Zakharov solution (Annenkov & Shrira Reference Annenkov and Shrira2006; Korotkevich Reference Korotkevich2008) or spectral peak downshift (Annenkov & Shrira Reference Annenkov and Shrira2006). The time interval of our simulation,
$80T_p$
, is simply too short for these quartet-based processes to become relevant.
One may instead think that triad interactions resulting in the spectral evolution are quasi-resonant. Indeed, this is a question that often arises in the context of non-local triad interactions (Onorato et al. Reference Onorato, Osborne, Janssen and Resio2009; Korotkevich et al. Reference Korotkevich, Nazarenko, Pan and Shatah2024). In our case, a low-wavenumber mode
$\boldsymbol {k}_L$
can be excited by two high-wavenumber modes
$\boldsymbol {k}_1$
and
$\boldsymbol {k}_2$
(say
$\boldsymbol {k}_L+\boldsymbol {k}_1=\boldsymbol {k}_2$
), so the interaction can be non-local. Now consider the limit
$\boldsymbol {k}_L\rightarrow 0$
; then we have the frequency mismatch
$\Delta \omega \equiv \omega _L+\omega _1-\omega _2 \sim O(k_L^{1/2}) \rightarrow 0$
. Will this trigger energy transfer by triad quasi-resonances and invalidate the normal form transformation due to small/zero divisor in (2.7)? It turns out that this is not a problem since as
$k_L\rightarrow 0$
, the numerator (2.8) of the transformation coefficient satisfies
$V^{[2]}_{L12}\sim O(k_L^{3/4})\rightarrow 0$
faster than the denominator, so the full term approaches zero instead of blowing up. Physically, this means that as
$k_L\rightarrow 0$
, the triad interaction coefficient approaches zero sufficiently fast so that there is no energy transfer by quasi-resonances. We note that the situation is different in shallow-water gravity waves as discussed in Onorato et al. (Reference Onorato, Osborne, Janssen and Resio2009). One can show that in the case of shallow-water waves, the denominator is
$\Delta \omega \sim O(k_L)$
, and the numerator approaches zero following
$O(k_L^{1/2})$
, leading to a blow-up of the whole term. This means that for shallow water, the normal form transformation should be applied with more caution, and the quasi-resonant triad interactions can indeed be relevant.
The absence of triad quasi-resonances is consistent with figure 4, which plots the evolution of modal energy
$E(k_r)$
with
$k_r=3$
at three nonlinearity levels. We see that the modal energy rises from 0 to a stationary value in a very short linear time scale
$O(T_p)$
, which is a feature of non-resonant interactions instead of quasi-resonant interactions. Does this mean that the energy level at low wavenumbers is sustained by bound modes (according to the reasoning for non-resonant interactions)? To answer this question, we examine the wavenumber–frequency spectrum
$S_{\eta }(k,\omega )$
, defined as
\begin{equation} S_\eta (k,\omega )=\int _0^{2\pi } |\hat {\eta }(\boldsymbol {k},\omega )|^2 k\, {\rm d}\theta , \end{equation}
where
$\hat {\eta }(\boldsymbol {k},\omega )$
is the spatiotemporal Fourier transform of
$\eta (\boldsymbol {x},t)$
:
\begin{equation} \hat {\eta }(\boldsymbol {k},\omega )=\frac {1}{4\pi ^2T_L}\iiint _{[0,T_L]\times [0,2\pi ]\times [0,2\pi ]} \eta (\boldsymbol {x},t)\,h_T(t)\,{\rm e}^{-\mathrm {i}(\boldsymbol {k}\boldsymbol{\cdot} \boldsymbol {x}-\omega t)}\,{\rm d}\boldsymbol {x}\,{\rm d}t, \end{equation}
with
$h_T(t)$
the Tukey window (Bloomfield Reference Bloomfield2004) of length
$T_L$
. Figure 5(a) shows
$S_{\eta }(k,\omega )$
with a zoom-in view of the low-wavenumber region. A more specific view is also provided in figure 5(b), which plots
$S_{\eta }(k=3,\omega )$
as a function of
$\omega$
, together with a vertical line denoting the free-mode frequency
$\omega =\sqrt {3}$
. Somewhat to our surprise, significant energy is located at the free mode, which is at least comparable to those away as bound modes. How do we reconcile the free-mode generation and non-resonant interactions? In § 4, we will show explicitly that non-resonant interactions indeed generate both bound and free modes, with the distribution between these two components solved analytically with validations.
Figure 4. Time evolution of
$E(k_r,t)$
with
$\epsilon =0.0629$
(blue),
$\epsilon =0.0944$
(green) and
$\epsilon =0.1259$
(red) at
$k_r=3$
.
Figure 5. (a) Wavenumber–frequency spectrum
$S_\eta (k,\omega )$
in log scale with the dispersion relation marked by a red line, obtained with
$\epsilon =0.0315$
and
$T_L=40T_p$
. A zoom-in view at small wavenumbers is shown to illustrate the portion of interest. (b) Plot of
$S_\eta (k=3,\omega )$
as a function of
$\omega$
(blue), with the corresponding free-mode frequency
$\omega _f(k=3)$
marked by a black dashed line.
4. Analytical study on triad non-resonant interactions with validations
We consider a simple model of a wave field consisting of two modes as the initial condition. In physical space, the initial field is described by
\begin{align} \eta (\boldsymbol {x})=\eta _1(\boldsymbol {x})+\eta _2(\boldsymbol {x}),\end{align}
\begin{align} \psi (\boldsymbol {x})=\psi _1(\boldsymbol {x})+\psi _2(\boldsymbol {x}),\end{align}
with
\begin{align} \eta _i(\boldsymbol {x})=\tilde {A}_i\exp {(\mathrm {i}\boldsymbol {k}_i\cdot \boldsymbol {x})}+\tilde {A}_i^*\exp {(-\mathrm {i}\boldsymbol {k}_i\cdot \boldsymbol {x})},\end{align}
\begin{align} \psi _i(\boldsymbol {x})=-\frac {\mathrm {i}}{\omega _i} [\tilde {A}_i\exp {(\mathrm {i}\boldsymbol {k}_i\cdot \boldsymbol {x})}-\tilde {A}_i^*\exp {(-\mathrm {i}\boldsymbol {k}_i\cdot \boldsymbol {x})]},\end{align}
where
$\tilde {A}_i\in \mathbb {C}$
,
$\boldsymbol {k}_i=(k_{xi},k_{yi})\in \mathbb {R}^2$
and
$\omega _i=k_i^{1/2}$
, with
$k_i=|\boldsymbol {k}_i|$
for
$i=1,2$
. This model allows two triad non-resonant interactions to occur, generating new modes at wavenumbers
\begin{align} \boldsymbol {k}_{1+2}=\boldsymbol {k}_1+\boldsymbol {k}_2,\end{align}
\begin{align} \boldsymbol {k}_{1-2}=\boldsymbol {k}_1-\boldsymbol {k}_2,\end{align}
respectively. Our goal is to show and elucidate the energy distribution between free and bound modes in this simpler problem. With the solution to this problem obtained, the spectral behaviour at low wavenumbers observed in figures 4 and 5 can be simply understood as the result of a collection of such non-resonant interactions.
We first provide an illustrative HOS simulation by setting
$\boldsymbol {k}_1=(25,0)$
,
$\boldsymbol {k}_2=(20,0)$
,
$|\tilde {A}_1|=4\times 10^{-4}$
,
$|\tilde {A}_2|=5\times 10^{-4}$
, with initial phases
$\arg {(\tilde {A}_i)}$
assigned as random values within
$[0,2\pi )$
. The simulation is conducted up to quadratic nonlinearity to avoid potential higher-order effects. Our interest is in the wavenumbers
$\boldsymbol {k}_{1+2}=(45,0)$
and
$\boldsymbol {k}_{1-2}=(5,0)$
. Figure 6 plots
$S_\eta (k,\omega )$
(computed over
$[0,30T_0]$
, where
$T_0=2\pi k_0^{-1/2}$
is the fundamental wave period) as a function of
$\omega$
for
$k=k_{1+2}=45$
(figure 6
a) and
$k=k_{1-2}=5$
(figure 6
b). We see that for each wavenumber of interest, there are two peaks located at the free-wave frequency (
$\omega _f=\omega _{1+2}=6.71$
for
$k_{1+2}$
, and
$\omega _f=\omega _{1-2}=2.24$
for
$k_{1-2}$
) and bound-wave frequency (
$\omega _b=\omega _1+\omega _2=9.47$
for
$k_{1+2}$
, and
$\omega _b=\omega _1-\omega _2=0.53$
for
$k_{1-2}$
). This is clear evidence that non-resonant triad interactions can distribute comparable amounts of energy between the free and bound modes. Defining
$S_b=S_{\eta }(k_{1\pm 2},\omega _b)$
and
$S_f=S_{\eta }(k_{1\pm 2},\omega _f)$
, we can calculate that
$S_b/S_f=0.63$
and
$0.1$
for
$k_{1+2}$
and
$k_{1-2}$
in this case.
Figure 6. Plots of
$S_\eta (k,\omega )$
as functions of
$\omega$
at (a)
$k=k_{1+2}=45$
and (b)
$k=k_{1-2}=5$
, with
$T_L=30T_0$
. The corresponding free-mode frequencies
$\omega _f$
are marked by black dashed lines, and bound-mode frequencies
$\omega _b$
are marked by red dashed lines. Insets: same plots in semi-log scale.
We next develop an analytical solution to describe quantitatively the distribution of energy, i.e. to obtain the value of
$S_b/S_f$
, for a general configuration of
$\boldsymbol {k}_1$
and
$\boldsymbol {k}_2$
. For this purpose, we perform a perturbation analysis based on (2.1) and (2.2). As detailed in Appendix A, we set the first-order linear solution to match the linear wave field corresponding to the initial conditions (4.1) and (4.2). Then we seek a second-order solution to describe the energy generation at free mode
$S_f$
and bound mode
$S_b$
. A critical condition applied at second order is a quiescent wave field for wavenumbers
$\boldsymbol {k}_{1+2}$
and
$\boldsymbol {k}_{1-2}$
at
$t=0$
, which allows the calculation of distribution of energy between the two frequencies at each wavenumber. The final solution yields the ratios between
$S_b$
and
$S_f$
at
$\boldsymbol {k}_{1+2}$
and
$\boldsymbol {k}_{1-2}$
respectively:
\begin{align} (S_b/S_f)_{1+2}=\frac {|\Gamma _1+\Gamma _2|^2}{|\Gamma _1|^2+|\Gamma _2|^2},\end{align}
\begin{align} (S_b/S_f)_{1-2}=\frac {|\Lambda _1+\Lambda _2|^2}{|\Lambda _1|^2+|\Lambda _2|^2},\end{align}
where
\begin{align} \Gamma _1=\frac {\alpha _{1,2}+\alpha _{2,1}+\beta _{1,2}+\beta _{2,1}}{\omega _1+\omega _2-\omega _{1+2}},\end{align}
\begin{align} \Gamma _2=\frac {\alpha _{1,2}+\alpha _{2,1}-\beta _{1,2}-\beta _{2,1}}{\omega _1+\omega _2+\omega _{1+2}},\end{align}
\begin{align} \Lambda _1=\frac {\alpha _{1,-2}-\alpha _{-2,1}-\beta _{1,-2}-\beta _{-2,1}}{\omega _1-\omega _2-\omega _{1-2}},\end{align}
\begin{align} \Lambda _2=\frac {\alpha _{1,-2}-\alpha _{-2,1}+\beta _{1,-2}+\beta _{-2,1}}{\omega _1-\omega _2+\omega _{1-2}},\end{align}
with
\begin{align} \alpha _{i,j}=\frac {1}{\omega _i\omega _{i+j}}(\boldsymbol {k}_i\cdot \boldsymbol {k}_{i+j}-k_ik_{i+j}),\end{align}
\begin{align} \beta _{i,j}=\frac {1}{2\omega _i\omega _j}(\boldsymbol {k}_i\cdot \boldsymbol {k}_j+k_ik_j).\end{align}
If we further consider the one-dimensional (1-D) case where
$\boldsymbol {k}_1=(k_{x1},0)$
and
$\boldsymbol {k}_2=(k_{x2},0)$
are along the
$x$
direction, then with some manipulations, the expressions for
$S_b/S_f$
are reduced to
\begin{equation} (S_b/S_f)_{1+2}=2\bigg /\left [\left (\frac {\omega _1+\omega _2}{\omega _{1+2}}\right )^{\mu }+1\right ], \end{equation}
\begin{equation} (S_b/S_f)_{1-2}=2\bigg /\left [\left (\frac {\omega _{1-2}}{\omega _1-\omega _2}\right )^{\mu }+1\right ], \end{equation}
where
\begin{equation} \mu =2\frac {\boldsymbol {k}_1\cdot \boldsymbol {k}_2}{k_1k_2}. \end{equation}
Physical intuition on the coexistence of bound and free modes can also be established as a reflection of the analytical derivation: this phenomenon is physically similar to a frictionless pendulum subject to an external forcing at non-natural frequency. Both the homogeneous solution at natural frequency (i.e. free mode) and the inhomogeneous solution at forcing frequency (i.e. bound mode) can be excited, with the former persisting under the frictionless condition (which is usually not satisfied for vibration systems, but is indeed true for our water wave problems, at least for long waves not significantly affected by physical dissipation mechanisms such as those from breaking and viscosity).
We then conduct numerical validations for the analytical solution (4.7) and (4.8) via the HOS method. Two cases are considered: (1) a 1-D case where we fix
$\boldsymbol {k}_1=(25,0)$
and vary
$\boldsymbol {k}_2=(k_{x2},0)$
; (2) a 2-D case where we fix
$k_1=20$
,
$k_2=25$
, and vary angle
$\theta _{1,2}$
between them. The numerical solutions of
$S_b/S_f$
are plotted in figures 7(a,b) as functions of
$k_{x2}$
and
$\theta _{1,2}$
, respectively, comparing against the analytical solutions (4.7) and (4.8) for 2-D cases, and (4.15) and (4.16) for 1-D cases. We see a very good match overall. The minor deviation at a few points may be attributed to the finiteness of wave amplitude in the simulation, and a finite time window in data processing that cannot exactly resolve
$\omega _b$
and
$\omega _f$
.
Finally, we notice that the analytical solution of
$S_b/S_f$
is bounded in
$[0, 2]$
in figure 7. This is, in fact, a general condition for
$S_b/S_f$
that can be seen easily in the formulae (4.7) and (4.8) (simply because
$0\leqslant (u+v)^2/(u^2+v^2)\leqslant 2$
for any
$u,v\in \mathbb {R}$
). This indicates that it is possible for the free-wave energy to be much larger than the bound-wave energy, but not vice versa. In other words, the free waves, instead of bound waves, should always be considered as a major generation from the non-resonant triad interactions.
Figure 7. Plots of
$S_b/S_f$
as functions of (a)
$k_{x2}$
in the 1-D case, with
$\boldsymbol {k}_1=(25,0)$
,
$\boldsymbol {k}_2=(k_{x2},0)$
and
$T_L=225T_0$
, and (b)
$\theta _{1,2}$
in the 2-D case, with
$k_1=20$
,
$k_2=25$
and
$T_L=20T_0$
. For
$k=k_{1+2}$
, analytical solutions are shown by blue lines, and numerical data are shown by circles. For
$k=k_{1-2}$
, analytical solutions are shown by red lines, and numerical data are shown by triangles.
5. Conclusions and discussions
In this paper, we present a study on the role of triad interactions in the spectral evolution of surface gravity waves in deep water. A decomposition technique is developed for Euler equations, which allows us to quantify the contributions from quadratic and cubic terms in spectral evolution. We find that the contribution from quadratic terms is comparable to that from cubic terms at most wavenumbers. This is consistent with the normal form transformation applied in theory, which converts the quadratic to cubic terms so that both effects can be observed on the quartet kinetic time scale for all wavenumbers. On the other hand, the quadratic terms dominates the spectral evolution at small wavenumbers with low initial energy levels. This is shown to be from non-resonant triad interactions that are found to distribute energy between free modes and bound modes. To understand this discovery of energy distribution, we start from a simple model with two existing modes and solve analytically for the generated modes. Our analytical solution to describe the distribution is validated by numerical simulations for a range of configurations, thus providing sound interpretation of our observation in the spectral evolution.
The results in this paper should not be considered as a challenge to the gravity-wave kinetic equation. The normal form transformation is indeed valid (see the discussion for non-local triads in § 3.2), although the rigorous mathematical justification exists only for one dimension (Berti, Feola & Pusateri Reference Berti, Feola and Pusateri2023). However, some new understanding under the current framework should be noted. First, the contribution from triad interactions to spectral evolution is very significant, instead of null, which may be expected by some researchers due to their non-resonant nature. Such contributions can come from resonant interactions connected by triads, which is comparable to cubic-term resonant interactions, or simply non-resonant triad interactions. In the latter case, the non-resonant triad interactions can fill in the low-energy portion of the spectrum on a fast time scale, generating both free waves and bound waves. Only in a longer time scale does the spectral evolution become dominated by quartet interactions as described by the wave kinetic equation.
Acknowledgements.
We thank Professor J. Shatah, Professor M. Onorato, Professor S. Nazarenko, Professor F. Pusateri and Dr M. Shavit for their valuable discussions and insights during the development of this work.
Funding.
We acknowledge the Simons Foundation for the funding support for this research.
Declaration of interests.
The authors report no conflict of interest.
Appendix A. Analytical solution of the ratio
$S_b/S_f$
in § 4
To compute the ratio between the bound- and free-mode energies
$S_b/S_f$
at a given wavenumber
$\boldsymbol {k}$
, we perform a perturbation analysis in this appendix. This can be done by solving either the governing (2.1) and (2.2) directly or the equation for the complex wave amplitude
$a_{\boldsymbol {k}}$
followed by the calculation of
$\hat {\eta }_{\boldsymbol {k}}$
. We have checked that both approaches are equivalent and give the same analytical result. Here, we follow the former way, starting from the transformation of (2.1) and (2.2) into the wavenumber domain with a 2-D spatial Fourier transform. The key procedure of this process is to express the surface vertical velocity
$\phi _z$
in terms of
$\eta$
and
$\psi$
in the wavenumber domain, which is based on the small steepness assumption (i.e.
$k\eta \sim O(\epsilon )\ll 1$
). The details of this process can be found in Mei et al. (Reference Mei, Stiassnie and Yue2005)and Nazarenko & Lukaschuk (Reference Nazarenko and Lukaschuk2016) (see also Pan (Reference Pan2017) in the context of capillary waves). For our purpose of investigating triad interactions, we consider the equations truncated up to
$O(\epsilon ^2)$
:
\begin{equation} \frac {\partial \hat {\eta }_{\boldsymbol {k}}}{\partial t}=k\hat {\psi }_{\boldsymbol {k}}+\sum _{\boldsymbol {k}_1,\boldsymbol {k}_2} (\boldsymbol {k}_1\cdot \boldsymbol {k}-k_1k)\hat {\psi }_{\boldsymbol {k}_1}\hat {\eta }_{\boldsymbol {k}_2}\delta _{0-1-2}, \end{equation}
\begin{equation} \frac {\partial \hat {\psi }_{\boldsymbol {k}}}{\partial t}=-\hat {\eta }_{\boldsymbol {k}}+\frac {1}{2}\sum _{\boldsymbol {k}_1,\boldsymbol {k}_2} (\boldsymbol {k}_1\cdot \boldsymbol {k}_2+k_1k_2)\hat {\psi }_{\boldsymbol {k}_1}\hat {\psi }_{\boldsymbol {k}_2}\delta _{0-1-2}, \end{equation}
where
$\boldsymbol {k}_1,\boldsymbol {k}_2\in \mathbb {Z}_L^2=2\pi \mathbb {Z}^2/L$
, with
$L$
being the size of the periodic square domain. We set the initial conditions as in (4.1)–(4.4), expressed in the wavenumber domain as
\begin{equation} \hat {\eta }_{\boldsymbol {k}}(t=0)=\tilde {A}_1\delta _{0-1}+\tilde {A}_2\delta _{0-2}+\tilde {A}^*_1\delta _{0+1}+\tilde {A}^*_2\delta _{0+2}, \end{equation}
\begin{equation} \hat {\psi }_{\boldsymbol {k}}(t=0)=-\frac {\mathrm {i}}{\omega _k}(\tilde {A}_1\delta _{0-1}+\tilde {A}_2\delta _{0-2}-\tilde {A}^*_1\delta _{0+1}-\tilde {A}^*_2\delta _{0+2}). \end{equation}
We note that the Kronecker delta functions on the right-hand side mean that
$\hat {\eta }_{\boldsymbol {k}}$
and
$\hat {\psi }_{\boldsymbol {k}}$
take non-zero values only for
$\boldsymbol {k}=\boldsymbol {k}_1,\boldsymbol {k}_2,-\boldsymbol {k}_1, -\boldsymbol {k}_2$
, with the whole expression guaranteeing the two quantities to be real in the physical domain.
Our goal is to perform a perturbation analysis on (A1) and (A2), with the first-order solution set as the wave field corresponding to initial conditions (A3) and (A4). Then we seek the second-order solution on wavenumbers
$\boldsymbol {k}_{1+2}$
and
$\boldsymbol {k}_{1-2}$
, which identifies
$S_b/S_f$
as our quantity of interest. For this purpose, we first write the solutions of (A1) and (A2) as perturbation series in the wave steepness
$\epsilon$
:
\begin{equation} \hat {\eta }_{\boldsymbol {k}}=\hat {\eta }_{\boldsymbol {k}}^{(1)}+\hat {\eta }_{\boldsymbol {k}}^{(2)}+O(\epsilon ^3), \end{equation}
\begin{equation} \hat {\psi }_{\boldsymbol {k}}=\hat {\psi }_{\boldsymbol {k}}^{(1)}+\hat {\psi }_{\boldsymbol {k}}^{(2)}+O(\epsilon ^3), \end{equation}
where
$\hat {\eta }_{\boldsymbol {k}}^{(m)},\hat {\psi }_{\boldsymbol {k}}^{(m)}\sim O(\epsilon ^m)$
(
$m=1,2,\ldots$
). By substituting (A5) and (A6) into (A1) and (A2), we collect terms to obtain a set of equations at each order. The first-order linear equation reads
\begin{equation} \frac {\partial \hat {\eta }_{\boldsymbol {k}}^{(1)}}{\partial t}=k\hat {\psi }^{(1)}_{\boldsymbol {k}}, \end{equation}
\begin{equation} \frac {\partial \hat {\psi }_{\boldsymbol {k}}^{(1)}}{\partial t}=-\hat {\eta }^{(1)}_{\boldsymbol {k}}. \end{equation}
The solution to (A7) and (A8) is a combination of linear waves propagating in the domain. Here, we set this solution to correspond to the initial wave field (A3) and (A4), written as
\begin{align} \hat {\eta }_{\boldsymbol {k}}^{(1)}=(\tilde {A}_1 \delta _{0-1} + \tilde {A}_2 \delta _{0-2})\,{\rm e}^{-\mathrm {i}\omega _kt} + (\tilde {A}^*_1 \delta _{0+1} + \tilde {A}^*_2 \delta _{0+2})\,{\rm e}^{\mathrm {i}\omega _kt}, \end{align}
\begin{align} \hat {\psi }_{\boldsymbol {k}}^{(1)}=-\frac {\mathrm {i}}{\omega _k} [(\tilde {A}_1 \delta _{0-1} + \tilde {A}_2 \delta _{0-2})\,{\rm e}^{-\mathrm {i}\omega _kt} - (\tilde {A}^*_1 \delta _{0+1} + \tilde {A}^*_2 \delta _{0+2})\,{\rm e}^{\mathrm {i}\omega _kt}]. \end{align}
Now we consider the second-order equations
\begin{equation} \frac {\partial \hat {\eta }_{\boldsymbol {k}}^{(2)}}{\partial t}=k\hat {\psi }^{(2)}_{\boldsymbol {k}}+\sum _{\boldsymbol {k}_1,\boldsymbol {k}_2} (\boldsymbol {k}_1\cdot \boldsymbol {k}-k_1k)\hat {\psi }^{(1)}_{\boldsymbol {k}_1}\hat {\eta }^{(1)}_{\boldsymbol {k}_2}\delta _{0-1-2}, \end{equation}
\begin{equation} \frac {\partial \hat {\psi }_{\boldsymbol {k}}^{(2)}}{\partial t}=-\hat {\eta }^{(2)}_{\boldsymbol {k}}+\frac {1}{2}\sum _{\boldsymbol {k}_1,\boldsymbol {k}_2} (\boldsymbol {k}_1\cdot \boldsymbol {k}_2+k_1k_2)\hat {\psi }^{(1)}_{\boldsymbol {k}_1}\hat {\psi }^{(1)}_{\boldsymbol {k}_2}\delta _{0-1-2}. \end{equation}
With the linear solutions
$\hat {\eta }_{\boldsymbol {k}}^{(1)}$
and
$\hat {\psi }_{\boldsymbol {k}}^{(1)}$
available, (A11) and (A12) are solved as a system of non-homogeneous linear differential equations. The corresponding second-order solutions at
$\boldsymbol {k}=\boldsymbol {k}_{1+2}$
and
$\boldsymbol {k}=\boldsymbol {k}_{1-2}$
can be expressed as
\begin{equation} \hat {\eta }_{\boldsymbol {k}}^{(2)}=\mathrm {i}\omega _k(C_1\,{\rm e}^{-\mathrm {i}\omega _kt}-C_2\,{\rm e}^{\mathrm {i}\omega _kt}+N_1-N_2), \end{equation}
\begin{equation} \hat {\psi }_{\boldsymbol {k}}^{(2)}=C_1\,{\rm e}^{-\mathrm {i}\omega _kt}+C_2\,{\rm e}^{\mathrm {i}\omega _kt}+N_1+N_2, \end{equation}
where
$N_1$
,
$N_2$
have the forms
\begin{equation} N_1= \left \{ \begin{aligned} &\frac {1}{2} (-\alpha _{1,2}-\alpha _{2,1}-\beta _{1,2}-\beta _{2,1})\frac {\mathrm {i}\tilde {A}_1\tilde {A}_2}{\omega _1+\omega _2-\omega _k}\,{\rm e}^{-\mathrm {i}(\omega _1+\omega _2)t},& \boldsymbol {k}=\boldsymbol {k}_{1+2},\\ &\frac {1}{2} (-\alpha _{1,-2}-\alpha _{-2,1}+\beta _{1,-2}+\beta _{-2,1})\frac {\mathrm {i}\tilde {A}_1\tilde {A}^*_2}{\omega _1-\omega _2-\omega _k}\,{\rm e}^{-\mathrm {i}(\omega _1-\omega _2)t},& \boldsymbol {k}=\boldsymbol {k}_{1-2}, \end{aligned} \right . \end{equation}
\begin{equation} N_2= \left \{ \begin{aligned} &\frac {1}{2} (\alpha _{1,2}+\alpha _{2,1}-\beta _{1,2}-\beta _{2,1})\frac {\mathrm {i}\tilde {A}_1\tilde {A}_2}{\omega _1+\omega _2+\omega _k}\,{\rm e}^{-\mathrm {i}(\omega _1+\omega _2)t},& \boldsymbol {k}=\boldsymbol {k}_{1+2},\\ &\frac {1}{2} (\alpha _{1,-2}+\alpha _{-2,1}+\beta _{1,-2}+\beta _{-2,1})\frac {\mathrm {i}\tilde {A}_1\tilde {A}^*_2}{\omega _1-\omega _2+\omega _k}e^{-\mathrm {i}(\omega _1-\omega _2)t},& \boldsymbol {k}=\boldsymbol {k}_{1-2}, \end{aligned} \right . \end{equation}
with
\begin{equation} \alpha _{i,j}=\frac {1}{\omega _i\omega _{i+j}}(\boldsymbol {k}_i\cdot \boldsymbol {k}_{i+j}-k_ik_{i+j}), \end{equation}
\begin{equation} \beta _{i,j}=\frac {1}{2\omega _i\omega _j}(\boldsymbol {k}_i\cdot \boldsymbol {k}_j+k_ik_j). \end{equation}
According to the initial conditions (A3) and (A4), the wave field is quiescent except for modes
$\boldsymbol {k}_1$
and
$\boldsymbol {k}_2$
. We therefore have the condition at
$\boldsymbol {k}=\boldsymbol {k}_{1+2}$
and
$\boldsymbol {k}=\boldsymbol {k}_{1-2}$
,
\begin{equation} \hat {\eta }_{\boldsymbol {k}}^{(2)}(t=0)=0, \end{equation}
\begin{equation} \hat {\psi }_{\boldsymbol {k}}^{(2)}(t=0)=0, \end{equation}
from which the coefficients
$C_1$
and
$C_2$
in (A13) and (A14) are solved as
\begin{equation} C_1=-N_1(t=0), \end{equation}
\begin{equation} C_2=-N_2(t=0). \end{equation}
With the second-order solutions obtained, we can now express
$\hat {\eta }_{\boldsymbol {k}}$
at
$\boldsymbol {k}_{1+2}$
and
$\boldsymbol {k}_{1-2}$
in compact forms. For
$\boldsymbol {k}_{1+2}$
, we have
\begin{equation} \hat {\eta }_{1+2}=-\frac {\omega _{1+2}\tilde {A}_1\tilde {A}_2}{2}[\Gamma _1\,{\rm e}^{-\mathrm {i}\omega _{1+2}t}+\Gamma _2\,{\rm e}^{\mathrm {i}\omega _{1+2}t}-(\Gamma _1+\Gamma _2)\,{\rm e}^{-\mathrm {i}(\omega _1+\omega _2)t}], \end{equation}
where
$\Gamma _i$
(
$i=1,2$
) are coefficients that can be evaluated explicitly from (A15), (A16), (A21) and (A22). On the right-hand side of (A23), the first two terms correspond to free waves with frequency
$\omega _{1+2}$
(propagating in opposite directions), while the third term corresponds to bound waves with frequency
$\omega _1+\omega _2$
. Therefore the ratio between energy in bound and free modes is expressed as
\begin{equation} (S_b/S_f)_{1+2}=\frac {|\Gamma _1+\Gamma _2|^2}{|\Gamma _1|^2+|\Gamma _2|^2}, \end{equation}
where
\begin{equation} \Gamma _1=\frac {\alpha _{1,2}+\alpha _{2,1}+\beta _{1,2}+\beta _{2,1}}{\omega _1+\omega _2-\omega _{1+2}}, \end{equation}
\begin{equation} \Gamma _2=\frac {\alpha _{1,2}+\alpha _{2,1}-\beta _{1,2}-\beta _{2,1}}{\omega _1+\omega _2+\omega _{1+2}}. \end{equation}
For
$\boldsymbol {k}_{1-2}$
, similar procedures are applied, yielding
\begin{equation} \hat {\eta }_{1-2}=-\frac {\omega _{1-2}\tilde {A}_1\tilde {A}^*_2}{2}[\Lambda _1\,{\rm e}^{-\mathrm {i}\omega _{1-2}t}+\Lambda _2\,{\rm e}^{\mathrm {i}\omega _{1-2}t}-(\Lambda _1+\Lambda _2)\,{\rm e}^{-\mathrm {i}(\omega _1+\omega _2)t}], \end{equation}
which gives the ratio
\begin{equation} (S_b/S_f)_{1-2}=\frac {|\Lambda _1+\Lambda _2|^2}{|\Lambda _1|^2+|\Lambda _2|^2}, \end{equation}
where
\begin{equation} \Lambda _1=\frac {\alpha _{1,-2}-\alpha _{-2,1}-\beta _{1,-2}-\beta _{-2,1}}{\omega _1-\omega _2-\omega _{1-2}}, \end{equation}
\begin{equation} \Lambda _2=\frac {\alpha _{1,-2}-\alpha _{-2,1}+\beta _{1,-2}+\beta _{-2,1}}{\omega _1-\omega _2+\omega _{1-2}}. \end{equation}
Finally, we consider the results in the 1-D case where
$\boldsymbol {k}_1=(k_{x1},0)$
and
$\boldsymbol {k}_2=(k_{x2},0)$
are along the
$x$
direction. Then with some manipulations, the expressions for
$S_b/S_f$
are reduced to
\begin{equation} (S_b/S_f)_{1+2}=2\bigg /\left [\left (\frac {\omega _1+\omega _2}{\omega _{1+2}}\right )^{\mu }+1\right ], \end{equation}
\begin{equation} (S_b/S_f)_{1-2}=2\bigg /\left [\left (\frac {\omega _{1-2}}{\omega _1-\omega _2}\right )^{\mu }+1\right ], \end{equation}
where
\begin{equation} \mu =2\frac {\boldsymbol {k}_1\cdot \boldsymbol {k}_2}{k_1k_2}. \end{equation}
Open access
