Stability of internal gravity wave modes: from triad resonance to broadband instability

Abstract

A theoretical study is made of the stability of propagating internal gravity wave modes along a horizontal stratified fluid layer bounded by rigid walls. The analysis is based on the Floquet eigenvalue problem for infinitesimal perturbations to a wave mode of small amplitude. The appropriate instability mechanism hinges on how the perturbation spatial scale relative to the basic-state wavelength, controlled by a parameter $\mu$, compares to the basic-state amplitude parameter, $\epsilon \ll 1$. For $\mu ={O}(1)$, the onset of instability arises due to perturbations that form resonant triads with the underlying wave mode. For short-scale perturbations such that $\mu \ll 1$ but $\alpha =\mu /\epsilon \gg 1$, this triad resonance instability reduces to the familiar parametric subharmonic instability (PSI), where triads comprise fine-scale perturbations with half the basic-wave frequency. However, as $\mu$ is further decreased holding $\epsilon$ fixed, higher-frequency perturbations than these two subharmonics come into play, and when $\alpha ={O}(1)$ Floquet modes feature broadband spectrum. This broadening phenomenon is a manifestation of the advection of small-scale perturbations by the basic-wave velocity field. By working with a set of ‘streamline coordinates’ in the frame of the basic wave, this advection can be ‘factored out’. Importantly, when $\alpha ={O}(1)$ PSI is replaced by a novel, multi-mode resonance mechanism which has a stabilising effect that provides an inviscid short-scale cut-off to PSI. The theoretical predictions are supported by numerical results from solving the Floquet eigenvalue problem for a mode-1 basic state.

1. Introduction

The original motivation for the present work comes from a recent asymptotic treatment of small-scale instabilities of finite-width internal gravity wave beams in an unbounded uniformly stratified fluid (Fan & Akylas 2021). The focus of the earlier study was on validating the approximate models for parametric subharmonic instability (PSI) of internal wave beams proposed by Karimi & Akylas (2014, 2017). PSI is a particular case of triad resonance instability (TRI), where the unstable perturbations which form resonant triads with the basic wave state are fine-scale disturbances at half the basic-wave frequency. This small-scale instability mechanism has been studied widely for sinusoidal plane waves (e.g. Staquet & Sommeria 2002) and more recently for wave beams (e.g. Dauxois et al. 2018) owing to its potential significance in the dissipation of the oceanic internal tide (Hibiya, Nagasawa & Niwa 2002; MacKinnon & Winters 2005).

As wave beams are time-periodic states, Fan & Akylas (2021) used Floquet-type normal mode analysis. The associated eigenvalue problem was solved asymptotically in the limit where PSI may arise (namely, for a small-amplitude uniform wave beam subject to fine-scale perturbations under nearly inviscid conditions), but without assuming that perturbations at half the beam frequency are the dominant components of the unstable Floquet modes. Apart from assessing the validity of the earlier PSI models, this asymptotic treatment also revealed a short-scale instability that, unlike PSI, involves a broadband spectrum of frequency components. Importantly, this novel instability mechanism can affect wave beams that are not susceptible to PSI.

In view of the findings of Fan & Akylas (2021) for wave beams in an unbounded fluid, it is natural to ask whether a similar small-scale instability would apply to internal wave modes propagating along a waveguide, such as the ocean thermocline. This question is addressed here for propagating gravity wave modes in the simplest waveguide configuration of a horizontal stratified fluid layer bounded by rigid walls. Background rotation also is ignored.

In contrast to plane waves in an unbounded stratified fluid, there are only few prior studies devoted to the stability of internal wave modes. Thorpe (1966) first showed that wave modes in a continuously stratified fluid layer with a rigid bottom and a fixed or free upper surface, can form resonant triads. The conditions for such triads require that the horizontal wavevectors and the frequencies of the participating modes sum up to zero. Davis & Acrivos (1967) presented experimental evidence and theoretical confirmation of TRI for a propagating mode-1 wave in a thin stratified layer separating two homogeneous fluid layers. In a follow-up laboratory experiment, Martin, Simmons & Wunsch (1972) investigated the TRI of a mode-3 wave in a uniformly stratified fluid layer (constant background buoyancy frequency) bounded by rigid walls. Their observations of unstable disturbances generally are consistent with the theoretically predicted modes forming resonant triads with the mode-3 wave.

In more recent related work, Joubaud et al. (2012) reported the first experimental measurement of TRI growth rates for a mode-1 wave propagating along a uniformly stratified fluid tank. Varma & Mathur (2017) examined theoretically resonant triads that comprise two modes with the same frequency, in a stratified layer bounded by rigid walls and also including background rotation. Their results recover the triad resonance conditions of Thorpe (1966) and confirm that resonant triad interactions are more likely to occur in non-uniform background stratification than a uniformly stratified fluid. Sutherland & Jefferson (2020) explored the stability of mode-1 waves in a stratified layer with rigid bottom and top and in the presence of background rotation, via numerical simulations of an initial-value problem using small-amplitude noise as initial perturbation. These simulations suggest that PSI is the dominant instability for uniform background stratification but this may not be the case for other stratifications. Finally, Young, Tsang & Balmforth (2008) developed a theory for near-inertial PSI, where the frequencies of the subharmonic perturbations are assumed to be close to the inertial frequency, based on the approximate equations of Young & Jelloul (1997). This theory predicts strong instability for a mode-1 wave under conditions representative of the oceanic internal tide.

The present stability analysis is based on the eigenvalue problem that governs Floquet-mode perturbations in the moving frame where the basic wave mode is steady in time and spatially periodic along the horizontal. It follows from this problem that the onset of instability in the limit of small basic-state amplitude parameter ($\epsilon \ll 1$), arises due to TRI. The triad resonance conditions are consistent with Thorpe (1966), and the associated ${O}(\epsilon )$ growth rate is computed for general background stratification from a certain $2\times 2$ eigenvalue problem.

Attention then is focused on small-scale instabilities that involve perturbations of high modal number and short wavelength relative to the basic wave mode, assuming for simplicity constant background buoyancy frequency. The appropriate small-scale instability mechanism hinges on how the perturbation scale, controlled by a parameter $\mu \ll 1$, compares to the basic-state amplitude $\epsilon$: for $\epsilon \ll \mu \ll 1$ ($\alpha =\mu /\epsilon \gg 1$) TRI reduces to PSI; however, as $\mu$ is further decreased holding $\epsilon$ fixed, higher-frequency perturbations than the two subharmonics at half the basic-wave frequency come into play, and when $\alpha ={O}(1)$ Floquet modes feature broadband spectrum.

Similar to Fan & Akylas (2021), this broadening phenomenon is a result of the advection of small-scale perturbations by the basic-state velocity field and can be ‘factored out’ by working with a set of ‘streamline coordinates’ in the frame of the basic wave. We find that when $\alpha ={O}(1)$ PSI is replaced by a multi-mode resonance mechanism, which has a stabilising effect and provides a short-scale cut-off to PSI. An important factor in this stabilisation is the Lagrangian mean flow due to the ‘Stokes drift’ of the basic wave mode (Thorpe 1968). The theoretical predictions are supported by numerical results from solving the Floquet eigenvalue problem for a mode-1 basic state. Furthermore, estimates of instability growth rates based on dimensional scales representative of the oceanic internal tide, suggest that PSI and the inviscid cut-off discussed here could be relevant in the field.

It appears that the present asymptotic analysis of small-scale instabilities of internal wave modes can be extended to allow for non-uniform stratification and background rotation. However, a wave mode with nearly twice the inertial frequency, where near-inertial PSI becomes relevant, would require special treatment.

2. Floquet stability problem

Consider an inviscid, continuously stratified, horizontal fluid layer of uniform depth bounded by rigid walls. This configuration supports a countable infinity of horizontally propagating internal gravity wave modes (e.g. Yih 1979). The focus here is on the stability of these modes to infinitesimal two-dimensional perturbations. We use dimensionless variables with $\lambda _*/2{\rm \pi}$ as the length scale and $1/N_*$ as the time scale. Here, $\lambda _*$ denotes the wavelength in the horizontal ($x$) direction of the basic mode and $N_*$ is a characteristic value of the background buoyancy frequency, which generally is a function of the vertical coordinate $z$ pointing upwards (antiparallel to gravity). The fluid is assumed to be incompressible and the Boussinesq approximation will be made. Incompressibility is satisfied automatically by working with a streamfunction $\varPsi (x,z,t)$ such that $\varPsi _z$ and $-\varPsi _x$ are the horizontal and vertical velocity components, respectively.

It is convenient for the stability analysis to make the basic wave steady by moving along $x$ with the wave speed $c>0$. In this reference frame, the (linearised) equations governing the perturbation streamfunction $\psi (x,z,t)$ and density $\rho (x,z,t)$ read

(2.1a)\begin{gather} \frac{{\rm D}}{{\rm D}t}\nabla^2\psi - \rho_x + J(\nabla^2\bar{\varPsi},\psi) = 0, \end{gather}

(2.1b)\begin{gather}\frac{{\rm D}}{{\rm D}t}\rho + N^2\psi_x + J(\bar{\rho},\psi) = 0. \end{gather}

Here, $\bar {\varPsi }(x,z)$ is the basic-wave streamfunction (i.e. $\varPsi = \bar {\varPsi }+\psi$), $\bar {\rho }(x,z)$ is the basic-wave density, $N(z)$ is the background buoyancy frequency profile and

(2.2)\begin{equation} \frac{{\rm D}}{{\rm D}t}\equiv\frac{\partial}{\partial t} - c\frac{\partial}{\partial x} + J({\cdot},\bar{\varPsi}),\end{equation}

denotes the (linearised) advective time derivative, where $J(a,b)=a_xb_z-a_zb_x$ is the Jacobian. Furthermore, $\psi$ obeys

