2.1. Formulation in the physical plane

Consider the periodic and irrotational plane motion of steady gravity waves progressing to the left in permanent form with constant speed $c$ at the interface between two unbounded homogeneous fluids with the different densities $\rho _1$ and $\rho _2$ ($\rho _1 < \rho _2$), as shown in figure 1(a), where subscripts $1$ and $2$ refer to the upper and lower fluid domains, respectively. It is convenient to formulate the problem in the flow domains for one period with wavelength $\lambda$ of the steady waves, namely the two domains surrounded by $\mathrm {A}_j \mathrm {A} \mathrm {C} \mathrm {B} \mathrm {B}_j$ ($\kern0.09em j = 1, 2$) in figure 1(a), in the frame of reference moving to the left with the waves. The steady wave profiles are assumed to be symmetric with respect to the vertical line passing through the wave crest. It is also assumed that the fluids are incompressible and inviscid, and that the fluid motion in each layer is two-dimensional and irrotational in the vertical cross-section $(x, y)$ along the propagation direction of the waves. Surface tension at the interface is neglected. The origin is placed such that the wave profile $y = \tilde {Y}(x,t)$ satisfies the zero mean level condition

(2.1)\begin{equation} \displaystyle{\int_0^{\lambda}} \tilde{Y}(x,t) \, \mathrm{d}\kern0.06em x = 0, \end{equation}

where $t$ denotes the time.

For irrotational plane flows, we can introduce the complex coordinate $z = x + \mathrm {i} y$ and the complex velocity potential $f_j = \phi _j + \mathrm {i} \psi _j$ ($\kern0.09em j = 1, 2$), where $\phi _j$ and $\psi _j$ denote the velocity potential and the streamfunction, respectively, for the $\kern0.09em j$th layer. When the physical variables are non-dimensionalized, with respect to $c$ and $k$ (${=}2{\rm \pi} /\lambda$), as

(2.2)\begin{equation} z_\ast = k z , \quad \tilde{Y}_\ast = k \tilde{Y} , \quad t_\ast = c k t , \quad f_{j \ast} = \displaystyle{\frac{f_j}{c/k}} , \quad p_{j \ast} = \displaystyle{\frac{p_j}{\rho_1 c^{2}}} \quad (\kern0.09em j = 1,2) , \end{equation}

the following dimensionless parameters appear in the problem:

(2.3a,b)\begin{equation} h_\ast = k h , \quad c_\ast = \displaystyle{\frac{c}{\sqrt{g/k}}} , \end{equation}

where $p_j$ denote the pressure, $h$ the crest-to-trough wave height and $g$ the gravitational acceleration. Note that $h_\ast = kh = 2{\rm \pi} h/\lambda$ is the wave steepness. The asterisks for dimensionless variables will be omitted hereinafter for brevity.

The complex velocity potentials $f_j = f_j(z,t)$ ($\kern0.09em j=1,2$) are both analytic in each layer. The boundary conditions at the interface $y = \tilde {Y}(x,t)$ are given, from the kinematic conditions, by

(2.4)\begin{equation} \tilde{Y}_t + \phi_{j, x} \tilde{Y}_x = \phi_{j, y} \quad (\kern0.09em j = 1,2)\quad \mbox{at} \ y = \tilde{Y}(x,t), \end{equation}

and, from the dynamic condition $p_1 = p_2$ at the interface $y = \tilde {Y}(x,t)$ with Bernoulli's equation for each layer, by

(2.5)\begin{align} &\phi_{1, t} + \displaystyle{\frac{1}{2}} ( {\phi_{1, x}}^{2} + {\phi_{1, y}}^{2} ) + \displaystyle{\frac{1}{c^{2}}} \tilde{Y} - B(t) \nonumber\\ &\quad = \displaystyle{\frac{\rho_2}{\rho_1}} \left\{ \phi_{2, t} + \displaystyle{\frac{1}{2}} ( {\phi_{2, x}}^{2} + {\phi_{2, y}}^{2} ) + \displaystyle{\frac{1}{c^{2}}} \tilde{Y} \right\} \quad \mbox{at} \ y = \tilde{Y}(x,t) , \end{align}

where an arbitrary function $B(t)$ can be absorbed into the velocity potentials and the subscripts denote the partial derivatives with respect to $x$, $y$ and $t$ as

(2.6)\begin{align} \tilde{Y}_{x} := \displaystyle{\frac{\partial\tilde{Y}}{\partial x}}, \quad \tilde{Y}_{t} := \displaystyle{\frac{\partial\tilde{Y}}{\partial t}}, \quad \phi_{j, x} := \displaystyle{\frac{\partial\phi_j}{\partial x}}, \quad \phi_{j, y} := \displaystyle{\frac{\partial\phi_j}{\partial y}}\quad \mbox{and} \quad \phi_{j, t} := \displaystyle{\frac{\partial\phi_j}{\partial t}}\quad (\kern0.09em j = 1, 2). \end{align}

To avoid any confusion, a comma is introduced to denote a derivative of an indexed function with respect to the variable following the comma.

It should be remarked that, as a velocity jump across the interface is allowed under the inviscid assumption, the velocity potential is discontinuous across the interface while the streamfunction is continuous so that

(2.7)\begin{equation} \psi_1(x,y,t) = \psi_2(x,y,t)\quad \mbox{at} \ y = \tilde{Y}(x,t). \end{equation}

In addition, the conditions that the flow is uniform far above and below the interface ($y \to \pm \infty$) yield

(2.8)\begin{equation} w_1 := \displaystyle{\frac{\partial f_1}{\partial z}} \to 1 \quad \mbox{as} \ y \to \infty \quad \mbox{and} \quad w_2 := \displaystyle{\frac{\partial f_2}{\partial z}} \to 1 \quad \mbox{as} \ y \to -\infty , \end{equation}

where $w_j = u_j - \mathrm {i} v_j$ ($\kern0.09em j = 1, 2$) denote the complex velocity.

2.2. Unsteady conformal mapping of the flow domains

As described in § 1, in order to overcome numerical difficulties for the two-dimensional linear stability analysis of finite-amplitude interfacial waves, we generalize the unsteady conformal mapping technique used for the stability of surface capillary and gravity waves (Tiron & Choi Reference Tiron and Choi2012; Murashige & Choi Reference Murashige and Choi2020) to the two-layer problem. Here, we map the upper flow domain ($y > \tilde {Y}(x,t$)) and the lower flow domain ($y < \tilde {Y}(x,t$)) into the upper half of the $\zeta _1$-plane ($\zeta _1 = \xi _1 + \mathrm {i} \eta _1$ with $\eta _1 > 0$) and the lower half of the $\zeta _2$-plane ($\zeta _2 = \xi _2 + \mathrm {i} \eta _2$ with $\eta _2 < 0$), respectively, as shown in figures 1(b) and 1(c). The time-varying interface $y = \tilde {Y}(x,t)$ is always mapped onto the real axes $\eta _j = 0$ in the $\zeta _j$-plane ($\kern0.09em j = 1, 2$) although the corresponding physical locations on the interface are different even when $\xi _1=\xi _2$.

The complex coordinates $z_j=z(\zeta _j, t) = x_j(\xi _j, \eta _j, t) + \mathrm {i} y_j(\xi _j, \eta _j, t)$ and the complex velocity potentials $f_j(\zeta _j, t) = \phi _j(\xi _j, \eta _j, t) + \mathrm {i} \psi _j(\xi _j, \eta _j, t)$ are analytic in the $\zeta _j$-plane ($\kern0.09em j = 1, 2$). We write $x_j$, $y_j$, $\phi _j$ and $\psi _j$ at the interface $\eta _j = 0$ ($\kern0.09em j = 1, 2$) as

(2.9)\begin{equation} \left.\begin{array}{ll@{}} \tilde{x}_j(\xi_j, t) = x_j(\xi_j, \eta_j = 0, t) , & \tilde{y}_j(\xi_j, t) = y_j(\xi_j, \eta_j = 0, t) \\ \tilde{\phi}_j(\xi_j, t) = \phi_j(\xi_j, \eta_j = 0, t), & \tilde{\psi}_j(\xi_j, t) = \psi_j(\xi_j, \eta_j = 0, t) \end{array}\right\}\quad (\kern0.09em j = 1, 2). \end{equation}

Note that, at the interface, $\xi _j$ ($\kern0.09em j=1,2$) range from $-{\rm \pi}$ to ${\rm \pi}$ for one wavelength, as shown in figures 1(b) and 1(c). We also introduce the following shorthand notation for partial derivatives with respect to $\xi _j$ and $t$, respectively, as

