Center for Geophysics and Planetary Physics, Department of Mathematics, Embry-
Riddle Aeronautical University, Daytona Beach, FL, U.S.A.
emails<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Interaction of Tollmien-Schlichting Waves in the Air with the Sea Surface

Shahrdad G. Sajjadi and Harihar Khanal

###### Abstract

Tollmien-Schlichting waves, Hydrodynamic stability, Air-sea interaction, Two-
fluid interface, Solitons Linear stability of fully developed flows of air
over water is carried out in order to study non-linear effects in the
generation of solitons by wind. A linear stability analysis of the basic flow
is made and the conditions at which solitons first begin to grow is
determined. Then, following [10], the non-linear stability of the flow is
examined and the quintic non-linear Schrödinger equation is derived for the
amplitude of disturbances. The coefficients of the non-linear Schrödinger
equation are calculated from the eigenvalue problem which determines the
stability of air-water interface.

An asymptotic and a numerical stability analysis is carried out to determine
the neutrally stable flow conditions for air-sea interface. Four different
profiles are considered for the airflow blowing over the surface of the sea,
namely, plane Couette flow (pCf), plane Poiseuille flow (pPf), laminar and
turbulent boundary layer (L,TBL) profiles. For each of the above cases the
shear flow counterpart in the water is assumed to be a pPf.

A nonlinear stability analysis results in the nonlinear Schrödinger equation

$\displaystyle\displaystyle\frac{\partial
A}{\partial\tau}-a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}=\displaystyle\frac{d_{1}}{d_{1r}}A+\kappa
A|A|^{2}+\varpi A|A|^{4},$

where $\xi$ and $\tau$ are local variables, and the amplitude of the surface
wave is proportional to $A$. The complex constants $a_{2}$, $d_{1}$, $\kappa$
and $\varpi$ are evaluated from the linear stability of the two-fluid
interface. The profile for the initial condition considered here is that of
the Stokes wave

$\displaystyle\eta=a\cos kx+\textstyle\frac{1}{2}a^{2}k\cos 2kx,$

where $a$ is the amplitude and $k$ is the wavenumber of the surface Stokes
wave.

It is shown that the above amplitude equation produces ‘snake’ solitons [9]
for pCf, pPf and LBL profiles, with striking similarities. On the other hand,
for TBL we observe a very violent surface motion. For cases of pCf and LBL
remarkable similarity is observed with observations made at sea.

We conclude that the effect of nonlinearity in the airflow over the sea
surface is much larger than nonlinear interactions in the water, and hence it
is not possible to decouple the motion in the air and the water for finite
amplitude wind-wave interactions, particularly in the case of wind-generated
solitons in shallow waters.

##  1 Introduction

The wind-generated solitons which moves nondispersively on the thin wind-
driven drift layer of wind-excited surface wave waves are important because
they are products of the processes under which surface waves are generated by
wind. When the local atmospheric conditions are such that significant bubble
densities are created, principally by breaking waves, the resulting bubble
layer will dominate and mask these solitons, in underwater acoustic
scattering. A qualitative mixed linear-nonlinear model of ocean-wave
generation is briefly described, a key element of which is the generation of
solitons by wind. These solitons are created, and destroyed on moving gravity
waves, caused by waves being generated by wind along with accompanying
capillary components. The gravity-capillary waves are dispersive, in contrast
the solitons are nondispersive, that is to say, they move at constant speed,
on the thin moving wind drift layer generated by the upper part of the water
surface.

In this paper, we will examine the stability of laminar and turbulent flows of
air over water surface where viscous effects are taken into account in both
fluids. The flow that will be studied here is the fully developed flow of air
over water. In contrast to [2], we shall restrict our analysis to only the
case when the motion is driven solely by a constant pressure gradient.

A linear stability analysis will be performed to determine the neutrally
stable flow conditions. This approach is different to [4], or indeed other
similar theories such as [1], and most experimental work, where the emphasis
is on determining wave growth rates rather than the conditions at which waves
will first start to grow. Following [10], a non-linear analysis will be
performed which will yield the non-linear quintic Schrödinger equation

$\displaystyle\displaystyle\frac{\partial
A}{\partial\tau}-a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}=\displaystyle\frac{d_{1}}{d_{1r}}A+\kappa
A|A|^{2}+\varpi A|A|^{4},$ (1)

governing the wave amplitude. In equation (1) $A$ is proportional to amplitude
of the non-linear surface wave. The scaled length and time, $\xi$ and $\tau$,
and the complex constants $a_{2},d_{1}$ and $\kappa$ and $\varpi$ will be
defined in the proceeding analysis. The constants $a_{2}$ and $d_{1}$ are
properties of the flow which is obtained from the linear stability theory,
while the effect of non-linear interactions is determined by $\kappa$ and
$\varpi$.

It is to be noted that although equation (1) is strictly applicable to
velocity profiles of laminar flow, there are experimental and theoretical
suggestions that equations similar to that of (1) are also applicable to wind-
wave generation situations. For example, the experimental work [6] suggests
that wind-waves at fixed fetch and under steady wind conditions have
properties that are similar to non-linear Stokes wavetrains. Indeed, the wind-
waves given by the equilibrium amplitude solutions of (1) (for which
$\kappa_{r}<0$ and $\varpi=0$) are similar to Stokes wavetrains

$\displaystyle\eta_{w}=a\cos\left\\{k_{c}x-k_{c}c_{r}t\left[\textstyle\frac{1}{2}\alpha(ak_{c})^{2}+\beta(ak_{c})^{2}\right]\right\\}$
(2)

where $k_{c}$ is the critical wavenumber, $\alpha$ and $\beta$ are the
dimensionless energy transfer parameters [13], and $c$ is the wave phase
speed. It has been shown [11] that such wavetrains with non-linear effects
lead to changes in the phase speed of the wave and provides the sideband
instability mechanism, similar to that of Benjamin-Fier, which will eventually
result to the break-up of the wave-train.

The wind-generated solitons which moves nondispersively on the thin wind-
driven drift layer of wind-excited surface wave waves are important because
they are products of the processes under which surface waves are generated by
wind. When the local atmospheric conditions are such that significant bubble
densities are created, principally by breaking waves, the resulting bubble
layer will dominate and mask these solitons, in underwater acoustic
scattering. A qualitative mixed linear-nonlinear model of ocean-wave
generation is briefly described, a key element of which is the generation of
solitons by wind. These solitons are created, and destroyed on moving gravity
waves, caused by waves being generated by wind along with accompanying
capillary components. The gravity-capillary waves are dispersive, in contrast
the solitons are nondispersive, that is to say, they move at constant speed,
on the thin moving wind drift layer generated by the upper part of the water
surface.

In this paper, we will examine the stability of laminar and turbulent flows of
air over water surface where viscous effects are taken into account in both
fluids. The flow that will be studied here is the fully developed flow of air
over water. In contrast to [2], we shall restrict our analysis to only the
case when the motion is driven solely by a constant pressure gradient.

##  2 Formulation of the problem

We consider the steady, flow of two super-imposed fluid layers confined
between a stationary horizontal parallel plate, placed at $y=-D_{\ell}$, and
open to atmosphere above at a distance $y=D_{u}\gg 1$. We assume the flow of
the upper fluid is turbulent, satisfying Prandtl’s mixing-length model, and
that of the lower one is a fully developed laminar flow. The velocity vector
$\mbox{\boldmath$U$}=[U_{b}(y),0,0]$, for the upper fluid, and
$\mbox{\boldmath$u$}=[u_{b}(y),0,0]$, for the lower fluid,111Hereafter, the
upper and the lower symbols (Greek or Roman) denote the upper (air) and the
lower (water) fluids, respectively. Also the subscripts $a$ and $w$ denote air
and water, respectively. both satisfy the Navier-Stokes equations. In the
present case we shall confine ourselves to the case when the motion of the
lower fluid is driven solely by a constant pressure gradient. Thus, the
equations of motion and their corresponding boundary conditions reduce to

$\displaystyle\left.\begin{array}[]{lll}U_{b}\rightarrow
U_{\infty},&&\mbox{as}\hskip 7.11317pty\rightarrow D_{u}\\\
l^{2}\left(\displaystyle\frac{dU_{b}}{dy}\right)^{2}=U_{*}^{2};\,\,l=\kappa(y+y_{0}),&&\mbox{in
upper fluid}\\\ u_{b}=U_{b},\hskip
7.11317pt\mu_{w}\displaystyle\frac{du_{b}}{dy}=\mu_{a}\displaystyle\frac{dU_{b}}{dy},&&\mbox{on}\hskip
7.11317pty=0\\\ \displaystyle\frac{G}{\rho_{w}}+\nu_{w}\displaystyle\frac{{\rm
d}^{2}u_{b}}{{\rm d}y^{2}}=0,&&\mbox{in lower fluid}\\\
u_{b}=0,&&\mbox{on}\hskip 7.11317pty=-D_{\ell}\\\ \end{array}\right\\}$ (8)

where $-G=\displaystyle\frac{dp}{dx}.$ The solutions of (8) are given by

$\displaystyle
U_{b}(y)=U_{1}\ln\left(\displaystyle\frac{y+y_{0}}{y_{0}}\right)+U_{0}$ (9)
$\displaystyle
u_{b}(y)=\displaystyle\frac{G}{2\mu_{w}}(D^{2}_{\ell}-y^{2})+\displaystyle\frac{\mu_{a}}{\mu_{w}}\displaystyle\frac{U_{1}}{y_{0}}(D_{\ell}+y)$
(10)

where

$\displaystyle
U_{0}=\displaystyle\frac{D_{\ell}}{\mu_{w}}\left(\displaystyle\frac{U_{1}\mu_{a}}{y_{0}}+\displaystyle\frac{GD_{\ell}}{2}\right),\hskip
14.22636ptU_{1}=\displaystyle\frac{U_{*}}{\kappa}$ (11)

$U_{*}$ is the friction velocity of the upper fluid, $\kappa$ is the von
Kármán’s constant, and $y_{0}$ is an effective roughness length.

Having obtained the velocity fields governing the basic flow, we shall now
formulate the non-dimensional equations that govern linear instability. The
motion is governed by the three-dimensional Navier-Stokes equations in each
layer of fluid. The boundary conditions are the no-slip conditions on the
plates and continuity of stress and velocity across the unknown interface
$y=\zeta(x,z,t)$, which will be determined from the proceeding non-linear
analysis.

We will introduce the following curvilinear coordinates [5]

$\displaystyle\left.\begin{array}[]{l}\chi=x-\displaystyle\frac{a\cosh k(y\mp
d)}{\sinh kd}\sin k(x-ct)\\\ \\\ \eta=y-\displaystyle\frac{a\sinh k(y\mp
d)}{\sinh kd}\cos k(x-ct)\\\ \\\ \vartheta=z\end{array}\right\\}$ (17)

where in (17), the minus sign is for the upper media and the plus sign for the
lower media. The Jacobian of the transformation is

$\displaystyle J=1-\displaystyle\frac{2ak}{\sinh kd}\cosh k(\eta\mp d)\cosh
k(\chi-ct)+O(ak)^{2}$ (18)

Note that, in the linear stability theory we shall be restricted to
infinitesimal disturbances for which the Jacobian (18) reduce to

$\displaystyle J=1+O(ak)$

Thus, within this order of approximation, the governing equations are
symbolically the same as that of Cartesian coordinates. We remark that under
the transformation (17), the lines $y=d$ and $y=-d$ will map to lines $\eta=1$
and $\eta=-d^{-1}$, respectively, and the unknown interface $y=\zeta$ will
then be at $\eta=0$.

The boundary conditions are:
on $\eta=1$ no-slip conditions $\hskip 34.14322ptU=V=W=0$,
on $\eta=0$ (the unknown interface);
(i) the continuty of velocities $\hskip
25.6073pt\mbox{\boldmath$u$}=\mbox{\boldmath$U$}$
(ii) the kinematic conditions

$\displaystyle\left.\begin{array}[]{l}\displaystyle\frac{\partial\zeta}{\partial
t}+U\displaystyle\frac{\partial\zeta}{\partial\chi}+V\displaystyle\frac{\partial\zeta}{\partial\eta}=W\\\
\\\ \displaystyle\frac{\partial\zeta}{\partial
t}+u\displaystyle\frac{\partial\zeta}{\partial\chi}+v\displaystyle\frac{\partial\zeta}{\partial\eta}=w\\\
\end{array}\right\\}$ (22)

(iii) the continuity of tangential stresses

$\displaystyle\left.\begin{array}[]{l}\displaystyle\frac{\partial
U}{\partial\eta}+\displaystyle\frac{\partial
V}{\partial\chi}=\mu^{-1}\left(\displaystyle\frac{\partial
u}{\partial\eta}+\displaystyle\frac{\partial v}{\partial\chi}\right)\\\ \\\
\displaystyle\frac{\partial W}{\partial\eta}+\displaystyle\frac{\partial
V}{\partial\vartheta}=\mu^{-1}\left(\displaystyle\frac{\partial
w}{\partial\eta}+\displaystyle\frac{\partial v}{\partial\vartheta}\right)\\\
\\\ \displaystyle\frac{\partial
U}{\partial\vartheta}+\displaystyle\frac{\partial
W}{\partial\chi}=\mu^{-1}\left(\displaystyle\frac{\partial
u}{\partial\vartheta}+\displaystyle\frac{\partial w}{\partial\chi}\right)\\\
\end{array}\right\\}$ (28)

(iv) the dynamic condition

$\displaystyle-\rho P+\displaystyle\frac{2\rho}{R}\displaystyle\frac{\partial
V}{\partial\eta}+\rho{\mathscr{F}}\zeta=-\rho
p+\displaystyle\frac{2}{R\nu}\displaystyle\frac{\partial
v}{\partial\eta}-{\mathscr{T}}\nabla_{\parallel}^{2}\zeta+{\mathscr{F}}\zeta$
(29)

where
$\nabla_{\parallel}^{2}\equiv\partial^{2}/\partial\chi^{2}+\partial^{2}/\partial\eta^{2}$,
and
$d=D_{u}/D_{\ell},\mu=\mu_{a}/\mu_{w},\rho=\rho_{a}/\rho_{w},\nu=\nu_{a}/\nu_{w}.$
Also, on $\eta=-d^{-1}$ there is no-slip conditions $u=v=w=0$.

In (29) $R=U_{\infty}d/\nu_{a}$ is the Reynolds number,
${\mathscr{F}}=gd/U^{2}_{\infty}$ is the Froude number and
${\mathscr{T}}=S/\rho_{w}U_{\infty}^{2}d$ is the Weber number, where
$U_{\infty}$ is a reference velocity, $S$ is the dimensional surface tension
and $g$ is the acceleration due to gravity.

We next perturb the basic flow by infinitesimal disturbances, i.e.

$\displaystyle\mbox{\boldmath$U$}=[U_{b}(\eta),0,0]+\delta[U(\eta),V(\eta),W(\eta)]e^{i(k\chi+\ell\vartheta-
kct)}+{\rm c.c.}$
$\displaystyle\mbox{\boldmath$u$}=[u_{b}(\eta),0,0]+\delta[u(\eta),v(\eta),w(\eta)]e^{i(k\chi+\ell\vartheta-
kct)}+{\rm c.c.}$ $\displaystyle
P=\displaystyle\frac{\chi{U}^{\prime\prime}_{b}}{\rho\mu R}+\delta
P(\eta)e^{i(k\chi+\ell\vartheta-kct)}+{\rm c.c.}$ $\displaystyle
p=\displaystyle\frac{\chi{u}^{\prime\prime}_{b}}{\rho\mu R}+\delta
p(\eta)e^{i(k\chi+\ell\vartheta-kct)}+{\rm c.c.}$ $\displaystyle\zeta=\delta
ae^{i(k\chi+\ell\vartheta-kct)}+{\rm c.c.},\,\,\,|\delta|\ll 1$

where c.c. denotes a complex conjugate, then the linearized equations of
motion may be cast as

$\displaystyle\left.\begin{array}[]{lcl}U=V=W=0&&\mbox{on $\eta=1$}\\\ \\\
ik(U_{b}-c)U+VU_{b}^{\prime}=-ikP+R^{-1}[U^{\prime\prime}-(k^{2}+\ell^{2})U]&&\\\
\\\ ik(U_{b}-c)V=-P^{\prime}+R^{-1}[V^{\prime\prime}-(k^{2}+\ell^{2})V]&&\\\
\\\ ik(U_{b}-c)W=-i\ell P+R^{-1}[W^{\prime\prime}-(k^{2}+\ell^{2})W]&&\\\ \\\
ikU+V^{\prime}+i\ell W=0&&\\\ \\\
u+au^{\prime}_{b}=U+aU_{b}^{\prime},\,\,\,v=V,\,\,w=W&&\mbox{on $\eta=0$}\\\
\\\ ik(u_{b}-c)a=v,\,\,\,ik(U_{b}-c)a=V&&\mbox{on $\eta=0$}\\\ \\\
u^{\prime}+ikv=\mu(U^{\prime}+ikV),\,\,\,w^{\prime}+i\ell
v=\mu(W^{\prime}+i\ell V)&&\mbox{on $\eta=0$}\\\ \\\ \rho
p-{\mathscr{F}}a-2v^{\prime}/\nu
R-{\mathscr{T}}a(k^{2}+\ell^{2})=\rho\left(P-{\mathscr{F}}a-2V^{\prime}/R\right)&&\mbox{on
$\eta=0$}\\\ \\\ u=v=w=0&&\mbox{on $\eta=d^{-1}$}\\\ \\\
ik(u_{b}-c)u+vu_{b}^{\prime}=-ik\rho p+(\nu
R)^{-1}[u^{\prime\prime}-(k^{2}+\ell^{2})u]&&\\\ \\\ ik(u_{b}-c)v=-\rho
P^{\prime}+(\nu R)^{-1}[v^{\prime\prime}-(k^{2}+\ell^{2})v]&&\\\ \\\
ik(u_{b}-c)w=-i\ell\rho P+(\nu
R)^{-1}[w^{\prime\prime}-(k^{2}+\ell^{2})w]&&\\\ \\\ iku+v^{\prime}+i\ell
w=0&&\\\ \end{array}\right\\}$ (57)