(2.3)\begin{equation} \psi_x=0\quad (z=0,H), \end{equation}

on the channel walls ($z=0,H$), where $H$ is the dimensionless fluid depth.

Before proceeding to the stability analysis, we specify $\{\bar {\varPsi },\bar {\rho };c\}$. Our choice for basic state is a finite-amplitude progressive wave mode with (normalised) wavelength $2{\rm \pi}$ along $x$ and phase speed $c$, which generally is a function of the wave amplitude (Yih 1974). Such finite-amplitude waves of permanent form are steady solutions of the nonlinear stratified flow equations and need to be computed separately, by amplitude expansions (Thorpe 1968; Yih 1974) or numerically.

Our interest here is on the onset of instability, which occurs when the amplitude of the basic state is small. Accordingly, we introduce the small-amplitude parameter $\epsilon$,

(2.4)\begin{equation} \epsilon = \frac{2{\rm \pi} U_*}{\lambda_*N_*}\ll1, \end{equation}

where $U_*$ is a characteristic velocity of the basic wave. In this limit, $\{\bar {\varPsi },\bar {\rho };c\}$ may be approximated as Yih (1974)

(2.5a)\begin{gather} \bar{\varPsi} = \epsilon\bar{f}_n(z)\cos x + {O}(\epsilon^2), \end{gather}

(2.5b)\begin{gather}\bar{\rho} = \epsilon\frac{N^2}{c_n}\bar{f}_n(z)\cos x + {O}(\epsilon^2), \end{gather}

(2.5c)\begin{gather}c = c_n +{O}(\epsilon^2). \end{gather}

Here, the wave speed $c_n$ and mode shape $\bar {f}_n(z)$ are the eigenvalue and eigenfunction corresponding to a certain eigensolution ($n=1,2,\ldots$) of the problem

(2.6a)\begin{gather} \frac{\mathrm{d}^2\bar{f}_n}{\mathrm{d}z^2}+ \left(\frac{N^2}{c_n^2}-1\right)\bar{f}_n = 0, \end{gather}

(2.6b)\begin{gather}\bar{f}_n = 0 \quad(z=0,H). \end{gather}

It should be noted that for uniform background stratification ($N=1$), the problem (2.6) has the simple closed-form solution

(2.7a,b)\begin{equation} \bar{f}_n=\frac{H}{n{\rm \pi}}\sin\frac{n{\rm \pi}}{H}z; \quad c_n=\frac{H}{(n^2{\rm \pi}^2+H^2)^{1/2}}\quad(n=1,2,\ldots), \end{equation}

where $\bar {f}_n$ has been normalised such that $U_*$ in (2.4) is the peak horizontal velocity of the mode-$n$ wave. Moreover, when $N=1$ the assumed background flow conditions conform to ‘Long's model’ for steady stratified flow (Dubreil-Jacotin 1932; Long 1953), so the small-amplitude basic state (2.5) with $\bar {f}_n$ and $c_n$ given by (2.7a,b) also satisfies the full nonlinear equations of motion. Accordingly, in this instance it is permissible to use this basic state in a stability analysis for any $\epsilon$ below the threshold for overturning $\epsilon _c=c_n$.

Returning now to the governing equations (2.1), we follow the procedure used in the stability analysis of a Stokes surface gravity wave (e.g. McLean 1982): as the basic state (2.5) is steady in $t$ and $2{\rm \pi}$-periodic in $x$, we look for Floquet-mode solutions in the form

(2.8)\begin{equation} (\psi,\rho) = \mathrm{e}^{-\mathrm{i}\sigma t}\,\mathrm{e}^{\mathrm{i} px}\sum_{m={-}\infty}^{\infty}(Q_m(z),R_m(z))\,\mathrm{e}^{\mathrm{i}mx}. \end{equation}

Here, in keeping with a temporal stability analysis, $p$ is a prescribed real wavenumber and $\sigma$ is a possibly complex frequency to be determined along with the Fourier coefficients $Q_m$ and $R_m$. It should be noted that, without loss of generality, $p$ can be restricted in the range $0\leq p<1$; however, here we find it convenient to treat $p$ as a free parameter $(-\infty < p<\infty )$.

Upon substituting (2.8) in (2.1), we first eliminate the density coefficients $R_m(z)$ and then simplify the equations for the streamfunction coefficients $Q_m(z)$ by using the ${O}(1)$ balance of terms to eliminate third-order $z$-derivatives at ${O}(\epsilon )$. Finally, $Q_m(z)$ are governed by the infinite equation system ($-\infty < m<\infty$)

(2.9)\begin{align} &\frac{\mathrm{d}^2Q_m}{\mathrm{d}z^2}+(p+m)^2 \left(\frac{N^2}{\varOmega_m^2}-1\right)Q_m\nonumber\\ &\quad =\frac{\epsilon}{2\varOmega_m}\left(K_m^+Q_{m+1}+K_m^-Q_{m-1}+D_m^+ \frac{\mathrm{d}Q_{m+1}}{\mathrm{d}z}+D_m^- \frac{\mathrm{d}Q_{m-1}}{\mathrm{d}z}\right)+ {O}(\epsilon^2), \end{align}

where

(2.10a)\begin{align} K_m^\pm & = (p+m\pm1)\left\{\vphantom{\left(\frac{p+m}{\varOmega_m\varOmega_{m\pm1}}+ \frac{p+m\pm1}{\varOmega_{m\pm1}^2}\mp\frac{1}{c_n^2}\mp \frac{p+m}{c_n\varOmega_m}\right)}N^2\frac{\mathrm{d}\bar{f}_n}{\mathrm{d}z} \left[\frac{1}{c_n^2}+\frac{p+m}{c_n\varOmega_m}- \frac{(p+m\pm1)}{\varOmega_{m\pm1}} \left(\frac{p+m}{\varOmega_m}+\frac{p+m\pm1}{\varOmega_{m\pm1}} \right)\right]\right.\nonumber\\ &\quad \mp\left.(N^2)_z\bar{f}_n\left(\frac{p+m}{\varOmega_m\varOmega_{m\pm1}}+ \frac{p+m\pm1}{\varOmega_{m\pm1}^2}\mp\frac{1}{c_n^2}\mp \frac{p+m}{c_n\varOmega_m}\right)\right\}, \end{align}

(2.10b)\begin{align} D_m^\pm & ={\mp} N^2\bar{f}_n\left(\frac{(p+m)(p+m\pm1)}{\varOmega_m \varOmega_{m\pm1}}+\frac{(p+m\pm1)^2}{\varOmega_{m\pm1}^2} -\frac{1}{c_n^2}-\frac{p+m}{c_n\varOmega_m}\right), \end{align}

with

(2.11)\begin{equation} \varOmega_m=\sigma+(p+m)c_n.\end{equation}

Furthermore, in view of (2.3), $Q_m(z)$ satisfy the boundary conditions

(2.12)\begin{equation} Q_m = 0\quad(z=0,H).\end{equation}

The equation system (2.9) along with the boundary conditions (2.12) constitute an eigenvalue problem for $\sigma =\sigma _\mathrm {r}+\mathrm {i}\sigma _{\mathrm {i}}$. Given that the governing equations (2.1) are real, eigenvalues appear in complex conjugate pairs and $\sigma _{\mathrm {i}}\neq 0$ is sufficient for instability. The following discussion focuses on solving the eigenvalue problem (2.9) and (2.12) for $0<\epsilon \ll 1$ and understanding the various instability mechanisms in this limit.

3. Triad resonance instability

3.1. Resonant triads

As expected, if the basic wave is absent ($\epsilon =0$), the eigenvalue problem (2.9) and (2.12) recovers the free propagating modes in the fluid layer. Specifically, as the coefficients $Q_m$ are entirely uncoupled when $\epsilon =0$, we write

(3.1)\begin{equation} Q_m = q_{m,l}(z;p+m)\quad(l=1,2,\ldots), \end{equation}

for any given $-\infty < m<\infty$, with the rest of the $Q$ being zero. Here, $q_{m,l}$ denote the eigenfunctions of the problem

(3.2a)\begin{gather} \frac{\mathrm{d}^2q_{m,l}}{\mathrm{d}z^2}+ \left(\frac{N^2}{c_{m,l}^2}-(p+m)^2\right)q_{m,l} = 0, \end{gather}

(3.2b)\begin{gather}q_{m,l} = 0 \quad(z=0,H), \end{gather}

and $c_{m,l}^2$ are the corresponding eigenvalues, with

(3.3)\begin{equation} c_{m,l}^2=\frac{\varOmega_m^2}{(p+m)^2}\quad(l=1,2,\ldots). \end{equation}

It should be noted that (3.3) are the dispersion relations of (the countable infinity of) free propagating modes $(l=1,2,\ldots )$, and $\varOmega _m=\pm (p+m)c_{m,l}$ (with $c_{m,l}>0$) are the frequencies (in the rest frame) of these modes at the wavenumber $p+m$. Hence, in view of (2.11), for $\epsilon =0$, the stability eigenvalues $\sigma$ are simply the (Doppler-shifted) frequencies of these waves in the frame moving with $c_n$,

(3.4)\begin{equation} \sigma=\sigma^\pm{\equiv} -(p+m)c_n\pm(p+m)c_{m,l}\quad(l=1,2,\ldots).\end{equation}

Next, we inquire into how the interaction with the underlying wave affects $\sigma$ in the small-amplitude limit ($0<\epsilon \ll 1$). According to (2.9), to leading order in $\epsilon$, each $Q_m$ is coupled to its nearest neighbours $Q_{m\pm 1}$ only. As this coupling is weak, it is natural to attempt to solve the eigenvalue problem (2.9) and (2.12) approximately in an iterative manner, starting from the known solution for $\epsilon =0$ (cf. (3.1)–(3.3)). Specifically, for any given $-\infty < m<\infty$ and $l=1,2,\ldots$, one may anticipate that