(2.10)\begin{equation} x_{j,j} := \displaystyle{\frac{\partial x_j}{\partial \xi_j}}\quad \mbox{and} \quad x_{j,t} := \displaystyle{\frac{\partial x_j}{\partial t}} \quad (\kern0.09em j = 1, 2) . \end{equation}

Then, using the chain rule, it can be shown that the derivatives in the $z$-plane can be transformed to those in the $\zeta _j$-plane ($\kern0.09em j=1,2$) as

(2.11)$$\begin{gather} \phi_{j, x} = \displaystyle{\frac{1}{J_j}} \displaystyle{\frac{\partial(\phi_j, y_j)}{\partial(\xi_j, \eta_j)}}, \quad \phi_{j, y} = \displaystyle{\frac{1}{J_j}} \displaystyle{\frac{\partial(x_j, \phi_j)}{\partial(\xi_j, \eta_j)}}, \quad \phi_{j, t} = \displaystyle{\frac{1}{J_j}} \displaystyle{\frac{\partial(x_j, y_j, \phi_j)}{\partial(\xi_j, \eta_j, t)}} \quad (\kern0.09em j = 1, 2), \end{gather}$$

(2.12)$$\begin{gather}\tilde{Y}_x = \displaystyle{\frac{\tilde{y}_{j, j}}{\tilde{x}_{j, j}}} \quad \mbox{and} \quad \tilde{Y}_t = \displaystyle{\frac{1}{\tilde{x}_{j, j}}} \displaystyle{\frac{\partial(\tilde{x}_j, \tilde{y}_j)}{\partial(\xi_j, t)}} \quad (\kern0.09em j = 1, 2) , \end{gather}$$

where $J_j$ ($\kern0.09em j = 1, 2$) are defined by

(2.13)\begin{equation} J_j := {x_{j,j}}{}^{2} + {y_{j,j}}{}^{2} \quad (\kern0.09em j = 1, 2) . \end{equation}

In order to match at the interface $\eta _j = 0$ ($\kern0.09em j=1,2$) the upper- and lower-layer solutions, we introduce a subsidiary real variable $\hat {\xi }$ along the interface and write $\xi _1$ and $\xi _2$ at the interface as

(2.14)\begin{equation} \xi_1(\hat{\xi}, t) = \hat{\xi} + \gamma(\hat{\xi},t) \quad \mbox{and} \quad \xi_2(\hat{\xi}, t) = \hat{\xi} - \gamma(\hat{\xi},t) \quad \mbox{for} \ -{\rm \pi} \leqslant \hat{\xi} \leqslant {\rm \pi}, \end{equation}

with

(2.15)\begin{equation} \xi_j(\hat{\xi} ={\pm}{\rm \pi}, t) ={\pm}{\rm \pi} \quad (\kern0.09em j = 1, 2) \quad \mbox{and} \quad \gamma(\hat{\xi} ={\pm}{\rm \pi}, t) = 0 \quad \mbox{for} \ t \geqslant 0 , \end{equation}

where $\gamma (\hat {\xi },t)$ is an unknown real function. Note that similar parametric representations of the interface have been previously adopted for the numerical computation of steady interfacial waves (Saffman & Yuen Reference Saffman and Yuen1982; Părău & Dias Reference Părău and Dias2001) and unsteady interfacial waves (Baker, Meiron & Orszag Reference Baker, Meiron and Orszag1982; Grue et al. Reference Grue, Friis, Palm and Rusås1997).

Using $\hat {\xi }$ and $\gamma$ defined by (2.14), the conditions that the wave profiles in the $\zeta _1$- and $\zeta _2$-planes coincide are given by

(2.16)\begin{equation} \left.\begin{array}{l@{}} \tilde{x}_1(\xi_1=\xi_1(\hat{\xi},t), t) = \tilde{x}_2(\xi_2=\xi_2(\hat{\xi},t), t)\\ \tilde{y}_1(\xi_1=\xi_1(\hat{\xi},t), t) = \tilde{y}_2(\xi_2=\xi_2(\hat{\xi},t), t) \end{array} \right\} . \end{equation}

The conditions in (2.16) are referred to as the contact conditions at the interface in this work. Also the continuity condition (2.7) of the streamfunctions at the interface can be expressed as

(2.17)\begin{equation} \tilde{\psi}_1(\xi_1=\xi_1(\hat{\xi},t), t) = \tilde{\psi}_2(\xi_2=\xi_2(\hat{\xi},t), t) . \end{equation}

Substituting (2.11) and (2.12) with the Cauchy–Riemann relations into (2.4) and (2.5), we obtain the kinematic conditions and the dynamic condition at the interface $\eta _j = 0$ in the $\zeta _j$-plane ($\kern0.09em j=1,2$), respectively, as

(2.18)\begin{equation} \tilde{x}_{j,j} \tilde{y}_{j,t} - \tilde{y}_{j,j} \tilde{x}_{j,t} ={-}\tilde{\psi}_{j,j} \quad (\kern0.09em j = 1, 2) , \end{equation}

and

(2.19)\begin{align} &\tilde{\phi}_{2,t} - \displaystyle{\frac{1}{J_2}} (\tilde{x}_{2,2} \tilde{x}_{2,t} + \tilde{y}_{2,2} \tilde{y}_{2,t}) \tilde{\phi}_{2,2} + \displaystyle{\frac{1}{2}} \displaystyle{\frac{1}{J_2}} ( {\tilde{\phi}_{2,2}}{}^{2} - {\tilde{\psi}_{2,2}}{}^{2} ) + \displaystyle{\frac{1}{c^{2}}} \tilde{y}_2- B(t)\nonumber\\ &\quad = \displaystyle{\frac{\rho_1}{\rho_2}} \left\{ \tilde{\phi}_{1,t} - \displaystyle{\frac{1}{J_1}} (\tilde{x}_{1,1} \tilde{x}_{1,t} + \tilde{y}_{1,1} \tilde{y}_{1,t}) \tilde{\phi}_{1,1} + \displaystyle{\frac{1}{2}} \displaystyle{\frac{1}{J_1}} ( {\tilde{\phi}_{1,1}}{}^{2} - {\tilde{\psi}_{1,1}}{}^{2} ) + \displaystyle{\frac{1}{c^{2}}} \tilde{y}_1 \right\} , \end{align}

where the notation for the derivatives in (2.10) has been used. Here, note that $\xi _j = \xi _j(\hat {\xi }, t)$ at the interface $\eta _j = 0$ ($\kern0.09em j = 1, 2$) are given by (2.14). The conditions at infinities (2.8) are satisfied with

(2.20)\begin{equation} \left.\begin{array}{l@{}} z_1 \to \zeta_1\quad \mbox{and} \quad f_1 \to \zeta_1\quad \mbox{as} \ \eta_1 \to \infty\\ z_2 \to \zeta_2\quad \mbox{and} \quad f_2 \to \zeta_2\quad \mbox{as} \ \eta_2 \to -\infty \end{array} \right\}. \end{equation}

The analyticity of $z_j = x_j + \mathrm {i} y_j$ and $f_j = \phi _j + \mathrm {i} \psi _j$ ($\kern0.09em j = 1,2$) with (2.20) yields the relations between the real and imaginary parts at the interface as

(2.21a–d)\begin{align} \tilde{x}_1 - \xi_1 = \mathcal{H}_1[\tilde{y}_1] , \quad \tilde{\psi}_1 ={-}\mathcal{H}_1[\tilde{\phi}_1 - \xi_1] , \quad \tilde{x}_2 - \xi_2 ={-}\mathcal{H}_2[\tilde{y}_2] , \quad \tilde{\psi}_2 = \mathcal{H}_2[\tilde{\phi}_2 - \xi_2] , \end{align}

where $\mathcal {H}_j[F]$ is the Hilbert transform for a real-valued function $F = F(\xi _j)$ in the $\zeta _j$-plane ($\kern0.09em j= 1, 2$) defined by

(2.22)\begin{equation} \mathcal{H}_j[F](\xi_j) := \displaystyle{\frac{1}{\rm \pi}} \mathrm{P.V.} \displaystyle{\int_{-\infty}^{\infty}} \displaystyle{\frac{F(\xi_j^{\prime})}{\xi_j^{\prime} - \xi_j}} \, \mathrm{d}\xi_j^{\prime} \quad (\kern0.09em j = 1, 2), \end{equation}