Substituting in equations of motion and their corresponding boundary
conditions, making use of the Squire transformation [14]
$m^{2}=k^{2}+\ell^{2},\hskip 7.11317pt\hat{U}=m^{-1}(kU+\ell W),\hskip
7.11317pt\hat{V}=V,\hskip 7.11317pt\hat{c}=c,\hat{P}=(m/k)P,\hskip
7.11317pt\hat{a}=(k/m)a$, etc. (dropping the hat notation) and introducing the
stream functions $\Phi(\eta)$ and $\phi(\eta)$ such that

$\displaystyle\left.\begin{array}[]{lll}U={\Phi}^{\prime},&&u={\phi}^{\prime}\\\
&&\\\ V=-ik\Phi,&&v=-ik\phi\end{array}\right\\}$ (61)

we obtain, after linearization and some manipulations,

$\displaystyle\left.\begin{array}[]{lll}\Phi={\Phi}^{\prime}=0,&&\mbox{on}\hskip
7.11317pt\eta=1\\\ \\\
ik(U_{b}-c)({\Phi}^{\prime\prime}-k^{2}\Phi)-ik{U}^{\prime\prime}_{b}\Phi=&&\\\
&&\\\ \hskip
85.35826ptR^{-1}(\Phi^{iv}-2k^{2}{\Phi}^{\prime\prime}+k^{4}\Phi)&&\mbox{in
upper fluid}\\\ \\\ \left.\begin{array}[]{l}\phi=\Phi\\\ \\\
{\phi}^{\prime}-\displaystyle\frac{{u}^{\prime}_{b}\phi}{\sigma}={\Phi}^{\prime}-\displaystyle\frac{{U}^{\prime}_{b}\Phi}{\sigma}\\\
\\\ {\phi}^{\prime\prime}+k^{2}\phi=\mu({\Phi}^{\prime\prime}+k^{2}\Phi)\\\
\\\ \displaystyle\frac{1}{\nu
R}({\phi}^{\prime\prime\prime}-3k^{2}{\phi}^{\prime})+ik(\phi
u_{b}-{\phi}^{\prime}\sigma)+i\displaystyle\frac{k\phi}{\sigma}\displaystyle\frac{({\sf
F}+k^{2}{\sf T})}{R^{2}}=\\\ \\\
\displaystyle\frac{\rho}{R}({\Phi}^{\prime\prime\prime}-3k^{2}{\Phi}^{\prime})+\rho
ik(\Phi U_{b}-{\Phi}^{\prime}\sigma)+\rho
i\displaystyle\frac{k\Phi}{\sigma}\displaystyle\frac{{\sf F}}{R^{2}}\\\
\end{array}\right]&&\mbox{on}\hskip 7.11317pt\eta=0\\\ \\\
ik(u_{b}-c)({\phi}^{\prime\prime}-k^{2}\phi)-ik{u}^{\prime\prime}_{b}\phi=&&\\\
&&\\\ \hskip 85.35826pt(\nu
R)^{-1}(\phi^{iv}-2k^{2}{\phi}^{\prime\prime}+k^{4}\phi)&&\mbox{in lower
fluid}\\\ \\\ \phi={\phi}^{\prime}=0,&&\mbox{on}\hskip
7.11317pt\eta=-d^{-1}\\\ \end{array}\right\\}$ (84)

where $\sigma=u_{b}(0)-c=U_{b}(0)-c$. Note that, we have redefined the Froude
and Weber number to be ${\sf F}=gD_{u}^{3}/\nu_{a}^{2}$ and ${\sf
T}=SD_{u}/\rho_{w}\nu_{a}^{2}$, respectively.

Systems similar to that of (84) have been derived [7,15–17]. However, as was
pointed out by Blennerhassett [2], in the derivations [7,17] the relevance of
${\mathscr{F}}$ and ${\mathscr{T}}$ to the stability of two-dimensional waves
was not noted. He further remarks that in earlier derivations [15,16], for the
stability of superposed fluid layers, the terms $u_{b}^{\prime}\phi/\sigma$
and $U_{b}^{\prime}\Phi/\sigma$ was incorrectly omitted from the second
interface condition in the system (84).

We remark that in general the solution of (84), and determining the location
of the neutral stability curve in $(R,k)$–plane is a computational task. The
full details of the numerical methods is given in section 6 below.

##  3 The long-wave limit

Following [2], we consider the stability of the basic flow in the limit of
long wavelength disturbances and expand the eigenfunctions and the eigenvalue
in regular perturbation series in $k$ [7]. Thus we may write

$\displaystyle\left.\begin{array}[]{l}c=c_{0}+kc_{1}+...\\\ \\\
\phi=\phi_{0}+k\phi_{1}+...\\\ \\\ \Phi=\Phi_{0}+k\Phi_{1}+...\\\
\end{array}\right\\}$ (90)

as $k\rightarrow 0$ for a fix value of $R,{\sf F},{\sf T}$. Substituting (90)
in (84) we obtain up to $O(1)$ the following homogeneous system

$\displaystyle\left.\begin{array}[]{lll}\Phi_{0}={\Phi}^{\prime}_{0}=0,&&\mbox{on}\hskip
7.11317pt\eta=1\\\ \\\ 0=R^{-1}\Phi_{0}^{iv}&&\mbox{in upper fluid}\\\ \\\ \\\
\left.\begin{array}[]{l}\phi_{0}=\Phi_{0}\\\ \\\
{\phi}^{\prime}_{0}-\displaystyle\frac{{u}^{\prime}_{b}\phi_{0}}{\sigma}={\Phi}^{\prime}_{0}-\displaystyle\frac{{U}^{\prime}_{b}\Phi_{0}}{\sigma},\hskip
7.11317pt\sigma=u_{b}c_{0}\\\ \\\
{\phi}^{\prime\prime}_{0}=\mu{\Phi}^{\prime\prime}_{0}\\\ \\\ (\nu
R)^{-1}{\phi}^{\prime\prime\prime}_{0}=\rho
R^{-1}{\Phi}^{\prime\prime\prime}_{0}\\\ \end{array}\right]&&\mbox{on}\hskip
7.11317pt\eta=0\\\ \\\ \\\ 0=(\nu R)^{-1}\Phi_{0}^{iv}&&\mbox{in lower
fluid}\\\ \\\ \phi_{0}={\phi}^{\prime}_{0}=0,&&\mbox{on}\hskip
7.11317pt\eta=-d^{-1}\\\ \end{array}\right\\}$ (109)

In the case of pPf, the system (109) has a eigensolution

$\Phi_{0}=(\eta-1)^{2}\left(\eta\displaystyle\frac{d}{2\mu}\displaystyle\frac{\mu-d^{2}}{1+d}+1\right)$

and

$\phi_{0}=(\eta+d^{-1})^{2}\left(\eta\displaystyle\frac{d}{2}\displaystyle\frac{\mu-d^{2}}{1+d}+d^{2}\right).$

The eigenvalue $c_{0}$ is given by

$\displaystyle c_{0}=u_{b}(0)+\displaystyle\frac{2\mu
d(1+d)({U}^{\prime}_{b}-{u}^{\prime}_{b})}{4\mu
d(1+d)^{2}+(\mu-d^{2})^{2}}\hskip 14.22636pt\mbox{provided $\mu\neq 1$}$ (110)

Similarly, in the limit of slow basic flow, when the Reynolds number for the
basic flow is small, we may expand $c\sim R^{-1}c_{0}+c_{1}+...$,
$\phi\sim\phi_{0}+R\phi_{1}+...$ and $\Phi\sim\Phi_{0}+R\Phi_{1}+...$ as
$R\rightarrow 0$ for fixed ${\mathscr{F}}$ and ${\mathscr{T}}$. Thus we obtain
a system similar to (109) whose solution is

$\displaystyle\Phi_{0}=A[\cosh
k(\eta-1)-\cosh\lambda(\eta-1)]+B[\sinh\lambda(\eta-1)-\lambda k^{-1}\sinh
k(\eta-1)]$ $\displaystyle\phi_{0}=a[\cosh
k(\eta+d^{-1})-\cosh\gamma(\eta+d^{-1})]+b[\sinh\gamma(\eta+d^{-1})-\gamma
k^{-1}\sinh k(\eta+d^{-1})]$

where $\lambda^{2}=k(k-ic_{0})$ and $\gamma^{2}=k(k-i\nu c_{0})$, and the
eigenvalues $c_{0}$ and the constants $A,B,a$ and $b$ are determined by
applying the interface conditions at $\eta=0$. This yields a homogeneous
$4\times 4$ set of equations whose non-trivial solutions can only be obtained
for certain values of $c_{0}$.

We next consider $O(k)$. To this order the system (84) reduce to

$\displaystyle\left.\begin{array}[]{lll}\Phi^{(1)}={\Phi^{(1)}}^{\prime}=0,&&\mbox{on}\hskip
7.11317pt\eta=1\\\ \\\
i(U_{b}-c_{0}){\Phi}^{\prime\prime}_{0}-i{U}^{\prime\prime}_{b}\Phi_{0}=R^{-1}\Phi^{iv}_{1}&&\mbox{in
upper fluid}\\\ &&\\\ \left.\begin{array}[]{l}\phi_{0}=\Phi_{0}\\\ \\\
{\phi}^{\prime}_{1}-\displaystyle\frac{{u}^{\prime}_{b}}{\sigma}\left(\phi_{1}+\displaystyle\frac{c_{1}}{\sigma}\phi_{0}\right)=\\\
\\\ \hskip
78.24507pt{\Phi}^{\prime}_{1}-\displaystyle\frac{{U}^{\prime}_{b}}{\sigma}\left(\Phi_{1}+\displaystyle\frac{c_{1}}{\sigma}\Phi_{0}\right)\\\
\\\ {\phi}^{\prime\prime}_{1}=\mu{\Phi}^{\prime\prime}_{1}\\\ \\\ (\nu
R)^{-1}{\phi}^{\prime\prime\prime}_{1}-i\left({\phi}^{\prime}_{1}\sigma-\phi_{0}{u}^{\prime}_{b}\right)+\displaystyle\frac{i\phi_{0}}{\sigma}\displaystyle\frac{{\sf
F}}{R^{2}}=\\\ \\\ \hskip 28.45274pt\rho
R^{-1}{\Phi}^{\prime\prime\prime}_{1}-\rho
i\left({\Phi}^{\prime}_{1}\sigma-\Phi_{0}{U}^{\prime}_{b}\right)+\\\ \\\
\hskip 28.45274pt\displaystyle\frac{\rho i\Phi_{0}{\sf F}}{\sigma R^{2}}\\\
\\\ \end{array}\right]&&\mbox{on}\hskip 7.11317pt\eta=0\\\ \\\
i(u_{b}-c_{0}){\phi}^{\prime\prime}_{0}-i{u}^{\prime\prime}_{b}\phi_{0}=(\nu
R)^{-1}\phi^{iv}_{1}&&\mbox{in lower fluid}\\\ \\\
\phi_{1}={\phi}^{\prime}_{1}=0,&&\mbox{on}\hskip 7.11317pt\eta=-d^{-1}\\\
\end{array}\right\\}$ (134)

Since $\phi^{(0)},\Phi^{(0)}$ and $c^{(0)}$ are real, we can easily see that
$c^{(1)}$ is purely imaginary. Note that the homogeneous part (134) is the
same as (109). Thus, in order that the inhomogeneous system to have a solution
we must satisfy certain solvability criterion. We introduce the adjoint
functions [6], $\Psi$ and $\psi$, and write down the system which is adjoint
to (134)

$\displaystyle\left.\begin{array}[]{lll}\Psi={\Psi}^{\prime}=0,&&\mbox{on}\hskip
7.11317pt\eta=1\\\ \\\ 0=R^{-1}\Psi^{iv}&&\mbox{in upper fluid}\\\ \\\
\left.\begin{array}[]{l}\psi=\Psi,\hskip
14.22636pt{\psi}^{\prime}={\Psi}^{\prime},\hskip
14.22636pt{\psi}^{\prime\prime}=\mu{\Psi}^{\prime\prime}\\\ \\\
{\psi}^{\prime\prime\prime}-\displaystyle\frac{{u}^{\prime}_{b}}{\sigma}{\psi}^{\prime}=\mu\left({\Psi}^{\prime\prime\prime}-\displaystyle\frac{{U}^{\prime}_{b}}{\sigma}{\Psi}^{\prime}\right)\\\
\\\ \end{array}\right]&&\mbox{on}\hskip 7.11317pt\eta=0\\\ \\\ 0=(\nu
R)^{-1}\psi^{iv}&&\mbox{in lower fluid}\\\ \\\
\psi={\psi}^{\prime}=0,&&\mbox{on}\hskip 7.11317pt\eta=-d^{-1}\\\
\end{array}\right\\}$ (148)

with the adjoint condition 222Stuart’s solvability criterion states that if
$L$ is a differential operator and $N$ is its adjoint operator, then if
$L\Phi=f$ and its adjoint is given by $N\Psi=0$ the solvability condition
states
$\int_{a}^{b}\Psi L\Phi\,d\eta=\int_{a}^{b}\Phi\underbrace{N\Psi}_{=0}\,d\eta$
which implies $\int_{a}^{b}\Psi f\,d\eta=0.$

$\displaystyle(\nu
R)^{-1}\displaystyle\int_{-d^{-1}}^{0}\psi\phi^{iv}_{0}\,d\eta+\rho
R^{-1}\displaystyle\int_{0}^{1}\Psi\Phi^{iv}_{0}\,d\eta=(\nu
R)^{-1}\displaystyle\int_{-d^{-1}}^{0}\phi_{0}\psi^{iv}\,d\eta+\rho
R^{-1}\displaystyle\int_{0}^{1}\Phi_{0}\Psi^{iv}\,d\eta=0$

The eigenvalue $c_{0}$ for the adjoint system (148) is the same as that of
(109), while the adjoint functions for pPf are

$\displaystyle\Psi=(\eta-1)^{2}\left[\eta\displaystyle\frac{\mu-d^{2}+4d(1+d)}{2(\mu+d)}+1\right]$
$\displaystyle\psi=(\eta+d^{-1})^{2}\left[\eta\displaystyle\frac{d^{2}(\mu-d^{2})-4\mu
d^{2}(1+d)}{2(\mu+d)}+d^{2}\right]$

When the solvability condition is applied to (109) we obtain the following
expression for $c_{1}$:

$\displaystyle\displaystyle\frac{c_{1}\phi_{0}}{\sigma^{2}}(\nu
R)^{-1}{\psi}^{\prime\prime}({u}^{\prime}_{b}-{U}^{\prime}_{b})$
$\displaystyle=$
$\displaystyle-i\psi\left[(1-\rho)({\phi}^{\prime}_{0}\sigma-\phi_{0}{u}^{\prime}_{b})-\displaystyle\frac{\phi_{0}{\sf
F}(1-\rho)}{\sigma R^{2}}\right]$ (149) $\displaystyle+$ $\displaystyle
i\displaystyle\int_{-d^{-1}}^{0}\psi\left[(u_{b}-c_{0}){\phi}^{\prime\prime}_{0}-{u}^{\prime\prime}_{b}\phi_{0}\right]\,d\eta$
$\displaystyle+$ $\displaystyle\rho
i\displaystyle\int_{0}^{1}\Psi\left[(U_{b}-c_{0}){\Phi}^{\prime\prime}_{0}-{U}^{\prime\prime}_{b}\Phi_{0}\right]\,d\eta$

To $O(k)$, the condition for neutral stability is that $c_{1}=0$. Hence, from
(149) we obtain the following relationship between $k$ and $R$:

$\displaystyle\displaystyle\frac{{\sf
F}}{R^{2}}=\sigma({\phi}^{\prime}_{0}\sigma-{u}^{\prime}_{b})$
$\displaystyle-$
$\displaystyle\sigma\left\\{\underbrace{\displaystyle\frac{1}{1-\rho}\displaystyle\int_{-d^{-1}}^{0}\psi\left[(u_{b}-c_{0}){\phi}^{\prime\prime}_{0}-{u}^{\prime\prime}_{b}\phi_{0}\right]\,d\eta}_{I}\right.$
(150) $\displaystyle+$
$\displaystyle\left.\underbrace{\rho\displaystyle\int_{0}^{1}\Psi\left[(U_{b}-c_{0}){\Phi}^{\prime\prime}_{0}-{U}^{\prime\prime}_{b}\Phi_{0}\right]\,d\eta}_{J}\right\\}$

Equation (150) is a quadratic in $\sigma$ which may be expressed as

$\displaystyle{\phi}^{\prime}_{0}\sigma^{2}-({u}^{\prime}_{b}+I+J)\sigma-\displaystyle\frac{{\sf
F}}{R^{2}}=0$

with its solutions given by

$\displaystyle\sigma=\displaystyle\frac{({u}^{\prime}_{b}+I+J)\pm\sqrt{({u}^{\prime}_{b}+I+J)^{2}+4{\phi}^{\prime}_{0}{\sf
F}/R^{2}}}{2{\phi}^{\prime}_{0}}\equiv\gamma(R)$

Since $\sigma=u_{b}(0)-c$, we obtain the following compact relation between
the wavenumber $k$ and the Reynolds number $R$, i.e.

$\displaystyle k=\displaystyle\frac{ng}{[u_{b}(0)-\gamma(R)]^{2}}$ (151)

##  4 The short-wave limit

We shall now focus our attention to conditions under which short-wavelengths
disturbances may grow along the interface. We shall pose the following
restrictions on the problem

$\displaystyle\left.\begin{array}[]{lclccclcl}{\rm(a)}&&d\ll
1,&&&&{\rm(c)}&&kd^{3}\gg 1\\\ \\\ {\rm(b)}&&k\ll 1,&&&&{\rm(d)}&&kd\ll 1\\\
\end{array}\right\\}$ (155)