(3.5a,b)\begin{equation} Q_m=q_{m,l}(z;p+m)+{O}(\epsilon^2), \quad Q_{m\pm1}=\epsilon \hat{q}_{m\pm1,l}(z)+{O}(\epsilon^2), \end{equation}

with the rest of the $Q$ being smaller than ${O}(\epsilon )$ and

(3.6)\begin{equation} \sigma = \sigma^\pm{+}{O}(\epsilon^2). \end{equation}

Here, $q_{m,l}$ and $\sigma ^\pm$ are the free-mode eigenfunctions and (real) frequencies defined in (3.2) and (3.4), respectively, and the correction terms $\hat {q}_{m\pm 1,l}$ satisfy the forced equations

(3.7a)\begin{equation} \frac{\mathrm{d}^2\hat{q}_{m\pm1,l}}{\mathrm{d}z^2}+(p+m\pm1)^2 \left(\frac{N^2}{\varOmega_{m\pm1}^2}-1\right)\hat{q}_{m\pm1,l} = \frac{1}{2\varOmega_{m\pm1}}\left(K_{m\pm1}^\mp q_{m,l} +D_{m\pm1}^\mp \frac{\mathrm{d}q_{m,l}}{\mathrm{d}z}\right),\end{equation}

subject to the boundary conditions

(3.7b)\begin{equation} \hat{q}_{m\pm1,l} = 0 \quad(z=0,H).\end{equation}

According to (3.5a,b), to leading order, the interaction with the basic wave induces the nearest two neighbours of $Q_m$ to ${O}(\epsilon )$. Furthermore, as indicated by (3.6), no instability is predicted at ${O}(\epsilon )$.

It is important to note, however, that the above (naive) approximation procedure breaks down if $\varOmega _{m\pm 1}^2/(p+m\pm 1)^2$ in (3.7a) happens to coincide with a free-mode eigenvalue (i.e. an eigenvalue $c^2_{m\pm 1,l}$ of the problem (3.2)). Under this resonance condition, which as discussed in the following may be interpreted as two free modes forming a resonant triad with the underlying wave, the forced problems (3.7) generally cannot be solved, and (3.5a,b)–(3.6) need to be revised. Moreover, in this instance it turns out that $\sigma _\mathrm {i}={O}(\epsilon )$, so triad resonances are associated with the onset of instability in the limit $\epsilon \ll 1$.

To analyse this TRI, without loss of generality ($-\infty < p<\infty$ is a free parameter), suppose that $(m=0,l)$ and $(m=1,l+r)$ are resonant free modes, where $l$ and $l+r$ are positive integers. Then, in view of (2.11) and (3.3), the corresponding eigenvalues $c_{0,l}$ and $c_{1,l+r}$ must satisfy

(3.8a)\begin{gather} (\sigma+pc_n)^2=p^2c_{0,l}^2, \end{gather}

(3.8b)\begin{gather}(\sigma+(p+1)c_n)^2=(p+1)^2c_{1,l+r}^2. \end{gather}

Hence,

(3.9)\begin{equation} \pm pc_{0,l}={-}c_n\pm(p+1)c_{1,l+r}, \end{equation}

where the $\pm$ signs above can be chosen independently. Therefore, the wavenumbers $k_0=p$ and $k_1=p+1$, along with the corresponding frequencies (in the rest frame) $\omega _0=\pm pc_{0,l}$ and $\omega _1=\pm (p+1)c_{1,l+r}$ of free modes that satisfy the resonance conditions (3.8), are linked via

(3.10a,b)\begin{equation} k_1-k_0=1,\quad \omega_1\pm\omega_0=c_n. \end{equation}

This confirms that such modes form a resonant triad with the basic wave as the latter has (normalised) wavenumber $1$ and frequency (in the rest frame) $c_n$. Alternatively, in the moving frame where the basic wave has zero frequency, the two resonant free modes have the same frequency, $\sigma$, and the frequency condition in (3.10a,b) is met trivially.

3.2. TRI eigenvalue problem

For given $l$ and $r$, conditions (3.8) determine specific (real) $p=p_c$ and $\sigma =\sigma _c$, say, at which the triad conditions (3.10a,b) are met. (Equation (3.8) may admit multiple such solutions.) Close to these critical values, we write

(3.11a,b)\begin{equation} p=p_c+\hat{p}\epsilon,\quad \sigma=\sigma_c+\lambda\epsilon, \end{equation}

where $\hat {p}$ is a real ${O}(1)$ wavenumber detuning and $\lambda$ is a possibly complex eigenvalue perturbation. In this neighbourhood, we seek solutions of the eigenvalue problem (2.9) and (2.12) in the form

(3.12a)\begin{gather} Q_0 = A_0q_{0,l}(z;p_c)+\epsilon\hat{q}_{0,l}(z) +{O}(\epsilon^2), \end{gather}

(3.12b)\begin{gather}Q_1 = A_1q_{1,l+r}(z;p_c+1)+\epsilon\hat{q}_{1,l+r}(z) +{O}(\epsilon^2), \end{gather}

with $Q_{-1},Q_2={O}(\epsilon )$ and the rest of the $Q$ $\mathscr {o}(\epsilon )$. Here, $q_{0,l}$ and $q_{1,l+r}$ are the eigenfunctions corresponding to the resonant eigenvalues $c_{0,l}$ and $c_{1,l+r}$ of the problem (3.2), under the normalisation

(3.13)\begin{equation} \int_{0}^{H}N^2q_{0,l}^2\, \mathrm{d}z=\int_{0}^{H}N^2q_{1,l+r}^2\, \mathrm{d}z=1,\end{equation}

and $A_0,A_1$ are constants.

Upon substituting (3.12) along with (3.11a,b) in (2.9) and (2.12), the ${O}(1)$ balance of terms is satisfied automatically. Next, the ${O}(\epsilon )$ corrections to $Q_0$ and $Q_1$ in (3.12) are to be found by solving the forced problems

(3.14a)\begin{gather} \frac{\mathrm{d}^2\hat{q}_{0,l}}{\mathrm{d}z^2}+ \left(\frac{N^2}{c_{0,l}^2}-k_0^2\right)\hat{q}_{0,l} = \mathcal{R}_{0,l}, \end{gather}

(3.14b)\begin{gather}\hat{q}_{0,l} = 0 \quad(z=0,H), \end{gather}

and

(3.15a)\begin{gather} \frac{\mathrm{d}^2\hat{q}_{1,l+r}}{\mathrm{d}z^2}+ \left(\frac{N^2}{c_{1,l+r}^2}-k_1^2\right)\hat{q}_{1,l+r} = \mathcal{R}_{1,l+r}, \end{gather}

(3.15b)\begin{gather}\hat{q}_{1,l+r} = 0 \quad(z=0,H). \end{gather}

Here,

(3.16a)\begin{align} \mathcal{R}_{0,l} &=2A_0q_{0,l}\left\{k_0\hat{p}-N^2 \frac{k_0^2}{\omega_0^2}\left(\frac{\hat{p}}{k_0}- \frac{\lambda+c_n\hat{p}}{\omega_0}\right)\right\}\nonumber\\ &\quad+\frac{A_1}{2\omega_0}\left(K_0^+q_{1,l+r}+D_0^+ \frac{\mathrm{d}q_{1,l+r}}{\mathrm{d}z}\right), \end{align}

(3.16b)\begin{align} \mathcal{R}_{1,l+r} &= 2A_1q_{1,l+r}\left\{k_1\hat{p}-N^2 \frac{k_1^2}{\omega_1^2}\left(\frac{\hat{p}}{k_1}-\frac{\lambda +c_n\hat{p}}{\omega_1}\right)\right\}\nonumber\\ &\quad+\frac{A_0}{2\omega_1}\left(K_1^-q_{0,l}+D_1^- \frac{\mathrm{d}q_{0,l}}{\mathrm{d}z}\right), \end{align}

where $k_0=p_c, \omega _0=\sigma _c+c_np_c$ and $k_1=p_c+1$, $\omega _1=\sigma _c+c_n(p_c+1)$ are the resonant triad wavenumbers and frequencies (in the rest frame), and the constants $K_0^+,K_1^-,D_0^+$ and $D_1^-$ are evaluated using (2.10) at $p=p_c$ and $\sigma =\sigma _c$.

Now, similar to the forced problems (3.7), we ask whether the forced problems (3.14) and (3.15) can be solved, given that the corresponding homogeneous problems have non-trivial solutions, namely the eigenfunctions $q_{0,l}(z;p_c)$ and $q_{1,l+r}(z;p_c+1)$, respectively. It turns out that, for (3.14) and (3.15) to be solvable, the forcing terms $\mathcal {R}_{0,l}$ and $\mathcal {R}_{1,l+r}$ must be orthogonal to these homogeneous solutions:

(3.17a,b)\begin{equation} \int_{0}^H\mathcal{R}_{0,l}q_{0,l}\, \mathrm{d}z = 0,\quad \int_{0}^H\mathcal{R}_{1,l+r}q_{1,l+r}\, \mathrm{d}z = 0. \end{equation}

The above solvability conditions are a particular instance of the Fredholm alternative (e.g. Haberman 2012). Here, they are obtained by multiplying both sides of (3.14a) and (3.15a) with $q_{0,l}(z;p_c)$ and $q_{1,l+r}(z;p_c+1)$, respectively, and integrating in $z$ from $0$ to $H$. After two integrations by parts and using the boundary conditions (3.14b) and (3.15b), it follows that the left-hand sides of these equations vanish; thus, the right-hand sides must do so as well, implying (3.17a,b).

Inserting the forcing terms (3.16) in the solvability conditions (3.17a,b) yields the following $2\times 2$ eigenvalue problem for $\lambda$