with $\mathrm {P.V.}$ denoting Cauchy's principal value. Thus the governing equations (2.16), (2.17), (2.18) ($\kern0.09em j = 1$ or $2$) and (2.19) determine the five unknown variables $\tilde {y}_j$, $\tilde {\phi }_j$ ($\kern0.09em j = 1, 2$) and $\gamma$.

4.1. Linearization around steady wave solutions in the $\zeta _1$- and $\zeta _2$-planes

For linear stability analysis, we superimpose time-dependent small-amplitude disturbances $z_j^{(1)}(\zeta _j, t)$, $f_j^{(1)}(\zeta _j, t)$ and $\gamma ^{(1)}(\hat {\xi },t)$ on the steady wave solutions $z_j^{(0)}(\zeta _j)$, $f_j^{(0)}(\zeta _j)$ and $\gamma ^{(0)}(\hat {\xi })$, respectively, as

(4.1)\begin{align} \left.\begin{array}{l@{}} z_j(\zeta_j, t) = z_j^{(0)}(\zeta_j) + z_j^{(1)}(\zeta_j, t) \\ f_j(\zeta_j, t) = f_j^{(0)}(\zeta_j) + f_j^{(1)}(\zeta_j, t) \end{array} \right\} \quad (\kern0.09em j = 1, 2) \quad \mbox{and} \quad \gamma(\hat{\xi}, t) = \gamma^{(0)}(\hat{\xi}) + \gamma^{(1)}(\hat{\xi},t) , \end{align}

and linearize the governing equations around the steady wave solutions. Note that $z_j^{(1)}(\zeta _j, t)$ and $f_j^{(1)}(\zeta _j, t)$ are both analytic in the $\zeta _j$-plane ($\kern0.09em j = 1, 2$). Similarly to (2.9), we write the disturbances at the interface as

(4.2)\begin{equation} \left.\begin{array}{ll@{}} \tilde{x}_j^{(1)}(\xi_j, t) = x_j^{(1)}(\xi_j, \eta_j = 0, t), & \tilde{y}_j^{(1)}(\xi_j, t) = y_j^{(1)}(\xi_j, \eta_j = 0, t)\\ \tilde{\phi}_j^{(1)}(\xi_j, t) = \phi_j^{(1)}(\xi_j, \eta_j = 0, t), & \tilde{\psi}_j^{(1)}(\xi_j, t) = \psi_j^{(1)}(\xi_j, \eta_j = 0, t) \end{array}\right\}\quad (\kern0.09em j=1,2). \end{equation}

We look for solutions of the linearized equations in the form of

(4.3)\begin{equation} \left.\begin{array}{ll@{}} \tilde{x}_j^{(1)}(\xi_j, t) = \mathrm{e}^{\sigma t} \check{x}_j^{(1)}(\xi_j), & \tilde{y}_j^{(1)}(\xi_j, t) = \mathrm{e}^{\sigma t} \check{y}_j^{(1)}(\xi_j)\\ \tilde{\phi}_j^{(1)}(\xi_j, t) = \mathrm{e}^{\sigma t} \check{\phi}_j^{(1)}(\xi_j), & \tilde{\psi}_j^{(1)}(\xi_j, t) = \mathrm{e}^{\sigma t} \check{\psi}_j^{(1)}(\xi_j) \end{array}\right\}\quad (\kern0.09em j = 1, 2), \end{equation}

and

(4.4)\begin{equation} \gamma^{(1)}(\hat{\xi}, t) = \mathrm{e}^{\sigma t} \check{\gamma}^{(1)}(\hat{\xi}) , \end{equation}

where $\sigma = \sigma _{{r}} + \mathrm {i} \sigma _{{i}} \in \mathbb {C}$. The steady interfacial waves are linearly unstable if and only if there exists a positive real part $\sigma _{{r}} > 0$. Notice that this condition is equivalent to the existence of any real part as both $\sigma$ and $-\bar {\sigma }$ are eigenvalues to be discussed later, where $\bar {\sigma }$ denotes the complex conjugate of $\sigma$.

Substituting these into the governing equations (2.18) for $\kern0.09em j = 1$, (2.19), (2.16) and (2.17), and expanding them around the steady wave solutions and $\xi _j = \xi _j^{(0)}(\hat {\xi })$ ($\kern0.09em j = 1, 2$), we can obtain the linearized equations for $\check {x}_j^{(1)}$, $\check {y}_j^{(1)}$, $\check {\phi }_j^{(1)}$, $\check {\psi }_j^{(1)}$ ($\kern0.09em j=1,2$) and $\check {\gamma }^{(1)}$ as follows:

(4.5)\begin{gather} \sigma (\tilde{x}_{1, 1}^{(0)} \check{y}_1^{(1)} - \tilde{y}_{1, 1}^{(0)} \check{x}_1^{(1)}) ={-}\check{\psi}_{1, 1}^{(1)} , \end{gather}

(4.6)\begin{align} &\sigma \left\{\check{\phi}_2^{(1)} - \displaystyle{\frac{1}{J_2^{(0)}}} (\tilde{x}_{2, 2}^{(0)} \check{x}_2^{(1)} + \tilde{y}_{2, 2}^{(0)} \check{y}_2^{(1)}) \right\} + \displaystyle{\frac{1}{J_2^{(0)}}} \check{\phi}_{2, 2}^{(1)} - \displaystyle{\frac{1}{(J_2^{(0)})^{2}}} ( \tilde{x}_{2, 2}^{(0)} \check{x}_{2, 2}^{(1)} + \tilde{y}_{2, 2}^{(0)} \check{y}_{2, 2}^{(1)} )\nonumber\\ &\qquad + \displaystyle{\frac{1}{c^{2}}} \check{y}_2^{(1)} + \left\{ \displaystyle{\frac{\frac{1}{2} J_{2, 2}^{(0)}}{(J_2^{(0)})^{2}}} - \displaystyle{\frac{1}{c^{2}}} \tilde{y}_{2, 2}^{(0)} \right\} \check{\gamma}^{(1)}\nonumber\\ &\quad = \displaystyle{\frac{\rho_1}{\rho_2}} \left[ \sigma \left\{ \check{\phi}_1^{(1)} - \displaystyle{\frac{1}{J_1^{(0)}}} (\tilde{x}_{1, 1}^{(0)} \check{x}_1^{(1)} + \tilde{y}_{1, 1}^{(0)} \check{y}_1^{(1)}) \right\} + \displaystyle{\frac{1}{J_1^{(0)}}} \check{\phi}_{1, 1}^{(1)} - \displaystyle{\frac{1}{(J_1^{(0)})^{2}}} ( \tilde{x}_{1, 1}^{(0)} \check{x}_{1, 1}^{(1)} + \tilde{y}_{1, 1}^{(0)} \check{y}_{1,1}^{(1)} ) \right. \nonumber\\ &\qquad \left.+ \displaystyle{\frac{1}{c^{2}}} \check{y}_1^{(1)} - \left\{ \displaystyle{\frac{\frac{1}{2} J_{1, 1}^{(0)}}{(J_1^{(0)})^{2}}} - \displaystyle{\frac{1}{c^{2}}} \tilde{y}_{1, 1}^{(0)} \right\} \check{\gamma}^{(1)} \right] , \end{align}

(4.7)\begin{gather} \check{x}_1^{(1)} - \check{x}_2^{(1)}={-} ( \tilde{x}_{1, 1}^{(0)} + \tilde{x}_{2, 2}^{(0)} ) \check{\gamma}^{(1)} , \end{gather}

(4.8)\begin{gather} \check{y}_1^{(1)} - \check{y}_2^{(1)} ={-} ( \tilde{y}_{1, 1}^{(0)} + \tilde{y}_{2, 2}^{(0)} )\check{\gamma}^{(1)} , \end{gather}

(4.9)\begin{gather} \check{\psi}_1^{(1)}(\xi_1 = \xi_1^{(0)}(\hat{\xi})) = \check{\psi}_2^{(1)}(\xi_2 = \xi_2^{(0)}(\hat{\xi})) , \end{gather}

where $\xi _j = \xi _j^{(0)}(\hat {\xi })$ ($\kern0.09em j = 1, 2$) at the interface are given by (3.2). In addition, eliminating $\gamma ^{(1)}$ from the contact conditions (4.7) and (4.8), we obtain