To study this, we will adopt an asymptotic analysis for pPf. However, before
we commence with our analysis, it is convenient to express the interfacial
boundary conditions (84) in more compact mathematical form. We do this by
defining

$\displaystyle\left.\begin{array}[]{l}\Theta=U_{b}+\displaystyle\frac{i\rho
k}{\sigma}\displaystyle\frac{\sf F}{R^{2}}\\\ \\\ \Xi=\nu
u_{b}+\displaystyle\frac{1}{\sigma}({\sf F}+k^{2}{\sf
S})\left(\displaystyle\frac{\nu}{R^{2}}+1\right)\\\ \end{array}\right\\}$
(159)

Thus, the second to the fourth interfacial conditions (84) may be written as

$\displaystyle\left.\begin{array}[]{l}\Phi^{\prime}-\phi^{\prime}=\Phi(\Theta-\Xi)/\sigma\\\
\\\ \Phi^{\prime\prime}+k^{2}\Phi=\mu(\phi^{\prime\prime}+k^{2}\phi)\\\ \\\
-ikR(\sigma\Phi^{\prime}+\Xi\Phi)-\Phi^{\prime\prime\prime}+3k^{2}\Phi^{\prime}\\\
\\\ \hskip
28.45274ptikR\rho(\sigma\phi^{\prime}+\Xi\phi)+\mu(\phi^{\prime\prime\prime}-3k^{2}\phi^{\prime})\\\
\\\ \hskip 28.45274pt=ikR({\sf F}+k^{2}{\sf
S})\Phi/\sigma\end{array}\right\\}\mbox{on $\eta=0$}$ (169)

We shall also scale the coordinates $(\chi,\eta)$ with the wavenumber $k$ such
that

$\displaystyle\theta=k\chi,\hskip 28.45274pt\varphi=k\eta$ (170)

and assume the following asymptotic expansions for $c,\phi$ and $\Phi$:

$\displaystyle\left.\begin{array}[]{l}c\sim
c_{0}+k^{-1}c_{1}+k^{-2}c_{2}+...\\\ \\\
\phi\sim\phi_{0}+k^{-1}\phi_{1}+k^{-2}\phi_{2}+...\\\ \\\
\Phi\sim\Phi_{0}+k^{-1}\Phi_{1}+k^{-2}\Phi_{2}+...\\\ \end{array}\right\\}$
(176)

Neglecting terms of $O(d^{-1})$ by invoking the restriction (a) in (24),
replacing the independent variable by $\hat{\varphi}=-\varphi$ and defining
${\mathscr{Y}}=\hat{\varphi}/d$ $(0\leq{\mathscr{Y}}\leq 1)$, the
eigenfunction for the lower fluid will satisfy the following Orr-Sommerfeld
equation

$\displaystyle\phi^{iv}-2(kd)^{2}\phi^{\prime\prime}+(kd)^{4}\phi=\displaystyle\frac{iR\rho
d^{3}}{\mu^{2}}k\left\\{\left[{\mathscr{Y}}-{\mathscr{Y}}^{2}+(1-c)\mu
d^{-1}\right]\left[\phi^{\prime\prime}-(kd)^{2}\phi\right]+2\phi\right\\}$
(177)

in which we anticipate $c\sim 1-O(k^{1/3})$.

Since the effective Reynolds number, $R\rho d^{3}/\mu^{2}$, is large, we may
neglect the viscous terms in (177) and thus reduce it to the Rayleigh equation

$\displaystyle{\mathscr{Y}}(1-{\mathscr{Y}})\phi^{\prime\prime}+2\phi=0$ (178)

Equation (178) has a regular solution

$\displaystyle\phi_{2}^{o}=K{\mathscr{Y}}(1-{\mathscr{Y}})$ (179)

and a solution with logarithmic branch point at ${\mathscr{Y}}=0$ and
${\mathscr{Y}}=1$, namely

$\displaystyle\phi_{1}^{o}=\Lambda\left\\{{\mathscr{Y}}(1-{\mathscr{Y}})\log{\mathscr{Y}}-\log(1-{\mathscr{Y}})+{\mathscr{Y}}^{2}-\textstyle\frac{1}{2}\right\\}$
(180)

where the superscript $o$ in (179) and (180) denotes the outer solution.

Now to recover the boundary layer near the interface, and following the usual
practice in the theory of hydrodynamic stability, we introduce
$\varepsilon=(iR\rho kd^{3}/\mu^{2})^{-1/3}$ and let ${\sf
Y}={\mathscr{Y}}/\varepsilon$. Thus, (177) becomes

$\displaystyle\phi^{iv}-2(\varepsilon kd)^{2}\phi^{\prime\prime}+(\varepsilon
kd)^{4}\phi=\left\\{\left[{\sf Y}(1-\varepsilon{\sf Y})+(1-c)\mu/\varepsilon
d\right]\left[\phi^{\prime\prime}-(\varepsilon
kd)^{2}\phi\right]+2\varepsilon\phi\right\\}\hskip 11.38092pt$ (181)

Expanding $\phi$ in powers of $\varepsilon$ we obtain, at zero order in
$\varepsilon$

$\displaystyle\phi_{0}^{iv}-{\sf Y}\phi_{0}^{\prime\prime}=0$ (182)

and at $O(\varepsilon)$ we get

$\displaystyle\phi_{1}^{iv}-{\sf Y}\phi_{1}^{\prime\prime}=-{\sf
Y}^{2}\phi_{0}^{\prime\prime}+2\phi_{0}$ (183)

The solution of (182) is

$\displaystyle\phi^{I}=m_{1}+m_{2}{\sf Y}+m_{3}{\mathscr{A}}_{1}({\sf
Y},2)+m_{4}{\mathscr{A}}_{2}({\sf Y},2)$ (184)

In (184), ${\mathscr{A}}_{k}({\sf Y},2),\,\,k=1,2$, denote generalized Airy
function

$\displaystyle{\mathscr{A}}_{k}({\sf Y},2)=\displaystyle\frac{1}{2\pi
i}\int_{{\mathscr{L}}_{k}}t^{-2}\exp\left({\sf
Y}t-\textstyle\frac{1}{3}t^{3}\right)\,dt$ (185)

where ${\mathscr{L}}_{k}$ is the appropriate path of integrals shown in figure
1.

Figure 1: The path ${\mathscr{L}}_{0}$, ${\mathscr{L}}_{1}$,
${\mathscr{L}}_{2}$ and ${\mathscr{L}}_{3}$ of integral (185).

For ${\sf Y}\gg 0$ the generalized Airy function ${\mathscr{A}}_{2}({\sf
Y},2)$ is exponentially large. Thus we set the constant $m_{4}=0$ in (184).
The constants $m_{1}$ and $m_{2}$ are found by matching (184) with the
inviscid solutions (179) and (180). However, to do this we require the next-
order solutions to the constant and linear parts in (184).

From (184) we see that

$\displaystyle\phi_{1}^{I}=m_{1}\left\\{1+2\varepsilon[n_{1}+n_{2}{\sf
Y}+B_{3}({\sf Y},2,1)]\right\\}$ (186)
$\displaystyle\phi_{2}^{I}=m_{2}\sf{Y}(1-{\sf Y})$ (187)

where the superscript $I$ represent the inner solution and
${\mathscr{B}}_{3}({\sf Y},2,1)$ is a generalized Airy function given by

$\displaystyle{\mathscr{B}}_{3}({\sf Y},2,1)=\displaystyle\frac{1}{2\pi
i}\int_{\infty\exp[4\pi i/3]}^{(0+)}t^{-2}\ln t\exp\left({\sf
Y}t-\textstyle\frac{1}{3}t^{3}\right)\,dt$ (188)

Now matching $\phi_{2}^{I}$ with $\phi_{2}^{o}$, and $\phi_{1}^{I}$ with
$\phi_{1}^{o}$, we obtain, respectively

$\displaystyle
m_{1}=-\Lambda/2,\,\,\,\,m_{2}=K\varepsilon,\,\,\,\,n_{1}=0,\,\,\,n_{2}=-\log\varepsilon-[\psi(2)+2\pi
i]$ (189)

where $\psi(2)$ is the psi function.

To evaluate the above integral we note that

$\displaystyle{\mathscr{A}}_{1}({\sf Y},2)=-\displaystyle\frac{1}{2}\left[{\sf
Y}^{2}\int_{\infty}^{\sf Y}\mbox{Ai}(t)\,dt-{\sf Y}\mbox{Ai}^{\prime}({\sf
Y})-\mbox{Ai}({\sf Y})\right]$

where Ai is the usual Airy function. Using the asymptotic representation of
the generalized Airy function

${\mathscr{A}}_{1}(z,p)\sim\displaystyle\frac{1}{2\sqrt{\pi}}z^{-(2p+3)/4}\exp\left(-\textstyle\frac{2}{3}z^{3/2}\right)+...$

and the following relationship

${\mathscr{A}}_{2}(z,p)=e^{-2\pi i/3}{\mathscr{A}}_{1}\left(ze^{2\pi
i/3},p\right)$

it can be shown that

$\displaystyle\phi_{1}^{I}=-\varepsilon\Lambda\left\\{\displaystyle\frac{1}{2\varepsilon}-[\log\varepsilon+\psi(2)+2\pi
i]{\sf Y}-\displaystyle\frac{{\sf Y}^{1/4}}{\sqrt{\pi}}\exp\left[-i\left(2{\sf
Y}^{1/2}-\textstyle\frac{3}{4}\pi\right)\right]\right.$
$\displaystyle\left.-\displaystyle\frac{{\sf
Y}^{-5/4}}{\sqrt{\pi}}\exp\left(-\textstyle\frac{2}{3}{\sf
Y}^{3/2}\right)\right\\}$

Note that, solutions of this kind are used as an initial guess for the full
numerical integration.

##  5 Nonlinear stability

The governing equations and boundary conditions for the two-dimensional motion
of two superposed fluids in curvilinear coordinates are

$\displaystyle\left.\begin{array}[]{c}U=V=0\hskip 8.5359pt{\rm on}\hskip
8.5359pt\eta=1\\\ \\\ \displaystyle\frac{\partial U}{\partial
t}+J^{\tiny\frac{1}{2}}\left\\{U\displaystyle\frac{\partial}{\partial\chi}\left(J^{-\tiny\frac{1}{2}}U\right)+V\displaystyle\frac{\partial}{\partial\eta}\left(J^{-\tiny\frac{1}{2}}U\right)\right\\}+\textstyle\frac{1}{2}J^{-1}\displaystyle\frac{\partial
J}{\partial\chi}(U^{2}+V^{2})=\\\ \\\ \hskip
28.45274pt-\displaystyle\frac{\partial
P}{\partial\chi}+R^{-1}\left\\{\displaystyle\frac{\partial
J}{\partial\eta}\left[\displaystyle\frac{\partial}{\partial\eta}\left(J^{-\tiny\frac{1}{2}}U\right)-\displaystyle\frac{\partial}{\partial\chi}\left(J^{-\tiny\frac{1}{2}}V\right)\right]+J\nabla^{2}\left(J^{-\tiny\frac{1}{2}}U\right)\right\\}\\\
\\\ \displaystyle\frac{\partial V}{\partial
t}+J^{\tiny\frac{1}{2}}\left\\{U\displaystyle\frac{\partial}{\partial\chi}\left(J^{-\tiny\frac{1}{2}}V\right)+V\displaystyle\frac{\partial}{\partial\eta}\left(J^{-\tiny\frac{1}{2}}V\right)\right\\}+\textstyle\frac{1}{2}J^{-1}\displaystyle\frac{\partial
J}{\partial\eta}(U^{2}+V^{2})=\\\ \\\ \hskip
28.45274pt-\displaystyle\frac{\partial
P}{\partial\eta}+R^{-1}\left\\{\displaystyle\frac{\partial
J}{\partial\chi}\left[\displaystyle\frac{\partial}{\partial\eta}\left(J^{-\tiny\frac{1}{2}}U\right)-\displaystyle\frac{\partial}{\partial\chi}\left(J^{-\tiny\frac{1}{2}}V\right)\right]-J\nabla^{2}\left(J^{-\tiny\frac{1}{2}}V\right)\right\\}\\\
\\\
\displaystyle\frac{\partial}{\partial\chi}\left(J^{-\tiny\frac{1}{2}}U\right)+\displaystyle\frac{\partial}{\partial\eta}\left(J^{-\tiny\frac{1}{2}}V\right)=0\\\
\\\ \begin{array}[]{clc}\left.\begin{array}[]{l}\hskip 85.35826ptu=U,\hskip
7.11317ptv=V,\hskip 7.11317pt\displaystyle\frac{\partial\zeta}{\partial
t}+J^{\tiny\frac{1}{2}}u\displaystyle\frac{\partial\zeta}{\partial\chi}=v,\hskip
7.11317pt\displaystyle\frac{\partial\zeta}{\partial
t}+J^{\tiny\frac{1}{2}}U\displaystyle\frac{\partial\zeta}{\partial\chi}=V\\\
\\\ \hskip 28.45274pt2J\zeta_{\chi}\displaystyle\frac{\partial
v}{\partial\eta}+\textstyle\frac{1}{2}J^{\tiny\frac{1}{2}}\left(1-J\zeta_{\chi}^{2}\right)\left(\displaystyle\frac{\partial
u}{\partial\eta}+\displaystyle\frac{\partial
v}{\partial\chi}\right)=\mu\left[2J\zeta_{\chi}\displaystyle\frac{\partial
V}{\partial\eta}+\textstyle\frac{1}{2}\left(1-\zeta_{\chi}^{2}\right)\left(\displaystyle\frac{\partial
U}{\partial\eta}+\displaystyle\frac{\partial V}{\partial\chi}\right)\right]\\\
\\\ \rho
p-\displaystyle\frac{\mathscr{F}}{R^{2}}\zeta-\displaystyle\frac{2J^{\tiny\frac{1}{2}}}{\nu
R\left(1-J\zeta_{\chi}^{2}\right)}\left[\left(1-J\zeta_{\chi}^{2}\right)\displaystyle\frac{\partial
v}{\partial\eta}-J^{\tiny\frac{1}{2}}\zeta_{\chi}\left(\displaystyle\frac{\partial
u}{\partial\eta}+\displaystyle\frac{\partial
v}{\partial\chi}\right)\right]+\displaystyle\frac{J\mathscr{T}\left(J^{-\tiny\frac{1}{2}}J^{\tiny\frac{1}{2}}_{\chi}\zeta_{\chi}+\zeta_{\chi\chi}\right)}{R^{2}\left(1+J\zeta_{\chi}^{2}\right)^{\tiny\frac{3}{2}}}\\\
\hskip
42.67912pt=\rho\left\\{P-\displaystyle\frac{\mathscr{F}}{R^{2}}\zeta-\displaystyle\frac{2J^{\tiny\frac{1}{2}}}{\left(1+J\zeta_{\chi}^{2}\right)}\left[\left(1-J\zeta_{\chi}^{2}\right)\displaystyle\frac{\partial
V}{\partial\eta}-J^{\tiny\frac{1}{2}}\zeta_{\chi}\left(\displaystyle\frac{\partial
U}{\partial\eta}+\displaystyle\frac{\partial
V}{\partial\chi}\right)\right]\right\\}\\\
\end{array}\right]&&\\!\\!\\!\\!\mbox{on $\eta=0$}\end{array}\\\ \\\
\displaystyle\frac{\partial u}{\partial
t}+J^{\tiny\frac{1}{2}}\left\\{u\displaystyle\frac{\partial}{\partial\chi}\left(J^{-\tiny\frac{1}{2}}u\right)+v\displaystyle\frac{\partial}{\partial\eta}\left(J^{-\tiny\frac{1}{2}}v\right)\right\\}+\textstyle\frac{1}{2}J^{-1}\displaystyle\frac{\partial
J}{\partial\chi}(u^{2}+v^{2})=\\\ \\\ \hskip
28.45274pt-\rho\displaystyle\frac{\partial p}{\partial\chi}+(\nu
R)^{-1}\left\\{\displaystyle\frac{\partial
J}{\partial\eta}\left[\displaystyle\frac{\partial}{\partial\eta}\left(J^{-\tiny\frac{1}{2}}u\right)-\displaystyle\frac{\partial}{\partial\chi}\left(J^{-\tiny\frac{1}{2}}v\right)\right]+J\nabla^{2}\left(J^{-\tiny\frac{1}{2}}u\right)\right\\}\\\
\\\ \displaystyle\frac{\partial v}{\partial
t}+J^{\tiny\frac{1}{2}}\left\\{u\displaystyle\frac{\partial}{\partial\chi}\left(J^{-\tiny\frac{1}{2}}v\right)+v\displaystyle\frac{\partial}{\partial\eta}\left(J^{-\tiny\frac{1}{2}}v\right)\right\\}+\textstyle\frac{1}{2}J^{-1}\displaystyle\frac{\partial
J}{\partial\eta}(u^{2}+v^{2})=\\\ \\\ \hskip
28.45274pt-\rho\displaystyle\frac{\partial p}{\partial\eta}+(\nu
R)^{-1}\left\\{\displaystyle\frac{\partial
J}{\partial\chi}\left[\displaystyle\frac{\partial}{\partial\eta}\left(J^{-\tiny\frac{1}{2}}u\right)-\displaystyle\frac{\partial}{\partial\chi}\left(J^{-\tiny\frac{1}{2}}v\right)\right]-J\nabla^{2}\left(J^{-\tiny\frac{1}{2}}v\right)\right\\}\\\
\\\
\displaystyle\frac{\partial}{\partial\chi}\left(J^{-\tiny\frac{1}{2}}u\right)+\displaystyle\frac{\partial}{\partial\eta}\left(J^{-\tiny\frac{1}{2}}v\right)=0\\\
\\\ u=v=0\hskip 8.5359pt{\rm on}\hskip 8.5359pt\eta=d^{-1}\\\
\end{array}\right\\}\,\,\,\,\,\,\,$ (222)