(3.18a)\begin{gather} (\lambda-c_g^0\hat{p})A_0 = E_1A_1, \end{gather}

(3.18b)\begin{gather}(\lambda-c_g^1\hat{p})A_1 = E_2A_0. \end{gather}

Here, using (2.10), the interaction coefficients $E_1$ and $E_2$ can be brought to the form

(3.19a)\begin{align} E_1 &=\frac{\omega_0^2}{4k_0^2}\left\{(k_1I_1+I_3) \left(\frac{k_1^2}{\omega_1^2}+\frac{k_0k_1}{\omega_0\omega_1}- \frac{k_0}{c_n\omega_0}-\frac{1}{c_n^2}\right)\right.\nonumber\\ &\quad\left.+k_1I_2\left(\frac{k_0}{\omega_0\omega_1}+ \frac{k_1}{\omega_1^2}-\frac{k_0}{c_n\omega_0}- \frac{1}{c_n^2}\right)\right\},\end{align}

(3.19b)\begin{align} E_2 &= \frac{\omega_1^2}{4k_1^2}\left\{(k_0I_1-I_4) \left(\frac{k_0^2}{\omega_0^2}+\frac{k_0k_1}{\omega_0\omega_1}- \frac{k_1}{c_n\omega_1}-\frac{1}{c_n^2}\right)\right.\nonumber\\ &\quad\left.-k_0I_2\left(\frac{k_1}{\omega_0\omega_1}+ \frac{k_0}{\omega_0^2}+\frac{k_1}{c_n\omega_1}+ \frac{1}{c_n^2}\right)\right\}, \end{align}

where the constants $I_1,\ldots,I_4$ are given by

(3.20a)$$\begin{gather} I_1 = \int_{0}^HN^2\bar{f}_n^\prime q_{0,l}q_{1,l+r}\, \mathrm{d}z, \end{gather}$$

(3.20b)$$\begin{gather}I_2 = \int_{0}^H(N^2)^\prime\bar{f}_n q_{0,l}q_{1,l+r}\, \mathrm{d}z, \end{gather}$$

(3.20c)$$\begin{gather}I_3 = \int_{0}^HN^2\bar{f}_n q_{0,l}q_{1,l+r}^\prime\,\mathrm{d}z, \end{gather}$$

(3.20d)$$\begin{gather}I_4 = \int_{0}^HN^2\bar{f}_nq_{0,l}^\prime q_{1,l+r}\, \mathrm{d}z ={-}(I_1+I_2+I_3), \end{gather}$$

with prime denoting derivative with respect to $z$. Finally, the constants $c_g^0$ and $c_g^1$ in (3.18) are associated with the wavenumber detuning in (3.11a,b). From (3.17a,b), also making use of (3.13), these constants can be expressed as

(3.21a)\begin{gather} c_g^0 ={-}c_n+\frac{\omega_0}{k_0}\left(1-\omega_0^2 \int_{0}^Hq_{0,l}^2\, \mathrm{d}z\right), \end{gather}

(3.21b)\begin{gather}c_g^1 ={-}c_n+\frac{\omega_1}{k_1}\left(1-\omega_1^2 \int_{0}^Hq_{1,l+r}^2\, \mathrm{d}z\right), \end{gather}

and they represent the group velocities (in the frame moving with the basic wave) of the modes that form a resonant triad with the basic wave. (Ignoring their interaction with the underlying wave, these modes would be free propagating waves, so a wavenumber shift $\hat {p}\epsilon$ would cause a frequency shift $c_g\hat {p}\epsilon$ in (3.11a,b); i.e. $\lambda =c_g^0\hat {p},c_g^1\hat {p}$, consistent with (3.18) for $E_1=E_2=0$.)

Based on the eigenvalue problem (3.18), instability ($\lambda =\lambda _\mathrm {r}+\mathrm {i}\lambda _\mathrm {i}$ complex) requires $E_1E_2<0$. Moreover, under this condition, instability is present within the ${O}(\epsilon )$ wavenumber window $p=p_c+\hat {p}\epsilon$ specified by

(3.22)\begin{equation} \hat{p}^2 <{-}\frac{4E_1E_2}{(c_g^1-c_g^0)^2}, \end{equation}

with the maximum growth rate

(3.23)\begin{equation} \sigma_{\mathrm{i}}\mid_{max}=\epsilon\lambda_\mathrm{i} \mid_{max}=\epsilon({-}E_1E_2)^{1/2}, \end{equation}

realised at $\hat {p}=0$ $(p=p_c)$.

3.3. The case $N=1$

In the case of uniform background stratification ($N=1$), the eigenfunctions $q_{m,l}$ of the eigenvalue problem (3.2) are sines that depend on the modal number $l$ but not on the wavenumber $p+m$,

(3.24)\begin{equation} q_{m,l} = C\sin\frac{l{\rm \pi}}{H}z\quad(l=1,2,\ldots), \end{equation}

where $C$ is a normalisation constant, and the free-mode dispersion relations (3.3) read

(3.25)\begin{equation} \varOmega_m^2 = \frac{(p+m)^2}{(p+m)^2+(l{\rm \pi}/H)^2}\quad(l=1,2,\ldots). \end{equation}

Thus, the eigenfunctions of possible resonant modes $(m=0,l)$ and $(m=1,l+r)$, under the normalisation (3.13), take the form

(3.26a,b)\begin{equation} q_{0,l} = \left(\frac{2}{H}\right)^{1/2}\sin \frac{l{\rm \pi}}{H}z, \quad q_{1,l+r} = \left(\frac{2}{H}\right)^{1/2} \sin\frac{(l+r){\rm \pi}}{H}z. \end{equation}

Next, based on (3.19) and (3.20), we compute the interaction coefficients $E_1,E_2$ in the TRI eigenvalue problem (3.18). These are linear combinations of $I_1,\ldots,I_4$ and, according to (3.20b), $I_2\equiv 0$ when $N$ is constant. Moreover, using (3.26a,b) and the basic-wave mode (2.7a,b), it follows from (3.20) that the rest of the $I$ vanish as well, unless $r=\pm n$. Thus, in the case of uniform background stratification, TRI requires that perturbations, apart from the resonant triad conditions (3.10a,b), also satisfy

(3.27)\begin{equation} l_1-l_0={\pm} n, \end{equation}

where $n$ is the basic-wave modal number (cf. (2.7a,b)) and $l_0,l_1$ denote the modal numbers of the perturbations. This condition is reminiscent of that satisfied along the vertical by the wavevectors of resonant triads in the TRI of propagating plane waves in an unbounded fluid (e.g. Mied 1976). Here, however, the basic state as well as the perturbations are standing waves in the vertical; moreover, in contrast to the triad conditions (3.10a,b), the constraint (3.27) applies only when $N$ is constant. The fact that uniform background stratification limits possible resonant triad interactions of wave modes was also noted by Varma & Mathur (2017).

To be specific, we satisfy (3.27) by taking $l_0=l$ and $l_1=l+n$, where $l=1,2,\ldots$. Then, from (3.19), making also use of (2.7a,b), (3.20) and (3.26a,b), we find that

(3.28a)\begin{gather} E_1 = \frac{1}{8n}(nk_0-l)\frac{\omega_0^2}{k_0^2} \left\{\frac{k_1^2}{\omega_1^2}+\frac{k_0k_1}{\omega_0\omega_1}- \frac{k_0}{c_n\omega_0}-\frac{1}{c_n^2}\right\}, \end{gather}

(3.28b)\begin{gather}E_2 = \frac{1}{8n}(nk_0-l)\frac{\omega_1^2}{k_1^2} \left\{\frac{k_0^2}{\omega_0^2}+\frac{k_0k_1}{\omega_0\omega_1}- \frac{k_1}{c_n\omega_1}-\frac{1}{c_n^2}\right\}. \end{gather}

Here, in keeping with (3.10a,b), $k_0=p_c$, $k_1=p_c+1$, $\omega _0=\sigma _c+c_nk_0$ and $\omega _1=\sigma _c+c_nk_1$ are the triad wavenumbers and frequencies, where $p=p_c$ and $\sigma =\sigma _c$ are obtained from the resonance conditions (3.8). In view of (3.25), for $N=1$ these conditions take the form

(3.29a)\begin{gather} \omega_0^2=\frac{p_c^2}{p_c^2+(l{\rm \pi}/H)^2}, \end{gather}

(3.29b)\begin{gather}\omega_1^2=\frac{(p_c+1)^2}{(p_c+1)^2+((l+n){\rm \pi}/H)^2}. \end{gather}

Expressions (3.28) agree with Martin et al. (1972) after converting to their non-dimensional variables.

3.4. Comparison with Joubaud et al. (2012)

The laboratory experiments of Joubaud et al. (2012) employed a wave generator at one end of a uniformly stratified fluid tank to excite monochromatic mode-1 waves which eventually became unstable due to TRI as they propagated along the tank. For each basic wave, Joubaud et al. (2012) verified experimentally that the unstable disturbances satisfied the triad resonance conditions and also measured the instability growth rate. Furthermore, they compared the observed TRI growth rates with theoretical estimates based on the TRI of a sinusoidal plane wave in an unbounded fluid.