(4.10)\begin{equation} ( \tilde{y}_{1, 1}^{(0)} + \tilde{y}_{2, 2}^{(0)} ) \check{x}_1^{(1)} - ( \tilde{x}_{1, 1}^{(0)} + \tilde{x}_{2, 2}^{(0)} ) \check{y}_1^{(1)} = ( \tilde{y}_{1, 1}^{(0)} + \tilde{y}_{2, 2}^{(0)} ) \check{x}_2^{(1)} - ( \tilde{x}_{1, 1}^{(0)} + \tilde{x}_{2, 2}^{(0)} ) \check{y}_2^{(1)} . \end{equation}

For later convenience, we use (4.10) instead of (4.8). Five equations (4.5), (4.6), (4.7), (4.9) and (4.10) determine the linear stability of steady interfacial waves.

Following Floquet theory, we can write the general solutions of (4.5), (4.6), (4.7), (4.9) and (4.10) in the form of

(4.11)\begin{equation} \left.\begin{array}{ll} \check{x}_j^{(1)}(\xi_j) = \mathrm{e}^{\mathrm{i} p \xi_j} \displaystyle{\sum_{m={-}\infty}^{\infty}} \alpha_{j m}^{(1)} \mathrm{e}^{\mathrm{i} m \xi_j} , & \check{y}_j^{(1)}(\xi_j) = \mathrm{e}^{\mathrm{i} p \xi_j} \displaystyle{\sum_{m={-}\infty}^{\infty}} a_{j m}^{(1)} \mathrm{e}^{\mathrm{i} m \xi_j} \\ \check{\phi}_j^{(1)}(\xi_j) = \mathrm{e}^{\mathrm{i} p \xi_j} \displaystyle{\sum_{m={-}\infty}^{\infty}} b_{j m}^{(1)} \mathrm{e}^{\mathrm{i} m \xi_j} , & \check{\psi}_j^{(1)}(\xi_j) = \mathrm{e}^{\mathrm{i} p \xi_j} \displaystyle{\sum_{m={-}\infty}^{\infty}} \beta_{j m}^{(1)} \mathrm{e}^{\mathrm{i} m \xi_j} \end{array} \right\} \quad (\kern0.09em j = 1, 2) , \end{equation}

and

(4.12)\begin{equation} \check{\gamma}^{(1)}(\hat{\xi}) = \mathrm{e}^{\mathrm{i} p \hat{\xi}} \displaystyle{\sum_{m={-}\infty}^{\infty}} c_m^{(1)} \mathrm{e}^{\mathrm{i} m \hat{\xi}} , \end{equation}

where the exponent $p$ is assumed to be real so that the solutions are bounded for $-\infty < \xi _j < \infty$ ($\kern0.09em j = 1, 2$). Since the linearized equations (4.5), (4.6), (4.7), (4.9) and (4.10) are invariant under the transformation $\xi _j \to -\xi _j$, $\check {\phi }_j^{(1)} \to -\check {\phi }_j^{(1)}$, $\check {\psi }_j^{(1)} \to -\check {\psi }_j^{(1)}$ ($\kern0.09em j = 1, 2$), $\hat {\xi } \to -\hat {\xi }$, $t \to -t$ and $\gamma ^{(1)} \to -\gamma ^{(1)}$, there is a degeneracy in the choice of $p$, similarly to the case of surface waves (for e.g. see Tiron & Choi (Reference Tiron and Choi2012, § 2.4) and Murashige & Choi (Reference Murashige and Choi2020, § 4.3)). In particular, if $\sigma$ is an eigenvalue for $p$, then $-\bar {\sigma }$ is another eigenvalue for the same value of $p$, and $\bar {\sigma }$ and $-\sigma$ are the eigenvalues for $1-p$. Thus the range of $p$ can be restricted to $0 \leqslant p \leqslant 1/2$. In addition, from the analyticity given by (2.21a–d) and the following relations for the Hilbert transform

(4.13) \begin{equation} \mathcal{H}_j[ \mathrm{e}^{\mathrm{i} \nu \xi_j} ] = \mathrm{i}\, \mathrm{sgn}(\nu)\, \mathrm{e}^{\mathrm{i} \nu \xi_j} \quad (\kern0.09em j = 1, 2)\ \mbox{with} \ \mathrm{sgn}(\nu) = \left\{ \begin{array}{@{}ll} +1 & (\nu > 0) \\ \hphantom{0\,}0 & (\nu = 0) \\ -1 & (\nu < 0) \end{array} \right. , \end{equation}

the coefficients $\alpha _{j m}^{(1)}$ and $\beta _{j m}^{(1)}$ ($\kern0.09em j = 1, 2$) in (4.11) can be related to $a_{j m}^{(1)}$ and $b_{j m}^{(1)}$, respectively, by

(4.14) \begin{equation} \left.\begin{array}{ll@{}} \alpha_{1 m}^{(1)} = \mathrm{i}\, \mathrm{sgn}(p + m)\, a_{1 m}^{(1)} , & \beta_{1 m}^{(1)} ={-}\mathrm{i}\, \mathrm{sgn}(p + m)\, b_{1 m}^{(1)}\\ \alpha_{2 m}^{(1)} ={-}\mathrm{i}\, \mathrm{sgn}(p + m)\, a_{2 m}^{(1)} , & \beta_{2 m}^{(1)} = \mathrm{i}\, \mathrm{sgn}(p + m)\, b_{2 m}^{(1)} \end{array} \right\} . \end{equation}

Substituting these into the linearized equations (4.5), (4.6), (4.7), (4.9) and (4.10), we obtain

(4.15)\begin{gather} \sigma \displaystyle{\sum_{m ={-}\infty}^{\infty}} A_{11, m}(\hat{\xi}) a_{1 m}^{(1)} = \displaystyle{\sum_{m ={-}\infty}^{\infty}} B_{11, m}(\hat{\xi}) b_{1 m}^{(1)} , \end{gather}

(4.16)\begin{align} &\sigma \displaystyle{\sum_{m ={-}\infty}^{\infty}} \{ A_{21, m}(\hat{\xi}) a_{1 m}^{(1)} + A_{22, m}(\hat{\xi}) a_{2 m}^{(1)} + B_{21, m}(\hat{\xi}) b_{1 m}^{(1)} + B_{22, m}(\hat{\xi}) b_{2 m}^{(1)} \} \nonumber\\ &\quad = \displaystyle{\sum_{m ={-}\infty}^{\infty}} \{ A_{31, m}(\hat{\xi}) a_{1 m}^{(1)} + A_{32, m}(\hat{\xi}) a_{2 m}^{(1)} + B_{31, m}(\hat{\xi}) b_{1 m}^{(1)} + B_{32, m}(\hat{\xi}) b_{2 m}^{(1)} + C_{3, m}(\hat{\xi}) c_{m}^{(1)} \}, \end{align}

(4.17)\begin{gather} \displaystyle{\sum_{m ={-}\infty}^{\infty}} A_{41, m}(\hat{\xi}) a_{1 m}^{(1)} +\displaystyle{\sum_{m ={-}\infty}^{\infty}} A_{42, m}(\hat{\xi}) a_{2 m}^{(1)} = \displaystyle{\sum_{m ={-}\infty}^{\infty}} C_{4, m}(\hat{\xi}) c_{m}^{(1)} , \end{gather}

(4.18a,b)\begin{gather} \displaystyle{\sum_{m ={-}\infty}^{\infty}} A_{51, m}(\hat{\xi}) a_{1 m}^{(1)} \!=\! \displaystyle{\sum_{m ={-}\infty}^{\infty}} A_{52, m}(\hat{\xi}) a_{2 m}^{(1)} \quad \mbox{and} \quad \displaystyle{\sum_{m ={-}\infty}^{\infty}} B_{51, m}(\hat{\xi}) b_{1 m}^{(1)} \!=\! \displaystyle{\sum_{m ={-}\infty}^{\infty}} B_{52, m}(\hat{\xi}) b_{2 m}^{(1)} , \end{gather}

where $A_{\ast, m}(\hat {\xi })$, $B_{\ast, m}(\hat {\xi })$ and $C_{\ast, m}(\hat {\xi })$ are $2{\rm \pi}$-periodic functions determined by the steady wave solutions and summarized in Appendix B. Furthermore, from periodicity, the coefficient functions $A_{\ast, m}(\hat {\xi })$, $B_{\ast, m}(\hat {\xi })$ and $C_{\ast, m}(\hat {\xi })$ can be expanded in the form of Fourier series as