where the upper and lower case symbols refer to fluids above and below the
interface, respectively, and $\rho=\rho_{a}/\rho_{w}$; $\mu=\mu_{a}/\mu_{w}$;
$\nu=\nu_{a}/\nu_{w}$ and $d=D_{a}/D_{w}$, being the ratio of the height above
and below the interface. Here, $\rho$, $\mu$, and $\nu$ are dimensionless
density, dynamic viscosity and kinematic viscosity, respectively, with
understanding that the suffices $a$ and $w$ refer to air and water. Also,
${\mathscr{F}}=gD_{a}/U_{\infty}$,
${\mathscr{T}}=S/\rho_{w}U_{\infty}^{2}D_{a}$ and $R=U_{\infty}D_{a}/\nu_{a}$
are the Froude, Weber and Reynolds numbers respectively, with $g$ denoting the
acceleration due to gravity, $S$ the coefficient of surface tension,
$U_{\infty}$ the free stream speed of the air, and $\eta$ is the interface
profile which is the Stokes wave.

For the disturbance near theocratical conditions we expand the frequency and
eigenfunction as a Taylor series. Thus

$-ikc=-ik_{c}c-r+ia_{1}(k-k_{c})-a_{2}(k-k_{c})^{2}+d_{1}(R-R_{c})+...$

Now, since $d(kc)/dk|_{k_{c}}=-a_{1}$, then $-a_{1}$ is the group velocity of
the disturbance; $a_{2}$ is related to the curvature of the nose of the
neutral curve and $d_{1}$ is the exponential growth of the wave for $R>R_{c}$.
Likewise we expand $\Phi$ in the form

$\Phi={\it\Phi}_{1}+(k-k_{c}){\it\Phi}_{10}+(R-R_{c}){\it\Phi}_{11}+(k-k_{c})^{2}{\it\Phi}_{12}+...$

with a similar expansion for $\phi$ in the lower fluid.

It is well known that the basic flow becomes unstable to travelling wave
disturbances when $R=R_{c}$ [3,8]. To examine the growth of these Tollmien-
Schlichting waves in the vicinity of $R_{c}$ we introduce the small parameter

$\displaystyle\varepsilon=d_{1r}|R-R_{c}|\hskip 8.5359pt{\rm as}\hskip
8.5359ptR\rightarrow R_{c}$ (223)

and scaled variable

$\displaystyle\tau=\varepsilon t,\hskip
14.22636pt\xi=\varepsilon^{\tiny\frac{1}{2}}(\chi+a_{1}t)$ (224)

The continuity equations in (222) are satisfied by the streamfunctions
$\Lambda$ and $\lambda$ for the upper and lower fluids, respectively, thus

$\displaystyle
U=\displaystyle\frac{\partial\Lambda}{\partial\eta},\,\,\,\,V=-\displaystyle\frac{\partial\Lambda}{\partial\chi}\hskip
8.5359pt{\rm and}\hskip
8.5359ptu=\displaystyle\frac{\partial\lambda}{\partial\eta},\,\,\,\,v=-\displaystyle\frac{\partial\lambda}{\partial\chi}$
(225)

Introducing the expansions

$\displaystyle\Lambda=\Phi_{0}(\eta,\xi,\tau)+[\Phi_{1}(\eta,\xi,\tau){\sf
E}+\Phi_{2}(\eta,\xi,\tau){\sf E}^{2}+{\rm c.c.}]+...$ (226)
$\displaystyle\lambda=\phi_{0}(\eta,\xi,\tau)+[\phi_{1}(\eta,\xi,\tau){\sf
E}+\phi_{2}(\eta,\xi,\tau){\sf E}^{2}+{\rm c.c.}]+...$ (227)

where c.c. denotes complex conjugate, and ${\sf E}=\exp[ik_{c}(\chi-c_{r}t)].$
We further expand the functions in (227) as

$\displaystyle\phi_{0}(\eta,\xi,\tau)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}u_{b}\,d\eta+\varepsilon\phi_{02}(\eta,\xi,\tau)+\varepsilon^{\tiny\frac{3}{2}}\phi_{03}(\eta,\xi,\tau)+\varepsilon^{2}\phi_{04}(\eta,\xi,\tau)+O\left(\varepsilon^{\tiny\frac{5}{2}}\right)$
$\displaystyle\phi_{1}(\eta,\xi,\tau)$ $\displaystyle=$
$\displaystyle\varepsilon^{\tiny\frac{1}{2}}\phi_{11}(\eta,\xi,\tau)+\varepsilon\phi_{12}(\eta,\xi,\tau)+\varepsilon^{\tiny\frac{3}{2}}\phi_{13}(\eta,\xi,\tau)+\varepsilon^{2}\phi_{14}(\eta,\xi,\tau)$
$\displaystyle+$
$\displaystyle\varepsilon^{\tiny\frac{5}{2}}\phi_{15}(\eta,\xi,\tau)+...$
$\displaystyle\phi_{2}(\eta,\xi,\tau)$ $\displaystyle=$
$\displaystyle\varepsilon\phi_{22}(\eta,\xi,\tau)+\varepsilon^{\tiny\frac{3}{2}}\phi_{23}(y\eta,\xi,\tau)+\varepsilon^{2}\phi_{24}(\eta,\xi,\tau)+O\left(\varepsilon^{\tiny\frac{5}{2}}\right)$

The pressure and interface displacement are expanded as

$\displaystyle p=p_{0}(\chi,\eta,\xi,\tau)+[p_{1}(\eta,\xi,\tau){\sf
E}+p_{2}(\eta,\xi,\tau){\sf E}^{2}+{\rm c.c.}]+...$ (228)
$\displaystyle\zeta=\theta_{0}(\xi,\tau)+[\theta_{1}(\xi,\tau){\sf
E}+\theta_{2}(\xi,\tau)\sf{E}^{2}+{\rm c.c.}]+...$ (229)

where

$\displaystyle p_{0}(\chi,\eta,\xi,\tau)$ $\displaystyle=$
$\displaystyle\displaystyle\frac{u_{b}^{\prime\prime}\chi}{\mu
R_{c}}+\varepsilon^{\tiny\frac{1}{2}}p_{01}(\xi,\tau)+\varepsilon\left\\{p_{02}(\eta,\xi,\tau)-\displaystyle\frac{u_{b}^{\prime\prime}\chi}{\mu
R_{c}^{2}d_{1r}}\right\\}+...+O\left(\varepsilon^{\tiny\frac{3}{2}}\right)$
$\displaystyle\theta_{0}(\xi,\tau)$ $\displaystyle=$
$\displaystyle\varepsilon\theta_{02}(\xi,\tau)+\varepsilon^{\tiny\frac{3}{2}}\theta_{03}(\xi,\tau)+\varepsilon^{2}\theta_{04}(\xi,\tau)+O\left(\varepsilon^{\tiny\frac{5}{2}}\right)$

Substituting in (222) and equating powers of
$O(\varepsilon^{\tiny\frac{1}{2}}{\sf E})$ we obtain the coupled Orr-
Sommerfeld problem

$\displaystyle\left.\begin{array}[]{c}{\it\Phi}_{11}={{\it\Phi}}^{\prime}_{11}=0,\hskip
8.5359pt{\rm on}\hskip 8.5359pt\eta=1\\\ \\\ \mbox{in upper fluid}\\\ \\\
\mathscr{L}(k_{c},{\it\Phi}_{11})=ik_{c}(U_{b}-c_{r})({{\it\Phi}}^{\prime\prime}_{11}-k_{c}^{2}{\it\Phi}_{11})-ik_{c}{U}^{\prime\prime}_{b}{\it\Phi}_{11}-R_{c}^{-1}({\it\Phi}_{11}^{iv}-2k_{c}^{2}{\it\Phi}^{\prime\prime}_{11}+k_{c}^{4}{\it\Phi}_{11})=0\\\
\\\ \mbox{on $\eta=0$}\\\ \\\ \begin{array}[]{l}\phi_{11}={\it\Phi}_{11}\\\
\\\
{\phi}^{\prime}_{11}-\displaystyle\frac{{u}^{\prime}_{b}\phi_{11}}{\sigma_{r}}={{\it\Phi}}^{\prime}_{11}-\displaystyle\frac{{U}^{\prime}_{b}{\it\Phi}_{11}}{\sigma_{r}}\\\
\\\
{\phi}^{\prime\prime}_{11}+k_{c}^{2}\phi_{11}=\mu({{\it\Phi}_{11}}^{\prime\prime}+k_{c}^{2}{\it\Phi}_{11})\\\
\\\ {\sf q}(k_{c},\phi_{11})=\displaystyle\frac{1}{\nu
R_{c}}({\phi}^{\prime\prime\prime}_{11}-3k_{c}^{2}{\phi}^{\prime}_{11})+ik_{c}(\phi_{11}u^{\prime}_{b}-{\phi}^{\prime}_{11}\sigma_{r})+i\displaystyle\frac{k_{c}\phi_{11}}{\sigma_{r}}\displaystyle\frac{(\mathscr{F}+k_{c}^{2}\mathscr{T})}{R_{c}^{2}}=\\\
\hskip
28.45274pt\displaystyle\frac{\rho}{R_{c}}({\it\Phi}^{\prime\prime\prime}_{11}-3k_{c}^{2}{\it\Phi}^{\prime}_{11})+\rho
ik_{c}({\it\Phi}_{11}U_{b}^{\prime}-{\it\Phi}^{\prime}_{11}\sigma_{r})+\rho
i\displaystyle\frac{k_{c}{\it\Phi}_{11}}{\sigma_{r}}\displaystyle\frac{{\mathscr{F}}}{R_{c}^{2}}=\rho\mathscr{Q}(k_{c},{\it\Phi}_{11})\\\
\end{array}\\\ \\\ \mbox{in lower fluid}\\\ \\\
{\mathfrak{l}}(k_{c},{\it\phi}_{11})=ik(u_{b}-c_{r})({\phi}^{\prime\prime}_{11}-k_{c}^{2}\phi_{11})-ik_{c}{u}^{\prime\prime}_{b}\phi_{11}-(\nu
R_{c})^{-1}(\phi_{11}^{iv}-2k_{c}^{2}{\phi}^{\prime\prime}_{11}+k_{c}^{4}\phi_{11})=0\\\
\\\ \phi_{11}={\phi}^{\prime}_{11}=0,\hskip 8.5359pt{\rm on}\hskip
8.5359pt\eta=-d^{-1}\\\ \end{array}\right\\}$ (253) (254)

Note that the system (254) is the same as (84) with $R=R_{c}$, $k=k_{c}$ and
$c=c_{r}$. The solution of (254) may be written as

$\displaystyle\left.\begin{array}[]{l}{\it\Phi}_{11}=A(\tau,\xi){\it\Psi}_{1}(\eta)\\\
\\\ {\it\phi}_{11}=A(\tau,\xi){\it\psi}_{1}(\eta)\end{array}\right\\}$ (258)

where $A(\tau,\xi)$ is the amplitude function to be determined.

The distortion of the mean flow occurs at $O(\varepsilon{\sf E}^{0})$, which
is

$\displaystyle\left.\begin{array}[]{c}{\it\Phi}_{02}={{\it\Phi}}^{\prime}_{02}=0,\hskip
8.5359pt{\rm on}\hskip 8.5359pt\eta=1\\\ \\\
R_{c}^{-1}{\it\Phi}_{02}^{iv}=ik_{c}|A|^{2}({\it\Psi}^{\prime\prime}_{1}\thinspace{\it\Psi}_{1}-{\it\Psi}_{1}\thinspace{\it\Psi}^{\prime}_{1})^{\prime\prime}\hskip
8.5359pt\mbox{in upper fluid}\\\ \\\ \mbox{on $\eta=0$}\\\ \\\
\begin{array}[]{l}{\phi}_{02}+\psi_{1}^{\prime}\thinspace{\zeta}_{1}|A|^{2}={\it\Phi}_{02}+{\it\Psi}_{1}^{\prime}\thinspace{\zeta}_{1}|A|^{2}\\\
\\\
\phi_{02}^{\prime\prime}-\displaystyle\frac{u_{b}^{\prime}\phi_{02}}{u_{b}+a_{1}}+|A|^{2}\left\\{u_{b}^{\prime\prime}|\zeta_{1}|^{2}-\displaystyle\frac{u_{b}^{\prime}}{u_{b}+a_{1}}(\thinspace{u}_{11}\zeta_{1}+\psi_{1}^{\prime}\thinspace{\zeta}_{1})+[\zeta_{1}\thinspace{\psi}_{1}+{\rm
c.c.}]\right\\}=\\\ \\\ \hskip
14.22636pt{\it\Phi}_{02}^{\prime\prime}-\displaystyle\frac{U_{b}^{\prime}{\it\Phi}_{02}}{U_{b}+a_{1}}+|A|^{2}\left\\{U_{b}^{\prime\prime}|\zeta_{1}|^{2}-\displaystyle\frac{U_{b}^{\prime}}{U_{b}+a_{1}}(\thinspace{U}_{11}\zeta_{1}+{\it\Psi}_{1}^{\prime}\thinspace{\zeta}_{1})+[\zeta_{1}\thinspace{\it\Psi}_{1}+{\rm
c.c.}]\right\\}\\\ \\\ \phi_{02}^{\prime\prime}-\nu
R_{c}ik_{c}|A|^{2}(\thinspace{\psi}_{1}\psi_{1}^{\prime}-\psi_{1}\thinspace{\psi}_{1}^{\prime})=\mu\left[{\it\Phi}_{02}^{\prime\prime}-R_{c}ik_{c}|A|^{2}(\thinspace{\it\Psi}_{1}{\it\Psi}_{1}^{\prime}-{\it\Psi}_{1}\thinspace{\it\Psi}_{1}^{\prime})\right]\\\
\\\ (\nu
R_{c})^{-1}{\it\phi}^{\prime\prime\prime}_{02}-ik_{c}|A|^{2}(\thinspace{\psi}_{1}\psi_{1}^{\prime}-\psi_{1}\thinspace{\psi}_{1}^{\prime})^{\prime}=\rho\left[R_{c}^{-1}{\it\Phi}^{\prime\prime\prime}_{02}-ik_{c}|A|^{2}(\thinspace{\it\Psi}_{1}{\it\Psi}_{1}^{\prime}-{\it\Psi}_{1}\thinspace{\it\Psi}_{1}^{\prime})^{\prime}\right]\\\
\end{array}\\\ \\\ (\nu
R_{c})^{-1}{\it\phi}_{02}^{iv}=ik_{c}|A|^{2}({\it\psi}^{\prime\prime}_{1}\thinspace{\it\psi}_{1}-{\it\psi}_{1}\thinspace{\it\psi}^{\prime}_{1})^{\prime\prime}\hskip
8.5359pt\mbox{in lower fluid}\\\ \\\
{\it\phi}_{02}={{\it\phi}}^{\prime}_{02}=0,\hskip 8.5359pt{\rm on}\hskip
8.5359pt\eta=-d^{-1}\\\ \end{array}\right\\}$ (279) (280)

where

$\displaystyle\left.\begin{array}[]{l}\zeta_{1}=-\psi_{1}(0)/\sigma_{r}=-{\it\Psi}_{1}(0)/\sigma_{r}\\\
\\\ u_{11}=\psi_{1}^{\prime}(0)+\zeta_{1}u_{b}^{\prime}(0)\\\ \\\
U_{11}={\it\Psi}_{1}^{\prime}(0)+\zeta_{1}U_{b}^{\prime}(0)\\\ \\\
ik_{c}\rho\mbox{\boldmath$p$}_{11}^{\prime}=(k_{c}^{2}/\nu
R_{c})[\psi_{1}^{\prime\prime}(0)-k_{c}^{2}\psi_{1}(0)]-ik_{c}^{3}\sigma_{r}\psi_{1}(0)\\\
\\\
ik_{c}\mbox{\boldmath$P$}_{11}^{\prime}=(k_{c}^{2}/R_{c})[{\it\Psi}_{1}^{\prime\prime}(0)-k_{c}^{2}{\it\Psi}_{1}(0)]-ik_{c}^{3}\sigma_{r}{\it\Psi}_{1}(0)\\\
\end{array}\right\\}$ (290)

and a tilde denotes complex conjugate. The system (254) has unique solution of
the form

$\displaystyle\left.\begin{array}[]{l}{\it\Phi}_{02}={\mathfrak{F}}(\eta)|A|^{2}\\\
\\\ {\it\phi}_{02}={\mathfrak{f}}(\eta)|A|^{2}\\\ \end{array}\right\\}$ (294)

where

$\displaystyle\begin{array}[]{lcl}{\mathfrak{f}}(\eta)=\displaystyle\int_{-d^{-1}}^{\eta}s(y)\,dy+(\eta+d^{-1})^{2}(b\eta+c),&&{\mathfrak{F}}(\eta)=\displaystyle\int_{1}^{y}S(y)\,dy+(\eta-1)^{2}(B\eta+C)\\\
\\\ s(\eta)=ik_{c}\nu
R_{c}\displaystyle\int_{-d^{-1}}^{\eta}(\thinspace{\psi}_{1}\psi_{1}^{\prime}-\psi_{1}\thinspace{\psi}_{1}^{\prime})\,dy,&&S(\eta)=ik_{c}R_{c}\displaystyle\int_{1}^{\eta}(\thinspace{\it\Psi}_{1}{\it\Psi}_{1}^{\prime}-{\it\Psi}_{1}\thinspace{\it\Psi}_{1}^{\prime})\,dy\\\
\end{array}$ (298)

The constants $B,C,b$ and $c$ are found by applying the interface conditions.

At $O(\varepsilon{\sf{E}}^{2})$ we obtain the system governing the second
harmonic, namely