Here, we make a brief comparison of these observations with the theoretically predicted TRI for mode-1 ($n=1$) waves in a uniformly stratified fluid ($N=1$). Specifically, we focus on the basic wave corresponding to $c_n=0.95$, $H=9.2$ and $\epsilon =0.14$ (in our dimensionless variables), for which Joubaud et al. (2012) report the strongest TRI. In this instance, the resonance conditions (3.29) for $l=9$ yield $k_0=1.3$, $\omega _0=-0.39$, $k_1=2.3$ and $\omega _1=0.56$. This resonant triad is a good approximation to the frequencies $\omega _0=-0.38$, $\omega _1=0.57$ as well as the wavenumbers of the observed unstable disturbances in figures 1 and 2 of Joubaud et al. (2012). The growth rate found from (3.23) and (3.28) for this triad is $\sigma _\mathrm {i}|_{{theor}}=6.8\times 10^{-2}$ whereas the measured growth rate is $\sigma _\mathrm {i}|_{{exp}}\approx 5.3\times 10^{-2}$. This fair agreement seems reasonable given that the theory does not account for viscous damping so $\sigma _\mathrm {i}|_{theor}$ is expected to overpredict $\sigma _\mathrm {i}|_{exp}$.

4. Short-scale disturbances

We now focus on small-scale instabilities that involve disturbances with high modal number $(l\gg n)$ and large wavenumber $(p\gg 1)$ relative to the basic wave mode. In this limit, although the eigenvalue problem (3.2) generally cannot be solved exactly by analytical means, it is possible to compute the eigenfunctions (3.1) and dispersion relations (3.3) of free modes via the WKB approximation. For simplicity, however, here and in the rest of the paper, we assume uniform background stratification $(N=1)$, where exact expressions are available (cf. (3.24)–(3.25)).

4.1. Parametric subharmonic instability

To analyse short-scale instabilities, we introduce a parameter $\mu$ that controls the perturbation vertical length scale and also we scale the horizontal wavenumber $p$ in sympathy with $1/\mu$,

(4.1a,b)\begin{equation} \mu=\frac{H}{{\rm \pi} l},\quad p=\frac{\kappa}{\mu}, \end{equation}

where $\kappa ={O}(1)$ is a re-scaled wavenumber. Thus, conditions (3.29) for determining $p=p_c$ and $\sigma =\sigma _c$ at which the onset of TRI occurs for given $l$, take the form

(4.2a)\begin{gather} \kappa_c^2\left(\frac{1}{\omega_0^2}-1\right)=1, \end{gather}

(4.2b)\begin{gather}(\kappa_c+\mu)^2\left(\frac{1}{\omega_1^2}-1\right)= \left(1+\frac{n{\rm \pi}}{H}\mu\right)^2. \end{gather}

These re-scaled conditions specify $\kappa _c=\mu p_c$ and $\omega _0=\sigma _c+c_np_c$ ($\omega _1=\omega _0+c_n$), for given $\mu$.

In the short-scale limit of interest here, (4.2) are solved by expanding in $\mu \ll 1$

(4.3a)\begin{equation} \kappa_c=\kappa_0+\mu{\rm \Delta}\kappa+\cdots,\quad \omega_0={-}\frac{c_n}{2}+\mu{\rm \Delta}\omega_0+\cdots,\end{equation}

where

(4.3b)\begin{equation} \kappa_0={\pm}\frac{c_n}{(4-c_n^2)^{1/2}},\quad {\rm \Delta} \kappa=\frac{1}{2}\left(\frac{n{\rm \pi}}{H}\kappa_0-1\right),\quad {\rm \Delta}\omega_0={-}\frac{c_n^3}{8\kappa_0^3}{\rm \Delta}\kappa.\end{equation}

Therefore, in this limit, TRI involves two short-scale modes with frequencies (in the rest frame) half the basic-wave frequency: $\omega _0=-c_n/2$ and $\omega _1= c_n/2$. This is the hallmark of the widely studied PSI of sinusoidal plane internal waves and plane wave beams in an unbounded uniformly stratified fluid (Staquet & Sommeria 2002; Dauxois et al. 2018).

Using (4.3), we may compute asymptotically the interaction coefficients $E_1, E_2$ in (3.28) of the TRI eigenvalue problem (3.18), in the PSI regime. Specifically,

(4.4a)\begin{gather} E_1\sim{-}\frac{1}{16(1-{c_n}^2)^{1/2}} \{(1-c_n^2)^{3/2}\pm(2c_n^2+1)(1-c_n^2/4)^{1/2}\}, \end{gather}

(4.4b)\begin{gather}E_2\sim{-}E_1, \end{gather}

where the $\pm$ sign corresponds to $\kappa _0=\pm c_n/(4-c_n^2)^{1/2}$ in (4.3b). This confirms that $E_1E_2<0$ so, in view of (3.23), PSI is always possible. Furthermore, numerical results (see § 6) indicate that PSI (for the $+$ sign in (4.4), which provides a higher growth rate) is the dominant resonant triad instability.

4.2. Beyond PSI

Recent asymptotic analysis of the Floquet stability eigenvalue problem for internal wave beams (Fan & Akylas 2021) pointed out that, as the length scale of the perturbation is decreased (for small but fixed beam amplitude), Floquet modes become ‘broadband’: they develop higher-frequency components than the two subharmonics at half the basic wave frequency which are dominant in PSI. This broadening of the frequency spectrum had been noted in earlier numerical work (Onuki & Tanaka 2019) and was attributed to the advection of the perturbation by the underlying wave beam. By adopting a frame riding with the wave beam, Fan & Akylas (2021) were able to ‘factor out’ this advection effect and reveal a novel small-scale instability mechanism, distinct from PSI.

Motivated by these findings, we now return to expansion (3.12) and examine the behaviour of the ${O}(\epsilon )$ Fourier coefficients $Q_{-1}$ and $Q_2$ in the short-scale limit ($\mu \ll 1$). It should be noted that, because $\omega _0\sim -c_n/2$ and $\omega _1\sim c_n/2$ in this limit according to (4.3a), these Fourier coefficients are associated with the $\pm 3c_n/2$ frequency components (in the rest frame) of the Floquet mode (2.8).

Specifically, from (2.9) and (2.12) combined with (3.12a), $Q_{-1}$ satisfies the forced equation

(4.5a)\begin{equation} \frac{\mathrm{d}^2Q_{{-}1}}{\mathrm{d}z^2}+(k_0-1)^2 \left(\frac{1}{\varOmega_{{-}1}^2}-1\right)Q_{{-}1}= \epsilon\mathcal{R}_{{-}1}, \end{equation}

subject to

(4.5b)\begin{equation} Q_{{-}1}=0\quad(z=0,H), \end{equation}

where

(4.6)\begin{equation} \mathcal{R}_{{-}1} = \frac{A_0}{2\varOmega_{{-}1}}\left(K_{{-}1}^+q_{0,l}+D_{{-}1}^+ \frac{\mathrm{d}q_{0,l}}{\mathrm{d}z}\right). \end{equation}

In general, the boundary-value problem (4.5) is solvable as $\varOmega _{-1}=\omega _0-c_n$ does not match the frequency of a free mode at the wavenumber $k_{-1}\equiv k_0-1$; i.e. $k_{-1},\varOmega _{-1}$ do not satisfy the dispersion relation (3.25) for any (integer) modal number $l$. Rather than determining the detailed solution, however, here it suffices to look at the asymptotic behaviour of $Q_{-1}$ for $\mu \ll 1$. Briefly, upon combining (2.7a,b), (2.10) and (3.26a,b) with (4.1), (4.3) and $\varOmega _{-1}\sim -3c_n/2$, we find from (4.6)

(4.7)\begin{equation} \mathcal{R}_{{-}1}\sim\frac{16}{9}\frac{A_0}{c_n^3} \left(\frac{2}{H}\right)^{1/2}\frac{\kappa_0^2}{\mu^3} \left\{\kappa_0\cos\frac{n{\rm \pi}}{H}z\sin\frac{z}{\mu}+ \frac{H}{n{\rm \pi}}\sin\frac{n{\rm \pi}}{H}z\cos\frac{z}{\mu}\right\}. \end{equation}

Therefore, the solution of problem (4.5), in the limit $\mu \ll 1$, schematically, takes the form

(4.8)\begin{equation} Q_{{-}1}\sim \frac{\epsilon}{\mu}\left\{A_{{-}1}^+\sin \left(\frac{z}{\mu}+\frac{n{\rm \pi}}{H}z\right)+A_{{-}1}^-\sin \left(\frac{z}{\mu}-\frac{n{\rm \pi}}{H}z\right)\right\}, \end{equation}

where $A^\pm _{-1}$ are certain ${O}(1)$ constants. Thus, $Q_{-1}={O}(\epsilon /\mu )$ and, by a similar procedure, it can be deduced that $Q_2={O}(\epsilon /\mu )$ as well.

The fact that $Q_{-1},Q_2={O}(\epsilon /\mu )$ in the joint limit $\epsilon,\mu \ll 1$ suggests that the coupling of the two resonant free modes in (3.12) with the basic wave, actually is ${O}(\epsilon /\mu )$. Hence, the assumption of weak coupling, which enables these modes to form a resonant triad with the basic wave, is valid when $\mu \gg \epsilon$ only. If this condition is violated (as will be the case for sufficiently fine-scale perturbations), all Fourier coefficients in the Floquet mode (2.8), formally, are expected to be equally important and the disturbance frequency spectrum would be broadband. A similar situation was encountered in the Floquet stability analysis of internal wave beams by Fan & Akylas (2021). The treatment of broadband instability of internal wave modes below follows along the lines of this earlier study.

5. Broadband instability

5.1. Streamline coordinates

Returning to the governing equations (2.1), the dominant coupling of the perturbations to the basic wave in the limit $\epsilon,\mu \ll 1$ derives from the Jacobian term in (2.2) which accounts for the advection due to the underlying wave velocity field. This effect can be ‘factored out’ by working with a new set of coordinates, $(x,\,z)\to (\xi,\,\zeta )$, defined by

(5.1a,b)\begin{equation} \xi=x+\frac{1}{c}\int_{}^{x}\bar{\varPsi}_z\,\mathrm{d} x^\prime,\quad\zeta=z-\frac{\bar{\varPsi}}{c}. \end{equation}