(4.19a–c)\begin{align} A_{{\ast}, m}(\hat{\xi}) = \displaystyle{\sum_{k={-}\infty}^{\infty}} A_{{\ast}, km} \, \mathrm{e}^{\mathrm{i} k \hat{\xi}} , \quad B_{{\ast}, m}(\hat{\xi}) = \displaystyle{\sum_{k={-}\infty}^{\infty}} B_{{\ast}, km} \, \mathrm{e}^{\mathrm{i} k \hat{\xi}} , \quad C_{{\ast}, m}(\hat{\xi}) = \displaystyle{\sum_{k={-}\infty}^{\infty}} C_{{\ast}, km} \, \mathrm{e}^{\mathrm{i} k \hat{\xi}} . \end{align}

Then, we can transform (4.15), (4.16), (4.17) and (4.18a,b) into

(4.20)\begin{gather} \sigma \displaystyle{\sum_{m ={-}\infty}^{\infty}} A_{11, km} a_{1 m}^{(1)} = \displaystyle{\sum_{m ={-}\infty}^{\infty}} B_{11, km} b_{1 m}^{(1)} , \end{gather}

(4.21) \begin{align} &\sigma \displaystyle{\sum_{m ={-}\infty}^{\infty}} \{ A_{21, km} a_{1 m}^{(1)} + A_{22, km} a_{2 m}^{(1)} + B_{21, km} b_{1 m}^{(1)} + B_{22, km} b_{2 m}^{(1)} \}\nonumber\\ &\quad= \displaystyle{\sum_{m ={-}\infty}^{\infty}} \{ A_{31, km} a_{1 m}^{(1)} + A_{32, km} a_{2 m}^{(1)} + B_{31, km} b_{1 m}^{(1)} + B_{32, km}b_{2 m}^{(1)} + C_{3, km} c_{m}^{(1)} \} , \end{align}

(4.22)\begin{gather} \displaystyle{\sum_{m ={-}\infty}^{\infty}} A_{41, km} a_{1 m}^{(1)} + \displaystyle{\sum_{m ={-}\infty}^{\infty}} A_{42, km} a_{2 m}^{(1)} = \displaystyle{\sum_{m ={-}\infty}^{\infty}} C_{4, km} c_{m}^{(1)} , \end{gather}

(4.23a,b)\begin{gather} \displaystyle{\sum_{m ={-}\infty}^{\infty}} A_{51, km} a_{1 m}^{(1)} = \displaystyle{\sum_{m ={-}\infty}^{\infty}} A_{52, km} a_{2 m}^{(1)} \quad \mbox{and} \quad \displaystyle{\sum_{m ={-}\infty}^{\infty}} B_{51, km} b_{1 m}^{(1)} = \displaystyle{\sum_{m ={-}\infty}^{\infty}} B_{52, km} b_{2 m}^{(1)} , \end{gather}

for all $k \in \mathbb {Z}$. These equations (4.20), (4.21), (4.22) and (4.23a,b) determine $\sigma$, $a_{j m}^{(1)}$ and $b_{j m}^{(1)}$ ($\kern0.09em j = 1, 2$). In particular, in the limit of wave steepness $h \to 0$, we can analytically obtain the expression of $\sigma$, in terms of the exponent $p$ and the mode number $m$, as

(4.24)\begin{equation} \sigma \to \sigma_{p,m}^{{\pm}} := \mathrm{i} \{ -(p + m) \pm \sqrt{|p + m|} \} \quad \mbox{as} \ h \to 0 . \end{equation}

Here, note that this limit $\sigma _{p,m}^{\pm }$ is independent of the density ratio $\rho _1/\rho _2$.

4.2. The method of computation

First, we can numerically determine the Fourier coefficients $A_{\ast, km}$, $B_{\ast, km}$ and $C_{\ast, km}$ ($k, m = -N, \ldots, N-1$) in (4.19a–c) using the steady wave solutions at $\hat {\xi } = \hat {\xi }_\ell = \ell {\rm \pi}/N$   ($\ell = -N, \ldots, N$) obtained in § 3 and the fast Fourier transform. Here, $N$ is the total number of terms of the truncated Fourier series (A5) for the steady wave solutions. Next, we substitute these Fourier coefficients into (4.19a–c), (4.20), (4.21), (4.22) and (4.23a,b), and truncate the resulting infinite series as

(4.25)\begin{equation} \displaystyle{\sum_{k ={-}\infty}^{\infty}} \sim \displaystyle{\sum_{k ={-}M}^{M}} \quad \mbox{and} \quad \displaystyle{\sum_{m ={-}\infty}^{\infty}} \sim \displaystyle{\sum_{m ={-}M}^{M}} \quad \mbox{with} \ M < N . \end{equation}

To reduce numerical errors, large-wavenumber components are filtered out by setting $M < N$. Then (4.20), (4.21), (4.22) and (4.23a,b) can be rewritten, respectively, as

(4.26)\begin{gather} \sigma \boldsymbol{\mathsf{A}}_{11} \boldsymbol{a}_{1}^{(1)} = \boldsymbol{\mathsf{B}}_{11} \boldsymbol{b}_{1}^{(1)} , \end{gather}

(4.27)\begin{align} &\sigma ( \boldsymbol{\mathsf{A}}_{21} \boldsymbol{a}_{1}^{(1)} + \boldsymbol{\mathsf{A}}_{22} \boldsymbol{a}_{2}^{(1)} + \boldsymbol{\mathsf{B}}_{21} \boldsymbol{b}_{1}^{(1)} + \boldsymbol{\mathsf{B}}_{22} \boldsymbol{b}_{2}^{(1)} ) \nonumber\\ &\quad = \boldsymbol{\mathsf{A}}_{31} \boldsymbol{a}_{1}^{(1)} + \boldsymbol{\mathsf{A}}_{32} \boldsymbol{a}_{2}^{(1)} + \boldsymbol{\mathsf{B}}_{31} \boldsymbol{b}_{1}^{(1)} + \boldsymbol{\mathsf{B}}_{32} \boldsymbol{b}_{2}^{(1)} + \boldsymbol{\mathsf{C}}_{3} \boldsymbol{c}^{(1)} , \end{align}

(4.28)\begin{gather} \boldsymbol{\mathsf{A}}_{41} \boldsymbol{a}_{1}^{(1)} + \boldsymbol{\mathsf{A}}_{42} \boldsymbol{a}_{2}^{(1)} = \boldsymbol{\mathsf{C}}_{4} \boldsymbol{c}^{(1)} , \end{gather}

(4.29a,b)\begin{gather} \boldsymbol{\mathsf{A}}_{51} \boldsymbol{a}_{1}^{(1)} = \boldsymbol{\mathsf{A}}_{52} \boldsymbol{a}_{2}^{(1)} \quad \mbox{and} \quad \boldsymbol{\mathsf{B}}_{51} \boldsymbol{b}_{1}^{(1)} = \boldsymbol{\mathsf{B}}_{52} \boldsymbol{b}_{2}^{(1)} , \end{gather}

where $\boldsymbol{\mathsf{A}}_{\ast } = (A_{\ast, km})_{k,m=-M}^{M}$, $\boldsymbol{\mathsf{B}}_{\ast } = (B_{\ast, km})_{k,m=-M}^{M}$ and $\boldsymbol{\mathsf{C}}_{\ast } = (C_{\ast, km})_{k,m=-M}^{M}$ are $(2M+1) \times (2M+1)$ matrices, and $\boldsymbol {a}_{j}^{(1)} \!=\! (a_{j m}^{(1)})_{m=-M}^{M}$, $\boldsymbol {b}_{j}^{(1)} \!=\! (b_{j m}^{(1)})_{m=-M}^{M}$ and $\boldsymbol {c}^{(1)} = (c_{m}^{(1)})_{m=-M}^{M}$ ($\kern0.09em j = 1, 2$) are vectors with $2M+1$ components. Furthermore, from (4.28) and (4.29a,b), we can represent $\boldsymbol {a}_{2}^{(1)}$, $\boldsymbol {b}_{2}^{(1)}$ and $\boldsymbol {c}^{(1)}$ by $\boldsymbol {a}_{1}^{(1)}$ and $\boldsymbol {b}_{1}^{(1)}$ as

(4.30a–c)\begin{equation} \boldsymbol{a}_{2}^{(1)} = \boldsymbol{\mathsf{A}}_{53} \boldsymbol{a}_{1}^{(1)} , \quad \boldsymbol{b}_{2}^{(1)} = \boldsymbol{\mathsf{B}}_{53} \boldsymbol{b}_{1}^{(1)} , \quad \boldsymbol{c}^{(1)} = \boldsymbol{\mathsf{C}}_5 ~ \boldsymbol{a}_{1}^{(1)} , \end{equation}