$\displaystyle\left.\begin{array}[]{c}{\it\Phi}_{22}={{\it\Phi}}^{\prime}_{22}=0,\hskip
8.5359pt{\rm on}\hskip 8.5359pt\eta=1\\\ \\\
{\mathscr{L}}(2k_{c},{\it\Phi}_{22})=ik_{c}A^{2}({\it\Psi}_{1}{\it\Psi}_{1}^{\prime\prime\prime}-{\it\Psi}_{1}^{\prime}{\it\Psi}^{\prime}_{1})\hskip
8.5359pt\mbox{in upper fluid}\\\ \\\ \mbox{on $\eta=0$}\\\ \\\
\begin{array}[]{l}{\phi}_{22}+\textstyle\frac{1}{2}{\zeta}_{1}\psi_{1}^{\prime}A^{2}={\it\Phi}_{22}+\textstyle\frac{1}{2}{\zeta}_{1}{\it\Psi}_{1}^{\prime}A^{2}\\\
\\\
\phi_{22}^{\prime}-\displaystyle\frac{u_{b}^{\prime}\phi_{22}}{\sigma_{r}}+A^{2}\left[\zeta_{1}\psi_{1}^{\prime\prime}-\textstyle\frac{1}{2}u_{b}^{\prime}\zeta_{1}(u_{11}+\psi_{1}^{\prime})+\textstyle\frac{1}{2}u_{b}^{\prime\prime}\zeta_{1}^{2}\right]=\\\
\\\ \hskip
14.22636pt{\it\Phi}_{22}^{\prime}-\displaystyle\frac{U_{b}^{\prime}{\it\Phi}_{22}}{\sigma_{r}}+A^{2}\left[\zeta_{1}{\it\Psi}_{1}^{\prime\prime}-\textstyle\frac{1}{2}U_{b}^{\prime}\zeta_{1}(U_{11}+{\it\Psi}_{1}^{\prime})+\textstyle\frac{1}{2}U_{b}^{\prime\prime}\zeta_{1}^{2}\right]\\\
\\\
\phi_{22}^{\prime\prime}+4k_{c}^{2}\phi_{22}+A^{2}\zeta_{1}\left[\psi_{1}^{\prime}(8k_{c}^{2}+\nu
R_{c}ik_{c}\sigma_{r})-\nu
R_{c}ik\psi_{1}\left(u_{b}^{\prime}+\displaystyle\frac{\mathscr{F}+k_{c}^{2}\mathscr{T}}{\sigma_{r}R_{c}^{2}}\right)\right]=\\\
\hskip
17.07182pt\mu\left\\{{\it\Phi}_{22}^{\prime\prime}+4k_{c}^{2}{\it\Phi}_{22}+A^{2}\zeta_{1}\left[{\it\Psi}_{1}^{\prime}(8k_{c}^{2}+R_{c}ik_{c}\sigma_{r})-R_{c}ik{\it\Psi}_{1}\left(U_{b}^{\prime}+\displaystyle\frac{\mathscr{F}}{\sigma_{r}R_{c}^{2}}\right)\right]\right\\}\\\
\\\ {\sf
q}(2k_{c},\phi_{22})+A^{2}\left[ik_{c}(\psi_{1}\psi_{1}^{\prime\prime}-\psi_{1}^{\prime}\psi_{1}^{\prime})+\displaystyle\frac{ik_{c}\zeta_{1}(u_{11}+\psi_{1}^{\prime})(\mathscr{F}+4k_{c}^{2}\mathscr{T})}{\sigma_{r}R_{c}^{2}}\right.\\\
\hskip
17.07182pt\left.+\zeta_{1}\left(2ik_{c}\rho\mbox{\boldmath$p$}^{\prime}_{11}-\displaystyle\frac{4k_{c}^{2}\psi_{1}^{\prime\prime}}{\nu
R_{c}}\right)\right]\\\ \\\ \hskip
17.07182pt=\rho\left\\{{\mathscr{Q}}(2k_{c},{\it\Phi}_{22})+A^{2}\left[ik_{c}({\it\Psi}_{1}{\it\Psi}_{1}^{\prime\prime}-{\it\Psi}_{1}^{\prime}{\it\Psi}_{1}^{\prime})+\displaystyle\frac{ik_{c}\zeta_{1}(U_{11}+{\it\Psi}_{1}^{\prime})\mathscr{F}}{\sigma_{r}R_{c}^{2}}\right.\right.\\\
\hskip
17.07182pt\left.\left.+\zeta_{1}\left(2ik_{c}\mbox{\boldmath$P$}^{\prime}_{11}-\displaystyle\frac{4k_{c}^{2}\psi_{1}^{\prime\prime}}{R_{c}}\right)\right]\right\\}\\\
\end{array}\\\ \\\
{\mathfrak{l}}(2k_{c},{\it\phi}_{22})=ik_{c}A^{2}({\it\psi}_{1}{\it\psi}_{1}^{\prime\prime\prime}-{\it\psi}_{1}^{\prime}{\it\psi}^{\prime}_{1})\hskip
8.5359pt\mbox{in lower fluid}\\\ \\\
{\it\phi}_{22}={{\it\phi}}^{\prime}_{22}=0,\hskip 8.5359pt\mbox{on
$\eta=-d^{-1}$}\\\ \end{array}\right\\}$ (324) (325)

The solution of (325) may be expressed as

$\displaystyle\left.\begin{array}[]{l}{\it\Phi}_{22}=A^{2}{\it\Psi}_{2}\\\ \\\
{\it\phi}_{22}=A^{2}{\it\psi}_{2}\\\ \end{array}\right\\}$ (329)

with the assumption that the solution is unique.

We emphasize that the adjoint equation in the upper and he lower fluid layer
is the adjoint of the Orr-Sommerfeld equation for a homogeneous fluid. The
adjoint system possesses the following relationship

$\int_{-d^{-1}}^{0}\psi{\mathfrak{l}}(\phi)\,d\eta+\rho\int_{0}^{1}\Psi{\mathscr{L}}(\Phi)\,d\eta=\int_{-d^{-1}}^{0}\phi{\mathfrak{l}}^{+}(\psi)\,d\eta+\rho\int_{0}^{1}\Phi{\mathscr{L}}^{+}(\Psi)\,d\eta$

where the operators ${\mathscr{L}}$ and ${\mathfrak{l}}$ are the same as that
given in (254), and

$\displaystyle{\mathscr{L}}^{+}(\Psi)=ik_{c}(U_{b}-c_{r})(\Psi^{\prime\prime}-k_{c}^{2}\Psi)+2ik_{c}U_{b}^{\prime}\Psi^{\prime}-R_{c}^{-1}(\Psi^{iv}-2k_{c}^{2}\Psi^{\prime\prime}+k_{c}^{4}\Psi)$
$\displaystyle{\mathfrak{l}}^{+}(\psi)=ik_{c}(u_{b}-c_{r})(\psi^{\prime\prime}-k_{c}^{2}\psi)+2ik_{c}u_{b}^{\prime}\psi^{\prime}-(\nu
R_{c})^{-1}(\psi^{iv}-2k_{c}^{2}\psi^{\prime\prime}+k_{c}^{4}\psi)$

The system of equations for $O(\varepsilon{\sf E})$ is given by

$\displaystyle\left.\begin{array}[]{c}{\it\Phi}_{12}={{\it\Phi}}^{\prime}_{12}=0,\hskip
8.5359pt\mbox{on $\eta=1$}\\\ \\\
{\mathscr{L}}(k_{c},{\it\Phi}_{12})=\displaystyle\frac{\partial
A}{\partial\xi}\left\\{{\it\Psi}_{1}[U_{b}^{\prime\prime}+2k_{c}^{2}(U_{b}-c_{r})]-[(U_{b}+a_{1})-4ik_{c}/R_{c}]({\it\Psi}_{1}^{\prime\prime}-k_{c}^{2}{\it\Psi}_{1})\right\\},\hskip
8.5359pt\mbox{in UF}\\\ \\\ \mbox{on $\eta=0$}\\\ \\\
\begin{array}[]{l}{\phi}_{12}={\it\Phi}_{12}\\\ \\\
\phi_{12}^{\prime}-\displaystyle\frac{u_{b}^{\prime}\phi_{12}}{\sigma_{r}}-i\displaystyle\frac{\partial
A}{\partial\xi}\displaystyle\frac{u_{b}(a_{1}+c_{r})}{k_{c}\sigma_{r}^{2}}\psi_{1}={\it\Phi}_{12}^{\prime}-\displaystyle\frac{U_{b}^{\prime}{\it\Phi}_{12}}{\sigma_{r}}-i\displaystyle\frac{\partial
A}{\partial\xi}\displaystyle\frac{U_{b}(a_{1}+c_{r})}{k_{c}\sigma_{r}^{2}}{\it\Psi}_{1}\\\
\\\ {\sf q}(k_{c},\phi_{12})+\displaystyle\frac{\partial
A}{\partial\xi}\left\\{\displaystyle\frac{6ik_{c}}{\nu
R_{c}}-[\psi_{1}^{\prime}(u_{b}-a_{1})-\psi_{1}u_{b}^{\prime}]\right.\\\
\hskip
14.22636pt\left.-\psi_{1}\displaystyle\frac{\mathscr{F}(-u_{b}+2c_{r}+a_{1})-k_{c}^{2}\mathscr{T}(3u_{b}-4c_{r}-a_{1})}{R_{c}^{2}\sigma_{r}^{2}}\right\\}\\\
\hskip
14.22636pt={\mathscr{Q}}(k_{c},{\it\Phi}_{12})+\displaystyle\frac{\partial
A}{\partial\xi}\left\\{\displaystyle\frac{6ik_{c}}{\nu
R_{c}}-[{\it\Psi}_{1}^{\prime}(U_{b}-a_{1})-{\it\Psi}_{1}U_{b}^{\prime}]-{\it\Psi}_{1}\displaystyle\frac{\mathscr{F}(-U_{b}+2c_{r}+a_{1})}{R_{c}^{2}\sigma_{r}^{2}}\right\\}\\\
\end{array}\\\ \\\
{\mathfrak{l}}(k_{c},{\it\phi}_{12})=\displaystyle\frac{\partial
A}{\partial\xi}\left\\{{\it\psi}_{1}[u_{b}^{\prime\prime}+2k_{c}^{2}(u_{b}-c_{r})]-[(u_{b}+a_{1})-4ik_{c}/R_{c}]({\it\psi}_{1}^{\prime\prime}-k_{c}^{2}{\it\psi}_{1})\right\\},\hskip
8.5359pt\mbox{in LF}\\\ \\\ {\it\phi}_{12}={{\it\phi}}^{\prime}_{12}=0,\hskip
8.5359pt\mbox{on $\eta=-d^{-1}$}\\\ \end{array}\right\\}$ (348) (349)

where UF and LF stands for upper and lower fluid, respectively. The solution
of the $O(\varepsilon{\sf E})$ system is

$\displaystyle\left.\begin{array}[]{l}{\it\Phi}_{12}=-i\displaystyle\frac{\partial
A}{\partial\xi}{\it\Psi}_{10}+A_{2}{\it\Psi}_{1}\\\ \\\
{\it\phi}_{12}=-i\displaystyle\frac{\partial
A}{\partial\xi}{\it\psi}_{10}+A_{2}{\it\psi}_{1}\\\ \end{array}\right\\}$
(353)

where $A_{2}$ is and arbitrary function of $\tau$ and $\xi$, and
${\it\Psi}_{10}$ and ${\it\psi}_{10}$ are the solutions of the following
system

$\displaystyle\left.\begin{array}[]{c}\\\
{\it\Psi}_{10}={{\it\Psi}}^{\prime}_{10}=0,\hskip 8.5359pt\mbox{on
$\eta=1$}\\\ \\\
{\mathscr{L}^{+}}(k_{c},{\it\Psi}_{10})=i{\it\Psi}_{1}[U_{b}^{\prime\prime}+2k_{c}^{2}(U_{b}-c_{r})]-[4ik_{c}/R_{c}+i(U_{b}+a_{1})]({\it\Psi}_{1}^{\prime\prime}-k_{c}^{2}{\it\Psi}_{1}),\hskip
8.5359pt\mbox{in UF}\\\ \\\ \mbox{on $\eta=0$}\\\
\begin{array}[]{l}{\psi}_{10}={\it\Psi}_{10}\\\ \\\
\psi_{10}^{\prime}-\displaystyle\frac{u_{b}^{\prime}\psi_{10}}{\sigma_{r}}+\displaystyle\frac{u_{b}(a_{1}+c_{r})}{k_{c}\sigma_{r}^{2}}\psi_{1}={\it\Psi}_{10}^{\prime}-\displaystyle\frac{U_{b}^{\prime}{\it\Psi}_{10}}{\sigma_{r}}+\displaystyle\frac{U_{b}(a_{1}+c_{r})}{k_{c}\sigma_{r}^{2}}{\it\Psi}_{1}\\\
\\\
\psi_{10}^{\prime\prime}+k_{c}^{2}\psi_{10}+2k_{c}\psi_{1}={\it\Psi}_{10}^{\prime\prime}+k_{c}^{2}{\it\Psi}_{10}+2k_{c}{\it\Psi_{1}}\\\
\\\ {\sf q}(k_{c},\phi_{10})-\displaystyle\frac{6k_{c}}{\nu
R_{c}}-i[\psi_{1}^{\prime}(u_{b}-a_{1})-\psi_{1}u_{b}^{\prime}]-i\psi_{1}\displaystyle\frac{\mathscr{F}(-u_{b}+2c_{r}+a_{1})-k_{c}^{2}\mathscr{T}(3u_{b}-4c_{r}-a_{1})}{R_{c}^{2}\sigma_{r}^{2}}\\\
\\\ \hskip
14.22636pt={\mathscr{Q}}(k_{c},{\it\Psi}_{10})-\displaystyle\frac{6ik_{c}}{R_{c}}-i[{\it\Psi}_{1}^{\prime}(U_{b}-a_{1})-{\it\Psi}_{1}U_{b}^{\prime}]-i{\it\Psi}_{1}\displaystyle\frac{\mathscr{F}(-U_{b}+2c_{r}+a_{1})}{R_{c}^{2}\sigma_{r}^{2}}\\\
\end{array}\\\ \\\
{\mathfrak{l}}^{+}(k_{c},{\it\psi}_{10})=i{\it\psi}_{1}[u_{b}^{\prime\prime}+2k_{c}^{2}(u_{b}-c_{r})]-\left[\displaystyle\frac{4ik_{c}}{\nu
R_{c}}+i(u_{b}+a_{1})\right]({\it\psi}_{1}^{\prime\prime}-k_{c}^{2}{\it\psi}_{1}),\hskip
8.5359pt\mbox{in LF}\\\ \\\ {\it\psi}_{10}={{\it\psi}}^{\prime}_{10}=0,\hskip
8.5359pt\mbox{on $\eta=-d^{-1}$}\\\ \end{array}\right\\}$ (374) (375)

since the adjoint to the system (84), at critical conditions, satisfies

$\displaystyle\left.\begin{array}[]{lll}\Psi={\Psi}^{\prime}=0,&&\mbox{on}\hskip
7.11317pt\eta=1\\\ \\\ {\mathscr{L}}^{+}(\Psi)=0&&\mbox{in upper fluid}\\\ \\\
\left.\begin{array}[]{l}\psi=\Psi\\\ \\\ {\psi}^{\prime}={\Psi}^{\prime}\\\
\\\
{\psi}^{\prime\prime}+k_{c}^{2}\psi=\mu({\Psi}^{\prime\prime}+k_{c}^{2}\Psi)\\\
\\\ -\displaystyle\frac{{\psi}^{\prime\prime\prime}-3k_{c}^{2}\psi}{\nu
R_{c}}+\displaystyle\frac{u_{b}^{\prime}(\psi^{\prime\prime}+k_{c}^{2}\psi)}{\sigma_{r}\nu
R_{c}}+ik_{c}\sigma_{r}\psi^{\prime}-\displaystyle\frac{ik_{c}\psi({\sf
F}+k_{c}^{2}{\sf T})}{\sigma_{r}R_{c}^{2}}=\\\ \\\
\rho\left[-\displaystyle\frac{{\Psi}^{\prime\prime\prime}-3k_{c}^{2}\Psi}{R_{c}}+\displaystyle\frac{U_{b}^{\prime}(\Psi^{\prime\prime}+k_{c}^{2}\Psi)}{\sigma_{r}R_{c}}+ik_{c}\sigma_{r}\Psi^{\prime}-\displaystyle\frac{ik_{c}\Psi{\sf
F}}{\sigma_{r}R_{c}^{2}}\right]\\\ \end{array}\right]&&\mbox{on}\hskip
7.11317pt\eta=0\\\ \\\ {\mathfrak{l}}^{+}(\psi)=0&&\mbox{in lower fluid}\\\
\\\ \psi={\psi}^{\prime}=0,&&\mbox{on}\hskip 7.11317pt\eta=-d^{-1}\\\
\end{array}\right\\}$ (394) (395)

Next, at $O\left(\varepsilon^{\tiny\frac{3}{2}}{\sf E}\right)$, we obtain the
system of equations