It should be noted that the curves $\zeta =\mathrm {constant}$ coincide with the streamlines of the background steady flow $(-c +\bar {\varPsi }_z,-\bar {\varPsi }_x)$, so switching to these ‘streamline coordinates’ is analogous to the change of frame used by Fan & Akylas (2021).

For uniform background stratification, in particular, according to (2.5) and (2.7a,b),

(5.2)\begin{equation} \bar{\varPsi}=\epsilon\bar{\psi}\equiv\epsilon \frac{H}{n{\rm \pi}}\sin\frac{n{\rm \pi}}{H}z\cos x, \end{equation}

so the coordinates (5.1a,b) are given by

(5.3a,b)\begin{equation} \xi = x+\frac{\epsilon}{c_n}\cos\frac{n{\rm \pi}}{H}z\sin x, \quad \zeta=z-\frac{\epsilon}{c_n}\frac{H}{n{\rm \pi}}\sin\frac{n{\rm \pi}}{H}z\cos x. \end{equation}

We remark that the basic wave streamfunction $\bar \varPsi$ is $2{\rm \pi}$-periodic in the transformed horizontal coordinate $\xi$, a property that is utilised in the following Floquet stability analysis (see § 5.2). This holds because $\bar \varPsi$ in (5.2) does not involve a term uniform in $x$; i.e. there is no (Eulerian) horizontal mean flow (under more general flow conditions where such a mean flow may be present, the definition of $\xi$ in (5.1a,b) would need to be reconsidered).

Upon implementing the transformation (5.3a,b), the advective derivative (2.2) takes the form

(5.4)\begin{equation} \frac{{\rm D}}{{\rm D}t}\to\frac{\partial}{\partial t} - c_n\frac{\partial}{\partial\xi} +\frac{\epsilon^2}{c_n}\left(\bar{\psi}_z^2- \left(\frac{n{\rm \pi}}{H}\right)^2\bar{\psi}_x^2\right) \frac{\partial}{\partial\xi}. \end{equation}

Thus, the Jacobian term in (2.2) has been eliminated correct to ${O}(\epsilon )$. Furthermore, using (5.2), the ${O}(\epsilon ^2)$ residual is expressed as

(5.5)\begin{equation} \frac{\epsilon^2}{2c_n}\left(\cos2\frac{n{\rm \pi}}{H}\zeta+\cos2\xi\right) \frac{\partial}{\partial\xi}+{O}(\epsilon^3). \end{equation}

The first term represents the advection effect due to the ‘Stokes drift’ (Thorpe 1968),

(5.6)\begin{equation} \bar{U}^s=\frac{\epsilon^2}{2c_n}\cos2\frac{n{\rm \pi}}{H}\zeta, \end{equation}

which here coincides with the Lagrangian horizontal mean flow associated with the basic wave, because the Eulerian mean flow vanishes. This effect makes an important contribution to the eigenvalue problem governing broadband instability (see § 5.3). The second term in (5.5), by contrast, is relatively insignificant and could have been eliminated by modifying via an ${O}(\epsilon ^2)$ term the definition of $\xi$ in (5.3a,b).

In terms of $\xi$ and $\zeta$, the governing equations (2.1) now read

(5.7a)\begin{align} &\frac{{\rm D}}{{\rm D}t}\nabla^2\psi - \rho_\xi + \frac{\epsilon}{c_n}\left\{\bar{\psi}_z\left(\frac{1}{c_n}\psi-\rho \right)_\xi-\bar{\psi}_x\left(\frac{1}{c_n}\psi-\rho \right)_\zeta\right\}\nonumber\\ &\quad +\frac{\epsilon^2}{c_n^3}\left(\bar{\psi}_z^2- \left(\frac{n{\rm \pi}}{H}\right)^2\bar{\psi}_x^2\right)\psi_\xi = 0, \end{align}

(5.7b)\begin{align} &\frac{{\rm D}}{{\rm D}t}\rho + \psi_\xi -\frac{\epsilon^2}{c_n^2}\left(\bar{\psi}_z^2- \left(\frac{n{\rm \pi}}{H}\right)^2\bar{\psi}_x^2\right) \psi_\xi = 0, \end{align}

where

(5.8)\begin{align} &\nabla^2 \to \left(1+\frac{\epsilon}{c_n}\bar{\psi}_z\right)^2 \frac{\partial^2}{\partial\xi^2}+\left(1- \frac{\epsilon}{c_n}\bar{\psi}_z\right)^2\frac{\partial^2}{\partial \zeta^2}\nonumber\\ &\quad +2\frac{\epsilon}{c_n}\bar{\psi}_x \left(\left(\frac{n{\rm \pi}}{H}\right)^2-1- \frac{\epsilon}{c_n^3}\bar{\psi}_z\right) \frac{\partial^2}{\partial\xi\partial\zeta}+ \frac{\epsilon}{c_n^3}\left(\bar{\psi}\frac{\partial}{\partial\zeta} +\bar{\psi}_{xz}\frac{\partial}{\partial\xi}\right)\nonumber\\ &\quad +\frac{\epsilon^2}{c_n^2}\bar{\psi}_x^2 \left(\frac{\partial^2}{\partial\zeta^2}+\left(\frac{n{\rm \pi}}{H} \right)^4\frac{\partial^2}{\partial\xi^2}\right). \end{align}

In addition, as the channel walls $z=0,H$ correspond to the streamlines $\zeta =0,H$, the boundary conditions (2.3) translate into

(5.9)\begin{equation} \psi_\xi=0\quad(\zeta=0,H). \end{equation}

5.2. Floquet stability analysis

As the coefficients of the transformed equations (5.7) are steady in $t$ and $2{\rm \pi}$-periodic in $\xi$, we seek Floquet-mode solutions similar to (2.8),

(5.10)\begin{equation} (\psi,\rho) = \mathrm{e}^{-\mathrm{i}\sigma t}\mathrm{e}^{\mathrm{i} p\xi}\sum_{m={-}\infty}^{\infty}\left(\tilde{Q}_m(\zeta), \tilde{R}_m(\zeta)\right)\mathrm{e}^{\mathrm{i}m\xi}, \end{equation}

where the Fourier coefficients $\tilde {Q}_m, \tilde {R}_m$ $(-\infty < m<\infty )$ and the eigenvalue $\sigma$ are to be determined.

Here, our interest is on short-scale perturbations $(\mu \ll 1)$ in the regime

(5.11)\begin{equation} \alpha\equiv\frac{\mu}{\epsilon}={O}(1), \end{equation}

where, as argued in § 4.2, PSI is replaced by a broadband instability. To analyse this ‘distinguished limit’, we work with the scaled wavenumber $\kappa =p\mu ={O}(1)$ defined in (4.1a,b) and the ‘stretched’ coordinate

(5.12)\begin{equation} Z=\frac{\zeta}{\mu}. \end{equation}

Furthermore, we re-scale $\tilde {Q}_m\to \mu \tilde {Q}_m\ (Z,\zeta )$ so that

(5.13)\begin{equation} \frac{\mathrm{d}\tilde{Q}_m}{\mathrm{d}\zeta}\to \frac{\partial\tilde{Q}_m}{\partial Z}+\mu\frac{\partial\tilde{Q}_m}{\partial\zeta}. \end{equation}

It should be noted that, in view of the transformation (5.3a,b), $\exp (\mathrm {i}p\xi )$ in (5.10) involves all harmonics in $x$. Moreover, for $p={O}(1/\mu )$ and $\epsilon /\mu ={O}(1)$ these harmonics contribute at the same level. Thus, in the regime (5.11) the modes (5.10) are ‘broadband’ even though, as discussed below, the $m=0,1$ components are dominant in the Fourier series in $\xi$.

Now, we derive the equations governing $\tilde {R}_m,\tilde {Q}_m$ and $\sigma$ by substituting (5.10) in (5.7) and implementing the scalings (5.11)–(5.13). Specifically, also making use of (5.5), (5.7b) yields correct to ${O}(\epsilon )$

(5.14)\begin{align} \varOmega_m\tilde{R}_m&= \kappa\tilde{Q}_m+\epsilon \left\{\alpha m+\frac{\kappa^2}{2\alpha c_n\varOmega_m}\cos2 \frac{n{\rm \pi}}{H}\zeta\right\}\tilde{Q}_m\nonumber\\ &\quad + \epsilon\frac{\kappa^2}{4\alpha c_n} \left\{\frac{\tilde{Q}_{m+2}}{\varOmega_{m+2}}+ \frac{\tilde{Q}_{m-2}}{\varOmega_{m-2}}\right\}, \end{align}

where $\varOmega _m$ is given in (2.11). Next, using (5.5), (5.8) and upon eliminating $\tilde {R}_m$ via (5.14), we obtain from (5.7a) the following equation system for $\tilde {Q}_m$ $(-\infty < m<\infty )$ correct to ${O}(\epsilon )$

(5.15) \begin{align} &\left\{\left(\frac{\partial}{\partial Z}+\alpha\epsilon \frac{\partial}{\partial\zeta}\right)^2+ (\kappa+\alpha m\epsilon)^2\left(\frac{1}{\varOmega_m^2}-1 \right)\right\}\tilde{Q}_m\nonumber\\ &\quad +\epsilon\left\{\vphantom{\left(H^+_m\frac{\partial\tilde{Q}_{m+1}}{\partial Z}-H^-_m \frac{\partial\tilde{Q}_{m-1}}{\partial Z}\right)}\cos2\frac{n{\rm \pi}}{H}\zeta G_m \tilde{Q}_m-\cos\frac{n{\rm \pi}}{H}\zeta(G^+_m\tilde{Q}_{m+1} +G^-_m\tilde{Q}_{m-1})\right.\nonumber\\ &\quad\left.+\sin\frac{n{\rm \pi}}{H}\zeta \left(H^+_m\frac{\partial\tilde{Q}_{m+1}}{\partial Z}-H^-_m \frac{\partial\tilde{Q}_{m-1}}{\partial Z}\right)+ L^+_m\tilde{Q}_{m+2}+L^-_m\tilde{Q}_{m-2}\right\}=0, \end{align}