where

(4.31a–c)\begin{equation} \boldsymbol{\mathsf{A}}_{53} = {\boldsymbol{\mathsf{A}}_{52}}^{{-}1} \boldsymbol{\mathsf{A}}_{51}, \quad \boldsymbol{\mathsf{B}}_{53} = {\boldsymbol{\mathsf{B}}_{52}}^{{-}1} \boldsymbol{\mathsf{B}}_{51} , \quad \boldsymbol{\mathsf{C}}_5 = {\boldsymbol{\mathsf{C}}_{4}}^{{-}1}\ (\boldsymbol{\mathsf{A}}_{41} + \boldsymbol{\mathsf{A}}_{42} \boldsymbol{\mathsf{A}}_{53} ). \end{equation}

Substituting these into (4.26) and (4.27), we obtain

(4.32)\begin{equation} \sigma\left( \begin{array}{cc} \boldsymbol{\mathsf{A}}_{11} & \boldsymbol{\mathsf{0}} \\ \boldsymbol{\mathsf{A}}_{23} & \boldsymbol{\mathsf{B}}_{23} \end{array}\right)\left(\begin{array}{c}\boldsymbol{a}_{1}^{(1)} \\ \boldsymbol{b}_{1}^{(1)}\end{array}\right)=\left( \begin{array}{cc} \boldsymbol{\mathsf{0}} & \boldsymbol{\mathsf{B}}_{11}\\ \boldsymbol{\mathsf{A}}_{33} & \boldsymbol{\mathsf{B}}_{33} \end{array}\right)\left(\begin{array}{c} \boldsymbol{a}_{1}^{(1)}\\ \boldsymbol{b}_{1}^{(1)}\end{array} \right), \end{equation}

where

(4.33) \begin{equation} \left.\begin{array}{ll@{}} \boldsymbol{\mathsf{A}}_{23} = \boldsymbol{\mathsf{A}}_{21} + \boldsymbol{\mathsf{A}}_{22} \boldsymbol{\mathsf{A}}_{53}, & \boldsymbol{\mathsf{A}}_{33} = \boldsymbol{\mathsf{A}}_{31} + \boldsymbol{\mathsf{A}}_{32} \boldsymbol{\mathsf{A}}_{53} + \boldsymbol{\mathsf{C}}_{3} \boldsymbol{\mathsf{C}}_5\\ \boldsymbol{\mathsf{B}}_{23} = \boldsymbol{\mathsf{B}}_{21} + \boldsymbol{\mathsf{B}}_{22} \boldsymbol{\mathsf{B}}_{53} , & \boldsymbol{\mathsf{B}}_{33} = \boldsymbol{\mathsf{B}}_{31} + \boldsymbol{\mathsf{B}}_{32} \boldsymbol{\mathsf{B}}_{53} \end{array} \right\} . \end{equation}