$\displaystyle\left.\begin{array}[]{c}\\\
{\it\Phi}_{13}={{\it\Phi}}^{\prime}_{13}=0,\hskip 8.5359pt\mbox{on
$\eta=1$}\\\ \\\
{\mathscr{L}}(k_{c},{\it\Psi}_{13})=({\it\Psi}^{\prime}_{1}-k_{c}^{2}{\it\Psi}_{1})\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)-\displaystyle\frac{A({\it\Psi}_{1}^{iv}-2k_{c}^{2}{\it\Psi}^{\prime\prime}_{1}+k_{c}^{4}{\it\Psi}_{1})}{d_{1r}R_{c}^{2}}\\\
\hskip 14.22636pt-{\mathscr{M}}\left(i\displaystyle\frac{\partial
A_{2}}{\partial\xi}+\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}\right)+A|A|^{2}{\mathscr{G}}(\eta),\hskip
8.5359pt\mbox{in UF}\\\ \\\ \mbox{on $\eta=0$}\\\ \\\
\begin{array}[]{l}\phi_{13}-{\sf
w}A|A|^{2}={\it\Phi}_{13}-{\mathscr{W}}A|A|^{2},\\\ \\\
\phi_{13}^{\prime}-\displaystyle\frac{u_{b}^{\prime}\phi_{13}}{\sigma_{r}}+\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)\displaystyle\frac{iu_{b}^{\prime}\psi_{1}}{k_{c}\sigma_{r}^{2}}-{\sf
m}\left(i\displaystyle\frac{\partial
A_{2}}{\partial\xi}+\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}\right)\\\
\hskip
14.22636pt+A|A|^{2}\left[{\mathscr{U}}+\displaystyle\frac{u_{b}^{\prime}}{\sigma_{r}}(K_{l})\right]\\\
\hskip
14.22636pt={\it\Phi}_{13}^{\prime}-\displaystyle\frac{U_{b}^{\prime}{\it\Phi}_{13}}{\sigma_{r}}+\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)\displaystyle\frac{iU_{b}^{\prime}{\it\Psi}_{1}}{k_{c}\sigma_{r}^{2}}-{\mathscr{M}}\left(i\displaystyle\frac{\partial
A_{2}}{\partial\xi}+\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}\right)\\\
\hskip
14.22636pt+A|A|^{2}\left[{\mathscr{U}}+\displaystyle\frac{U_{b}^{\prime}}{\sigma_{r}}(K_{u})\right]\\\
\\\ \phi_{13}^{\prime\prime}+k_{c}^{2}\phi_{13}-{\sf
m}\left(i\displaystyle\frac{\partial
A_{2}}{\partial\xi}+\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}\right)+A|A|^{2}{\sf
s}\\\ \hskip
14.22636pt=\mu\left[{\it\Phi}_{13}^{\prime\prime}+k_{c}^{2}{\it\Phi}_{13}-{\mathscr{M}}\left(i\displaystyle\frac{\partial
A_{2}}{\partial\xi}+\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}\right)+A|A|^{2}{\mathscr{S}}\right],\\\
\\\ {\sf
q}(k_{c},\phi_{13})+\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)\left[\psi_{1}^{\prime}+\psi_{1}\displaystyle\frac{(\mathscr{F}+k_{c}^{2}\mathscr{T})}{\sigma_{r}R_{c}^{2}}\right]-{\sf
m}\left(\displaystyle\frac{\partial
A_{2}}{\partial\xi}+i\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}\right)\\\
\hskip
14.22636pt+\displaystyle\frac{Aik_{c}}{d_{1r}R_{c}}\left(u_{b}^{\prime}\psi_{1}-\sigma_{r}\psi_{1}^{\prime}-\psi_{1}\displaystyle\frac{\mathscr{F}+k_{c}^{2}\mathscr{T}}{\sigma_{r}R_{c}^{2}}\right)+{\sf
n}A|A|^{2},\\\ \hskip
14.22636pt=\rho\left\\{\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)\left[{\it\Psi}_{1}^{\prime}+{\it\Psi}_{1}\displaystyle\frac{\mathscr{F}}{\sigma_{r}R_{c}^{2}}\right]-{\mathscr{M}}\left(\displaystyle\frac{\partial
A_{2}}{\partial\xi}+i\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}\right)\right.\\\
\left.\hskip
14.22636pt+\displaystyle\frac{Aik_{c}}{d_{1r}R_{c}}\left(U_{b}^{\prime}{\it\Psi}_{1}-\sigma_{r}{\it\Psi}_{1}^{\prime}-{\it\Psi}_{1}\displaystyle\frac{\mathscr{F}}{\sigma_{r}R_{c}^{2}}\right)+{\mathscr{N}}A|A|^{2}\right\\},\\\
\\\ \end{array}\\\
{\mathfrak{l}}(k_{c},{\it\psi}_{13})=({\psi}^{\prime}_{1}-k_{c}^{2}\psi_{1})\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)-\displaystyle\frac{A(\psi_{1}^{iv}-2k_{c}^{2}{\psi}^{\prime\prime}_{1}+k_{c}^{4}\psi_{1})}{d_{1r}\nu
R_{c}^{2}}\\\ \hskip 14.22636pt-{\sf m}\left(i\displaystyle\frac{\partial
A_{2}}{\partial\xi}+\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}\right)+A|A|^{2}{\sf
g}(\eta),\hskip 8.5359pt\mbox{in LF}\\\ \\\
{\it\psi}_{13}={{\it\psi}}^{\prime}_{13}=0,\hskip 8.5359pt\mbox{on
$\eta=-d^{-1}$}\\\ \end{array}\right\\}$ (424) (425)

where

$\displaystyle{\mathscr{U}}=\zeta_{02}\psi_{1}^{\prime\prime}+{\mathfrak{f}}^{\prime\prime}\zeta_{1}+\zeta_{22}\thinspace{\psi}_{1}^{\prime\prime}+\thinspace{\zeta}_{1}\psi_{2}^{\prime}+\psi_{1}^{\prime\prime\prime}|\zeta_{1}|^{2}+u_{b}^{\prime\prime}(\zeta_{02}\zeta_{1}+\zeta_{22}\thinspace{\zeta}_{1})+\textstyle\frac{1}{2}\thinspace{\psi}_{1}^{\prime\prime\prime}\zeta_{1}^{2}$
$\displaystyle{\mathscr{U}}=\zeta_{02}{\it\Psi}_{1}^{\prime\prime}+{\mathfrak{F}}^{\prime\prime}\zeta_{1}+\zeta_{22}\thinspace{\it\Psi}_{1}^{\prime\prime}+\thinspace{\zeta}_{1}{\it\Psi}_{2}^{\prime}+{\it\Psi}_{1}^{\prime\prime\prime}|\zeta_{1}|^{2}+U_{b}^{\prime\prime}(\zeta_{02}\zeta_{1}+\zeta_{22}\thinspace{\zeta}_{1})+\textstyle\frac{1}{2}\thinspace{\it\Psi}_{1}^{\prime\prime\prime}\zeta_{1}^{2}$

Proceeding as above we obtain the following set of equations at
$O\left(\varepsilon^{\tiny\frac{5}{2}}{\sf E}^{2}\right)$.

$\displaystyle\left.\begin{array}[]{c}\\\
{\it\Phi}_{23}={{\it\Phi}}^{\prime}_{23}=0,\hskip 8.5359pt\mbox{on
$\eta=1$}\\\ \\\
{\mathscr{L}}(k_{c},{\it\Phi}_{23})=\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)({\it\Psi}^{\prime\prime}_{1}-k_{c}^{2}{\it\Psi}_{1})-\displaystyle\frac{A({\it\Psi}_{1}^{iv}-2k_{c}^{2}{\it\Psi}^{\prime\prime}_{1}+k_{c}^{4}{\it\Psi}_{1})}{d_{1r}R_{c}^{2}}\\\
\hskip 14.22636pt-i{\mathscr{\alpha}}_{u}\displaystyle\frac{\partial
A_{2}}{\partial\xi}-{\mathscr{\beta}}_{u}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}+A|A|^{2}{\mathscr{G}}(\eta)+A|A|^{4}{\mathscr{H}}(\eta),\hskip
8.5359pt\mbox{in UF}\\\ \\\ \mbox{on $\eta=0$}\\\ \\\
\begin{array}[]{l}\phi_{23}-{\sf v}A|A|^{2}-{\sf
w}A|A|^{4}={\it\Phi}_{23}-{\mathscr{V}}A|A|^{2}-{\mathscr{W}}A|A|^{4},\\\ \\\
\phi_{23}^{\prime}-\displaystyle\frac{u_{b}^{\prime}\phi_{23}}{\sigma_{r}}\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)\displaystyle\frac{iu_{b}^{\prime}\psi_{1}}{k_{c}\sigma_{r}^{2}}--i{\mathscr{\alpha}}_{l}\displaystyle\frac{\partial
A_{2}}{\partial\xi}-{\mathscr{\beta}}_{l}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}\\\
\hskip
14.22636pt+A|A|^{2}\left[{\mathscr{U}}+\displaystyle\frac{u_{b}^{\prime}}{\sigma_{r}}(K_{l})\right]+A|A|^{4}\left[{\mathscr{U}}+\displaystyle\frac{u_{b}^{\prime}}{\sigma_{r}}(J_{l})\right]\\\
\hskip
14.22636pt={\it\Phi}_{23}^{\prime}-\displaystyle\frac{U_{b}^{\prime}{\it\Phi}_{23}}{\sigma_{r}}\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)\displaystyle\frac{iU_{b}^{\prime}{\it\Psi}_{1}}{k_{c}\sigma_{r}^{2}}-i{\mathscr{\alpha}}_{u}\displaystyle\frac{\partial
A_{2}}{\partial\xi}\\\ \hskip
14.22636pt-{\mathscr{\beta}}_{u}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}+A|A|^{2}\left[{\mathscr{U}}+\displaystyle\frac{U_{b}^{\prime}}{\sigma_{r}}(K_{u})\right]+A|A|^{4}\left[{\mathscr{U}}+\displaystyle\frac{U_{b}^{\prime}}{\sigma_{r}}(J_{u})\right],\\\
\\\
\phi_{23}^{\prime\prime}+k_{c}^{2}\phi_{23}-i{\mathscr{\alpha}}_{l}\displaystyle\frac{\partial
A_{2}}{\partial\xi}+A|A|^{2}{\sf s}+A|A|^{4}{\sf r}\\\ \hskip
14.22636pt=\mu\left[{\it\Phi}_{23}^{\prime\prime}+k_{c}^{2}{\it\Phi}_{23}-i{\mathscr{\alpha}}_{u}\displaystyle\frac{\partial
A_{2}}{\partial\xi}+A|A|^{2}{\mathscr{S}}+A|A|^{4}{\mathscr{R}}\right],\\\ \\\
{\sf
q}(k_{c},\phi_{23})+\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)\left[\psi_{1}^{\prime}+\psi_{1}\displaystyle\frac{(\mathscr{F}+k_{c}^{2}\mathscr{T})}{\sigma_{r}^{2}R_{c}^{3}}\right]-i{\mathscr{\alpha}}_{l}\displaystyle\frac{\partial
A_{2}}{\partial\xi}-{\mathscr{\beta}}_{l}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}\\\
\hskip
14.22636pt+\displaystyle\frac{Aik_{c}}{d_{1r}R_{c}}\left(u_{b}^{\prime}\psi_{1}-\sigma_{r}\psi_{1}^{\prime}-\psi_{1}\displaystyle\frac{\mathscr{F}+k_{c}^{2}\mathscr{T}}{\sigma_{r}R_{c}^{2}}\right)+{\sf
n}A|A|^{2}+{\sf m}A|A|^{4}\\\ \hskip 14.22636pt=\rho\left\\{{\sf
Q}(k_{c},{\it\Phi}_{23})+\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)\left[{\it\Psi}_{1}^{\prime}+{\it\Psi}_{1}\displaystyle\frac{\mathscr{F}}{\sigma_{r}R_{c}^{2}}\right]-i{\mathscr{\alpha}}_{u}\displaystyle\frac{\partial
A_{2}}{\partial\xi}-{\mathscr{\beta}}_{u}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}\right.\\\
\left.\hskip
14.22636pt+\displaystyle\frac{Aik_{c}}{d_{1r}R_{c}}\left(U_{b}^{\prime}{\it\Psi}_{1}-\sigma_{r}{\it\Psi}_{1}^{\prime}-{\it\Psi}_{1}\displaystyle\frac{\mathscr{F}}{\sigma_{r}R_{c}^{2}}\right)+{\mathscr{N}}A|A|^{2}+{\mathscr{M}}A|A|^{4}\right\\},\\\
\\\ \end{array}\\\
{\mathfrak{l}}(k_{c},{\it\psi}_{23})=({\psi}^{\prime\prime}_{1}-k_{c}^{2}\psi_{1})\left(a_{2}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}-\displaystyle\frac{\partial
A}{\partial\tau}\right)-\displaystyle\frac{A(\psi_{1}^{iv}-2k_{c}^{2}{\psi}^{\prime\prime}_{1}+k_{c}^{4}\psi_{1})}{d_{1r}\nu
R_{c}^{2}}\\\ \hskip
14.22636pt-i{\mathscr{\alpha}}_{l}\displaystyle\frac{\partial
A_{2}}{\partial\xi}-{\mathscr{\beta}}_{l}\displaystyle\frac{\partial^{2}A}{\partial\xi^{2}}+A|A|^{2}{\sf
g}(\eta)+A|A|^{4}{\sf h}(\eta),\hskip 8.5359pt\mbox{in LF}\\\ \\\
{\it\phi}_{23}={{\it\phi}}^{\prime}_{23}=0,\hskip 8.5359pt\mbox{on
$\eta=-d^{-1}$}\\\ \end{array}\right\\}$ (454) (455)

where

$\displaystyle{\sf
w}=-\zeta_{02}\psi_{2}^{\prime}+\zeta_{22}\thinspace{\psi}_{2}^{\prime}-2\thinspace{\zeta}_{1}\psi_{3}^{\prime}-\psi_{2}^{\prime\prime}|\zeta_{1}|^{2}+\thinspace{\psi}_{2}^{\prime\prime}\zeta_{1}^{2}$
$\displaystyle{\mathscr{W}}=-\zeta_{02}{\it\Psi}_{2}^{\prime}+\zeta_{22}\thinspace{\it\Psi}_{2}^{\prime}-2\thinspace{\zeta}_{1}{\it\Psi}_{3}^{\prime}-{\it\Psi}_{2}^{\prime\prime}|\zeta_{1}|^{2}+\thinspace{\it\Psi}_{2}^{\prime\prime}\zeta_{1}^{2}$

and