where

(5.16a)\begin{gather} G_m= \frac{\kappa^3}{\alpha c_n\varOmega_m^3}, \end{gather}

(5.16b)\begin{gather}G_m^\pm{=} \frac{\kappa^2}{c_n}\left(2-\frac{1}{\varOmega_{m\pm1}^2}- \frac{1}{2\varOmega_m\varOmega_{m\pm1}}\right), \end{gather}

(5.16c)\begin{gather}H_m^\pm{=} \frac{\kappa}{c_n}\left(\frac{H}{n{\rm \pi}}\right) \left(\left(\frac{n{\rm \pi}}{H}\right)^2-1+\frac{1}{2\varOmega_m \varOmega_{m\pm1}}\right), \end{gather}

(5.16d)\begin{gather}L_m^\pm{=} \frac{\kappa^3}{4\alpha c_n\varOmega_m \varOmega_{m\pm2}}\left(\frac{1}{\varOmega_m}+ \frac{1}{\varOmega_{m\pm2}}\right). \end{gather}

5.3. Eigenvalue problem for $\alpha ={O}(1)$

Using the equation system (5.15), we now derive the stability eigenvalue problem appropriate to the asymptotic regime (5.11). To this end, $\kappa$ and $\varOmega _0$ are assumed to be in the vicinity of the critical values $\kappa =\kappa _c$ and $\varOmega _0=\omega _0$ in (4.3) where the triad resonance conditions (4.2) are met for $\mu \ll 1$. Accordingly, we write

(5.17)\begin{equation} \kappa = \kappa_0 + (\alpha{\rm \Delta}\kappa+s)\epsilon,\quad \varOmega_0={-}\frac{c_n}{2}+(\alpha{\rm \Delta}\omega_0+\lambda)\epsilon, \end{equation}

where $\kappa _0$, ${\rm \Delta} \kappa$ and ${\rm \Delta} \omega _0$ are given in (4.3b), $s={O}(1)$ is a real wavenumber detuning and $\lambda$ is a generally complex eigenvalue to be determined. It should be noted that, because $\mu =\alpha \epsilon$ according to (5.11), the ${O}(\alpha \epsilon )$ terms in (5.17) are the ${O}(\mu )$ corrections to $\kappa _c$ and $\omega _0$ in (4.3a). The eigenvalue $\lambda$ hinges on the detuning $s$ and, more importantly, the resonant interaction of perturbations with the basic wave, which can cause instability.

For $\kappa$, $\varOmega _0$ and $\varOmega _1=\varOmega _0+c_n$ in keeping with (5.17), the solution of (5.15) consistent with the boundary conditions (5.9) takes the form

(5.18a)\begin{gather} \tilde{Q}_0=\sum_{r={-}\infty}^\infty A_{0,r}\sin \left(Z+\frac{r{\rm \pi}}{H}\zeta\right)+{O}(\epsilon^2), \end{gather}

(5.18b)\begin{gather}\tilde{Q}_1=\sum_{r={-}\infty}^\infty A_{1,r}\sin \left(Z+\frac{r{\rm \pi}}{H}\zeta\right)+{O}(\epsilon^2), \end{gather}

with the rest of $\tilde {Q}_m(m\neq 0,1){O}(\epsilon )$ or smaller. The coefficients $A_{0,r}$ and $A_{1,r}$ $(-\infty < r<\infty )$ are determined by substituting (5.18) in (5.15) for $m=0,1$ and collecting terms proportional to $\sin (Z+({r{\rm \pi} }/{H})\zeta )$ correct to ${O}(\epsilon )$. Specifically, making also use of (5.16) and (5.17), we find

(5.19a)\begin{align} &\left\{\lambda\pm s(1-c_n^2/4)^{3/2}-\frac{\alpha}{2} \frac{r{\rm \pi}}{H}c_n(1-c_n^2/4)\right\}A_0,r \nonumber\\ &\quad = \tilde{E}^\pm A_{1,r+n}+\tilde{E}^\mp A_{1,r-n}\pm \frac{\tilde{D}}{\alpha}(A_{0,r-2n}+A_{0,r+2n}), \end{align}

(5.19b)\begin{align} &\left\{\lambda\mp s(1-c_n^2/4)^{3/2}+\frac{\alpha}{2} \frac{(r-n){\rm \pi}}{H}c_n(1-c_n^2/4)\right\}A_1,r \nonumber\\ &\quad ={-}\tilde{E}^\mp A_{0,r+n}-\tilde{E}^\pm A_{0,r-n}\pm \frac{\tilde{D}}{\alpha}(A_{1,r-2n}+A_{1,r+2n}). \end{align}

Here

(5.20a)\begin{gather} \tilde{E}^\pm{=}{-}\frac{1}{16(1-{c_n}^2)^{1/2}} \{(1-c_n^2)^{3/2}\pm(2c_n^2+1)(1-c_n^2/4)^{1/2}\}, \end{gather}

(5.20b)\begin{gather}\tilde{D}= \frac{1}{8(1-c_n^2/4)^{1/2}}, \end{gather}

where the upper (lower) sign in (5.19) and (5.20) corresponds to the positive (negative) value of $\kappa _0$ in (4.3b).

The equation system (5.19) is the desired stability eigenvalue problem for short-scale perturbations $(\mu \ll 1)$ in the broadband regime $\mu ={O}(\epsilon )$. This problem, in contrast to the $2\times 2$ system (3.18) that governs TRI, formally involves infinite number of mode amplitudes $A_{0,r}$ and $A_{1,r}$ ($-\infty < r<\infty$): when $\mu ={O}(\epsilon )$ all modes with wavenumber $k_0\sim \kappa _0/\mu (k_1=k_0+1)$ and modal number $(H/{\rm \pi} )/\mu +r$ are nearly resonant and participate in the interaction with the basic wave. Of particular note are the interaction terms proportional to $\tilde {D}/\alpha$ in the system (5.19); in view of (5.20b) and (4.3b), $\pm \tilde {D}/\alpha =k_0\epsilon /4c_n$ so these terms arise from the ‘Doppler shift’ $k_0\bar {U}^s$ of the perturbations by the Stokes drift $\bar {U}^s$ in (5.6).

Finally, as expected when $\epsilon \ll \mu \ll 1(\alpha \gg 1)$, the broadband instability eigenvalue problem (5.19) reduces to the PSI limit of the TRI eigenvalue problem (3.18). Specifically, in the limit $\alpha \gg 1$, $A_{0,0}$ and $A_{1,n}$ dominate the rest of the amplitudes, so (5.19) simplifies to

(5.21a)\begin{gather} \{\lambda\pm s(1-c_n^2/4)^{3/2}\}A_{0,0}=\tilde{E}^\pm A_{1,n}, \end{gather}

(5.21b)\begin{gather}\{\lambda\mp s(1-c_n^2/4)^{3/2}\}A_{1,n}={-}\tilde{E}^\pm A_{0,0}. \end{gather}

Returning to (5.20a) and noting that $\tilde {E}^\pm$ and $-\tilde {E}^\pm$ match the asymptotic expressions (4.4) for the TRI interaction coefficients $E_1$ and $E_2$, respectively, the $2\times 2$ eigenvalue problem (5.21) agrees with (3.18) in the PSI limit.

6. Numerical results

Here we compare the theoretical predictions for TRI, PSI and broadband instability with numerical results from solving the full Floquet eigenvalue problem for the $n=1$ wave mode in a uniformly stratified $(N=1)$ fluid layer. Having in mind the oceanic internal tide, we choose the horizontal length scale $\lambda _*/2{\rm \pi} =20$ km, the background buoyancy frequency $N_*=2 \times 10^{-3}\ {\rm s}^{-1}$ and the fluid depth as 4 km so the dimensionless depth $H=0.2$. It should be noted though that the assumption of constant $N$, made here for analytical convenience, is not realistic for the oceans and also that our analysis ignores the Earth's rotation.

The differential equation system (2.9) for the Floquet modes (2.8) was tackled by expanding $Q_m(z)$ in Fourier sine series in $0< z< H$ consistent with the boundary conditions (2.12)

(6.1)\begin{equation} Q_m(z) = \sum_{j=1}^\infty B_{m,j}\sin\frac{j{\rm \pi}}{H}z. \end{equation}

Thus, (2.9) and (2.12) reduce to an algebraic eigenvalue problem for $B_{m,j} (-\infty < m<\infty, j\geq 1)$ and $\sigma$. After truncating to a finite, but large enough to ensure convergence, number of Fourier modes, this problem was solved using standard MATLAB algorithms. The resolution used typically involved 10 modes in $x$ and 20 modes in $z$.

The results below are for the critical wavenumber $p=p_c$, computed from the triad resonance conditions (3.29) as a function of modal number $l$ (and $n=1$). It should be noted that (3.29) determine two solution branches $p_c(l)$ in which $p_c>0$ or $p_c<0$. Here, we report on $p_c>0$, which features a higher TRI growth rate (3.23). Furthermore, this choice of $p$ provides the dominant instability (highest growth rate $\sigma _\mathrm {i}$) for the values of the amplitude parameter $10^{-4}\leq \epsilon \leq 6\times 10^{-3}$ used in our computations. (This range of $\epsilon$ is below the threshold $\epsilon _c=6.35\times 10^{-2}$ for overturning of the mode-1 basic state.)