This is a generalized eigenvalue problem for the eigenvalue $\sigma$ and the eigenvector $(\boldsymbol {a}_{1}^{(1)}, \boldsymbol {b}_{1}^{(1)})$. When the density ratio $\rho _1/\rho _2$, the steady wave steepness $h$ and the exponent $p$ of disturbances are given, we can solve (4.32). In this work, we numerically solve (4.32) using computational routines for eigenvalue problems in LAPACK (http://www.netlib.org/lapack/).

4.3. The dominant wavenumber of disturbances and classification of eigenvalues

The generalized eigenvalue problem (4.32) produces $4M + 2$ sets of the eigenvalue $\sigma$ and the corresponding eigenvector $(\boldsymbol {a}_1^{(1)}, \boldsymbol {b}_1^{(1)})$. Here, as shown in (4.11) and (4.25), $\boldsymbol {a}_1^{(1)} = (a_{1 m}^{(1)})_{m=-M}^{M}$ and $\boldsymbol {b}_1^{(1)} = (b_{1 m}^{(1)})_{m=-M}^{M}$ are the truncated Fourier coefficient vectors for $\check {y}_1^{(1)}$ and $\check {\phi }_1^{(1)}$, respectively. For each eigenvalue $\sigma$, the disturbance has $2M+1$ Fourier components, from which the dominant mode number $\mu$ is defined as the value of $m$ at which the absolute value $|a_{1 m}^{(1)}|$ is the maximum for $-M \leqslant m \leqslant M$, namely

(4.34)\begin{equation} \left| a_{1 \mu}^{(1)} \right| = \displaystyle{\max_{{-}M \leqslant m \leqslant M}} \left| a_{1 m}^{(1)} \right|. \end{equation}

An alternative is to use $|b_{1 m}^{(1)}|$, but no difference is expected for the choice of $\mu$. Note that $\mu$ is an integer and $-M \leqslant \mu \leqslant M$. Using this definition, $4M + 2$ eigenvalues are classified into the $2M + 1$ dominant mode groups. Then, as $\check {y}_1^{(1)}(\xi _1)$ can be expressed as

(4.35)\begin{equation} \check {y}_1^{(1)}(\xi _1) = \sum _{m=-M}^{M} a_{1 m}^{(1)} \exp ({\mathrm {i} (m + p) \xi _1}), \end{equation}

the dominant wavenumber of the disturbance $\check {y}_1^{(1)}$ is $\mu + p$ and, if $| a_{1 \mu }^{(1)} |$ is much greater than $| a_{1 m}^{(1)} |$ for all $m$ ($m \neq \mu$), we may approximate $\check {y}_1^{(1)}(\xi _1)$ as $\check {y}_1^{(1)}(\xi _1) \sim a_{1 \mu }^{(1)} \exp ({\mathrm {i} (\mu + p) \xi _1})$.

Figure 5 shows the variation of $| a_{1 m}^{(1)} |$ with $m$ and the corresponding disturbance $\check {y}_1^{(1)}$ for two eigenvalues. The physical and numerical parameters for these solutions are $\rho _1/\rho _2 = 0.9$, $h = 0.5$, $p = 1/2$, $N = 128$ and $M=60$. One eigenvalue is for a stable mode ($\sigma _{{r}} \simeq 0.0$, $\sigma _{{i}} \simeq -0.257$) and the other is for an unstable mode ($\sigma _{{r}} \simeq 0.706$, $\sigma _{{i}} \simeq -24.7$). The results suggest that the dominant mode is well defined when its value is small, as shown in figure 5(a), but the bandwidth of $|a_{1 m}^{(1)}|$ increases with $|\mu |$. For large $|\mu |$, there is some ambiguity in the definition of $\mu$, as can be observed in figure 5(b). Furthermore, when $|\mu |$ is found large, the convergence of the Fourier series expansion for the corresponding disturbances in (4.11) and (4.12) becomes slow due to its truncation. Then $M$ needs to be further increased, but, unfortunately, the computational cost increases with $M$. Therefore, for a fixed value of $M$, this causes a limitation of the classification of eigenvalues in terms of $\mu$. Due to our limited computational resources, we choose $N = 128$ and $M = 60$, for which we focus on the modes of $|\mu | \leqslant 40$ for $h \leqslant 1$.

Figure 5.

Examples of the variation of $| a_{1 m}^{(1)} |$ with the mode number $m$ and the corresponding disturbance $\check {y}_1^{(1)}$ for $\rho _1/\rho _2 = 0.9$, $h = 0.5$ and $p = 1/2$: (a) stable case ($\sigma _{{r}} \simeq 0.0$, $\sigma _{{i}} \simeq -0.257$) for which $\mu = 1$; (b) unstable case ($\sigma _{{r}} \simeq 0.706$, $\sigma _{{i}} \simeq -24.7$) for which $\mu = 20$. (a1) and (b1) show $|a_{1 m}^{(1)}|$ for $\mu$ = 1 and 20, respectively. (a2) and (b2) show $\check{y}^{(1)}$ for $\mu = 1$ and 20, respectively. The computed results are obtained using the present method with $N = 128$ and $M = 60$. The dominant mode number $\mu$ and the disturbance $\check {y}_1^{(1)}$ are defined by (4.34) and (4.11), respectively. The red and blue lines in (a2) and (b2) represent the real and imaginary parts of $\check {y}_1^{(1)}$, respectively.

5.1. Wave-induced KH instability

As described in § 1, after assuming that the velocity jump near the crest is locally constant, Pullin & Grimshaw (Reference Pullin and Grimshaw1985, (4.1) on p. 330) proposed the approximate growth rate of the wave-induced KH instability for interfacial waves using the well-known result of the KH instability for two horizontal uniform currents of different speed. Based on their approximation, the positive real part, or the growth rate $\sigma _{{r}}^{(\mathrm {KH})}$ of the wave-induced KH instability, can be written as

(5.1)\begin{equation} \sigma_{{r}}^{(\mathrm{{KH}})} \sim \tilde{\sigma}_{{r}}^{(\mathrm{{KH}})} = \mathrm{Re}\left\{ \sqrt{ (\mu + p)^{2} \displaystyle{\frac{\rho_1/\rho_2}{(1 + \rho_1/\rho_2)^{2}}} (\Delta \tilde{q}_{{crest}})^{2} - |\mu + p| \left(\displaystyle{\frac{c_0}{\tilde{c}}}\right)^{2} } \right\} , \end{equation}

where $\mu + p$ is the dominant wavenumber of disturbances defined by (4.34), $\Delta \tilde {q}_{{crest}}$ is the tangential velocity jump at the crest of the steady interfacial wave, $\tilde {c}$ is the steady wave speed and $c_0$ is the linear wave speed given by (3.4). As $\Delta \tilde {q}_{{crest}}$ and $\tilde {c}$ depend on the steady wave steepness $h$, the approximate growth rate $\tilde {\sigma }_{{r}}^{(\mathrm {{KH}})}$ changes with $\rho _1/\rho _2$, $h$, $\mu$ and $p$.

In addition, (5.1) yields the critical mode number $\tilde {\mu }_{{c}}^{(\mathrm {{KH}})}$, at which the inside of the square root in (5.1) vanishes, as

(5.2)\begin{equation} \tilde{\mu}_{{c}}^{(\mathrm{{KH}})} := \displaystyle{\frac{(1 + \rho_1/\rho_2)^{2}}{\rho_1/\rho_2}} \displaystyle{\frac{1}{(\Delta \tilde{q}_{{crest}})^{2}}} \left(\displaystyle{\frac{c_0}{\tilde{c}}}\right)^{2} - p . \end{equation}

For $\mu < \tilde {\mu }_{{c}}^{(\mathrm {{KH}})}$ ($> \tilde {\mu }_{{c}}^{(\mathrm {{KH}})}$), the steady interfacial wave is stable (unstable). Figure 6 shows the variation of $\tilde {\mu }_{{c}}^{(\mathrm {{KH}})}$ with the steady wave steepness $h$ for the exponent $p = 1/2$ and different values of the density ratio $\rho _1/\rho _2$.

Figure 6.

Variation of the approximate critical mode number $\tilde {\mu }_{{c}}^{(\mathrm {{KH}})}$ given by (5.2) with the steady wave steepness $h$ for $p = 1/2$ and different values of the density ratio $\rho _1/\rho _2$. In computing $\tilde {\mu }_{{c}}^{(\mathrm {{KH}})}$, $\Delta \tilde {q}_{{crest}}$ and $\tilde {c}$ in (5.2) are estimated using both the fifth-order approximate solutions of Tsuji & Nagata (Reference Tsuji and Nagata1973) (solid) and the numerical solutions with $N=128$ (dashed).

Note that $\sigma _{{r}}^{(\mathrm {{KH}})}$ and $\tilde {\mu }_{{c}}^{(\mathrm {{KH}})}$ can be estimated using either the fifth-order approximate solutions (Tsuji & Nagata Reference Tsuji and Nagata1973), or the fully nonlinear computed results for $\Delta \tilde {q}_{{crest}}$ and $\tilde {c}$. As can be seen in figure 6, the two results show little difference, except for $\rho _1/\rho _2 = 0.1$. Therefore, in § 6, the fifth-order approximate solutions are used for $\Delta \tilde {q}_{{crest}}$ and $\tilde {c}$ as $\sigma _{{r}}^{(\mathrm {{KH}})}$ and $\tilde {\mu }_{{c}}^{(\mathrm {{KH}})}$ can be analytically estimated.

5.2. Modulational instability

Introducing the slow space scale $X = \epsilon (x^{\prime } + c_g t)$ and the slow time scale $\tau = \epsilon ^{2} t$ in the inertial frame $(x^{\prime }, y^{\prime }, t)$ with $x^{\prime } = x - c t$ and $y^{\prime } = y$, Grimshaw & Pullin (Reference Grimshaw and Pullin1985, (4.4), (4.8) and (4.9) on pp. 304–305) derived the nonlinear Schrödinger (NLS) equation for slowly modulated small-amplitude interfacial waves as

(5.3)\begin{equation} -\mathrm{i} A_\tau + \alpha A_{X X} + \beta |A|^{2} A = 0 , \end{equation}

with

(5.4a–c)\begin{equation} c_g = \displaystyle{\frac{1}{2}} c_0 , \quad \alpha ={-}\displaystyle{\frac{1}{8}} , \quad \beta ={-}\displaystyle{\frac{1 + (\rho_1/\rho_2)^{2}}{2(1 + \rho_1/\rho_2)^{2}}} , \end{equation}

where $c_g$ denotes the group velocity, $c_0$ is given by (3.4) and $A = A(X,\tau )$ is the complex wave amplitude of the leading-order wave train. Following Mei, Stiassnie & Yue (Reference Mei, Stiassnie and Yue2005, chap. 13), we can investigate the linear stability of the basic Stokes wave given by

(5.5)\begin{equation} A_0 = a_0 \mathrm{e}^{-\mathrm{i}\beta {a_0}^{2} \tau}\quad (a_0 \in \mathbb{R}) , \end{equation}

which is the solution of $-\mathrm {i} A_{0 \tau } + \beta |A_0|^{2} A_0 = 0$. In particular, the approximate growth rate $\tilde {\sigma }_{{r}}^{(\mathrm {{NLS}})}$ of disturbances for the modulational instability is given by

(5.6)\begin{equation} \tilde{\sigma}_{{r}}^{(\mathrm{{NLS}})} := \mathrm{Re}\left\{ \displaystyle{\frac{p}{8}} \sqrt{ 2 \displaystyle{\frac{1 + (\rho_1/\rho_2)^{2}}{(1 + \rho_1/\rho_2)^{2}}} h^{2} - p^{2} } \right\} . \end{equation}

This is supposed to be valid for small values of $h$ and $p$.

6.1. Wave-induced KH instability

Figure 7 shows the variation of the growth rate $\sigma _{{r}}$ with the dominant mode number $\mu$ for the density ratio $\rho _1/\rho _2 = 0.9$, the exponent $p = 1/2$ and different values of the wave steepness: (a) $h = 0.4$, (b) 0.5, (c) 0.6 and (d) 0.7. The computed results are compared with the approximate growth rate $\tilde {\sigma }_{{r}}^{(\mathrm {{{KH}}})}$ in (5.1) for the wave-induced KH instability. The computed results show that there exist unstable modes ($\sigma _{{r}} > 0$) for $|\mu | > \mu _{{c}}$, where $\mu _{{c}}$ is the critical mode number. This is analogous to the case of the well-known KH instability for two horizontal uniform currents of different speed. Since there are no background currents in this problem, the instability in figure 7 is the KH instability excited locally by the wave-induced tangential velocity jump. While its variation with $\mu$ is qualitatively similar to the computed results, the approximate growth rate $\tilde {\sigma }_{{r}}^{(\mathrm {{{KH}}})}$ in (5.1) is always overpredicted. Nevertheless, the critical mode number $\mu _{{c}}$ for the onset of the wave-induced KH instability is well approximated by $\tilde {\mu }_{{c}}^{(\mathrm {{KH}})}$ in (5.2). As the wave steepness increases, the critical mode number decreases, as can be seen in figure 7. For any wave steepness, as the growth rate increases with $\mu \ ({>}{\mu }_{{c}})$ and appears to be unbounded as $\mu \to \infty$, the inviscid initial value problem would be ill posed.

Figure 7.

Variation of the growth rate $\sigma _{{r}}$ with the dominant mode number $\mu$ of the corresponding disturbance for $\rho _1/\rho _2 = 0.9$, $p = 1/2$ and different values of the wave steepness $h$. The computed results (dots) using the present method with $N = 128$ and $M = 60$ are compared with the approximate growth rate (red solid lines) $\tilde {\sigma }_{{r}}^{(\mathrm {{{KH}}})}$ given by (5.1). The red dashed line represents the approximate critical mode number $\tilde {\mu }_{{c}}^{(\mathrm {{KH}})}$ given by (5.2) for the wave-induced KH instability.

Figure 8 shows the variation of the growth rate $\sigma _{{r}}$ $({\geqslant }0)$ with the wave steepness $h$ and the dominant mode number $\mu$ for the exponent $p = 1/2$ and different values of the density ratio: (a) $\rho _1/\rho _2 = 0.1$, (b) 0.3, (c) 0.5, (d) 0.7, (e) 0.9 and (f) 0.99. These three-dimensional plots allow us to observe the continuous change of $\sigma _{{r}}$ with both $h$ and $\mu$ simultaneously. The critical dominant wavenumber $\mu _{{c}}$ depends on the density ratio and the wave steepness, and the interfacial periodic steady waves are always unstable for $\mu > \mu _{{c}}$. We can observe that the value of the growth rate $\sigma _{{r}} > 0$ for the unstable range $\mu > \mu _{{c}}$ increases with the density ratio $\rho _1/\rho _2$. However, the computed results vary little for $\rho _1/\rho _2 \geqslant 0.9$. It is found that the estimated critical curve $\mu = \tilde {\mu }_{{c}}^{(\mathrm {{KH}})}(h)$ approximately agrees with the computed results for $\rho _1/\rho _2 \geqslant 0.3$ but not for $\rho _1/\rho _2 = 0.1$. This disagreement for $\rho _1/\rho _2 = 0.1$ may be due to the fact that the tangential velocity jump near the crest for $\rho _1/\rho _2 = 0.1$ changes more rapidly than that for $\rho _1/\rho _2 \geqslant 0.3$, as shown in figure 3, and the tangential velocity jump cannot be assumed to be constant near the crest of the steady waves. Therefore, the approximate local analysis breaks down for small values of $\rho _1/\rho _2$.

Figure 8.

Variation of the growth rate $\sigma _{{r}}$ with the steady wave steepness $h$ and the dominant mode number $\mu$ of disturbances for different values of the density ratio $\rho _1/\rho _2$ ($p = 1/2$). The computed results (blue dots) are obtained using the present method with $N = 128$, $M = 60$. The red solid line represents on the $(\mu, h)$-plane the approximate critical mode number $\mu = \tilde {\mu }_{{c}}^{(\mathrm {{KH}})}(h)$ defined by (5.2).

Figure 9 shows the variation of the real and imaginary parts of the eigenvalue $\sigma = \sigma _{{r}} + \mathrm {i} \sigma _{{i}}$ with the wave steepness $h$ for the density ratio $\rho _1/\rho _2 = 0.9$, the exponent $p = 1/2$ and the dominant mode number $|\mu | \leqslant 40$. The computed results of the imaginary part $\sigma _{{i}}$ demonstrate that some pairs of the eigenvalues collide successively at the wave steepness $h$ shown by the red dashed lines. When two eigenvalues collide, the corresponding growth rate $\sigma _{{r}}$ changes from $\sigma _{{r}} = 0$ to $\sigma _{{r}} \neq 0$, and the dominant mode number $\mu$ is equal to the critical mode number $\mu _{{c}}$. These successive collisions of the eigenvalues result in the wave-induced KH instability found in figures 7 and 8. This observation is consistent with the previous theoretical results that the instability arises when two eigenvalues of opposite Krein signature or opposite energy sign collide, for example, in MacKay & Saffman (Reference MacKay and Saffman1986) for surface waves and Benjamin & Bridges (Reference Benjamin and Bridges1997) for interfacial waves using Hamiltonian approaches.

Figure 9.

Variation of the real and imaginary parts of the eigenvalue $\sigma = \sigma _{{r}} + \mathrm {i} \sigma _{{i}}$ with the steady wave steepness $h$ for $\rho _1/\rho _2 = 0.9$ and $p = 1/2$: (a) $0.41 \leqslant h \leqslant 0.44$; (b) $0.62 \leqslant h \leqslant 0.7$. The results (dots) are computed with $N = 128$ and $M = 60$. At the wave steepness $h$ marked by the red dashed line, two eigenvalues collide and the growth rate $\sigma _{{r}}$ changes from $\sigma _{{r}} = 0$ to $\sigma _{{r}} \neq 0$.

In figures 7, 8 and 9, the value of the exponent $p$ is fixed to $p = 0.5$. Figure 10 shows the variation of the growth rate $\sigma _{{r}}$ with the wave steepness $h$ and the dominant mode number $\mu$ for different values of the exponent $p$ (the density ratio $\rho _1/\rho _2 = 0.9$). For relatively large values of $|\mu |$, the computed eigenvalue $\sigma$ almost remains unchanged with $p$, as shown in figure 10, because the dominant wavenumber $\mu + p$ with $0 \leqslant p \leqslant 1/2$ is nearly equal to $\mu$. On the other hand, for small values of $|\mu |$, namely for small-wavenumber (long-wavelength) disturbances, the eigenvalue $\sigma$ changes with $p$ and a different type of instability is observed, to be discussed in § 6.2.

Figure 10.

Variation of the growth rate $\sigma _{{r}}$ with the wave steepness $h$ and the dominant mode number $\mu$ for the density ratio $\rho _1/\rho _2 = 0.9$ and different values of the exponent $p$: (a) variation of $\sigma _{{r}}$ with $h$ and $\mu$; (b) variation of $\sigma _{{r}}$ with $\mu$. (a1) and (a2) show variation of $\sigma_r$ with h and $\mu$ for $p=0$ and 0.25, respectively. (b1) and (b2) show variation of $\sigma_r$ with $\mu$ for $h=0.5$ and 0.7, respectively. The results (dots) are computed with $N = 128$ and $M = 60$. The red solid line in (a) represents on the $(\mu, h)$-plane the approximate critical mode number $\mu =\tilde {\mu }_{{c}}^{(\mathrm {{KH}})}(h)$ defined by (5.2).

6.2. Small-wavenumber instability

Until now, we have discussed the instability associated with large values of the dominant mode number $\mu$, but the instability for which $\mu$ is small is also observed. Figure 11 shows the variation of the growth rate $\sigma _{{r}}$ with the exponent $p$ for $\mu = -1$, $1$ and $0$. Here, we fix the density ratio to $\rho _1/\rho _2 = 0.9$, but vary the wave steepness: (a) $h = 0.3$, (b) 0.5, (c) 0.7 and (d) 0.9. The computed results using the present method are compared with (i) the approximate growth rate $\tilde {\sigma }_{{r}}^{(\mathrm {{NLS}})}$ in (5.6) for the modulational instability and (ii) the computed results using the previous method (Yuen Reference Yuen1984) that was developed in the physical plane. It is found that the approximate growth rate $\tilde {\sigma }_{{r}}^{(\mathrm {{NLS}})}$ is valid only for small values of $h$ and $p$, as expected. It is also noticed that the computed results using the previous method agree with those using the present method only for $h \leqslant 0.7$, but not for $h = 0.9$. For $h > 0.7$, we could not obtain accurate numerical solutions using the previous method, as shown in figure 11(d).

Figure 11.

Variation of the growth rate $\sigma _{{r}}$ with the exponent $p$ for the density ratio $\rho _1/\rho _2 = 0.9$ and the dominant mode numbers $\mu = -1$, $0$ and $1$. The computed results using the present method with $N = 128$ and $M = 60$ (dots) are compared with those using the previous method by (Yuen Reference Yuen1984) (circles) and the approximate growth rate $\tilde{\sigma}_r^{({\rm NLS})}$ in (5.6) for the modulational instability (red lines).

It should be remarked that the value of the growth rate $\sigma _{{r}} > 0$ of small-wavenumber disturbances in figure 11 is much smaller than that of the wave-induced KH instability in § 6.1. Thus, it is difficult to identify the small-wavenumber instabilities with $\mu = -1$, $0$ and $1$, including the modulational instability in figures 7, 8, 9 and 10. It is interesting to note that the present method captures another type of instability for $h = 0.9$ and $0.21 < p < 0.25$, as shown in figure 11(d). This instability seems to be similar to the ‘nonlinear mode jumping’ phenomenon described by Yuen (Reference Yuen1984, figure 4(g) on p. 80) for $\rho _1/\rho _2 = 0.9$ and $h = 0.5$ with the background current jump $\Delta U = 0.5 (c_0/c)$. The source of this instability remains unknown and is a topic for future research.