$\displaystyle{\sf h}(\eta)$ $\displaystyle=$
$\displaystyle\psi_{3}^{\prime}(\thinspace{\psi}_{2}^{\prime\prime}-k_{c}^{2}\thinspace{\psi}_{2})+\thinspace{\psi}_{3}(\thinspace{\psi}_{2}^{\prime\prime\prime}-k_{c}^{2}\thinspace{\psi}_{2}^{\prime})-2\thinspace{\psi}_{2}^{\prime}(\psi_{3}^{\prime\prime}-4k_{c}^{2}\psi_{3})$
$\displaystyle-$
$\displaystyle\thinspace{\psi}_{2}(\psi_{3}^{\prime\prime\prime}-4k_{c}^{2}\psi_{3}^{\prime})-{\mathfrak{f}}^{\prime}(\psi_{2}^{\prime\prime}-k_{c}^{2}\psi_{2})+{\mathfrak{f}}^{\prime\prime\prime}\psi_{2}$
$\displaystyle{\mathscr{H}}(\eta)$ $\displaystyle=$
$\displaystyle{\it\Psi}_{3}^{\prime}(\thinspace{\it\Psi}_{2}^{\prime\prime}-k_{c}^{2}\thinspace{\it\Psi}_{2})+\thinspace{\it\Psi}_{3}(\thinspace{\it\Psi}_{2}^{\prime\prime\prime}-k_{c}^{2}\thinspace{\it\Psi}_{2}^{\prime})-2\thinspace{\it\Psi}_{2}^{\prime}({\it\Psi}_{3}^{\prime\prime}-4k_{c}^{2}{\it\Psi}_{3})$
$\displaystyle-$
$\displaystyle\thinspace{\it\Psi}_{2}({\it\Psi}_{3}^{\prime\prime\prime}-4k_{c}^{2}{\it\Psi}_{3}^{\prime})-{\mathfrak{F}}^{\prime}({\it\Psi}_{2}^{\prime\prime}-k_{c}^{2}{\it\Psi}_{2})+{\mathfrak{F}}^{\prime\prime\prime}{\it\Psi}_{2}$
$\displaystyle\zeta_{1}$ $\displaystyle=$
$\displaystyle-\psi_{1}/\sigma_{r}=-{\it\Psi}_{1}/\sigma_{r}$
$\displaystyle\zeta_{02}$ $\displaystyle=$
$\displaystyle-[{\mathfrak{f}}(0)+\thinspace{u}_{11}\zeta_{1}+\psi_{1}^{\prime}\thinspace{\zeta}_{1}]/(u_{b}+a_{1})$
$\displaystyle\zeta_{22}$ $\displaystyle=$
$\displaystyle-[\psi_{2}+\textstyle\frac{1}{2}\zeta_{1}(u_{11}+\psi_{1}^{\prime})]/\sigma_{r}$
$\displaystyle u_{11}$ $\displaystyle=$
$\displaystyle\psi_{1}^{\prime}+u_{b}^{\prime}\zeta_{1}$ $\displaystyle K_{l}$
$\displaystyle=$ $\displaystyle{\sf
v}-(u_{02}\zeta_{1}-u_{22}\thinspace{\zeta}_{1}+2\thinspace{u}_{11}\zeta_{22})$
$\displaystyle K_{u}$ $\displaystyle=$
$\displaystyle{\mathscr{V}}-(u_{02}\zeta_{1}-u_{22}\thinspace{\zeta}_{1}+2\thinspace{u}_{11}\zeta_{22})$
$\displaystyle{\sf v}$ $\displaystyle=$
$\displaystyle-\zeta_{02}\psi_{1}^{\prime}+\zeta_{22}\thinspace{\psi}_{1}^{\prime}-2\thinspace{\zeta}_{1}\psi_{2}^{\prime}-\psi_{1}^{\prime\prime}|\zeta_{1}|^{2}+\thinspace{\psi}_{1}^{\prime\prime}\zeta_{1}^{2}$
$\displaystyle{\mathscr{V}}$ $\displaystyle=$
$\displaystyle-\zeta_{02}{\it\Psi}_{1}^{\prime}+\zeta_{22}\thinspace{\it\Psi}_{1}^{\prime}-2\thinspace{\zeta}_{1}{\it\Psi}_{2}^{\prime}-{\it\Psi}_{1}^{\prime\prime}|\zeta_{1}|^{2}+\thinspace{\it\Psi}_{1}^{\prime\prime}\zeta_{1}^{2}$
$\displaystyle u_{02}$ $\displaystyle=$
$\displaystyle{\mathfrak{f}}^{\prime}(0)+u_{b}^{\prime}\zeta_{02}+(\zeta_{1}\thinspace{\psi}_{1}^{\prime\prime}+{\rm
c.c.})+u_{b}^{\prime\prime}|\zeta_{1}|^{2}$ $\displaystyle u_{22}$
$\displaystyle=$
$\displaystyle\psi_{2}^{\prime}+u_{b}^{\prime}\zeta_{22}+\zeta_{1}\psi_{1}^{\prime\prime}+\textstyle\frac{1}{2}u_{b}^{\prime\prime}\zeta_{1}^{2}$
$\displaystyle{\mathscr{\alpha}}_{u}$ $\displaystyle=$ $\displaystyle
i{\it\Psi}_{1}\left[{U}^{\prime\prime}_{b}+2k_{c}^{2}(U_{b}-c_{r})\right]-\left[\frac{4k_{c}}{R_{c}}+i(U_{b}+a_{1})\right]\left({\it\Psi}^{\prime\prime}_{1}-k_{c}^{2}{\it\Psi}_{1}\right)$
$\displaystyle{\mathscr{\alpha}}_{l}$ $\displaystyle=$ $\displaystyle
i\psi_{1}\left[{u}^{\prime\prime}_{b}+2k_{c}^{2}(u_{b}-c_{r})\right]-\left[\frac{4k_{c}}{\nu
R_{c}}+i(u_{b}+a_{1})\right]\left({\psi}^{\prime\prime}_{1}-k_{c}^{2}\psi_{1}\right)$
$\displaystyle{\mathscr{\beta}}_{u}$ $\displaystyle=$ $\displaystyle
i{\it\Psi}_{10}\left[{U}^{\prime\prime}_{b}+2k_{c}^{2}(U_{b}-c_{r})\right]-\left[\frac{4k_{c}}{R_{c}}+i(U_{b}+a_{1})\right]\left({\it\Psi}^{\prime\prime}_{10}-k_{c}^{2}{\it\Psi}_{10}\right)$
$\displaystyle+$ $\displaystyle
a_{2}\left({\it\Psi}^{\prime\prime}_{1}-k_{c}^{2}{\it\Psi}_{1}\right)+ik_{c}{\it\Psi}_{1}(3U_{b}-c_{r}+2a_{1})-\frac{2}{R_{c}}\left({\it\Psi}^{\prime\prime}_{1}-3k_{c}^{2}{\it\Psi}_{1}\right)$
$\displaystyle{\mathscr{\beta}}_{l}$ $\displaystyle=$ $\displaystyle
i\psi_{10}\left[{u}^{\prime\prime}_{b}+2k_{c}^{2}(u_{b}-c_{r})\right]-\left[\frac{4k_{c}}{\nu
R_{c}}+i(u_{b}+a_{1})\right]\left({\psi}^{\prime\prime}_{10}-k_{c}^{2}\psi_{10}\right)$
$\displaystyle+$ $\displaystyle
a_{2}\left({\psi}^{\prime\prime}_{1}-k_{c}^{2}\psi_{1}\right)+ik_{c}\psi_{1}(3u_{b}-c_{r}+2a_{1})-\frac{2}{\nu
R_{c}}\left({\psi}^{\prime\prime}_{1}-3k_{c}^{2}\psi_{1}\right)$
$\displaystyle{\sf m}$ $\displaystyle=$ $\displaystyle
ik_{c}(-{\mathfrak{f}}^{\prime}\psi_{2}^{\prime}-\thinspace{\psi}_{2}^{\prime}\psi_{3}^{\prime}+2\psi_{3}\thinspace{\psi}_{2}^{\prime\prime}-\psi_{3}^{\prime\prime}\thinspace{\psi}_{2})-ik_{c}K_{l}({\mathscr{F}}+k_{c}^{2}{\mathscr{T}})/\sigma_{r}R_{c}^{2}$
$\displaystyle+$ $\displaystyle
3Tik_{c}^{5}\zeta_{1}|\zeta_{1}|^{2}/2R_{c}^{2}+ik_{c}\rho{\sf p}_{13}$
$\displaystyle{\mathscr{M}}$ $\displaystyle=$ $\displaystyle
ik_{c}(-{\mathfrak{f}}^{\prime}{\it\Psi}_{2}^{\prime}-\thinspace{\it\Psi}_{2}^{\prime}{\it\Psi}_{3}^{\prime}+2{\it\Psi}_{3}\thinspace{\it\Psi}_{2}^{\prime\prime}-{\it\Psi}_{3}^{\prime\prime}\thinspace{\it\Psi}_{2})-ik_{c}K_{l}{\mathscr{F}}/\sigma_{r}R_{c}^{2}$
$\displaystyle+$ $\displaystyle
3Tik_{c}^{5}\zeta_{1}|\zeta_{1}|^{2}/2R_{c}^{2}+ik_{c}\rho{\mathscr{P}}_{13}$
$\displaystyle{\sf p}_{13}$ $\displaystyle=$ $\displaystyle
ik_{c}\rho\left[(\mbox{\boldmath$p$}_{02}^{\prime}-\mbox{\boldmath$p$}_{22}^{\prime})\zeta_{1}+(\zeta_{02}-\zeta_{22})\mbox{\boldmath$p$}_{11}^{\prime}+(|\zeta_{1}|^{2}-\textstyle\frac{1}{2}\zeta_{1}^{2})\mbox{\boldmath$p$}_{11}^{\prime\prime}\right]$
$\displaystyle{\mathscr{P}}_{13}$ $\displaystyle=$ $\displaystyle
ik_{c}\rho\left[(\mbox{\boldmath$P$}_{02}^{\prime}-\mbox{\boldmath$P$}_{22}^{\prime})\zeta_{1}+(\zeta_{02}-\zeta_{22})\mbox{\boldmath$P$}_{11}^{\prime}+(|\zeta_{1}|^{2}-\textstyle\frac{1}{2}\zeta_{1}^{2})\mbox{\boldmath$P$}_{11}^{\prime\prime}\right]$

In (425), the remaining coefficients are given by Blennerhassett [2].

This system of equations has a unique solution provided the inhomogeneous
terms are orthogonal to the adjoint function. Applying the solvability
condition yields the amplitude equation (1).

The complex coefficients $\kappa,a_{2}$, $d_{1}$ and $\varpi$ in equation (1)
are given by

$\displaystyle\kappa=\displaystyle\frac{u_{b}^{\prime}({\psi}^{\prime\prime}+k_{c}^{2}\psi)}{\sigma_{r}\nu
R_{c}}-\displaystyle\frac{{\psi}^{\prime\prime\prime}-3k_{c}^{2}{\psi}^{\prime}}{\nu
R_{c}}+ik_{c}\sigma_{r}{\psi}^{\prime}-\displaystyle\frac{ik_{c}\psi(\mathscr{F}+k_{c}^{2}\mathscr{T})}{\sigma_{r}R_{c}^{2}}$
$\displaystyle
a_{2}=-\displaystyle\frac{{\psi}^{\prime\prime}{+k_{c}^{2}\psi}}{\nu
R_{c}}(u_{b}^{\prime}-U_{b}^{\prime})(a_{1}+c_{r})\left(1-\displaystyle\frac{u_{b}+a_{1}}{k_{c}\sigma_{r}}\right)\displaystyle\frac{\phi}{k_{c}\sigma_{r}^{2}}+\displaystyle\frac{{\psi}^{\prime}}{\nu
R_{c}}(1-\mu)(2k_{c}+1)\phi$ $\displaystyle\hskip
8.5359pt+i\phi\displaystyle\frac{k_{c}\mathscr{T}(3u_{b}-5c_{r}-2a_{1})}{\sigma_{r}^{2}R_{c}^{2}}$
$\displaystyle\hskip 8.5359pt-\psi\left\\{i\left[\phi
u_{b}^{\prime}-{\phi}^{\prime}(u_{b}+a_{1})\right]-i\rho\left[\Phi
U_{b}^{\prime}-{\Phi}^{\prime}(U_{b}+a_{1})\right]-\displaystyle\frac{6k_{c}}{\nu
R_{c}}({\phi}^{\prime}-\mu{\Phi}^{\prime})-\displaystyle\frac{3}{\nu
R_{c}}({\phi}^{\prime}-\mu{\Phi}^{\prime})\right.$ $\displaystyle\hskip
8.5359pt\left.+i\phi\left[\displaystyle\frac{\mathscr{F}(u_{b}-2c_{r}-a_{1})(1-\rho)+k_{c}^{2}\mathscr{T}(3u_{b}-4c_{r}-a_{1})}{\sigma_{r}^{2}R_{c}^{2}}+(a_{1}+c_{r})^{2}\displaystyle\frac{\mathscr{F}(1-\rho)+k_{c}^{2}\mathscr{T}}{\sigma_{r}^{3}k_{c}R_{c}^{2}}\right]\right\\}$
$\displaystyle a_{1}=\displaystyle\frac{ic_{r}\phi}{\nu
k_{c}R_{c}\sigma_{r}^{2}}(u_{b}^{\prime}-U_{b}^{\prime})(\psi^{\prime\prime}+k_{c}^{2}\psi)-\displaystyle\frac{2ik_{c}\phi\psi^{\prime}(1-\mu)}{\nu
R_{c}}$ $\displaystyle\hskip 8.5359pt-\psi\left[\phi
u_{b}^{\prime}-\phi^{\prime}u_{b}-\rho(\Phi
U_{b}^{\prime}-\Phi^{\prime}U_{b})+\displaystyle\frac{6ik_{c}}{\nu
R_{c}}(\phi^{\prime}-\mu\Phi^{\prime})\right.$ $\displaystyle\hskip
8.5359pt\left.+\phi\displaystyle\frac{{\mathscr{F}}(u_{b}-2c_{r})(1-\rho)+k_{c}^{2}{\mathscr{T}}(3u_{b}-4c_{r})}{\sigma_{r}^{2}R_{c}^{2}}\right]$
$\displaystyle
d_{1}=\displaystyle\frac{ik_{c}\psi}{R_{c}}\left[\displaystyle\frac{\mathscr{F}(1-\rho)+k_{c}^{2}\mathscr{T}}{\sigma_{r}R_{c}^{2}}\phi+(1-\rho)({\phi}^{\prime}\sigma_{r}-\phi
u_{b}^{\prime})\right]$ $\displaystyle\hskip
8.5359pt-\displaystyle\frac{1}{\nu
R_{c}^{2}}\displaystyle\int_{-d^{-1}}^{0}\psi(\phi^{iv}-2k_{c}^{2}{\phi}^{\prime\prime}+k_{c}^{4}\phi)\,d\eta$
$\displaystyle\hskip
8.5359pt-\displaystyle\frac{\rho}{R_{c}^{2}}\displaystyle\int_{0}^{1}\Psi(\Phi^{iv}-2k_{c}^{2}{\Phi}^{\prime\prime}+k_{c}^{4}\Phi)\,d\eta$
$\displaystyle\varpi=\displaystyle\frac{-i\phi}{\nu
k_{c}R_{c}\sigma_{r}^{2}}(u^{\prime}_{b}-U^{\prime}_{b})(\psi^{\prime\prime}+k_{c}\psi)+\psi\left[\rho\Phi^{\prime}-\phi^{\prime}-\phi\displaystyle\frac{{\mathscr{F}}(1-\rho)+k_{c}^{2}{\mathscr{T}}}{\sigma_{r}^{2}R_{c}^{2}}\right]$
$\displaystyle\hskip
8.5359pt+\displaystyle\int_{-d^{-1}}^{0}\psi(\phi^{\prime\prime}-k_{c}^{2}\phi)\,d\eta+\rho\displaystyle\int_{0}^{1}\Psi(\Phi^{\prime\prime}-k_{c}^{2}\Phi)\,d\eta$

where $\Phi,\phi$ satisfies the Orr-Sommerfeld equation for the basic flow
with $\Psi,\psi$ satisfying their adjoints.

It is important to note that in the formulation presented here quantities such
as $\Phi$ and $\phi$, and their adjoints, are functions of $\eta$ alone, but,
by virtue of (17) they are functions of both spacial coordinates $x$, $y$ and
time $t$. Also, since the Jacobian of the transformation is $1+O(ak)$, the
resulting equations expressed in curvilinear coordinates are, to first order
in $ak$, symbollically identical to their Caretesian counterpart. Thus, the
solutions for the basic flow, in curvilinear coordinates, are the same as
those given by (10) and (489), in Cartesian coordinates, with $y$ simply
replaced with $\eta$. The error in this approximation is $O(ak)^{2}$ and may
be neglected since $ak\ll 1$.

##  6 Numerical Schemes

The complex cubic–quintic nonlinear Schrödinger equation (1) is solved
numerically, using a second-order central difference for the spacial part and
the third-order explicit Runge-Kutta method for the transient part, whose
complex coefficients are determined from the system of ordinary differential
equations (254), (329) and (329) which coupled via interface conditions which
are solved numerically using the finite difference scheme as described below.
First, we rewrite the system in a compact form as

$\displaystyle\left.\begin{array}[]{lcl}\Phi=\Phi^{\prime}=0&&{\rm
on}\,\,\eta=1,\\\ \\\
\Phi^{iv}+C(\eta)\Phi^{\prime\prime}+D(\eta)\Phi+2E(\eta)\Phi^{\prime}=F(\eta)&&{\rm
in}\,\,0<\eta<1,\\\ \\\
\sum_{p=0}^{3}FL_{np}\phi^{(p)}+FLNH_{n}=\sum_{p=0}^{3}FU_{np}\Phi^{(p)}+FUNH_{n},&&\eta=0,\,n=1,2,3,4,\\\
\\\
\phi^{iv}+c(\eta)\phi^{\prime\prime}+d(\eta)\phi+2e(\eta)\phi^{\prime}=f(\eta)&&{\rm
in}\,\,-d^{-1}<\eta<0,\\\ \\\ \phi=\phi^{\prime}=0&&{\rm
on}\,\,\eta=-d^{-1},\\\ \end{array}\right\\}$ (465) (466)

where $\Phi=({\it\Phi}_{11},{\it\Phi}_{22},\Psi)$,
$\phi=(\phi_{11},\phi_{22},\psi)$ and the definitions of the coefficients $c$,
$C$, $d$, $D$, $e$, $E$, $FL_{np}$, $FU_{np}$, $FLNH_{n}$, $FUNH_{n}$, $f$,
$F$ are those given by system of equation (254), (329) and (395),
respectively.

We approximate

$\displaystyle\left.\begin{array}[]{lcccl}\Phi=\left(1+\frac{1}{6}\delta^{2}-\frac{1}{720}\delta^{4}\right)G&&&&\mbox{for
the upper fluid}\\\ \\\
\phi=\left(1+\frac{1}{6}\delta^{2}-\frac{1}{720}\delta^{4}\right)g&&&&\mbox{for
the lower fluid}\\\ \end{array}\right\\}$ (470)

where $\delta$ is the usual central difference operator, i.e., $\delta
g_{i}=g_{i+1/2}-g_{i-1/2}$. Then

$\displaystyle h^{4}\phi^{iv}$ $\displaystyle=$
$\displaystyle\delta^{4}g+O\left(\delta^{10}\right)$ (471) $\displaystyle
h^{2}\phi^{\prime\prime}$ $\displaystyle=$
$\displaystyle(\delta^{2}+\frac{1}{12}\delta^{4})g+O\left(\delta^{6}\right)$
(472) $\displaystyle h\phi^{\prime}$ $\displaystyle=$ $\displaystyle\mu\delta
g+O\left(\delta^{5}\right)$ (473)

Here $h=\eta_{i}-\eta_{i-1}$ and $\mu$ is defined as $\mu
g_{i}=\frac{1}{2}\left(g_{i+1/2}+g_{i-1/2}\right)$.

For $\phi^{\prime\prime\prime}$ we use the following approximation

$\displaystyle
h^{3}\phi^{\prime\prime\prime}=\mu\delta^{3}\left(1+\frac{1}{6}\delta^{2}-\frac{1}{720}\delta^{4}\right)^{-1}\phi-\frac{1}{12}\mu\delta\left(h^{4}\phi^{(iv)}\right)+O\left(\delta^{7}\right).$
(474)

Using (474) and the appropriate equation of motion from (466), we have

$\displaystyle
h^{3}\phi^{\prime\prime\prime}=\mu\delta^{3}g+\frac{1}{12}h^{2}\mu\delta\left(c\delta^{2}g+h^{2}dg+2he\mu\delta
g-h^{2}f\right)+O\left(\delta^{7}\right).$ (475)

Hence we maintain $O\left(h^{4}\right)$ accuracy and a five point formula for
the third derivative in the interface. The similar expressions result for the
derivatives of $\Phi$. Substituting (471)–(473) in (466), we obtain the
following finite difference approximation.

$\displaystyle{\mathscr{A}}_{1}g_{i-2}+{\mathscr{A}}_{2}g_{i-1}+{\mathscr{A}}_{3}g_{i}+{\mathscr{A}}_{4}g_{i+1}+{\mathscr{A}}_{5}g_{i+2}=h^{4}f_{i},$
(476)

where

$\displaystyle{\mathscr{A}}_{1}={\mathscr{A}}_{5}=1+\frac{1}{12}h^{2}c_{i}$
(477)
$\displaystyle{\mathscr{A}}_{2}=-4+\frac{2}{3}h^{2}c_{i}+\frac{1}{6}h^{4}d_{i}-h^{3}e_{i}$
(478)
$\displaystyle{\mathscr{A}}_{3}=6-\frac{3}{2}h^{2}c_{i}+\frac{2}{3}h^{4}d_{i}$
(479) $\displaystyle{\mathscr{A}}_{4}={\mathscr{A}}_{2}+2h^{3}e_{i}$ (480)