Figure 1 compares the instability growth rates obtained from the TRI stability problem (3.18) for $3\leq l\leq 30$ (and $\hat {p}=0$) with those computed from the Floquet eigenvalue problem for the same $l$ and $\epsilon =10^{-4}$. TRI provides an excellent approximation near the onset of instability ($\epsilon \ll 1$), and for $l\gtrsim 10$ the TRI growth rate (3.23) is already very close to the PSI limit (4.4).

Figure 1. Comparison of TRI instability growth rate ($\times$) as a function of modal number $l$ with numerical results ($\circ$) from the full Floquet stability problem for $\epsilon =10^{-4}$. The dotted line indicates the PSI limit (4.4).

It was argued in § 5 that in the short-scale limit ($\mu =H/({\rm \pi} l)\ll 1$) PSI applies if $\alpha =\mu /\epsilon \gg 1$, but when $\alpha ={O}(1)$ it is replaced by broadband instability. As a check of this theoretical prediction, figure 2 plots as a function of $0<\alpha <5$ the growth rate predicted by the eigenvalue problem (5.19) (with $s=0$ and the upper sign which applies to $\kappa _0>0$) that pertains to PSI and broadband instability, together with numerical results computed from the Floquet problem for the same range of $\alpha$ and various $\epsilon$. Specifically, for $\epsilon =10^{-4},10^{-3}$ and $3\times 10^{-3}$, the computed growth rates are well approximated by PSI when $\alpha \gtrsim 2$, but for $\alpha$ less than about $1.5$ the growth rate exhibits a sharp drop and the instability is severely suppressed. Furthermore, for these $\epsilon$ the theoretical predictions based on (5.19) are in good quantitative agreement with the numerical computations. For the relatively larger value of $\epsilon =6\times 10^{-3}$, when $\alpha$ is decreased the growth rate behaves in a similar manner as for the smaller $\epsilon$, but there is only qualitative agreement between theoretical and numerical results.

Figure 2. Comparison of theoretical instability growth rate (—) based on the eigenvalue problem (5.19), as a function of $\alpha =\mu /\epsilon$, with numerical results from the full Floquet stability problem for $\epsilon =10^{-4}\ (\circ ),\,10^{-3}\ (\square ),\,3\times 10^{-3}\ (\diamond )$ and $6\times 10^{-3}\ (\triangle )$. The dotted line indicates the PSI limit (4.4).

The stabilisation of PSI in figure 2 is caused by the terms $\pm \tilde {D}/\alpha$ in the eigenvalue problem (5.19) that represent the effect of the Stokes drift of the basic wave (cf. (5.6)). This becomes apparent from the relative magnitudes of the interaction coefficients $\tilde {E}^\pm$ and $\tilde {D}$ in (5.19). Specifically, from (2.7a,b), the mode-1 speed $c_1=0.0635$ so, according to (5.20), $\tilde {E}^-=7.89\times 10^{-4}$ is much smaller than $\tilde {E}^+=-1.25\times 10^{-1}$ and $\tilde {D}=1.25\times 10^{-1}$. As a result, if the effect of the Stokes drift is ignored (by setting $\tilde {D}=0$ in (5.19)), $A_{0,r}$ and $A_{1,r+1}$ practically are coupled to each other only. The multi-mode resonance interaction thus degenerates to a set of (essentially uncoupled) resonant triads ($-\infty < r<\infty$) and the dominant instability arises for $r=0$, which recovers PSI.

For the dimensional scales chosen here, the (dimensionless) maximum instability growth rate $\sigma _{\mathrm {i}}\approx 0.12\epsilon$ in figure 2 translates to an e-folding time of roughly $5\times 10^{-2}/\epsilon$ days. As an example, taking $\epsilon =3\times 10^{-3}$, which corresponds to $U_*=12$ cm s$^{-1}$ for the peak horizontal velocity of the mode-1 basic state, this e-folding time is 16 days: about twice the estimate found by Young et al. (2008) for near-inertial PSI.

As discussed in § 5, the transition from PSI to broadband instability is associated with the broadening of the Floquet mode spectrum as the disturbance scale controlled by $\mu$ is decreased for given $\epsilon$. This is illustrated in figure 3, which shows the relative magnitudes of the Fourier coefficients $B_{m,j}$ in (6.1) for $\epsilon =10^{-3}$ and three values of $l=3$ ($\mu =2.12\times 10^{-2},\ \alpha =21.2$), $l=38$ ($\mu =1.68\times 10^{-3},\ \alpha =1.68$) and $l=100$ ($\mu =6.37\times 10^{-4},\ \alpha =0.637$). For $\alpha =21.2$ (figure 3a), $B_{0,l}$ and $B_{1,l+1}$ are clearly dominant, as expected in TRI. For $\alpha =1.68$ (figure 3b), $B_{0,l}$ and $B_{1,l+1}$ are still dominant in keeping with PSI, but the neighbouring Fourier coefficients $B_{-1,l+1}$ and $B_{2,l}$ associated with the frequency components $\pm 3c_n/2$ (in the rest frame) as well as $B_{0,l+2}$ and $B_{1,l-1}$ are starting to gain strength. Finally, for $\alpha =0.637$ (figure 3c), the PSI assumption is no longer valid as several Fourier coefficients are of comparable magnitude to $B_{0,l}$ and $B_{1,l+1}$; this transition to a ‘broadband’ spectrum is accompanied by a significantly reduced growth rate relative to PSI for $\alpha \lesssim 1$ (figure 2).

Figure 3. Relative magnitudes of Fourier coefficients $B_{m,j}$ in (6.1), normalised by the coefficient of largest magnitude, for $\epsilon =10^{-3}$: (a) $\alpha =21.2$ ($l=3$); (b) $\alpha =1.68$ ($l=38$); (c) $\alpha =0.637$ ($l=100$).

7. Concluding remarks

We have made a systematic stability analysis of internal gravity wave modes in a stratified fluid layer bounded by rigid walls. The temporal stability of Floquet modes is governed by an eigenvalue problem that involves an infinite system of differential equations subject to inviscid conditions on the walls. Examining this problem in the limit of small basic-state amplitude $(\epsilon \ll 1)$ shows that the onset of instability is triggered by perturbations that form resonant triads with the underlying wave mode, and the associated ${O}(\epsilon )$ growth rate is determined by the $2\times 2$ eigenvalue problem (3.18). In general, the resonant triad conditions (cf. (3.10a,b)) require that the (horizontal) wavevectors and the frequencies of the perturbations sum up to the wavevector and frequency of the underlying wave mode, consistent with Thorpe (1966). The case of uniform background stratification (constant $N$) is exceptional, as the modal numbers of the perturbations also need to satisfy the constraint (3.27) for TRI to be possible.

A particular case of TRI, where resonant triads comprise fine spatial-scale perturbations with half the basic wave frequency, is the so-called PSI. Owing to its potential geophysical significance, PSI has attracted considerable interest in the context of sinusoidal plane waves and finite-width beams in an unbounded, uniformly stratified fluid. In an effort to understand the role of PSI for propagating internal wave modes in a waveguide setting, we have studied the Floquet eigenvalue problem for a small-amplitude basic wave mode $(\epsilon \ll 1)$ subject to short-scale $(\mu \ll 1)$ disturbances, assuming for simplicity constant $N$ background stratification. Our analysis reveals that the nature of the instability mechanism in this joint limit hinges on the perturbation scale, controlled by $\mu$, relative to the basic-state amplitude $\epsilon$: PSI applies only when $\epsilon \ll \mu \ll 1$ ($\alpha = \mu /\epsilon \gg 1$); as $\mu$ is further decreased for fixed $\epsilon$, higher-frequency perturbations than the two subharmonics at half the basic-wave frequency come into play, and when $\alpha ={O}(1)$ Floquet modes feature broadband spectrum.

A similar situation was encountered in a recent Floquet stability analysis of finite-width wave beams (Fan & Akylas 2021), which confirmed an earlier claim (Onuki & Tanaka 2019) that the broadening of the Floquet-mode spectrum is due to the advection of the perturbation by the underlying wave beam. Furthermore, by riding on a frame moving with the beam velocity field, Fan & Akylas (2021) ‘factored out’ this advection effect and revealed a novel instability which features broadband frequency spectrum. Following an analogous approach, switching to the ‘streamline coordinates’ (5.1a,b) enabled us to factor out the advection due to the underlying-mode velocity field and obtain the eigenvalue problem (5.19) which pertains to the broadband regime $\alpha ={O}(1)$.

Unlike the broadband instability of a wave beam, which is of the resonant triad type after the advection effect has been removed, the instability mechanism found here for $\alpha ={O}(1)$ is a multi-mode resonance. This fundamental difference is reflected in the eigenvalue problem (5.19), which involves infinite number of mode amplitudes. In particular, the interaction terms in (5.19) that account for the effects of the ${O}(\epsilon ^2)$ Lagrangian mean flow due to the Stokes drift (5.6) of the basic mode, are responsible for the sharp drop of the instability growth rate below the PSI limit when $\alpha$ is less than about $1.5$ (figure 2).

Based on the results presented in § 6, the broadening of the Floquet-mode spectrum as $\alpha$ is decreased (figure 3) has a strong stabilising effect that provides a short-scale cut-off to PSI. A similar cut-off effect for $\mu \ll 1$ would be expected due to viscous dissipation, given that the viscous decay rate of internal waves is ${O}(\nu /\mu ^2)$ where $\nu$ is the inverse Reynolds number (see, e.g., Lighthill 1978). Although viscous effects would be dominant in a laboratory setting, the inviscid mechanism discussed here would prevail in a nearly inviscid environment where $\nu /\mu ^2\ll \epsilon$; i.e. $\nu \ll \epsilon ^3$ for $\alpha ={O}(1)$. Assuming a kinematic viscosity $\nu _*=10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$ and taking $\epsilon =10^{-3}$, this condition is met for the oceanic scales chosen in § 6.