where the subscript $i$ implies evaluation at $\eta=\eta_{i}$. Similar
expressions result for the approximation to $\Phi$. Thus, the finite
difference approximation to the coupled ODEs system with the boundary
conditions and the interface condition (466) leads to the matrix equation
${\bf{\mathscr{A}}x}={\bf b}$ with the following form:

$\left(\begin{array}[]{lcccccccccccccccccccccr}\ *&*&*&&&&&&&&&&&&&&&&&&&&\\\
\ *&*&*&*&&&&&&&&&&&&&&&&&&&\\\ \ *&*&*&*&*&&&&&&&&&&&&&&&&&&\\\
&*&*&*&*&*&&&&&&&&&&&&&&&&&\\\ &&&&\vdots&\vdots&&&&&&&&&&&&&&&\\\
&&&&&*&*&*&*&*&&&&&&&&&&&&&\\\ &&&&&&*&*&*&*&*&&&&&&&&&&&&\\\
&&&&&&&*&*&*&*&*&&&&&&&&&&&\\\ &&&&&&&*&*&*&*&*&*&*&*&*&*&&&&&&\\\
&&&&&&&*&*&*&*&*&*&*&*&*&*&&&&&&\\\ &&&&&&&*&*&*&*&*&*&*&*&*&*&&&&&&\\\
&&&&&&&*&*&*&*&*&*&*&*&*&*&&&&&&\\\ &&&&&&&&&&&&*&*&*&*&*&&&&&&\\\
&&&&&&&&&&&&&*&*&*&*&*&&&&&\\\ &&&&&&&&&&&&&&*&*&*&*&*&&&&\\\
&&&&&&&&&&&&&&&&&\vdots&\vdots&&\\\ &&&&&&&&&&&&&&&&&*&*&*&*&*&\\\
&&&&&&&&&&&&&&&&&&*&*&*&*&*\\\ &&&&&&&&&&&&&&&&&&&*&*&*&*\\\
&&&&&&&&&&&&&&&&&&&&*&*&*\\\
\end{array}\right)\left(\begin{array}[]{c}G_{NU}\\\ \cdot\\\ \cdot\\\ \cdot\\\
\cdot\\\ G_{2}\\\ G_{1}\\\ G_{0}\\\ G_{-1}\\\ G_{-2}\\\ G_{-3}\\\ g_{-3}\\\
g_{-2}\\\ g_{-1}\\\ g_{0}\\\ g_{1}\\\ g_{2}\\\ \cdot\\\ \cdot\\\ \cdot\\\
\cdot\\\ g_{NL}\end{array}\right)=\left(\begin{array}[]{c}F_{NU}\\\ \cdot\\\
\cdot\\\ \cdot\\\ \cdot\\\ F_{2}\\\ F_{1}\\\ F_{0}\\\ F_{-1}\\\ F_{-2}\\\
F_{-3}\\\ f_{-3}\\\ f_{-2}\\\ f_{-1}\\\ f_{0}\\\ f_{1}\\\ f_{2}\\\ \cdot\\\
\cdot\\\ \cdot\\\ \cdot\\\ f_{NL}\end{array}\right)$

where the entries of the matrix ${\bf\mathscr{A}}$ and unknown vector ${\bf
x}$ are complex numbers.

###  6.1 Complex cubic–quintic nonlinear Schrödinger equation

Let $h=(b-a)/n$ be the grid spacing and $x_{j}=a+jh,\,j=0,1,\cdots,n$ be the
grid points. Define $A_{j}(t)$ as an approximation to
$A(x_{j},t),\,j=0,1,\cdots,n$, ${\bf A}(t)=(A_{1}(t),\cdots,A_{n-1}(t))^{\sf
T}$ and ${\bf A}_{0}(x)=(A_{0}(x_{1}),\cdots,A_{0}(x_{n-1}))^{\sf T}$. Let
$F(A_{j}(t))=\displaystyle\frac{d_{1}}{d_{1r}}A_{j}(t)+\kappa
A_{j}(t)|A_{j}(t)|^{2}+\varpi A_{j}(t)|A_{j}(t)|^{4}$. Using the central
difference approximation to $\partial_{x}^{2}A(x_{j},t)$, we write the
semi–discretization of the IBVP for equation (1) as the following system of
ODEs

$\displaystyle\dot{\bf A}(t)={\mathscr{B}}{\bf A}(t)+{\bf F}({\bf
A}(t)),\,\forall t>0,$ (481) $\displaystyle{\bf A}(0)={\bf A}_{0},$ (482)

where ${\bf F}({\bf
A}(t))=\left(F(A_{1}(t)),\,\cdots,\,F(A_{n}(t))\right)^{\sf T}$. The upper dot
in (481) indicates derivative with respect to $t$ and ${\mathscr{B}}$ is the
finite difference matrix.

To construct an integration scheme to solve the ODE system (481)–(482), let
$t_{n+1}=t_{n}+\Delta t$, and let ${\bf A}^{n}$ denote the value of the
variable ${\bf A}$ at time $t_{n}$. Employing a low storage variant
third–order Runge–Kutta scheme [12], we write the fully discrete system as

$\begin{array}[]{lll}&{\bf Q}_{1}=\Delta t{\bf G}({\bf A}^{n}),&{\bf
A}_{1}={\bf A}^{n}+\frac{1}{3}{\bf Q}_{1},\\\ \\\ &{\bf Q}_{2}=\Delta t{\bf
G}({\bf A}_{1})-\frac{5}{9}{\bf Q}_{1},&{\bf A}_{2}={\bf
A}_{1}+\frac{15}{16}{\bf Q}_{2},\\\ \\\ &{\bf Q}_{3}=\Delta t{\bf G}({\bf
A}_{2})-\frac{153}{128}{\bf Q}_{2},&{\bf A}^{n+1}={\bf A}_{2}+\frac{8}{15}{\bf
Q}_{3},\\\ \end{array}$ (483)

where ${\bf G}({\bf A}^{n})={\mathscr{B}}{\bf A}^{n}+{\bf F}({\bf A}^{n})$.

The numerical scheme is parallelized for distributed memory clusters of
processors or heterogeneous networked computers using the MPI (Message Passing
Interface) library and implemented in FORTRAN.

##  7 Results

For the numerical simulation reported in this paper, we choose
${\mathscr{F}}=10^{3}k_{s}^{3}$ and ${\mathscr{T}}=10^{3}(1-\rho)k_{s}$ for
the present case of wind blowing over waves. Here the quantity $k_{s}$ is
defined as

$\displaystyle k_{s}=D_{u}\sqrt{\displaystyle\frac{(\rho_{w}-\rho_{a})g}{S}}$

with $\rho_{a}=1.225\times 10^{-3}\,{\rm g/cm^{3}}$, $\rho_{w}=0.9991\,{\rm
g/cm^{3}}$, $g=981\,{\rm cm/s^{2}}$, $S=73.5\,{\rm dyn/cm}$, and we choose
$D_{\ell}=D_{u}$ such that $d=1$. The kinematic viscosities are taken to be
$\nu_{a}=0.145\,{\rm cm^{2}/s}$ and $\nu_{w}=1.138\times 10^{-2}\,{\rm
cm^{2}/s}$.

The velocity profiles for the basic flow, for the case of turbulent boundary
layer, are those given by (8). The velocity profiles for pPf and pCf, for the
upper fluid, satisfy

$\displaystyle\left.\begin{array}[]{lll}U_{b}=U_{1},&&\mbox{on}\hskip
7.11317pty=D_{u}\\\
\displaystyle\frac{G}{\rho_{a}}+\nu_{w}\displaystyle\frac{{\rm
d}^{2}U_{b}}{{\rm d}y^{2}}=0,&&\mbox{in upper fluid}\\\ u_{b}=U_{b},\hskip
7.11317pt\mu_{w}\displaystyle\frac{du_{b}}{dy}=\mu_{a}\displaystyle\frac{dU_{b}}{dy},&&\mbox{on}\hskip
7.11317pty=0\\\ \displaystyle\frac{G}{\rho_{w}}+\nu_{w}\displaystyle\frac{{\rm
d}^{2}u_{b}}{{\rm d}y^{2}}=0,&&\mbox{in lower fluid}\\\
u_{b}=0,&&\mbox{on}\hskip 7.11317pty=-D_{\ell}\\\ \end{array}\right\\}$ (489)

The solution of this system is given by [2]

$\displaystyle U_{b}(y)=-{\sf
A}(y-1)\left[y+\displaystyle\frac{\mu(1+d)}{d(\mu+d)}\right]+{\sf B}(yd+\mu)$
$\displaystyle u_{b}(u)=-\displaystyle\frac{\mu}{2}{\sf
A}(yd+1)\left(y-\displaystyle\frac{1+d}{\mu+d}\right)+\mu{\sf B}(yd+1)$

where

$\displaystyle{\sf
A}=U_{p}\left[\displaystyle\frac{\mu+d^{2}}{6d^{2}(1+d)}+\displaystyle\frac{1}{2}\displaystyle\frac{\mu(1+d)}{d(\mu+d)}\right]^{-1}$
$\displaystyle{\sf B}=(1-U_{p})\displaystyle\frac{2(1+d)}{\mu+d(d+2\mu)}$

and

$\displaystyle
U_{p}=\displaystyle\frac{2(1+d)d}{\mu+d(d+2\mu)}\left[\displaystyle\frac{\mu+d^{2}}{6d^{2}(1+d)}+\displaystyle\frac{\mu}{2d}\displaystyle\frac{1+d}{\mu+d}\right]\left\\{1+\displaystyle\frac{\mu+d^{2}}{3d[\mu+d(d+2\mu)]}\right.$
$\displaystyle\left.+\displaystyle\frac{\mu(1+d)}{d(\mu+d)}\left[1+\displaystyle\frac{d(1+d)}{\mu+d(d+2\mu)}\right]\right\\}^{-1}$

Note that, for pPf $U_{p}=1$ which implies that he upper boundary is
stationary, and the flow is driven by the pressure gradient. In contrast, for
pCf $U_{p}=0$ and the motion is generated soley by the relative motion of the
boundaries.

Computations were performed on a Linux cluster (zeus.db.erau.edu: 256 Intel
Xeon 3.2GHz 1024 KB cache 4GB with Myrinet MX, GNU Linux) at Embry-Riddle
University. In all simulations presented here, we use the spacial domain of
$0\leq\xi\leq 3$ and temporal domain $0\leq\tau\leq 15$, the typical values of
the parameters $a=0.3$, $l=0.5$ for the amplitude and the wavelength of the
Stokes wave profile, respectively, and the other parameters of the model (1)
are given in Table 1 below.

Profile | $d_{1}$ | $a_{2}$ | $\kappa$ | $\varpi$
---|---|---|---|---
pCf | $(0.895\times 10^{-5},0.217\times 10^{-2})$ | $(0.895\times 10^{-5},0.217\times 10^{-2})$ | $(-29.7,237)$ | $(-5.94,47.4)$
pPf | $(0.404\times 10^{-5},0.742\times 10^{-3})$ | $(0.557\times 10^{-3},-0.17\times 10^{-1})$ | $(-9.95,912.3)$ | $(-2.49,228)$
LBL | $(0.554\times 10^{-5},0.275\times 10^{-2})$ | $(0.165\times 10^{-2},-0.159)$ | $(-0.221,20.61)$ | $(-0.0737,6.87)$
TBL | $(0.222\times 10^{-2},0.376\times 10^{-3})$ | $(0.042\times 10^{-3},0.11\times 10^{-1})$ | $(-19.2,1142)$ | $(-9.6,571)$
Table 1: Parameters for various models

Figure 2: Left: the surface profile for pCf in the air; Right: the surface
profile for pPf in the air.

Figure 3: Left: the surface profile for LBL in the air, Right: the surface
profile for TBL in the air.

Figure 4: Wind generated solitons in shallow water.

Figure 2 shows the surface structure for pCf (left) and pPf (right). For the
case of pPf we observe the classical ‘snake’ type solitons [9], whereas for
pCf we see a different structure similar to colliding solitons which, in time,
bifurcates to weaker soliton structure.

In figure 3 the surface structure for LBL (left) and TBL (right) are depicted.
As can be seen, the soliton structures for TBL resembles the chaotic solitons
[9]. This is, of course, not surprising as turbulence generates more complex
surface agitation compared to its laminar counterpart. We may term this as a
violent surface motion.

The most intriguing results obtained is that shown in figure 3 for LBL. This
results clearly confirms the generation of solitons by wind, which has
commonly been observed is shallow waters, c.f. figure 4.

##  8 Conclusions

In this paper we presented an asymptotic theory for linear stability of non-
linear air-sea interfaces by taking into account the motion in the air and the
water. The two motions are then matched at the free surface via the usual
kinematic, dynamic and free stress boundary conditions. The flow above water
is assumed to be turbulent, with a mean logarithmic velocity profile. On the
other hand, a plane Poiseuille flow is assumed for the motion in the water.
Also presented, is the linear stability calculations of two fluid interface
which suggest the flow become unstable for wavenumbers above that given by
equation (23). Finally, the growth of surface waves by the shear wind blowing
over is calculated. It is shown and extra term arises in the expression for
the dimensionless energy transfer parameter which is due to the coupling
between wind and surface waves.

In the second part of this paper [3], we use a non-linear hydrodynamic
stability [7] to derive the quintic non-linear Schrödinger equation (1) with
complex coefficients. We then use the above analysis to determine the complex
amplitude of the equation (1).

Using a nonlinear stability analysis of a coupled air-sea system is performed.
An amplitude equation, namely the quintic nonlinear Schrödinger equation with
complex coefficients is derived. The complex coefficients are obtained from
the linear stability of the same system in the neighborhood of the critical
Reynolds number. Four different profiles for the shear flow for the airflow
over the interface is considered, namely, plane Couette flow (pCf), plane
Poiseuille flow (pPf), laminar and turbulent boundary layer (L,TBL) profiles.
For each of the above cases the shear flow counterpart in the water is assumed
to be a pPf.

It is shown that the above amplitude equation produces ‘snake’ solitons for
pCf, pPf and LBL profiles, with striking similarities. On the other hand, for
TBL we observe a very violent surface motion. For cases of pCf and LBL
remarkable similarity is observed with observations made at sea.

We conclude that the effect of nonlinearity in the airflow over the sea
surface is much larger than nonlinear interactions in the water, and hence it
is not possible to decouple the motion in the air and the water for finite
amplitude wind-wave interactions. We remark, under the conditions shown in
this paper, the generation of solitons by wind is possible particularly in
shallow waters.

## References

  * (1)
  * (2) [1] S.E. Belcher and J.C.R. Hunt, Turbulent shear flow over slowly moving waves. J. Fluid Mech., 251, 109–148. 1993.
  * (3)
  * (4) [2] P.J. Blennerhassett, On the generation of waves by wind. Phil. Trans. Roy. Soc., A 298, 451-494, 1980.
  * (5)
  * (6) [3] S. Feldman, On the hydrodynamic stability of two viscous incompressible fluids in parallel uniform shearing motion. J. Fluid Mech., 2, 343-370, 1957.
  * (7)
  * (8) [4] H. Khanal and S.C. Mancas, Numerical simulations of five novel classes of dissipative solitons. CMMSE Proceedings, Vol II, 354-365, 2008.
  * (9)
  * (10) [5] H. Khanal and S.G. Sajjadi, Wind-generated solitons Part 2. Numerical simulations. CMMSE Proceedings, Vol II, 366-367, 2008.
  * (11)
  * (12) [6] B.M. Lake and H.C. Yuen, A new model for nonlinear wind waves. Part 1. Physical model and experimental evidence. J. Fluid Mech., 88, 33-62, 1978.
  * (13)
  * (14) [7] R.C. Lock, Hydodynamic stability of the flow in the laminar boundary layer between parallel streams, Proc. Camb. Phil. Soc., 50, 105-124, 1954
  * (15)
  * (16) [8] J.W. Miles, On the generation of surface waves by turbulent shear flows. J. Fluid Mech., 3, 185–204, 1957.
  * (17)
  * (18) [9] S.G. Sajjadi, Interaction of turbulence due to tropical cyclones with surface waves. Adv. Appl. Fluid Mech., 1, 101-145, 2007.
  * (19)
  * (20) [10] S.G. Sajjadi, On the growth of a fully non-linear Stokes wave by turbulent shear flow. Part 2. Rapud distortion theory. Math. Engng. Ind., 6, 247-260, 1998\.
  * (21)
  * (22) [11] S.G. Sajjadi, J. Wakefield, and A.J. Croft, On the growth of a fully non-linear Stokes wave by turbulent shear flow. Part 1. Eddy-viscosity model. Math. Engng. Ind., 6, 185-207, 1997.
  * (23)
  * (24) [12] K. Stewartson and J.T. Stuart, A non-linear instability theory for a wave system in plane Poiseuille flow. J. Fluid Mech., 48, 529-545, 1971.
  * (25)
  * (26) [13] J.T. Stuart and R.C. DiPrima, The Eckhaus and Benjamin-Feir resonance mechanisms. Proc. Roy. Soc., A 362, 27-41, 1978.
  * (27)
  * (28) [14] J.T. Stuart, On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. Part 1. The basic behaviour in plane Poiseuille flow. J. Fluid Mech., 9, 353-370, 1960.
  * (29)
  * (30) [15] G.R. Valenzuela, The growth of gravity-capillary waves in a coupled shear flow. J. Fluid Mech., 276, 229-250, 1976.
  * (31)
  * (32) [16] J. Williamson, Low-storage Runge-Kutta schemes, J. Comp. Phys., 35, 48-56, 1980.
  * (33)
  * (34) [17] C.S. Yih, 1967 Instability due to viscosity stratification. J. Fluid Mech., 27, 337-352. 1967.
  * (35)