A mathematical analysis of the Kakinuma model
for interfacial gravity waves.
Part I: Structures and well-posedness

Vincent Duchêne and Tatsuo Iguchi

###### Abstract

We consider a model, which we named the Kakinuma model, for interfacial
gravity waves. As is well-known, the full model for interfacial gravity waves
has a variational structure whose Lagrangian is an extension of Luke’s
Lagrangian for surface gravity waves, that is, water waves. The Kakinuma model
is a system of Euler–Lagrange equations for approximate Lagrangians, which are
obtained by approximating the velocity potentials in the Lagrangian for the
full model. In this paper, we first analyze the linear dispersion relation for
the Kakinuma model and show that the dispersion curves highly fit that of the
full model in the shallow water regime. We then analyze the linearized
equations around constant states and derive a stability condition, which is
satisfied for small initial data when the denser water is below the lighter
water. We show that the initial value problem is in fact well-posed locally in
time in Sobolev spaces under the stability condition, the non-cavitation
assumption and intrinsic compatibility conditions in spite of the fact that
the initial value problem for the full model does not have any stability
domain so that its initial value problem is ill-posed in Sobolev spaces.
Moreover, it is shown that the Kakinuma model enjoys a Hamiltonian structure
and has conservative quantities: mass, total energy, and in the case of the
flat bottom, momentum.

##  1 Introduction

We are concerned with the motion of interfacial gravity waves at the interface
between two layers of immiscible waters in a domain of the $(n+1)$-dimensional
Euclidean space in the rigid-lid case. Let $t$ be the time,
$\bm{x}=(x_{1},\ldots,x_{n})$ the horizontal spatial coordinates, and $z$ the
vertical spatial coordinate. We assume that the interface, the rigid-lid of
the upper layer, and the bottom of the lower layer are represented as
$z=\zeta(\bm{x},t)$, $z=h_{1}$, and $z=-h_{2}+b(\bm{x})$, respectively, where
$\zeta(\bm{x},t)$ is the elevation of the interface, $h_{1}$ and $h_{2}$ are
mean thicknesses of the upper and lower layers, and $b(\bm{x})$ represents the
bottom topography. The only external force applied to the system is the
constant and vertical gravity, and interfacial tension is neglected. Moreover,
we assume that the waters in the upper and the lower layers are both
incompressible and inviscid fluids with constant densities $\rho_{1}$ and
$\rho_{2}$, respectively, and that the flows are both irrotational. See Figure
1.1.

$\bm{x}$$z$$\zeta(\bm{x},t)$$h_{1}$$h_{2}$$\rho_{1}$$\rho_{2}$

Figure 1.1: Interfacial gravity waves

Then, the motion of the waters is described by the velocity potentials
$\Phi_{1}$ and $\Phi_{2}$ and the pressures $P_{1}$ and $P_{2}$ in the upper
and the lower layers, respectively, satisfying the basic equations in the
theory of fluid dynamics, which will be referred as the full model for
interfacial gravity waves throughout in this paper. As was shown by J. C. Luke
[20], the basic equations for the surface gravity waves, that is, the water
wave problem has a variational structure, whose Lagrangian is written in terms
of the surface elevation of the water and the velocity potential, and the
Lagrangian density is given by the vertical integral of the pressure in the
water region. The full model for interfacial gravity waves has also a
variational structure and the Lagrangian density
$\mathscr{L}(\Phi_{1},\Phi_{2},\zeta)$ is again given by the vertical integral
of the pressure in both water regions. T. Kakinuma [14, 15, 16] proposed a
model for interfacial gravity waves and applied his model to simulate
numerically the waves. To derive the model, he approximated the velocity
potentials $\Phi_{1}$ and $\Phi_{2}$ by

(1.1) $\Phi_{k}^{\rm
app}(\bm{x},z,t)=\sum_{i=0}^{N}Z_{k,i}(z;\tilde{h}_{k}(\bm{x}))\phi_{k,i}(\bm{x},t)$

for $k=1,2$, where $\\{Z_{1,i}\\}$ and $\\{Z_{2,i}\\}$ are appropriate
function systems in the vertical coordinate $z$ and may depend on
$\tilde{h}_{1}(\bm{x})$ and $\tilde{h}_{2}(\bm{x})$, respectively, which are
thickness of the upper and the lower layers in the rest state, whereas
$\bm{\phi}_{k}=(\phi_{k,0},\phi_{k,1},\ldots,\phi_{k,N})^{\rm T}$, $k=1,2$,
are unknown variables. Then, he derived an approximate Lagrangian density
$\mathscr{L}^{\rm
app}(\bm{\phi}_{1},\bm{\phi}_{2},\zeta)=\mathscr{L}(\Phi_{1}^{\rm
app},\Phi_{2}^{\rm app},\zeta)$ for unknowns
$(\bm{\phi}_{1},\bm{\phi}_{2},\zeta)$. The Kakinuma model is a corresponding
system of Euler–Lagrange equations for the approximated Lagrangian density
$\mathscr{L}^{\rm app}(\bm{\phi}_{1},\bm{\phi}_{2},\zeta)$. Different choices
of the function systems $\\{Z_{1,i}\\}$ and $\\{Z_{2,i}\\}$ give different
Kakinuma models and we have to carefully choose the function systems for the
Kakinuma model to provide good approximations for interfacial gravity waves.

The Kakinuma model is an extension to interfacial gravity waves of the so-
called Isobe–Kakinuma model for the surface gravity waves, that is, the water
waves. In the case of the surface gravity waves, the basic equations are known
to have a variational structure with Luke’s Lagrangian density
$\mathscr{L}_{\rm Luke}(\Phi,\zeta)$, where $\zeta$ is the surface elevation
and $\Phi$ is the velocity potential of the water. The Isobe–Kakinuma model is
a system of Euler–Lagrange equations for the approximated Lagrangian density
$\mathscr{L}^{\rm app}(\bm{\phi},\zeta)=\mathscr{L}_{\rm Luke}(\Phi^{\rm
app},\zeta)$, where $\Phi^{\rm app}$ is an approximate velocity potential

(1.2) $\Phi^{\rm
app}(\bm{x},z,t)=\sum_{i=0}^{N}Z_{i}(z;b(\bm{x}))\phi_{i}(\bm{x},t)$

and $\bm{\phi}=(\phi_{0},\phi_{1},\ldots,\phi_{N})^{\rm T}$ are unknown
variables. The model was first proposed by M. Isobe [12, 13] and then applied
by T. Kakinuma to simulate numerically the water waves. We note that a similar
model was derived by G. Klopman, B. van Groesen, and M. W. Dingemans [18], and
used to simulate the water waves. See also Ch. E. Papoutsellis and G. A.
Athanassoulis [25]. Recently, this model was analyzed from mathematical point
of view. One possible choice of the function system $\\{Z_{i}\\}$ is a set of
polynomials in $z$, for example, $Z_{i}(z;b(\bm{x}))=(z+h-b(\bm{x}))^{p_{i}}$
with integers $p_{i}$ satisfying $0=p_{0}<p_{1}<\cdots<p_{N}$. Under this
choice of the function system $\\{Z_{i}\\}$, the initial value problem to the
Isobe–Kakinuma model was analyzed by Y. Murakami and T. Iguchi [22] in a
special case and by R. Nemoto and T. Iguchi [23] in the general case. The
hypersurface $t=0$ in the space-time $\mathbf{R}^{n}\times\mathbf{R}$ is
characteristic for the Isobe–Kakinuma model, so that one needs to impose some
compatibility conditions on the initial data for the existence of the
solution. Under these compatibility conditions and a sign condition
$-\partial_{z}P^{\rm app}\geq c_{0}>0$ on the water surface, they showed the
well-posedness of the initial value problem locally in time, where $P^{\rm
app}$ is an approximate pressure in the Isobe–Kakinuma model calculated from
Bernoulli’s equation. Moreover, T. Iguchi [10, 11] showed that under the
choice of the function system

(1.3) $Z_{i}(z;b(\bm{x}))=\begin{cases}(z+h)^{2i}&\mbox{in the case of the
flat bottom},\\\ (z+h-b(\bm{x}))^{i}&\mbox{in the case of a variable
bottom},\end{cases}$

the Isobe–Kakinuma model is a higher order shallow water approximation for the
water wave problem in the strongly nonlinear regime. Furthermore, V. Duchêne
and T. Iguchi [7] showed that the Isobe–Kakinuma model also enjoys a
Hamiltonian structure analogous to the one exhibited by V. E. Zakharov [26] on
the full water wave problem. Our aim in the present paper is to extend these
results on the surface gravity waves to interfacial gravity waves.

In view of these results on the Isobe–Kakinuma model, in the present paper we
consider the Kakinuma model under the choice of the approximate velocity
potentials in (1.1) as

(1.4) $\begin{cases}\displaystyle\Phi_{1}^{\rm
app}(\bm{x},z,t)=\sum_{i=0}^{N}(-z+h_{1})^{2i}\phi_{1,i}(\bm{x},t),\\\\[10.76385pt]
\displaystyle\Phi_{2}^{\rm
app}(\bm{x},z,t)=\sum_{i=0}^{N^{*}}(z+h_{2}-b(\bm{x}))^{p_{i}}\phi_{2,i}(\bm{x},t),\end{cases}$

where $N,N^{*}$ and $p_{0},p_{1},\ldots,p_{N^{*}}$ are nonnegative integers
satisfying $0=p_{0}<p_{1}<\cdots<p_{N^{*}}$. In applications of the Kakinuma
model, it would be better to choose $N^{*}=N$ and $p_{i}=2i$ in the case of
the flat bottom, and $N^{*}=2N$ and $p_{i}=i$ in the case of a variable
bottom. In the case $N=N^{*}=0$, that is, if we choose the approximation
$\Phi_{k}^{\rm app}(\bm{x},z,t)=\phi_{k}(\bm{x},t)$ for $k=1,2$ the functions
independent of the vertical coordinate $z$, then the corresponding Kakinuma
model is reduced to the shallow water equations. In the case $N+N^{*}>0$, the
Kakinuma model is classified into a system of nonlinear dispersive equations.

It is well-known that in the case of the flat bottom $b=0$, the dispersion
relation of the linearized equations to the full model around the flow
$(\zeta,\Phi_{1},\Phi_{2})=(0,\bm{u}_{1}\cdot\bm{x},\bm{u}_{2}\cdot\bm{x})$
with constant horizontal velocities $\bm{u}_{1}$ and $\bm{u}_{2}$ is given by

$\displaystyle(\rho_{1}\coth(h_{1}|\bm{\xi}|)+\rho_{2}\coth(h_{2}|\bm{\xi}|))\omega^{2}$
$\displaystyle+2(\rho_{1}\bm{\xi}\cdot\bm{u}_{1}\coth(h_{1}|\bm{\xi}|)+\rho_{2}\bm{\xi}\cdot\bm{u}_{2}\coth(h_{2}|\bm{\xi}|))\omega$
$\displaystyle+\rho_{1}(\bm{\xi}\cdot\bm{u}_{1})^{2}\coth(h_{1}|\bm{\xi}|)+\rho_{2}(\bm{\xi}\cdot\bm{u}_{2})^{2}\coth(h_{2}|\bm{\xi}|)-(\rho_{2}-\rho_{1})g|\bm{\xi}|=0,$

where $\bm{\xi}\in\mathbf{R}^{n}$ is the wave vector, $\omega\in\mathbf{C}$
the angular frequency, and $g$ the gravitational constant. It is easy to see
that the roots $\omega$ of the above equation are always real for any wave
vector $\bm{\xi}\in\mathbf{R}^{n}$ if and only if $\bm{u}_{1}=\bm{u}_{2}$ and
$\rho_{2}\geq\rho_{1}$. Otherwise, the roots of the above equation have the
form $\omega=\omega_{r}(|\bm{\xi}|)\pm\mathrm{i}\omega_{i}(|\bm{\xi}|)$
satisfying $\omega_{i}(|\bm{\xi}|)\to+\infty$ as $|\bm{\xi}|\to+\infty$, which
leads to an instability of the interface. The instabilities in the case
$\rho_{2}>\rho_{1}$ and $\bm{u}_{1}\neq\bm{u}_{2}$ and in the case
$\rho_{2}<\rho_{1}$ and $\bm{u}_{1}=\bm{u}_{2}$ are known as the
Kelvin–Helmholtz and the Rayleigh–Taylor instabilities, respectively. For more
details, see for example P. G. Drazin and W. H. Reid [6]. In the following of
this paper, we are interested in the situation where

$(\rho_{2}-\rho_{1})g>0,$

that is, the denser water is below the lighter water. In the case
$\bm{u}_{1}=\bm{u}_{2}=\bm{0}$, the linear dispersion relation is written
simply as

$\omega^{2}=\frac{(\rho_{2}-\rho_{1})g|\bm{\xi}|}{\rho_{1}\coth(h_{1}|\bm{\xi}|)+\rho_{2}\coth(h_{2}|\bm{\xi}|)}.$

We denote the right-hand side by $\omega_{\mbox{\rm\tiny IW}}(\bm{\xi})^{2}$.
Then, the phase speed $c_{\mbox{\rm\tiny IW}}(\bm{\xi})$ of the plane wave
solution related to the wave vector $\bm{\xi}$ is given by

(1.5) $c_{\mbox{\rm\tiny IW}}(\bm{\xi})=\frac{\omega_{\mbox{\rm\tiny
IW}}(\bm{\xi})}{|\bm{\xi}|}=\pm\sqrt{\frac{(\rho_{2}-\rho_{1})g}{\rho_{1}|\bm{\xi}|\coth(h_{1}|\bm{\xi}|)+\rho_{2}|\bm{\xi}|\coth(h_{2}|\bm{\xi}|)}}.$

As a shallow water limit $h_{1}|\bm{\xi}|,h_{2}|\bm{\xi}|\to 0$, we have

(1.6) $c_{\mbox{\rm\tiny IW}}(\bm{\xi})\simeq c_{\mbox{\rm\tiny
SW}}=\pm\sqrt{\frac{(\rho_{2}-\rho_{1})gh_{1}h_{2}}{\rho_{1}h_{2}+\rho_{2}h_{1}}},$

where $c_{\mbox{\rm\tiny SW}}$ is the phase speed of infinitely long and small
interfacial gravity waves. In Section 3, we will analyze the linear dispersion
relation of the Kakinuma model and calculate the phase speed
$c_{\mbox{\rm\tiny K}}(\bm{\xi})$ of the plane wave solution related to the
wave vector $\bm{\xi}$. Under the choice $N^{*}=N$ and $p_{i}=2i$, or
$N^{*}=2N$ and $p_{i}=i$ in the approximation (1.4) of the velocity
potentials, it turns out that

(1.7) $|c_{\mbox{\rm\tiny IW}}(\bm{\xi})^{2}-c_{\mbox{\rm\tiny
K}}(\bm{\xi})^{2}|\lesssim(h_{1}|\bm{\xi}|+h_{2}|\bm{\xi}|)^{4N+2},$

which indicates that the Kakinuma model may be a good approximation of the
full model for interfacial gravity waves in the shallow water regime
$h_{1}|\bm{\xi}|,h_{2}|\bm{\xi}|\ll 1$. We note that Miyata–Choi–Camassa model
derived by M. Miyata [21] and W. Choi and R. Camassa [3] is a model for
interfacial gravity waves in the strongly nonlinear regime and can be regarded
as a generalization of the Green–Naghdi equations for water waves into a two-
layer system. Let $c_{\mbox{\rm\tiny MCC}}(\bm{\xi})$ be the phase speed of
the plane wave solution related to the wave vector $\bm{\xi}$ for the
linearized equations of the Miyata–Choi–Camassa model around the rest state.
Then, we have

$|c_{\mbox{\rm\tiny IW}}(\bm{\xi})^{2}-c_{\mbox{\rm\tiny
MCC}}(\bm{\xi})^{2}|\lesssim(h_{1}|\bm{\xi}|+h_{2}|\bm{\xi}|)^{4},$

so that the Kakinuma model gives a better approximation of the full model than
the Miyata–Choi–Camassa model in the shallow water regime, at least, at the
linear level. A rigorous analysis for the consistency of the Kakinuma model in
the shallow water regime will be analyzed in the subsequent paper V. Duchêne
and T. Iguchi [8]. On the other hand, in the deep water limit we have

$\lim_{h_{1}|\bm{\xi}|,h_{2}|\bm{\xi}|\to\infty}c_{\mbox{\rm\tiny
K}}(\bm{\xi})^{2}>0,$

which is not consistent with the limit of the full model

$\lim_{h_{1}|\bm{\xi}|,h_{2}|\bm{\xi}|\to\infty}c_{\mbox{\rm\tiny
IW}}(\bm{\xi})^{2}=0.$

We notice that the Miyata–Choi–Camassa model is only apparently consistent
with the full model in this deep water limit since

$\lim_{h_{1}|\bm{\xi}|,h_{2}|\bm{\xi}|\to\infty}c_{\mbox{\rm\tiny
MCC}}(\bm{\xi})^{2}=0,$

but we note also that

$\lim_{h_{1}|\bm{\xi}|,h_{2}|\bm{\xi}|\to\infty}\frac{c_{\mbox{\rm\tiny
IW}}(\bm{\xi})^{2}}{c_{\mbox{\rm\tiny MCC}}(\bm{\xi})^{2}}=\infty.$

We refer to V. Duchêne, S. Israwi, and R. Talhouk [9] for further discussion
and the derivation of modified Miyata–Choi–Camassa models having either the
same dispersion relation as the full model, or the same behavior as the
Kakinuma model in the deep water limit. As we discuss below, thanks to the
high-frequency behavior of the linearized equations, and contrarily to both
the full model and the Miyata–Choi–Camassa model, the Kakinuma model has a
non-trivial stability domain and, as a result, the initial value problem to
the Kakinuma model is well-posed locally in time in Sobolev spaces under
appropriate assumptions on the initial data.

As we have already seen, the roots $\omega\in\mathbf{C}$ of the dispersion
relation of the linearized equations of the full model around the rest state
are always real, so that the corresponding initial value problem is well-
posed. However, as for the nonlinear problem, even if the initial velocity is
continuous on the interface, a discontinuity of the velocity in the tangential
direction on the interface would be created instantaneously in general, so
that the Kelvin–Helmholtz instability appears locally in space. As a result,
the initial value problem for the full model turns out to be ill-posed. For
more details, we refer to T. Iguchi, N. Tanaka, and A. Tani [24]. See also V.
Kamotski and G. Lebeau [17] and D. Lannes [19]. In Section 4 we consider the
linearized equations of the Kakinuma model around an arbitrary flow. After
freezing the coefficients and neglecting lower order terms of the linearized
equations, we calculate the linear dispersion relation and derive a stability
condition, which is equivalent to

(1.8) $-\partial_{z}(P_{2}^{\rm app}-P_{1}^{\rm
app})-\frac{\rho_{1}\rho_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}|\nabla\Phi_{2}^{\rm
app}-\nabla\Phi_{1}^{\rm app}|^{2}\geq c_{0}>0$

on the interface, where $P_{1}^{\rm app}$ and $P_{2}^{\rm app}$ are
approximate pressures of the waters in the upper and the lower layers in the
Kakinuma model calculated from Bernoulli’s equations, $H_{1}$ and $H_{2}$ are
thickness of the upper and the lower layers, respectively, $\alpha_{1}$ is a
constant depending only on $N$, $\alpha_{2}$ is a constant determined from
$\\{p_{0},p_{1},\ldots,p_{N^{*}}\\}$, and
$\nabla=(\partial_{x_{1}},\ldots,\partial_{x_{n}})^{\rm T}$ is the nabla with
respect to the horizontal spatial coordinates $\bm{x}=(x_{1},\ldots,x_{n})$.
If $\rho_{1}=0$, then (1.8) coincides with the stability condition for the
Isobe–Kakinuma model for water waves derived by R. Nemoto and T. Iguchi [23].

As in the case of the Isobe–Kakinuma model, the hypersurface $t=0$ in the
space-time $\mathbf{R}^{n}\times\mathbf{R}$ is characteristic for the Kakinuma
model, so that one needs to impose some compatibility conditions on the
initial data for the existence of the solution. Under these compatibility
conditions, the non-cavitation assumption $H_{1}\geq c_{0}>0$and $H_{2}\geq
c_{0}>0$, and the stability condition (1.8), we will show in this paper that
the initial value problem to the Kakinuma model is well-posed locally in time
in Sobolev spaces. Here, we note that the coefficients $\alpha_{1}$ and
$\alpha_{2}$ in the stability condition (1.8) converge to $0$ as
$N,N^{*}\to\infty$, so that the domain of stability diminishes as $N$ and
$N^{*}$ grow. This fact is consistent with the aforementioned properties of
the full model.

As is well-known that the full model for interfacial gravity waves has a
conserved energy

(1.9) $\displaystyle\mathscr{E}$
$\displaystyle=\int\\!\\!\\!\int_{\Omega_{1}(t)}\frac{1}{2}\rho_{1}\bigl{(}|\nabla\Phi_{1}(\bm{x},z,t)|^{2}+(\partial_{z}\Phi_{1}(\bm{x},z,t))^{2}\bigr{)}{\rm
d}\bm{x}{\rm d}z$
$\displaystyle\quad\;+\int\\!\\!\\!\int_{\Omega_{2}(t)}\frac{1}{2}\rho_{2}\bigl{(}|\nabla\Phi_{2}(\bm{x},z,t)|^{2}+(\partial_{z}\Phi_{2}(\bm{x},z,t))^{2}\bigr{)}{\rm
d}\bm{x}{\rm d}z$
$\displaystyle\quad\;+\int_{\mathbf{R}^{n}}\frac{1}{2}(\rho_{2}-\rho_{1})g\zeta(\bm{x},t)^{2}{\rm
d}\bm{x},$

where $\Omega_{1}(t)$ and $\Omega_{2}(t)$ are the upper and the lower layers,
respectively. This is the total energy, that is, the sum of the kinetic
energies of the waters in the upper and the lower layers and the potential
energy due to the gravity. Moreover, T. B. Benjamin and T. J. Bridges [1]
found that the full model can be written in Hamilton’s canonical form

$\partial_{t}\zeta=\frac{\delta\mathscr{H}}{\delta\phi},\quad\partial_{t}\phi=-\frac{\delta\mathscr{H}}{\delta\zeta},$

where the canonical variable $\phi$ is defined by

(1.10)
$\phi(\bm{x},t)=\rho_{2}\Phi_{2}(\bm{x},\zeta(\bm{x},t),t)-\rho_{1}\Phi_{1}(\bm{x},\zeta(\bm{x},t),t)$

and the Hamiltonian $\mathscr{H}$ is the total energy $\mathscr{E}$ written in
terms of the canonical variables $(\zeta,\phi)$. Their result can be viewed as
a generalization into interfacial gravity waves of Zakharov’s Hamiltonian [26]
for water waves. For mathematical treatments of the Hamiltonian for
interfacial gravity waves, we refer to W. Craig and M. D. Groves [4] and W.
Craig, P. Guyenne, and H. Kalisch [5]. The Kakinuma model has also a conserved
energy $\mathscr{E}^{\rm K}$, which is the total energy given by (1.9) with
$\Phi_{1}$ and $\Phi_{2}$ replaced by $\Phi_{1}^{\rm app}$ and $\Phi_{2}^{\rm
app}$. Moreover, we will show that the Kakinuma model enjoys a Hamiltonian
structure with a Hamiltonian $\mathscr{H}^{\rm K}$ the total energy in terms
of canonical variables $\zeta$ and $\phi$, where $\phi$ is defined by (1.10)
with $\Phi_{1}$ and $\Phi_{2}$ replaced by $\Phi_{1}^{\rm app}$ and
$\Phi_{2}^{\rm app}$. This fact can be viewed as a generalization to the
Kakinuma model for interfacial gravity waves of a Hamiltonian structure of the
Isobe–Kakinuma model for water waves given by V. Duchêne and T. Iguchi [7].

The contents of this paper are as follows. In Section 2 we begin with
reviewing the full model for interfacial gravity waves and derive the Kakinuma
model. Then, we state one of the main results of this paper, that is, Theorem
2.1 about the well-posedness of the initial value problem to the Kakinuma
model locally in time. In Section 3 we analyze the linear dispersion relation
of the linearized equations of the Kakinuma model around the rest state in the
case of the flat bottom and show (1.7). In Section 4 we derive the stability
condition (1.8) by analyzing the linearized equations of the Kakinuma model
around an arbitrary flow. In Section 5 we derive an energy estimate for the
linearized equations with freezed coefficients and then transform the
equations into a standard positive symmetric system by introducing an
appropriate symmetrizer. In Section 6 we introduce several differential
operators related to the Kakinuma model and derive elliptic estimates for
these operators. In Section 7 we prove one of our main result, Theorem 2.1, by
using the method of parabolic regularization of the equations. In Section 8 we
prove another main result Theorem 8.4, which ensures that the Kakinuma model
enjoys a Hamiltonian structure. Finally, in Section 9 we derive conservation
laws of mass, momentum, and energy to the Kakinuma model together with the
corresponding flux functions.

Notation. We denote by $W^{m,p}(\mathbf{R}^{n})$ the $L^{p}$ Sobolev space of
order $m$ on $\mathbf{R}^{n}$ and $H^{m}=W^{m,2}(\mathbf{R}^{n})$. The norm of
a Banach space $B$ is denoted by $\|\cdot\|_{B}$. The $L^{2}$-inner product is
denoted by $(\cdot,\cdot)_{L^{2}}$. We put
$\partial_{t}=\frac{\partial}{\partial t}$,
$\partial_{j}=\partial_{x_{j}}=\frac{\partial}{\partial x_{j}}$, and
$\partial_{z}=\frac{\partial}{\partial z}$. $[P,Q]=PQ-QP$ denotes the
commutator and $[P;u,v]=P(uv)-(Pu)v-u(Pv)$ denotes the symmetric commutator.
For a matrix $A$ we denote by $A^{\rm T}$ the transpose of $A$. For a vector
$\mbox{\boldmath$\phi$}=(\phi_{0},\phi_{1},\ldots,\phi_{N})^{\rm T}$ we denote
the last $N$ components by
$\mbox{\boldmath$\phi$}^{\prime}=(\phi_{1},\ldots,\phi_{N})^{\rm T}$. We use
the notational convention $\frac{0}{0}=0$. We denote by
$C(a_{1},a_{2},\ldots)$ a positive constant depending on $a_{1},a_{2},\ldots$.
$f\lesssim g$ means that there exists a non-essential positive constant $C$
such that $f\leq Cg$ holds. $f\simeq g$ means that $f\lesssim g$ and
$g\lesssim f$ hold.

Acknowledgement
T. I. was partially supported by JSPS KAKENHI Grant Number JP17K18742 and
JP17H02856.

##  2 Kakinuma model and well-posedness

We begin with formulating mathematically the full model for interfacial
gravity waves. In what follows, the upper layer, the lower layer, the
interface, the rigid-lid of the upper layer, and the bottom of the lower
layer, at time $t$, are denoted by $\Omega_{1}(t)$, $\Omega_{2}(t)$,
$\Gamma(t)$, $\Sigma_{t}$, and $\Sigma_{b}$, respectively. Then, the motion of
the waters is described by the velocity potentials $\Phi_{1}$ and $\Phi_{2}$
and the pressures $P_{1}$ and $P_{2}$ in the upper and the lower layers
satisfying the equations of continuity

(2.1)
$\displaystyle\Delta\Phi_{1}+\partial_{z}^{2}\Phi_{1}=0\quad\mbox{in}\quad\Omega_{1}(t),$
(2.2)
$\displaystyle\Delta\Phi_{2}+\partial_{z}^{2}\Phi_{2}=0\quad\mbox{in}\quad\Omega_{2}(t),$

where $\Delta=\partial_{1}^{2}+\cdots+\partial_{n}^{2}$ is the Laplacian with
respect to the horizontal spatial coordinates $\bm{x}=(x_{1},\ldots,x_{n})$,
and Bernoulli’s equations

(2.3)
$\displaystyle\rho_{1}\left(\partial_{t}\Phi_{1}+\frac{1}{2}(|\nabla\Phi_{1}|^{2}+(\partial_{z}\Phi_{1})^{2})+gz\right)+P_{1}=0\quad\mbox{in}\quad\Omega_{1}(t),$
(2.4)
$\displaystyle\rho_{2}\left(\partial_{t}\Phi_{2}+\frac{1}{2}(|\nabla\Phi_{2}|^{2}+(\partial_{z}\Phi_{2})^{2})+gz\right)+P_{2}=0\quad\mbox{in}\quad\Omega_{2}(t).$

The dynamical boundary condition on the interface is given by

(2.5) $P_{1}=P_{2}\quad\mbox{on}\quad\Gamma(t).$

The kinematic boundary conditions on the interface, the rigid-lid, and the
bottom are given by

(2.6)
$\displaystyle\partial_{t}\zeta+\nabla\Phi_{1}\cdot\nabla\zeta-\partial_{z}\Phi_{1}=0\quad\mbox{on}\quad\Gamma(t),$
(2.7)
$\displaystyle\partial_{t}\zeta+\nabla\Phi_{2}\cdot\nabla\zeta-\partial_{z}\Phi_{2}=0\quad\mbox{on}\quad\Gamma(t),$
(2.8)
$\displaystyle\partial_{z}\Phi_{1}=0\makebox[78.79143pt]{}\quad\mbox{on}\quad\Sigma_{t},$
(2.9) $\displaystyle\nabla\Phi_{2}\cdot\nabla
b-\partial_{z}\Phi_{2}=0\makebox[27.12492pt]{}\quad\mbox{on}\quad\Sigma_{b}.$

These are the basic equations for interfacial gravity waves. We can remove the
pressures $P_{1}$ and $P_{2}$ from these basic equations. In fact, it follows
from Bernoulli’s equations (2.3)–(2.4) and the dynamical boundary condition
(2.5) that

(2.10)
$\displaystyle\rho_{1}\left(\partial_{t}\Phi_{1}+\frac{1}{2}(|\nabla\Phi_{1}|^{2}+(\partial_{z}\Phi_{1})^{2})+gz\right)$
$\displaystyle-\rho_{2}\left(\partial_{t}\Phi_{2}+\frac{1}{2}(|\nabla\Phi_{2}|^{2}+(\partial_{z}\Phi_{2})^{2})+gz\right)=0\quad\mbox{on}\quad\Gamma(t).$

Then, the basic equations consist of (2.1)–(2.2) and (2.6)–(2.10), and we can
regard Bernoulli’s equations (2.3)–(2.4) as the definition of the pressures
$P_{1}$ and $P_{2}$.

In the case of surface gravity waves, as was shown by J. C. Luke [20], the
basic equations have a variational structure and Luke’s Lagrangian density is
given by the vertical integral of the pressure $P-P_{\rm atm}$ in the water
region, where $P_{\rm atm}$ is a constant atmospheric pressure. Therefore, it
is natural to expect that even in the case of interfacial gravity waves the
vertical integral of the pressure in the water regions would give a Lagrangian
density $\mathscr{L}$, so that we first define $\mathscr{L}^{\rm pre}$ by

(2.11) $\mathscr{L}^{\rm
pre}=\int_{\zeta(\bm{x},t)}^{h_{1}}P_{1}(\bm{x},z,t){\rm
d}z+\int_{-h_{2}+b(\bm{x})}^{\zeta(\bm{x},t)}P_{2}(\bm{x},z,t){\rm d}z.$

By Bernoulli’s equations (2.3)–(2.4), this can be written in terms of the
velocity potentials $\Phi_{1}$, $\Phi_{2}$, and the elevation of the interface
$\zeta$ as

$\displaystyle\mathscr{L}^{\rm pre}$
$\displaystyle=-\rho_{1}\int_{\zeta}^{h_{1}}\left(\partial_{t}\Phi_{1}+\frac{1}{2}(|\nabla\Phi_{1}|^{2}+(\partial_{z}\Phi_{1})^{2})\right){\rm
d}z$
$\displaystyle\quad\;-\rho_{2}\int_{-h_{2}+b}^{\zeta}\left(\partial_{t}\Phi_{2}+\frac{1}{2}(|\nabla\Phi_{2}|^{2}+(\partial_{z}\Phi_{2})^{2})\right){\rm
d}z$
$\displaystyle\quad\;-\frac{1}{2}(\rho_{2}-\rho_{1})g\zeta^{2}-\frac{1}{2}\rho_{1}gh_{1}^{2}+\frac{1}{2}\rho_{2}g(-h_{2}+b)^{2}.$

The last two terms do not contribute to the calculus of variations of this
Lagrangian, so that we define the Lagrangian density
$\mathscr{L}(\Phi_{1},\Phi_{2},\zeta)$ by

(2.12) $\displaystyle\mathscr{L}(\Phi_{1},\Phi_{2},\zeta)$
$\displaystyle=-\rho_{1}\int_{\zeta}^{h_{1}}\left(\partial_{t}\Phi_{1}+\frac{1}{2}(|\nabla\Phi_{1}|^{2}+(\partial_{z}\Phi_{1})^{2})\right){\rm
d}z$
$\displaystyle\quad\;-\rho_{2}\int_{-h_{2}+b}^{\zeta}\left(\partial_{t}\Phi_{2}+\frac{1}{2}(|\nabla\Phi_{2}|^{2}+(\partial_{z}\Phi_{2})^{2})\right){\rm
d}z-\frac{1}{2}(\rho_{2}-\rho_{1})g\zeta^{2}$

and the action function $\mathscr{J}(\Phi_{1},\Phi_{2},\zeta)$ by

$\mathscr{J}(\Phi_{1},\Phi_{2},\zeta)=\int_{t_{0}}^{t_{1}}\\!\\!\\!\int_{\mathbf{R}^{n}}\mathscr{L}(\Phi_{1},\Phi_{2},\zeta){\rm
d}\bm{x}{\rm d}t.$

It is not difficult to check that the corresponding system of Euler–Lagrange
equations is exactly the same as the basic equations (2.1)–(2.2) and
(2.6)–(2.10) for interfacial gravity waves.

We proceed to derive the Kakinuma model for interfacial gravity waves. Let
$\Phi_{1}^{\rm app}$ and $\Phi_{2}^{\rm app}$ be approximate velocity
potentials defined by (1.4) and define an approximate Lagrangian density
$\mathscr{L}^{\rm app}(\bm{\phi}_{1},\bm{\phi}_{2},\zeta)$ for
$\bm{\phi}_{1}=(\phi_{1,0},\phi_{1,1},\ldots,\phi_{1,N})^{\rm T}$,
$\bm{\phi}_{2}=(\phi_{2,0},\phi_{2,1},\ldots,\phi_{2,N^{*}})^{\rm T}$, and
$\zeta$ by

(2.13) $\mathscr{L}^{\rm
app}(\bm{\phi}_{1},\bm{\phi}_{2},\zeta)=\mathscr{L}(\Phi_{1}^{\rm
app},\Phi_{2}^{\rm app},\zeta),$

which can be written explicitly as

$\displaystyle\mathscr{L}^{\rm app}$
$\displaystyle=\rho_{1}\left\\{\sum_{i=0}^{N}\frac{1}{2i+1}H_{1}^{2i+1}\partial_{t}\phi_{1,i}\right.$
$\displaystyle\qquad\left.+\frac{1}{2}\sum_{i,j=0}^{N}\left(\frac{1}{2(i+j)+1}H_{1}^{2(i+j)+1}\nabla\phi_{1,i}\cdot\nabla\phi_{1,j}+\frac{4ij}{2(i+j)-1}H_{1}^{2(i+j)-1}\phi_{1,i}\phi_{1,j}\right)\right\\}$
$\displaystyle\quad-\rho_{2}\left\\{\sum_{i=0}^{N^{*}}\frac{1}{p_{i}+1}H_{2}^{p_{i}+1}\partial_{t}\phi_{2,i}\right.$
$\displaystyle\qquad\quad+\frac{1}{2}\sum_{i,j=0}^{N^{*}}\left(\frac{1}{p_{i}+p_{j}+1}H_{2}^{p_{i}+p_{j}+1}\nabla\phi_{2,i}\cdot\nabla\phi_{2,j}-\frac{2p_{i}}{p_{i}+p_{j}}H_{2}^{p_{i}+p_{j}}\phi_{2,i}\nabla
b\cdot\nabla\phi_{2,j}\right.$
$\displaystyle\left.\phantom{\qquad\quad+\frac{1}{2}\sum_{i,j=0}^{N^{*}}\biggl{(}}\left.+\frac{p_{i}p_{j}}{p_{i}+p_{j}-1}H_{2}^{p_{i}+p_{j}-1}(1+|\nabla
b|^{2})\phi_{2,i}\phi_{2,j}\right)\right\\}$
$\displaystyle\quad-\frac{1}{2}(\rho_{2}-\rho_{1})g\zeta^{2},$

where $H_{1}$ and $H_{2}$ are thicknesses of the upper and the lower layers,
that is,

$H_{1}(\bm{x},t)=h_{1}-\zeta(\bm{x},t),\qquad
H_{2}(\bm{x},t)=h_{2}+\zeta(\bm{x},t)-b(\bm{x}).$

The corresponding system of Euler–Lagrange equations is the Kakinuma model,
which consists of the equations

(2.14)
$H_{1}^{2i}\partial_{t}\zeta-\sum_{j=0}^{N}\left\\{\nabla\cdot\left(\frac{1}{2(i+j)+1}H_{1}^{2(i+j)+1}\nabla\phi_{1,j}\right)-\frac{4ij}{2(i+j)-1}H_{1}^{2(i+j)-1}\phi_{1,j}\right\\}=0$

for $i=0,1,\ldots,N$,

(2.15) $\displaystyle
H_{2}^{p_{i}}\partial_{t}\zeta+\sum_{j=0}^{N^{*}}\left\\{\nabla\cdot\left(\frac{1}{p_{i}+p_{j}+1}H_{2}^{p_{i}+p_{j}+1}\nabla\phi_{2,j}-\frac{p_{j}}{p_{i}+p_{j}}H_{2}^{p_{i}+p_{j}}\phi_{2,j}\nabla
b\right)\right.$
$\displaystyle\phantom{H_{2}^{p_{i}}\partial_{t}\zeta+\sum_{j=0}^{N^{*}}}\left.+\frac{p_{i}}{p_{i}+p_{j}}H_{2}^{p_{i}+p_{j}}\nabla
b\cdot\nabla\phi_{2,j}-\frac{p_{i}p_{j}}{p_{i}+p_{j}-1}H_{2}^{p_{i}+p_{j}-1}(1+|\nabla
b|^{2})\phi_{2,j}\right\\}=0$

for $i=0,1,\ldots,N^{*}$, and

(2.16)
$\displaystyle\rho_{1}\left\\{\sum_{j=0}^{N}H_{1}^{2j}\partial_{t}\phi_{1,j}+g\zeta+\frac{1}{2}\left(\left|\sum_{j=0}^{N}H_{1}^{2j}\nabla\phi_{1,j}\right|^{2}+\left(\sum_{j=0}^{N}2jH_{1}^{2j-1}\phi_{1,j}\right)^{2}\right)\right\\}$
$\displaystyle\quad-\rho_{2}\left\\{\sum_{j=0}^{N^{*}}H_{2}^{p_{j}}\partial_{t}\phi_{2,j}+g\zeta\right.$
$\displaystyle\qquad\quad\left.+\frac{1}{2}\left(\left|\sum_{j=0}^{N^{*}}(H_{2}^{p_{j}}\nabla\phi_{2,j}-p_{j}H_{2}^{p_{j}-1}\phi_{2,j}\nabla
b)\right|^{2}+\left(\sum_{j=0}^{N^{*}}p_{j}H_{2}^{p_{j}-1}\phi_{2,j}\right)^{2}\right)\right\\}=0.$

Here and in what follows we use the notational convention $\frac{0}{0}=0$.
This system of equations is the Kakinuma model that we are going to consider
in this paper. We consider the initial value problem to the Kakinuma model
(2.14)–(2.16) under the initial condition

(2.17)
$(\zeta,\bm{\phi}_{1},\bm{\phi}_{2})=(\zeta_{(0)},\bm{\phi}_{1(0)},\bm{\phi}_{2(0)})\makebox[30.00005pt]{at}t=0.$

For notational convenience, we decompose $\bm{\phi}_{k}$ as
$\bm{\phi}_{k}=(\phi_{k,0},\bm{\phi}_{k}^{\prime})^{\rm T}$ for $k=1,2$ with
$\bm{\phi}_{1}^{\prime}=(\phi_{1,1},\ldots,\phi_{1,N})$ and
$\bm{\phi}_{2}^{\prime}=(\phi_{2,1},\ldots,\phi_{2,N^{*}})$. Accordingly, we
decompose the initial data $\bm{\phi}_{k(0)}$ as
$\bm{\phi}_{k(0)}=(\phi_{k,0(0)},\bm{\phi}_{k(0)}^{\prime})^{\rm T}$ for
$k=1,2$.

The hypersurface $t=0$ in the space-time $\mathbf{R}^{n}\times\mathbf{R}$ is
characteristic for the Kakinuma model (2.14)–(2.16), so that the initial value
problem (2.14)–(2.17) is not solvable in general. In fact, by eliminating the
time derivative $\partial_{t}\zeta$ from the equations, we see that if the
problem has a solution $(\zeta,\bm{\phi}_{1},\bm{\phi}_{2})$, then the
solution has to satisfy the $N+N^{*}+1$ relations

(2.18) $\displaystyle
H_{1}^{2i}\sum_{j=0}^{N}\nabla\cdot\left(\frac{1}{2j+1}H_{1}^{2j+1}\nabla\phi_{1,j}\right)$
$\displaystyle\quad-\sum_{j=0}^{N}\left\\{\nabla\cdot\left(\frac{1}{2(i+j)+1}H_{1}^{2(i+j)+1}\nabla\phi_{1,j}\right)-\frac{4ij}{2(i+j)-1}H_{1}^{2(i+j)-1}\phi_{1,j}\right\\}=0$

for $i=1,2,\ldots,N$,

(2.19) $\displaystyle
H_{2}^{p_{i}}\sum_{j=0}^{N^{*}}\nabla\cdot\left(\frac{1}{p_{j}+1}H_{2}^{p_{j}+1}\nabla\phi_{2,j}-\frac{p_{j}}{p_{j}}H_{2}^{p_{j}}\phi_{2,j}\nabla
b\right)$
$\displaystyle\quad-\sum_{j=0}^{N^{*}}\left\\{\nabla\cdot\left(\frac{1}{p_{i}+p_{j}+1}H_{2}^{p_{i}+p_{j}+1}\nabla\phi_{2,j}-\frac{p_{j}}{p_{i}+p_{j}}H_{2}^{p_{i}+p_{j}}\phi_{2,j}\nabla
b\right)\right.$
$\displaystyle\phantom{\quad+\sum_{j=0}^{N^{*}}\bigg{\\{}}\left.+\frac{p_{i}}{p_{i}+p_{j}}H_{2}^{p_{i}+p_{j}}\nabla
b\cdot\nabla\phi_{2,j}-\frac{p_{i}p_{j}}{p_{i}+p_{j}-1}H_{2}^{p_{i}+p_{j}-1}(1+|\nabla
b|^{2})\phi_{2,j}\right\\}=0$

for $i=1,2,\ldots,N^{*}$, and

(2.20)
$\displaystyle\sum_{j=0}^{N}\nabla\cdot\left(\frac{1}{2j+1}H_{1}^{2j+1}\nabla\phi_{1,j}\right)+\sum_{j=0}^{N^{*}}\nabla\cdot\left(\frac{1}{p_{j}+1}H_{2}^{p_{j}+1}\nabla\phi_{2,j}-\frac{p_{j}}{p_{j}}H_{2}^{p_{j}}\phi_{2,j}\nabla
b\right)=0.$

Therefore, as a necessary condition the initial date
$(\zeta_{(0)},\bm{\phi}_{1(0)},\bm{\phi}_{2(0)})$ and the bottom topography
$b$ have to satisfy the relation (2.18)–(2.20) for the existence of the
solution. These necessary conditions will be referred as the compatibility
conditions.

The following theorem is one of our main results in this paper, which
guarantees the well-posedness of the initial value problem to the Kakinuma
model (2.14)–(2.17) locally in time.

######  Theorem 2.1.

Let $g,\rho_{1},\rho_{2},h_{1},h_{2},c_{0},M_{0}$ be positive constants and
$m$ an integer such that $m>\frac{n}{2}+1$. There exists a time $T>0$ such
that for any initial data $(\zeta_{(0)},\bm{\phi}_{1(0)},\bm{\phi}_{2(0)})$
and bottom topography $b$ satisfying the compatibility conditions
(2.18)–(2.20), the stability condition (1.8), and

(2.21)
$\begin{cases}\|(\zeta_{(0)},\nabla\phi_{1,0(0)},\nabla\phi_{2,0(0)})\|_{H^{m}}+\|(\bm{\phi}_{1(0)}^{\prime},\bm{\phi}_{2(0)}^{\prime})\|_{H^{m+1}}+\|b\|_{W^{m+2,\infty}}\leq
M_{0},\\\ h_{1}-\zeta_{(0)}(\bm{x})\geq c_{0},\quad
h_{2}+\zeta_{(0)}(\bm{x})-b(\bm{x})\geq
c_{0}\makebox[30.00005pt]{for}\bm{x}\in\mathbf{R}^{n},\end{cases}$

the initial value problem (2.14)–(2.17) has a unique solution
$(\zeta,\bm{\phi}_{1},\bm{\phi}_{2})$ satisfying

$\begin{cases}\zeta,\nabla\phi_{1,0},\nabla\phi_{2,0}\in C([0,T];H^{m})\cap
C^{1}([0,T];H^{m-1}),\\\ \bm{\phi}_{1}^{\prime},\bm{\phi}_{2}^{\prime}\in
C([0,T];H^{m+1})\cap C^{1}([0,T];H^{m}).\end{cases}$

Moreover, the flow map is continuous.

######  Remark 2.2.

The term $\left(\partial_{z}(P_{2}^{\rm app}-P_{1}^{\rm
app})\right)|_{z=\zeta}$ in the stability condition (1.8) is explicitly given
in (4.4). It includes the terms $\partial_{t}\bm{\phi}_{k}(\bm{x},0)$ for
$k=1,2$. Although the hypersurface $t=0$ is characteristic for the Kakinuma
model, we can uniquely determine them in terms of the initial data and $b$.
For details, we refer to Remark 7.1. Under the condition
$(\rho_{2}-\rho_{1})g>0$ and if the initial data and the bottom topography are
suitably small, the stability condition (1.8) is automatically satisfied at
$t=0$.

######  Remark 2.3.

In the case $N=N^{*}=0$, that is, if we approximate the velocity potentials in
the Lagrangian by functions independent of the vertical spatial variable $z$
as $\Phi_{k}^{\rm app}(\bm{x},z,t)=\phi_{k}(\bm{x},t)$ for $k=1,2$, then the
Kakinuma model (2.14)–(2.16) is reduced to the nonlinear shallow water
equations

(2.22)
$\begin{cases}\partial_{t}\zeta-\nabla\cdot((h_{1}-\zeta)\nabla\phi_{1})=0,\\\
\partial_{t}\zeta+\nabla\cdot((h_{2}+\zeta-b)\nabla\phi_{2})=0,\\\
\displaystyle\rho_{1}\left(\partial_{t}\phi_{1}+g\zeta+\frac{1}{2}|\nabla\phi_{1}|^{2}\right)-\rho_{2}\left(\partial_{t}\phi_{2}+g\zeta+\frac{1}{2}|\nabla\phi_{2}|^{2}\right)=0.\end{cases}$

The compatibility conditions (2.18)–(2.20) are reduced to

$\nabla\cdot((h_{1}-\zeta)\nabla\phi_{1})+\nabla\cdot((h_{2}+\zeta-b)\nabla\phi_{2})=0$

and the stability condition (1.8) is reduced to

$g(\rho_{2}-\rho_{1})-\frac{\rho_{1}\rho_{2}}{\rho_{1}H_{2}+\rho_{2}H_{1}}|\nabla\phi_{2}-\nabla\phi_{1}|^{2}\geq
c_{0}>0.$

Therefore, we recover the conditions for the well-posedness in Sobolev spaces
of the initial value problem to the nonlinear shallow water equations (2.22)
proved by D. Bresch and M. Renardy [2].

######  Remark 2.4.

By analogy to the canonical variable (1.10) for interfacial gravity waves
introduced by T. B. Benjamin and T. J. Bridges [1], we introduce a canonical
variable for the Kakinuma model by

(2.23)
$\phi=\rho_{2}\sum_{j=0}^{N^{*}}H_{2}^{p_{j}}\phi_{2,j}-\rho_{1}\sum_{j=0}^{N}H_{1}^{2j}\phi_{1,j}.$

Given the initial data $(\zeta_{(0)},\phi_{(0)})$ for the canonical variables
$(\zeta,\phi)$ and the bottom topography $b$, the compatibility conditions
(2.18)–(2.20) and the relation (2.23) determine the initial data
$(\bm{\phi}_{1(0)},\bm{\phi}_{2(0)})$ for the Kakinuma model, which is unique
up to an additive constant of the form
$(\mathcal{C}\rho_{2},\mathcal{C}\rho_{1})$ to
$(\phi_{1,0(0)},\phi_{2,0(0)})$. In fact, we have the following proposition,
which is a simple corollary of Lemma 6.4 given in Section 6.

######  Proposition 2.5.

Let $\rho_{1},\rho_{2},h_{1},h_{2},c_{0},M_{0}$ be positive constants and $m$
an integer such that $m>\frac{n}{2}+1$. There exists a positive constant $C$
such that for any initial data $(\zeta_{(0)},\phi_{(0)})$ and bottom
topography $b$ satisfying

$\begin{cases}\|\zeta_{(0)}\|_{H^{m}}+\|b\|_{W^{m,\infty}}\leq
M_{0},\quad\|\nabla\phi_{(0)}\|_{H^{m-1}}<\infty,\\\
h_{1}-\zeta_{(0)}(\bm{x})\geq c_{0},\quad
h_{2}+\zeta_{(0)}(\bm{x})-b(\bm{x})\geq
c_{0}\makebox[30.00005pt]{for}\bm{x}\in\mathbf{R}^{n},\end{cases}$

the compatibility conditions (2.18)–(2.20) and the relation (2.23) determine
the initial data $(\bm{\phi}_{1(0)},\bm{\phi}_{2(0)})$ for the Kakinuma model,
uniquely up to an additive constant of the form
$(\mathcal{C}\rho_{2},\mathcal{C}\rho_{1})$ to
$(\phi_{1,0(0)},\phi_{2,0(0)})$. Moreover, we have

$\|(\nabla\phi_{1,0(0)},\nabla\phi_{2,0(0)})\|_{H^{m-1}}+\|(\bm{\phi}_{1(0)}^{\prime},\bm{\phi}_{2(0)}^{\prime})\|_{H^{m}}\leq
C\|\nabla\phi_{(0)}\|_{H^{m-1}}.$

Therefore, given the initial data $(\zeta_{(0)},\phi_{(0)})$, we infer initial
data for the Kakinuma model, which satisfy the compatibility conditions
(2.18)–(2.20).

##  3 Linear dispersion relation

In this section we consider the linearized equations of the Kakinuma model
(2.14)–(2.16) around the flow
$(\zeta,\bm{\phi}_{1},\bm{\phi}_{2})=(0,\bm{0},\bm{0})$ in the case of the
flat bottom. The linearized equations have the form

(3.1)
$\begin{cases}\displaystyle\partial_{t}\zeta-\sum_{j=0}^{N}\left(\frac{h_{1}^{2j+1}}{2(i+j)+1}\Delta\phi_{1,j}-\frac{4ij}{2(i+j)-1}h_{1}^{2j-1}\phi_{1,j}\right)=0\quad\mbox{for}\quad
i=0,1,\ldots,N,\\\
\displaystyle\partial_{t}\zeta+\sum_{j=0}^{N^{*}}\left(\frac{h_{2}^{p_{j}+1}}{p_{i}+p_{j}+1}\Delta\phi_{2,j}-\frac{p_{i}p_{j}}{p_{i}+p_{j}-1}h_{2}^{p_{j}-1}\phi_{2,j}\right)=0\quad\mbox{for}\quad
i=0,1,\ldots,N^{*},\\\
\displaystyle\rho_{1}\left(\sum_{j=0}^{N}h_{1}^{2j}\partial_{t}\phi_{1,j}+g\zeta\right)-\rho_{2}\left(\sum_{j=0}^{N^{*}}h_{2}^{p_{j}}\partial_{t}\phi_{2,j}+g\zeta\right)=0.\end{cases}$

Putting
$\bm{\psi}_{1}=(\phi_{1,0},h_{1}^{2}\phi_{1,1},\ldots,h_{1}^{2N}\phi_{1,N})^{\rm
T}$ and
$\bm{\psi}_{2}=(\phi_{2,0},h_{1}^{p_{1}}\phi_{2,1},\ldots,h_{1}^{p_{N^{*}}}\phi_{2,N^{*}})^{\rm
T}$, we can rewrite the above equations as the following simple matrix form

$\displaystyle\begin{pmatrix}0&-\rho_{1}\bm{1}^{\rm T}&\rho_{2}\bm{1}^{\rm
T}\\\ h_{1}\bm{1}&O&O\\\
-h_{2}\bm{1}&O&O\end{pmatrix}\partial_{t}\begin{pmatrix}\zeta\\\
\bm{\psi}_{1}\\\ \bm{\psi}_{2}\end{pmatrix}$
$\displaystyle\quad+\begin{pmatrix}(\rho_{2}-\rho_{1})g&\bm{0}^{\rm
T}&\bm{0}^{\rm T}\\\ \bm{0}&-h_{1}^{2}A_{1,0}\Delta+A_{1,1}&O\\\
\bm{0}&O&-h_{2}^{2}A_{2,0}\Delta+A_{2,1}\end{pmatrix}\begin{pmatrix}\zeta\\\
\bm{\psi}_{1}\\\ \bm{\psi}_{2}\end{pmatrix}=\bm{0},$

where ${\bm{1}}=(1,\ldots,1)^{\rm T}$ and matrices $A_{k,0}$ and $A_{k,1}$ for
$k=1,2$ are given by

$\displaystyle A_{1,0}=\biggl{(}\frac{1}{2(i+j)+1}\biggr{)}_{0\leq i,j\leq
N},\qquad A_{1,1}=\biggl{(}\frac{4ij}{2(i+j)-1}\biggr{)}_{0\leq i,j\leq N},$
$\displaystyle A_{2,0}=\biggl{(}\frac{1}{p_{i}+p_{j}+1}\biggr{)}_{0\leq
i,j\leq N^{*}},\qquad
A_{2,1}=\biggl{(}\frac{p_{i}p_{j}}{p_{i}+p_{j}-1}\biggr{)}_{0\leq i,j\leq
N^{*}}.$

Therefore, the linear dispersion relation is given by

$\det\begin{pmatrix}(\rho_{2}-\rho_{1})g&\mbox{\rm i}\rho_{1}\omega\bm{1}^{\rm
T}&-\mbox{\rm i}\rho_{2}\omega\bm{1}^{\rm T}\\\ -\mbox{\rm
i}h_{1}\omega\bm{1}&\mathcal{A}_{1}(h_{1}\bm{\xi})&O\\\ \mbox{\rm
i}h_{2}\omega\bm{1}&O&\mathcal{A}_{2}(h_{2}\bm{\xi})\end{pmatrix}=0,$

where $\bm{\xi}\in\mathbf{R}^{n}$ is the wave vector, $\omega\in\mathbf{C}$ is
the angular frequency, and
$\mathcal{A}_{k}(\bm{\xi})=|\bm{\xi}|^{2}A_{k,0}+A_{k,1}$ for $k=1,2$. We can
expand this dispersion relation as

(3.2)
$\left(\rho_{1}h_{1}\det\tilde{\mathcal{A}}_{1}(h_{1}\bm{\xi})\det\mathcal{A}_{2}(h_{2}\bm{\xi})+\rho_{2}h_{2}\det\tilde{\mathcal{A}}_{2}(h_{2}\bm{\xi})\det\mathcal{A}_{1}(h_{1}\bm{\xi})\right)\omega^{2}\\\
-(\rho_{2}-\rho_{1})g\det\mathcal{A}_{1}(h_{1}\bm{\xi})\det\mathcal{A}_{2}(h_{2}\bm{\xi})=0.$

Here and in what follows, we use the notation

$\tilde{\mathcal{A}}=\begin{pmatrix}0&{\bm{1}}^{\rm T}\\\
-{\bm{1}}&\mathcal{A}\end{pmatrix}$

for a matrix $\mathcal{A}$. Concerning the determinants appearing in the above
dispersion relation, we have the following proposition, which was proved by R.
Nemoto and T. Iguchi [23].

######  Proposition 3.1.

  1. 1.

For any $\bm{\xi}\in\mathbf{R}^{n}\setminus\\{{\bm{0}}\\}$, the symmetric
matrices $\mathcal{A}_{1}(\bm{\xi})$ and $\mathcal{A}_{2}(\bm{\xi})$ are
positive.

  2. 2.

There exists $c_{0}>0$ such that for any $\bm{\xi}\in\mathbf{R}^{n}$ we have
$\det\tilde{\mathcal{A}}_{k}(\bm{\xi})\geq c_{0}$ for $k=1,2$.

  3. 3.

$|\bm{\xi}|^{-2}\det\mathcal{A}_{1}(\bm{\xi})$ and
$|\bm{\xi}|^{-2}\det\mathcal{A}_{2}(\bm{\xi})$ are polynomials in
$|\bm{\xi}|^{2}$ of degree $N$ and $N^{*}$ and their coefficient of
$|\bm{\xi}|^{2N}$ and $|\bm{\xi}|^{2N^{*}}$ are $\det A_{1,0}$ and $\det
A_{2,0}$, respectively.

  4. 4.

$\det\tilde{\mathcal{A}}_{1}(\bm{\xi})$ and
$\det\tilde{\mathcal{A}}_{2}(\bm{\xi})$ are polynomials in $|\bm{\xi}|^{2}$ of
degree $N$ and $N^{*}$ and their coefficient of $|\bm{\xi}|^{2N}$ and
$|\bm{\xi}|^{2N^{*}}$ are $\det\tilde{A}_{1,0}$ and $\det\tilde{A}_{2,0}$,
respectively.

Thanks of this proposition and the dispersion relation (3.2), the linearized
system (3.1) is classified into the dispersive system in the case $N+N^{*}>0$,
so that the Kakinuma model (2.14)–(2.16) is a nonlinear dispersive system of
equations.

Therefore, we can define the phase speed $c_{\mbox{\rm\tiny K}}(\bm{\xi})$ of
the plane wave solution to (3.1) related to the wave vector
$\bm{\xi}\in\mathbf{R}^{n}$ by

(3.3) $c_{\mbox{\rm\tiny
K}}(\bm{\xi})^{2}=\frac{(\rho_{2}-\rho_{1})g|\bm{\xi}|^{-2}\det\mathcal{A}_{1}(h_{1}\bm{\xi})\det\mathcal{A}_{2}(h_{2}\bm{\xi})}{\rho_{1}h_{1}\det\tilde{\mathcal{A}}_{1}(h_{1}\bm{\xi})\det\mathcal{A}_{2}(h_{2}\bm{\xi})+\rho_{2}h_{2}\det\tilde{\mathcal{A}}_{2}(h_{2}\bm{\xi})\det\mathcal{A}_{1}(h_{1}\bm{\xi})}.$

It follows from Proposition 3.1 that

$\lim_{h_{1}|\bm{\xi}|,h_{2}|\bm{\xi}|\to\infty}c_{\mbox{\rm\tiny
K}}(\bm{\xi})^{2}=\frac{(\rho_{2}-\rho_{1})gh_{1}h_{2}\det A_{1,0}\det
A_{2,0}}{\rho_{1}h_{2}\det\tilde{A}_{1,0}\det
A_{2,0}+\rho_{2}h_{1}\det\tilde{A}_{2,0}\det A_{1,0}}>0,$

which is not consistent with the linear interfacial gravity waves

$\lim_{h_{1}|\bm{\xi}|,h_{2}|\bm{\xi}|\to\infty}c_{\mbox{\rm\tiny
IW}}(\bm{\xi})^{2}=0.$

However, as is shown by the following theorems the Kakinuma model gives a very
precise approximation in the shallow water regime
$h_{1}|\bm{\xi}|,h_{2}|\bm{\xi}|\ll 1$ under an appropriate choice of the
indices $p_{i}$ for $i=0,1,\ldots,N^{*}$.

######  Theorem 3.2.

If we choose $N^{*}=N$ and $p_{i}=2i$ for $i=0,1,\ldots,N^{*}$ or $N^{*}=2N$
and $p_{i}=i$ for $i=0,1,\ldots,N^{*}$, then for any
$\bm{\xi}\in\mathbf{R}^{n}$ and any $h_{1},h_{2},g>0$ we have

$\left|\left(\frac{c_{\mbox{\rm\tiny IW}}(\bm{\xi})}{c_{\mbox{\rm\tiny
SW}}}\right)^{2}-\left(\frac{c_{\mbox{\rm\tiny
K}}(\bm{\xi})}{c_{\mbox{\rm\tiny SW}}}\right)^{2}\right|\leq
C(h_{1}|\bm{\xi}|+h_{2}|\bm{\xi}|)^{4N+2},$

where $C$ is a positive constant depending only on $N$.

###### Proof..

The phase speeds $c_{\mbox{\rm\tiny IW}}(\bm{\xi})$ and $c_{\mbox{\rm\tiny
K}}(\bm{\xi})$ can be written in the form

$\left(\frac{c_{\mbox{\rm\tiny IW}}(\bm{\xi})}{c_{\mbox{\rm\tiny
SW}}}\right)^{2}=\frac{\displaystyle\frac{\tanh(h_{1}|\bm{\xi}|)}{h_{1}|\bm{\xi}|}\frac{\tanh(h_{2}|\bm{\xi}|)}{h_{2}|\bm{\xi}|}}{\displaystyle\theta\frac{\tanh(h_{1}|\bm{\xi}|)}{h_{1}|\bm{\xi}|}+(1-\theta)\frac{\tanh(h_{2}|\bm{\xi}|)}{h_{2}|\bm{\xi}|}}$

and

$\left(\frac{c_{\mbox{\rm\tiny K}}(\bm{\xi})}{c_{\mbox{\rm\tiny
SW}}}\right)^{2}=\frac{\displaystyle\frac{\det\mathcal{A}_{1}(h_{1}\bm{\xi})}{(h_{1}|\bm{\xi}|)^{2}\det\tilde{\mathcal{A}}_{1}(h_{1}\bm{\xi})}\frac{\det\mathcal{A}_{2}(h_{2}\bm{\xi})}{(h_{2}|\bm{\xi}|)^{2}\det\tilde{\mathcal{A}}_{2}(h_{2}\bm{\xi})}}{\displaystyle\theta\frac{\det\mathcal{A}_{1}(h_{1}\bm{\xi})}{(h_{1}|\bm{\xi}|)^{2}\det\tilde{\mathcal{A}}_{1}(h_{1}\bm{\xi})}+(1-\theta)\frac{\det\mathcal{A}_{2}(h_{2}\bm{\xi})}{(h_{2}|\bm{\xi}|)^{2}\det\tilde{\mathcal{A}}_{2}(h_{2}\bm{\xi})}},$

respectively, where
$\theta=\frac{\rho_{2}h_{1}}{\rho_{2}h_{1}+\rho_{1}h_{2}}\in(0,1)$. It has
been shown by R. Nemoto and T. Iguchi [23] that

$\left|\frac{\tanh|\bm{\xi}|}{|\bm{\xi}|}-\frac{\det\mathcal{A}_{k}(\bm{\xi})}{|\bm{\xi}|^{2}\det\tilde{\mathcal{A}}_{k}(\bm{\xi})}\right|\leq
C|\bm{\xi}|^{4N+2}$

for $k=1,2$, so that we obtain the desired inequality. ∎

##  4 Stability condition

In this section, we will derive the stability condition (1.8) by analyzing a
system of linearized equations to the Kakinuma model (2.14)–(2.16). We
linearize the Kakinuma model around an arbitrary flow
$(\zeta,{\bm{\phi}}_{1},{\bm{\phi}}_{2})$ and denote the variation by
$(\dot{\zeta},\dot{{\bm{\phi}}}_{1},\dot{{\bm{\phi}}}_{2})$. After neglecting
lower order terms, the linearized equations have the form

(4.1)
$\begin{cases}\displaystyle\partial_{t}\dot{\zeta}+{\bm{u}}_{1}\cdot\nabla\dot{\zeta}-\sum_{j=0}^{N}\frac{1}{2(i+j)+1}H_{1}^{2j+1}\Delta\dot{\phi}_{1,j}=0\quad\mbox{for}\quad
i=0,1,\ldots,N,\\\
\displaystyle\partial_{t}\dot{\zeta}+{\bm{u}}_{2}\cdot\nabla\dot{\zeta}+\sum_{j=0}^{N^{*}}\frac{1}{p_{i}+p_{j}+1}H_{2}^{p_{j}+1}\Delta\dot{\phi}_{2,j}=0\quad\mbox{for}\quad
i=0,1,\ldots,N^{*},\\\
\displaystyle\rho_{1}\sum_{j=0}^{N}H_{1}^{2j}(\partial_{t}\dot{\phi}_{1,j}+{\bm{u}}_{1}\cdot\nabla\dot{\phi}_{1,j})-\rho_{2}\sum_{j=0}^{N^{*}}H_{2}^{p_{j}}(\partial_{t}\dot{\phi}_{2,j}+{\bm{u}}_{2}\cdot\nabla\dot{\phi}_{2,j})-a\dot{\zeta}=0,\end{cases}$

where $H_{1}=h_{1}-\zeta$ and $H_{2}=h_{2}+\zeta-b$ are the thicknesses of the
layers,

(4.2) $\begin{cases}\displaystyle{\bm{u}}_{1}=(\nabla\Phi_{1}^{\rm
app})|_{z=\zeta}=\sum_{j=0}^{N}H_{1}^{2j}\nabla\phi_{1,j},\\\
\displaystyle{\bm{u}}_{2}=(\nabla\Phi_{2}^{\rm
app})|_{z=\zeta}=\sum_{j=0}^{N^{*}}(H_{2}^{p_{j}}\nabla\phi_{2,j}-p_{j}H_{2}^{p_{j}-1}\phi_{2,j}\nabla
b)\end{cases}$

are approximate horizontal velocities in the upper and the lower layers on the
interface,

(4.3) $\begin{cases}\displaystyle w_{1}=(\partial_{z}\Phi_{1}^{\rm
app})|_{z=\zeta}=-\sum_{j=0}^{N}2jH_{1}^{2j-1}\phi_{1,j},\\\ \displaystyle
w_{2}=(\partial_{z}\Phi_{2}^{\rm
app})|_{z=\zeta}=\sum_{j=0}^{N^{*}}p_{j}H_{2}^{p_{j}-1}\phi_{2,j}\end{cases}$

are approximate vertical velocities in the upper and the lower layers on the
interface, and

(4.4) $\displaystyle a$
$\displaystyle=\rho_{2}\left(\sum_{j=0}^{N^{*}}p_{j}H_{2}^{p_{j}-1}(\partial_{t}\phi_{2,j}+{\bm{u}}_{2}\cdot\nabla\phi_{2,j})+(w_{2}-{\bm{u}}_{2}\cdot\nabla
b)\sum_{j=0}^{N^{*}}p_{j}(p_{j}-1)H_{2}^{p_{j}-2}\phi_{2,j}+g\right)$
$\displaystyle\quad+\rho_{1}\left(\sum_{j=0}^{N}2jH_{1}^{2j-1}(\partial_{t}\phi_{1,j}+{\bm{u}}_{1}\cdot\nabla\phi_{1,j})-w_{1}\sum_{j=0}^{N}2j(2j-1)H^{2(j-1)}\phi_{1,j}-g\right)$
$\displaystyle=-\left(\partial_{z}(P_{2}^{\rm app}-P_{1}^{\rm
app})\right)|_{z=\zeta}.$

Here, $P_{1}^{\rm app}$ and $P_{2}^{\rm app}$ are approximate pressures in the
upper and the lower layers calculated from Bernoulli’s equations (2.3)–(2.4),
that is,

$P_{k}^{\rm app}=-\rho_{k}\left(\partial_{t}\Phi_{k}^{\rm
app}+\frac{1}{2}\left(|\nabla\Phi_{k}^{\rm
app}|^{2}+(\partial_{z}\Phi_{k}^{\rm app})^{2}\right)+gz\right)$

for $k=1,2$. Now, we freeze the coefficients in the linearized equations (4.1)
and put

(4.5)
$\begin{cases}\dot{{\bm{\psi}}}_{1}=(\dot{\phi}_{1,0},H_{1}^{2}\dot{\phi}_{1,1},\ldots,H_{1}^{2N}\dot{\phi}_{1,N})^{\rm
T},\\\
\dot{{\bm{\psi}}}_{2}=(\dot{\phi}_{2,0},H_{2}^{p_{1}}\dot{\phi}_{2,1},\ldots,H_{2}^{p_{N^{*}}}\dot{\phi}_{2,N^{*}})^{\rm
T}.\end{cases}$

Then, (4.1) can be written in the form

$\displaystyle\begin{pmatrix}0&-\rho_{1}{\bm{1}}^{\rm T}&\rho_{2}{\bm{1}}^{\rm
T}\\\ H_{1}{\bm{1}}&O&O\\\
-H_{2}{\bm{1}}&O&O\end{pmatrix}\partial_{t}\begin{pmatrix}\dot{\zeta}\\\
\dot{\bm{\psi}}_{1}\\\ \dot{\bm{\psi}}_{2}\end{pmatrix}$
$\displaystyle\quad+\begin{pmatrix}a&-\rho_{1}\bm{1}^{\rm
T}(\bm{u}_{1}\cdot\nabla)&\rho_{2}\bm{1}^{\rm T}(\bm{u}_{2}\cdot\nabla)\\\
H_{1}\bm{1}(\bm{u}_{1}\cdot\nabla)&-H_{1}^{2}A_{1,0}\Delta&O\\\
-H_{2}\bm{1}(\bm{u}_{2}\cdot\nabla)&O&-H_{2}^{2}A_{2,0}\Delta\end{pmatrix}\begin{pmatrix}\dot{\zeta}\\\
\dot{\bm{\psi}}_{1}\\\ \dot{\bm{\psi}}_{2}\end{pmatrix}=\bm{0}.$

Therefore, the linear dispersion relation for (4.1) is given by

$\det\begin{pmatrix}a&\mbox{\rm
i}\rho_{1}(\omega-\bm{u}_{1}\cdot\bm{\xi})\bm{1}^{\rm T}&-\mbox{\rm
i}\rho_{2}(\omega-\bm{u}_{2}\cdot\bm{\xi})\bm{1}^{\rm T}\\\ -\mbox{\rm
i}H_{1}(\omega-\bm{u}_{1}\cdot\bm{\xi})\bm{1}&(H_{1}|\bm{\xi}|)^{2}A_{1,0}&O\\\
\mbox{\rm
i}H_{2}(\omega-\bm{u}_{2}\cdot\bm{\xi})\bm{1}&O&(H_{2}|\bm{\xi}|)^{2}A_{2,0}\end{pmatrix}=0,$

where $\bm{\xi}\in\mathbf{R}^{n}$ is the wave vector and $\omega\in\mathbf{C}$
the angular frequency. The left-hand side can be expanded as

LHS $\displaystyle=\det\begin{pmatrix}a&\mbox{\rm
i}\rho_{1}(\omega-\bm{u}_{1}\cdot\bm{\xi})\bm{1}^{\rm T}&-\mbox{\rm
i}\rho_{2}(\omega-\bm{u}_{2}\cdot\bm{\xi})\bm{1}^{\rm T}\\\
\bm{0}&(H_{1}|\bm{\xi}|)^{2}A_{1,0}&O\\\
\bm{0}&O&(H_{2}|\bm{\xi}|)^{2}A_{2,0}\end{pmatrix}$
$\displaystyle\quad+\det\begin{pmatrix}0&\mbox{\rm
i}\rho_{1}(\omega-\bm{u}_{1}\cdot\bm{\xi})\bm{1}^{\rm T}&-\mbox{\rm
i}\rho_{2}(\omega-\bm{u}_{2}\cdot\bm{\xi})\bm{1}^{\rm T}\\\ -\mbox{\rm
i}H_{1}(\omega-\bm{u}_{1}\cdot\bm{\xi})\bm{1}&(H_{1}|\bm{\xi}|)^{2}A_{1,0}&O\\\
\mbox{\rm
i}H_{2}(\omega-\bm{u}_{2}\cdot\bm{\xi})\bm{1}&O&(H_{2}|\bm{\xi}|)^{2}A_{2,0}\end{pmatrix}$
$\displaystyle=a\det\left((H_{1}|\bm{\xi}|)^{2}A_{1,0}\right)\det\left((H_{2}|\bm{\xi}|)^{2}A_{2,0}\right)$
$\displaystyle\quad+\det\begin{pmatrix}0&\mbox{\rm
i}\rho_{1}(\omega-\bm{u}_{1}\cdot\bm{\xi})\bm{1}^{\rm T}\\\ -\mbox{\rm
i}H_{1}(\omega-\bm{u}_{1}\cdot\bm{\xi})\bm{1}&(H_{1}|\bm{\xi}|)^{2}A_{1,0}\end{pmatrix}\det\left((H_{2}|\bm{\xi}|)^{2}A_{2,0}\right)$
$\displaystyle\quad+\det\begin{pmatrix}0&-\mbox{\rm
i}\rho_{2}(\omega-\bm{u}_{2}\cdot\xi)\bm{1}^{\rm T}\\\ \mbox{\rm
i}H_{2}(\omega-\bm{u}_{2}\cdot\bm{\xi})\bm{1}&(H_{2}|\bm{\xi}|)^{2}A_{2,0}\end{pmatrix}\det\left((H_{1}|\bm{\xi}|)^{2}A_{1,0}\right)$
$\displaystyle=H_{1}^{2N+1}H_{2}^{2N^{*}+1}|\bm{\xi}|^{2(N+N^{*}+1)}\left\\{aH_{1}H_{2}|\bm{\xi}|^{2}\det
A_{1,0}\det A_{2,0}\right.$
$\displaystyle\qquad\left.-\rho_{1}H_{2}(\omega-\bm{u}_{1}\cdot\bm{\xi})^{2}\det\tilde{A}_{1,0}\det
A_{2,0}-\rho_{2}H_{1}(\omega-\bm{u}_{2}\cdot\bm{\xi})^{2}\det\tilde{A}_{2,0}\det
A_{1,0}\right\\},$

so that the linear dispersion relation is given simply as

(4.6)
$\frac{\rho_{1}}{H_{1}\alpha_{1}}(\omega-\bm{u}_{1}\cdot\bm{\xi})^{2}+\frac{\rho_{2}}{H_{2}\alpha_{2}}(\omega-\bm{u}_{2}\cdot\bm{\xi})^{2}-a|\bm{\xi}|^{2}=0,$

where

(4.7) $\alpha_{k}=\frac{\det
A_{k,0}}{\det\tilde{A}_{k,0}},\qquad\tilde{A}_{k,0}=\begin{pmatrix}0&\bm{1}^{\rm
T}\\\ -\bm{1}&A_{k,0}\end{pmatrix}$

for $k=1,2$. The discriminant of this quadratic equation in $\omega$ is

$\displaystyle\left(\frac{\rho_{1}}{H_{1}\alpha_{1}}\bm{u}_{1}\cdot\bm{\xi}+\frac{\rho_{2}}{H_{2}\alpha_{2}}\bm{u}_{2}\cdot\bm{\xi}\right)^{2}$
$\displaystyle\quad-\left(\frac{\rho_{1}}{H_{1}\alpha_{1}}+\frac{\rho_{2}}{H_{2}\alpha_{2}}\right)\left(\frac{\rho_{1}}{H_{1}\alpha_{1}}(\bm{u}_{1}\cdot\bm{\xi})^{2}+\frac{\rho_{2}}{H_{2}\alpha_{2}}(\bm{u}_{2}\cdot\bm{\xi})^{2}-a|\bm{\xi}|^{2}\right)$
$\displaystyle=\left(\frac{\rho_{1}}{H_{1}\alpha_{1}}+\frac{\rho_{2}}{H_{2}\alpha_{2}}\right)\left(a|\bm{\xi}|^{2}-\frac{\rho_{1}\rho_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}\left((\bm{u}_{2}-\bm{u}_{1})\cdot\bm{\xi}\right)^{2}\right).$

Therefore, the solutions $\omega$ to the dispersion relation (4.6) are real
for any wave vector $\bm{\xi}\in\mathbf{R}^{n}$ if and only if

$a-\frac{\rho_{1}\rho_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}|\bm{u}_{2}-\bm{u}_{1}|^{2}\geq
0.$

Otherwise, the roots of the linear dispersion relation (4.6) have the form
$\omega=\omega_{r}(\bm{\xi})\pm\mathrm{i}\omega_{i}(\bm{\xi})$ satisfying
$\omega_{i}(\bm{\xi})\to+\infty$ as $\bm{\xi}=(\bm{u}_{2}-\bm{u}_{1})\xi$ and
$\xi\to+\infty$, which leads to an instability of the problem. These
consideration leads us to the following stability condition

(4.8)
$a-\frac{\rho_{1}\rho_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}|\bm{u}_{2}-\bm{u}_{1}|^{2}\geq
c_{0}>0,$

which is equivalent to

$-\big{(}\partial_{z}(P_{2}^{\rm app}-P_{1}^{\rm
app})\big{)}|_{z=\zeta}-\frac{\rho_{1}\rho_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}\big{(}|\nabla\Phi_{2}^{\rm
app}-\nabla\Phi_{1}^{\rm app}|^{2}\big{)}|_{z=\zeta}\geq c_{0}.$

Here, we note that $\alpha_{1}$ and $\alpha_{2}$ are positive constants
depending only on $N$ and $\\{p_{0},p_{1},\ldots,p_{N^{*}}\\}$ and converge to
$0$ as $N,N^{*}\to\infty$. Therefore, as $N$ and $N^{*}$ go to infinity the
domain of stability diminishes.

##  5 Analysis of the linearized system

In this section, we still analyze the system of linearized equations (4.1)
with freezed coefficients. We first derive an energy estimate for solutions to
the linearized system by defining a suitable energy function, and then
transform the linearized system into a standard symmetric form, for which the
hypersurface $t=0$ in the space-time $\mathbf{R}^{n}\times\mathbf{R}$ is
noncharacteristic. These results motivate the subsequent analysis on the
nonlinear equations.

###  5.1 Energy estimate

With the notation (4.5), the linearized system (4.1) with freezed coefficients
can be written in a symmetric form as

(5.1)
$\mathscr{A}_{1}\partial_{t}\dot{\bm{U}}+\mathscr{A}_{0}\dot{\bm{U}}=\bm{0},$

where $\dot{\bm{U}}=(\dot{\zeta},\dot{\bm{\psi}}_{1},\dot{\bm{\psi}}_{2})^{\rm
T}$ and

$\displaystyle\mathscr{A}_{1}=\begin{pmatrix}0&-\rho_{1}{\bm{1}}^{\rm
T}&\rho_{2}{\bm{1}}^{\rm T}\\\ \rho_{1}{\bm{1}}&O&O\\\
-\rho_{2}{\bm{1}}&O&O\end{pmatrix},$
$\displaystyle\mathscr{A}_{0}=\begin{pmatrix}a&-\rho_{1}\bm{1}^{\rm
T}(\bm{u}_{1}\cdot\nabla)&\rho_{2}\bm{1}^{\rm T}(\bm{u}_{2}\cdot\nabla)\\\
\rho_{1}\bm{1}(\bm{u}_{1}\cdot\nabla)&-\rho_{1}H_{1}A_{1,0}\Delta&O\\\
-\rho_{2}\bm{1}(\bm{u}_{2}\cdot\nabla)&O&-\rho_{2}H_{2}A_{2,0}\Delta\end{pmatrix}.$

We note that $\mathscr{A}_{0}$ is symmetric in $L^{2}(\mathbf{R}^{n})$ whereas
$\mathscr{A}_{1}$ is skew-symmetric. Therefore, by taking $L^{2}$-inner
product of (5.1) with $\partial_{t}\dot{\bm{U}}$ we have

$\frac{\rm d}{{\rm
d}t}(\dot{{\bm{U}}},\mathscr{A}_{0}\dot{{\bm{U}}})_{L^{2}}=0$

for any regular solution $\dot{\bm{U}}$ to (5.1), so that
$(\dot{{\bm{U}}},\mathscr{A}_{0}\dot{{\bm{U}}})_{L^{2}}$ would give a
mathematical energy function to the linearized system (5.1) if we show the
positivity of the symmetric operator $\mathscr{A}_{0}$ in
$L^{2}(\mathbf{R}^{n})$. We proceed to check the positivity. For simplicity,
we consider first the case $N=N^{*}=0$ so that $A_{1,0}=A_{2,0}=1$. Then, we
see that

$(\dot{\bm{U}},\mathscr{A}_{0}\dot{\bm{U}})_{L^{2}}=\int_{\mathbf{R}^{n}}\begin{pmatrix}\dot{\zeta}\\\
\nabla\dot{\phi}_{1,0}\\\
\nabla\dot{\phi}_{2,0}\end{pmatrix}\cdot\begin{pmatrix}a&-\rho_{1}\bm{u}_{1}^{\rm
T}&\rho_{2}\bm{u}_{2}^{\rm T}\\\ -\rho_{1}\bm{u}_{1}&\rho_{1}H_{1}\mbox{\rm
Id}&O\\\ \rho_{2}\bm{u}_{2}&O&\rho_{2}H_{2}\mbox{\rm
Id}\end{pmatrix}\begin{pmatrix}\dot{\zeta}\\\ \nabla\dot{\phi}_{1,0}\\\
\nabla\dot{\phi}_{2,0}\end{pmatrix}{\rm d}\bm{x}.$

Therefore, it is sufficient to analyze the positivity of this
$(2n+1)\times(2n+1)$ matrix. The characteristic polynomial of this matrix is
given by

$\displaystyle 0$
$\displaystyle=\det\begin{pmatrix}\lambda-a&\rho_{1}\bm{u}_{1}^{\rm
T}&-\rho_{2}\bm{u}_{2}^{\rm T}\\\
\rho_{1}\bm{u}_{1}&(\lambda-\rho_{1}H_{1})\mbox{\rm Id}&O\\\
-\rho_{2}\bm{u}_{2}&O&(\lambda-\rho_{2}H_{2})\mbox{\rm Id}\end{pmatrix}$
$\displaystyle=(\lambda-a)(\lambda-\rho_{1}H_{1})^{n}(\lambda-\rho_{2}H_{2})^{n}$
$\displaystyle\quad-\rho_{1}^{2}|\bm{u}_{1}|^{2}(\lambda-\rho_{1}H_{1})^{n-1}(\lambda-\rho_{2}H_{2})^{n}-\rho_{2}^{2}|\bm{u}_{2}|^{2}(\lambda-\rho_{1}H_{1})^{n}(\lambda-\rho_{2}H_{2})^{n-1}$
$\displaystyle=(\lambda-\rho_{1}H_{1})^{n-1}(\lambda-\rho_{2}H_{2})^{n-1}\left\\{(\lambda-a)(\lambda-\rho_{1}H_{1})(\lambda-\rho_{2}H_{2})\right.$
$\displaystyle\quad\left.-\rho_{1}^{2}|\bm{u}_{1}|^{2}(\lambda-\rho_{2}H_{2})-\rho_{2}^{2}|\bm{u}_{2}|^{2}(\lambda-\rho_{1}H_{1})\right\\}.$

Therefore, the eigenvalues of the matrix are $\rho_{1}H_{1}$ and
$\rho_{2}H_{2}$ of multiplicity $n-1$ and
$\lambda_{1},\lambda_{2},\lambda_{3}$, which are the roots of the polynomial

$(\lambda-a)(\lambda-\rho_{1}H_{1})(\lambda-\rho_{2}H_{2})-\rho_{1}^{2}|\bm{u}_{1}|^{2}(\lambda-\rho_{2}H_{2})-\rho_{2}^{2}|\bm{u}_{2}|^{2}(\lambda-\rho_{1}H_{1})=0.$

Here, we see that

$\lambda_{1}\lambda_{2}\lambda_{3}=\rho_{1}\rho_{2}(aH_{1}H_{2}-\rho_{1}H_{2}|{\bm{u}}_{1}|^{2}-\rho_{2}H_{1}|{\bm{u}}_{2}|^{2}),$

which is not necessarily positive even if ${\bm{u}}_{1}={\bm{u}}_{2}$.
Therefore, for the positivity of the symmetric operator $\mathscr{A}_{0}$ we
need a smallness of the horizontal velocities ${\bm{u}}_{1}$ and
${\bm{u}}_{2}$. Such a condition is, of course, stronger restriction than the
stability condition (4.8). This means that
$(\dot{{\bm{U}}},\mathscr{A}_{0}\dot{{\bm{U}}})_{L^{2}}$ is not an optimal
energy function and we proceed to find out another one.

We are now considering the linearized system (5.1) with freezed coefficients,
that is,

(5.2)
$\begin{cases}H_{1}\bm{1}(\partial_{t}\dot{\zeta}+\bm{u}_{1}\cdot\nabla\dot{\zeta})-H_{1}^{2}A_{1,0}\Delta\dot{\bm{\psi}}_{1}=\bm{0},\\\
H_{2}\bm{1}(\partial_{t}\dot{\zeta}+\bm{u}_{2}\cdot\nabla\dot{\zeta})+H_{2}^{2}A_{2,0}\Delta\dot{\bm{\psi}}_{2}=\bm{0},\\\
\rho_{1}\bm{1}\cdot\left(\partial_{t}\dot{\bm{\psi}}_{1}+(\bm{u}_{1}\cdot\nabla)\dot{\bm{\psi}}_{1}\right)-\rho_{2}\bm{1}\cdot\left(\partial_{t}\dot{\bm{\psi}}_{2}+(\bm{u}_{2}\cdot\nabla)\dot{\bm{\psi}}_{2}\right)-a\dot{\zeta}=0.\end{cases}$

Applying $\Delta$ to the last equation in (5.2) we have

(5.3)
$\rho_{1}(A_{1,0})^{-1}\bm{1}\cdot(\partial_{t}+\bm{u}_{1}\cdot\nabla)A_{1,0}\Delta\dot{\bm{\psi}}_{1}-\rho_{2}(A_{2,0})^{-1}\bm{1}\cdot(\partial_{t}+\bm{u}_{2}\cdot\nabla)A_{2,0}\Delta\dot{\bm{\psi}}_{2}-a\Delta\dot{\zeta}=0.$

Plugging the first and the second equations in (5.2) into (5.3) to remove
$\dot{\bm{\psi}}_{1}$ and $\dot{\bm{\psi}}_{2}$, we obtain

$\left(\frac{\rho_{1}(A_{1,0})^{-1}\bm{1}\cdot\bm{1}}{H_{1}}(\partial_{t}+\bm{u}_{1}\cdot\nabla)^{2}+\frac{\rho_{2}(A_{2,0})^{-1}\bm{1}\cdot\bm{1}}{H_{2}}(\partial_{t}+{\bm{u}}_{2}\cdot\nabla)^{2}\right)\dot{\zeta}-a\Delta\dot{\zeta}=0.$

In view of the relation following from Cramer’s rule

$(A_{k,0})^{-1}\bm{1}\cdot\bm{1}=\frac{\det\tilde{A}_{k,0}}{\det
A_{k,0}}=\frac{1}{\alpha_{k}}$

for $k=1,2$, the above equation for $\dot{\zeta}$ can be written as

(5.4)
$\left(\frac{\rho_{1}}{H_{1}\alpha_{1}}+\frac{\rho_{2}}{H_{2}\alpha_{2}}\right)(\partial_{t}+\bm{u}\cdot\nabla)^{2}\dot{\zeta}-\left(a\Delta-\frac{\rho_{1}\rho_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}((\bm{u}_{2}-\bm{u}_{1})\cdot\nabla)^{2}\right)\dot{\zeta}=0,$

where $\bm{u}$ is an averaged horizontal velocity on the interface defined by

(5.5)
$\bm{u}=\frac{\rho_{1}H_{2}\alpha_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}\bm{u}_{1}+\frac{\rho_{2}H_{1}\alpha_{1}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}\bm{u}_{2}.$

Taking (5.4) into account, we consider the following constant coefficient
second order partial differential equation

(5.6)
$c_{1}(\partial_{t}+\bm{u}\cdot\nabla)^{2}\dot{\zeta}-\left(c_{2}\Delta-(\bm{v}\cdot\nabla)^{2}\right)\dot{\zeta}=0,$

where $c_{1}$ and $c_{2}$ are positive constants. By taking $L^{2}$-inner
product of (5.6) with $(\partial_{t}+{\bm{u}}\cdot\nabla)\dot{\zeta}$ and
using integration by parts, we see that

$\frac{\rm d}{{\rm
d}t}\left(c_{1}\|\partial_{t}\dot{\zeta}+{\bm{u}}\cdot\nabla\dot{\zeta}\|_{L^{2}}^{2}+c_{2}\|\nabla\dot{\zeta}\|_{L^{2}}^{2}-\|{\bm{v}}\cdot\nabla\dot{\zeta}\|_{L^{2}}^{2}\right)=0$

for any regular solution $\dot{\zeta}$ to (5.6). Here, we have

$c_{2}\|\nabla\dot{\zeta}\|_{L^{2}}^{2}-\|{\bm{v}}\cdot\nabla\dot{\zeta}\|_{L^{2}}^{2}=(\nabla\dot{\zeta},(c_{2}{\rm
Id}-{\bm{v}}\otimes{\bm{v}})\nabla\dot{\zeta})_{L^{2}}.$

The matrix $c_{2}{\rm Id}-{\bm{v}}\otimes{\bm{v}}$ is positive if and only if
$c_{2}-|{\bm{v}}|^{2}>0$. Under this assumption, we obtain an energy estimate
for the solutions to (5.6). Applying this consideration to (5.4), we see that
the positivity condition is exactly the same as the stability condition (4.8),
under which we can obtain an energy estimate for (5.4).

In view of this fact we rewrite the linearized system (5.1) with freezed
coefficients in the form

$\mathscr{A}_{1}(\partial_{t}+\bm{u}\cdot\nabla)\dot{\bm{U}}+\mathscr{A}_{0}^{\rm
mod}\dot{\bm{U}}=\bm{0},$

where

$\displaystyle\mathscr{A}_{0}^{\rm
mod}=\mathscr{A}_{0}-\mathscr{A}_{1}({\bm{u}}\cdot\nabla)$
$\displaystyle=\begin{pmatrix}a&\frac{\rho_{1}\rho_{2}H_{1}\alpha_{1}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}\bm{1}^{\rm
T}(\bm{v}\cdot\nabla)&\frac{\rho_{1}\rho_{2}H_{2}\alpha_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}\bm{1}^{\rm
T}(\bm{v}\cdot\nabla)\\\\[4.30554pt]
-\frac{\rho_{1}\rho_{2}H_{1}\alpha_{1}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}\bm{1}(\bm{v}\cdot\nabla)&-\rho_{1}H_{1}A_{1,0}\Delta&O\\\\[4.30554pt]
-\frac{\rho_{1}\rho_{2}H_{2}\alpha_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}\bm{1}(\bm{v}\cdot\nabla)&O&-\rho_{2}H_{2}A_{2,0}\Delta\end{pmatrix}$

and $\bm{v}=\bm{u}_{2}-\bm{u}_{1}$. By taking $L^{2}$-inner product of this
equation with $(\partial_{t}+{\bm{u}}\cdot\nabla)\dot{\bm{U}}$ and using
integration by parts, we see that

$\frac{\rm d}{{\rm d}t}(\mathscr{A}_{0}^{\rm
mod}\dot{\bm{U}},\dot{\bm{U}})_{L^{2}}=0$

for any regular solution to (5.1). We proceed to check the positivity of the
symmetric operator $\mathscr{A}_{0}^{\rm mod}$ in $L^{2}(\mathbf{R}^{n})$
under the stability condition (4.8). We see that

$\displaystyle(\mathscr{A}_{0}^{\rm mod}\dot{\bm{U}},\dot{\bm{U}})_{L^{2}}$
$\displaystyle=(a\dot{\zeta},\dot{\zeta})_{L^{2}}+\sum_{l=1}^{n}\left\\{(\rho_{1}H_{1}A_{1,0}\partial_{l}\dot{\bm{\psi}}_{1},\partial_{l}\dot{\bm{\psi}}_{1})_{L^{2}}+(\rho_{2}H_{2}A_{2,0}\partial_{l}\dot{\bm{\psi}}_{2},\partial_{l}\dot{\bm{\psi}}_{2})_{L^{2}}\right\\}$
$\displaystyle\quad+2(\frac{\rho_{1}\rho_{2}H_{1}\alpha_{1}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}(\bm{v}\cdot\nabla)(\bm{1}\cdot\dot{\bm{\psi}}_{1}),\dot{\zeta})_{L^{2}}$
$\displaystyle\quad+2(\frac{\rho_{1}\rho_{2}H_{2}\alpha_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}(\bm{v}\cdot\nabla)(\bm{1}\cdot\dot{\bm{\psi}}_{2}),\dot{\zeta})_{L^{2}}.$

On the other hand, the matrix $\tilde{A}_{k,0}$ is nonsingular and its inverse
matrix can be written as

$(\tilde{A}_{k,0})^{-1}=\begin{pmatrix}0&\bm{1}^{\rm T}\\\
-\bm{1}&A_{k,0}\end{pmatrix}^{-1}=\begin{pmatrix}q_{k,0}&(\bm{q}_{k,0})^{\rm
T}\\\ -\bm{q}_{k,0}&Q_{k,0}\end{pmatrix}$

with a symmetric matrix $Q_{k,0}$ for $k=1,2$. Moreover, $q_{k,0}=\frac{\det
A_{k,0}}{\det\tilde{A}_{k,0}}=\alpha_{k}$ is positive and $Q_{k,0}$ is
nonnegative. In fact, for any $\bm{\psi}$ putting $\zeta$ and $\bm{\phi}$ by
$\begin{pmatrix}\zeta\\\
\bm{\phi}\end{pmatrix}=(\tilde{A}_{k,0})^{-1}\begin{pmatrix}0\\\
\bm{\psi}\end{pmatrix}$ we have

$Q_{k,0}\bm{\psi}\cdot\bm{\psi}=\begin{pmatrix}q_{k,0}&(\bm{q}_{k,0})^{\rm
T}\\\ -\bm{q}_{k,0}&Q_{k,0}\end{pmatrix}\begin{pmatrix}0\\\
\bm{\psi}\end{pmatrix}\cdot\begin{pmatrix}0\\\
\bm{\psi}\end{pmatrix}=\begin{pmatrix}\zeta\\\
\bm{\phi}\end{pmatrix}\cdot\tilde{A}_{k,0}\begin{pmatrix}\zeta\\\
\bm{\phi}\end{pmatrix}=\bm{\phi}\cdot A_{k,0}\bm{\phi}\geq 0.$

We note that $Q_{k,0}$ is not positive because it has a zero eigenvalue with
an eigenvector $\bm{1}$. Now, for any $\bm{\phi}$, putting
$\eta=\bm{1}\cdot\bm{\phi}$ and $\bm{\psi}=A_{k,0}\bm{\phi}$, we have
$\tilde{A}_{k,0}\begin{pmatrix}0\\\
\bm{\phi}\end{pmatrix}=\begin{pmatrix}\eta\\\ \bm{\psi}\end{pmatrix}$ so that

$A_{k,0}\bm{\phi}\cdot\bm{\phi}=\tilde{A}_{k,0}\begin{pmatrix}0\\\
\bm{\phi}\end{pmatrix}\cdot\begin{pmatrix}0\\\
\bm{\phi}\end{pmatrix}=\begin{pmatrix}\eta\\\
\bm{\psi}\end{pmatrix}\cdot(\tilde{A}_{k,0})^{-1}\begin{pmatrix}\eta\\\
\bm{\psi}\end{pmatrix}=q_{k,0}\eta^{2}+Q_{k,0}\bm{\psi}\cdot\bm{\psi},$

from which we deduce the identity

(5.7)
$A_{k,0}\bm{\phi}\cdot\bm{\phi}=\alpha_{k}(\bm{1}\cdot\bm{\phi})^{2}+Q_{k,0}A_{k,0}\bm{\phi}\cdot
A_{k,0}\bm{\phi}.$

By using the decomposition (5.7) we see that

$\displaystyle(\mathscr{A}_{0}^{\rm mod}\dot{\bm{U}},\dot{\bm{U}})_{L^{2}}$
$\displaystyle=\sum_{l=1}^{n}\left\\{(\rho_{1}H_{1}Q_{1,0}A_{1,0}\partial_{l}\dot{\bm{\psi}}_{1},A_{1,0}\partial_{l}\dot{\bm{\psi}}_{1})_{L^{2}}+(\rho_{2}H_{2}Q_{2,0}A_{2,0}\partial_{l}\dot{\bm{\psi}}_{2},A_{2,0}\partial_{l}\dot{\bm{\psi}}_{2})_{L^{2}}\right\\}$
$\displaystyle\quad+\biggl{\\{}(a\dot{\zeta},\dot{\zeta})_{L^{2}}+(\rho_{1}H_{1}\alpha_{1}\nabla(\bm{1}\cdot\dot{\bm{\psi}}_{1}),\nabla(\bm{1}\cdot\dot{\bm{\psi}}_{1}))_{L^{2}}+(\rho_{2}H_{2}\alpha_{2}\nabla(\bm{1}\cdot\dot{\bm{\psi}}_{2}),\nabla(\bm{1}\cdot\dot{\bm{\psi}}_{2}))_{L^{2}}$
$\displaystyle\qquad+(\frac{2\rho_{1}\rho_{2}H_{1}\alpha_{1}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}(\bm{v}\cdot\nabla)(\bm{1}\cdot\dot{\bm{\psi}}_{1}),\dot{\zeta})_{L^{2}}+(\frac{2\rho_{1}\rho_{2}H_{2}\alpha_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}(\bm{v}\cdot\nabla)(\bm{1}\cdot\dot{\bm{\psi}}_{2}),\dot{\zeta})_{L^{2}}\biggr{\\}}$
$\displaystyle=:I_{1}+I_{2}.$

Here, $I_{1}\geq 0$ since $Q_{1,0}$ and $Q_{2,0}$ are nonnegative, and

$\displaystyle I_{2}$
$\displaystyle\geq\int_{\mathbf{R}^{n}}\biggl{\\{}a\dot{\zeta}^{2}+\rho_{1}H_{1}\alpha_{1}|\nabla(\bm{1}\cdot\dot{\bm{\psi}}_{1})|^{2}+\rho_{2}H_{2}\alpha_{2}|\nabla(\bm{1}\cdot\dot{\bm{\psi}}_{2})|^{2}$
$\displaystyle\makebox[40.00006pt]{}-\frac{2\rho_{1}\rho_{2}|\bm{v}|}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}\left(H_{1}\alpha_{1}|\nabla(\bm{1}\cdot\dot{\bm{\psi}}_{1})|+H_{2}\alpha_{2}|\nabla(\bm{1}\cdot\dot{\bm{\psi}}_{2})|\right)|\dot{\zeta}|\biggr{\\}}{\rm
d}\bm{x},$

so that it is sufficient to show the positivity of the matrix

$\mathfrak{A}_{0}:=\begin{pmatrix}a&-\frac{\rho_{1}\rho_{2}H_{1}\alpha_{1}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}|\bm{v}|&-\frac{\rho_{1}\rho_{2}H_{2}\alpha_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}|\bm{v}|\\\\[4.30554pt]
-\frac{\rho_{1}\rho_{2}H_{1}\alpha_{1}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}|\bm{v}|&\rho_{1}H_{1}\alpha_{1}&0\\\\[4.30554pt]
-\frac{\rho_{1}\rho_{2}H_{2}\alpha_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}|\bm{v}|&0&\rho_{2}H_{2}\alpha_{2}\end{pmatrix}.$

From Sylvester’s criterion and since $\rho_{k}H_{k}\alpha_{k}$ is positive for
$k=1,2$, the positivity of the matrix $\mathfrak{A}_{0}$ is equivalent to

$\displaystyle\det\mathfrak{A}_{0}$
$\displaystyle=a(\rho_{1}H_{1}\alpha_{1})(\rho_{2}H_{2}\alpha_{2})$
$\displaystyle\quad\;-\rho_{1}H_{1}\alpha_{1}\left(\frac{\rho_{1}\rho_{2}H_{2}\alpha_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}|\bm{v}|\right)^{2}-\rho_{2}H_{2}\alpha_{2}\left(\frac{\rho_{1}\rho_{2}H_{1}\alpha_{1}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}|\bm{v}|\right)^{2}$
$\displaystyle=(\rho_{1}H_{1}\alpha_{1})(\rho_{2}H_{2}\alpha_{2})\left(a-\frac{\rho_{1}\rho_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}}|\bm{v}|^{2}\right)>0.$

Since $\bm{v}=\bm{u}_{2}-\bm{u}_{1}$, under the stability condition (4.8) we
have the positivity of $\mathfrak{A}_{0}$, so that in view of (5.7) and the
positivity of the matrix $A_{k,0}$ for $k=1,2$ we finally obtain the
equivalence

$(\mathscr{A}_{0}^{\rm
mod}\dot{\bm{U}},\dot{\bm{U}})_{L^{2}}\simeq\|\dot{\zeta}\|_{L^{2}}^{2}+\|\nabla\dot{\bm{\phi}}_{1}\|_{L^{2}}^{2}+\|\nabla\dot{\bm{\phi}}_{2}\|_{L^{2}}^{2}.$

Therefore, $(\mathscr{A}_{0}^{\rm mod}\dot{\bm{U}},\dot{\bm{U}})_{L^{2}}$
would provide a useful mathematical energy function.

###  5.2 Symmetrization of the linearized equations

We still consider the linearized equations (4.1) with freezed coefficients.
However, for later use we define $\dot{\bm{\phi}}_{1}$ and
$\dot{\bm{\phi}}_{2}$ in place of (4.5) by

$\begin{cases}\dot{\bm{\phi}}_{1}=(\dot{\phi}_{1,0},\dot{\phi}_{1,1},\ldots,\dot{\phi}_{1,N})^{\rm
T},\\\
\dot{\bm{\phi}}_{2}=(\dot{\phi}_{2,0},\dot{\phi}_{2,1},\ldots,\dot{\phi}_{2,N^{*}})^{\rm
T}.\end{cases}$

Then, the linearized equations have the form

(5.8)
$\begin{cases}\bm{l}_{1}(H_{1})(\partial_{t}\dot{\zeta}+\bm{u}_{1}\cdot\nabla\dot{\zeta})-A_{1}(H_{1})\Delta\dot{\bm{\phi}}_{1}=\bm{0},\\\
-\bm{l}_{2}(H_{2})(\partial_{t}\dot{\zeta}+\bm{u}_{2}\cdot\nabla\dot{\zeta})-A_{2}(H_{2})\Delta\dot{\bm{\phi}}_{2}=\bm{0},\\\
-\rho_{1}\bm{l}_{1}(H_{1})\cdot(\partial_{t}\dot{\bm{\phi}}_{1}+(\bm{u}_{1}\cdot\nabla)\dot{\bm{\phi}}_{1})+\rho_{2}\bm{l}_{2}(H_{2})\cdot(\partial_{t}\dot{\bm{\phi}}_{2}+(\bm{u}_{2}\cdot\nabla)\dot{\bm{\phi}}_{2})+a\dot{\zeta}=0,\end{cases}$

where

(5.9) ${\bm{l}}_{1}(H_{1})=(1,H_{1}^{2},H_{1}^{4},\ldots,H_{1}^{2N})^{\rm
T},\quad{\bm{l}}_{2}(H_{2})=(1,H_{2}^{p_{1}},H_{2}^{p_{2}},\ldots,H_{2}^{p_{N^{*}}})^{\rm
T},$

and

(5.10) $\begin{cases}\displaystyle
A_{1}(H_{1})=\biggl{(}\frac{1}{2(i+j)+1}H_{1}^{2(i+j)+1}\biggr{)}_{0\leq
i,j\leq N},\\\ \displaystyle
A_{2}(H_{2})=\biggl{(}\frac{1}{p_{i}+p_{j}+1}H_{2}^{p_{i}+p_{j}+1}\biggr{)}_{0\leq
i,j\leq N^{*}}.\end{cases}$

In the following, we abbreviate simply ${\bm{l}}_{k}(H_{k})$ and
$A_{k}(H_{k})$ as ${\bm{l}}_{k}$ and $A_{k}$ for $k=1,2$. We are going to show
that the system can be transformed into a positive symmetric system of the
form

(5.11) $\mathscr{A}_{0}^{\rm
mod}\partial_{t}\dot{\bm{U}}+\mathscr{A}\dot{\bm{U}}=\bm{0},$

where $\dot{\bm{U}}=(\dot{\zeta},\dot{\bm{\phi}}_{1},\dot{\bm{\phi}}_{2})^{\rm
T}$, $\mathscr{A}_{0}^{\rm mod}$ is the positive operator defined in the
previous section with slight modification, and $\mathscr{A}$ is a skew-
symmetric operator in $L^{2}(\mathbf{R}^{n})$. As before, we put
$\bm{v}=\bm{u}_{2}-\bm{u}_{1}$ and define $\bm{u}$ by (5.5). Furthermore, we
introduce the notation

(5.12)
$\theta_{1}=\frac{\rho_{2}H_{1}\alpha_{1}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}},\qquad\theta_{2}=\frac{\rho_{1}H_{2}\alpha_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}},$

where $\alpha_{1}$ and $\alpha_{2}$ are positive constants defined by (4.7).
Then, we have $\bm{u}=\theta_{2}\bm{u}_{1}+\theta_{1}\bm{u}_{2}$ and
$\theta_{1}+\theta_{2}=1$. We can also express $\bm{u}_{1}$ and $\bm{u}_{2}$
in terms of $\bm{u}$ and $\bm{v}$ as

$\bm{u}_{1}=\bm{u}-\theta_{1}\bm{v},\qquad\bm{u}_{2}=\bm{u}+\theta_{2}\bm{v}.$

Applying $\Delta$ to the third equation in (5.8) and differentiating the first
and the second equations with respect to $t$, we obtain

$\displaystyle\begin{pmatrix}0&-\rho_{1}\bm{l}_{1}^{\rm
T}&\rho_{2}\bm{l}_{2}^{\rm T}\\\ -\rho_{1}\bm{l}_{1}&\rho_{1}A_{1}&O\\\
\rho_{2}\bm{l}_{2}&O&\rho_{2}A_{2}\end{pmatrix}\begin{pmatrix}\partial_{t}^{2}\dot{\zeta}\\\
\Delta\partial_{t}\dot{\bm{\phi}}_{1}\\\
\Delta\partial_{t}\dot{\bm{\phi}}_{2}\end{pmatrix}+\begin{pmatrix}0\\\
-\rho_{1}\bm{l}_{1}(\bm{u}_{1}\cdot\nabla)\\\
\rho_{2}\bm{l}_{2}(\bm{u}_{2}\cdot\nabla)\end{pmatrix}\partial_{t}\dot{\zeta}$
$\displaystyle\quad+\begin{pmatrix}a&-\rho_{1}\bm{l}_{1}^{\rm
T}(\bm{u}_{1}\cdot\nabla)&\rho_{2}\bm{l}_{2}^{\rm T}(\bm{u}_{2}\cdot\nabla)\\\
\bm{0}&O&O\\\ \bm{0}&O&O\end{pmatrix}\Delta\dot{\bm{U}}=\bm{0}.$

In view of this, we introduce a symmetric matrix as

(5.13) $\begin{pmatrix}q_{0}&\bm{q}_{1}^{\rm T}&\bm{q}_{2}^{\rm T}\\\
\bm{q}_{1}&Q_{11}&Q_{12}\\\
\bm{q}_{2}&Q_{21}&Q_{22}\end{pmatrix}=\begin{pmatrix}0&-\rho_{1}\bm{l}_{1}^{\rm
T}&\rho_{2}\bm{l}_{2}^{\rm T}\\\ -\rho_{1}\bm{l}_{1}&\rho_{1}A_{1}&O\\\
\rho_{2}\bm{l}_{2}&O&\rho_{2}A_{2}\end{pmatrix}^{-1},$

where $Q_{11}^{\rm T}=Q_{11}$, $Q_{22}^{\rm T}=Q_{22}$, and $Q_{12}^{\rm
T}=Q_{21}$. Moreover, we have

$\begin{cases}-\rho_{1}\bm{l}_{1}\cdot\bm{q}_{1}+\rho_{2}\bm{l}_{2}\cdot\bm{q}_{2}=1,\quad
A_{1}\bm{q}_{1}=q_{0}\bm{l}_{1},\quad A_{2}\bm{q}_{2}=-q_{0}\bm{l}_{2},\\\
\rho_{1}A_{1}Q_{11}=\mbox{\rm Id}+\rho_{1}\bm{l}_{1}\bm{q}_{1}^{\rm
T},\quad\rho_{2}A_{2}Q_{22}=\mbox{\rm Id}-\rho_{2}\bm{l}_{2}\bm{q}_{2}^{\rm
T},\\\ A_{1}Q_{12}=\bm{l}_{1}\bm{q}_{2}^{\rm T},\quad
A_{2}Q_{21}=-\bm{l}_{2}\bm{q}_{1}^{\rm T}\end{cases}$

and by Cramer’s rule,

$q_{0}=-\frac{H_{1}H_{2}\alpha_{1}\alpha_{2}}{\rho_{1}H_{2}\alpha_{2}+\rho_{2}H_{1}\alpha_{1}},\quad\bm{l}_{1}\cdot\bm{q}_{1}=\frac{-q_{0}}{H_{1}\alpha_{1}}=-\frac{\theta_{2}}{\rho_{1}},\quad\bm{l}_{2}\cdot\bm{q}_{2}=\frac{q_{0}}{H_{2}\alpha_{2}}=\frac{\theta_{1}}{\rho_{2}}.$

Using these notations, we have

$\displaystyle\begin{pmatrix}-\rho_{1}A_{1}\Delta\partial_{t}\dot{\bm{\phi}}_{1}\\\
-\rho_{2}A_{2}\Delta\partial_{t}\dot{\bm{\phi}}_{2}\end{pmatrix}+\begin{pmatrix}-\rho_{1}A_{1}&O\\\
O&-\rho_{2}A_{2}\end{pmatrix}\begin{pmatrix}\bm{q}_{1}&Q_{11}&Q_{12}\\\
\bm{q}_{2}&Q_{21}&Q_{22}\end{pmatrix}\begin{pmatrix}0\\\
-\rho_{1}\bm{l}_{1}(\bm{u}_{1}\cdot\nabla)\\\
\rho_{2}\bm{l}_{2}(\bm{u}_{2}\cdot\nabla)\end{pmatrix}\partial_{t}\dot{\zeta}$
$\displaystyle\quad+\begin{pmatrix}-\rho_{1}A_{1}&O\\\
O&-\rho_{2}A_{2}\end{pmatrix}\begin{pmatrix}\bm{q}_{1}&Q_{11}&Q_{12}\\\
\bm{q}_{2}&Q_{21}&Q_{22}\end{pmatrix}\begin{pmatrix}a&-\rho_{1}\bm{l}_{1}^{\rm
T}(\bm{u}_{1}\cdot\nabla)&\rho_{2}\bm{l}_{2}^{\rm T}(\bm{u}_{2}\cdot\nabla)\\\
\bm{0}&O&O\\\ \bm{0}&O&O\end{pmatrix}\Delta\dot{\bm{U}}=\bm{0}.$

Here, we see that

$\displaystyle\begin{pmatrix}-\rho_{1}A_{1}&O\\\
O&-\rho_{2}A_{2}\end{pmatrix}\begin{pmatrix}\bm{q}_{1}&Q_{11}&Q_{12}\\\
\bm{q}_{2}&Q_{21}&Q_{22}\end{pmatrix}\begin{pmatrix}0\\\
-\rho_{1}\bm{l}_{1}(\bm{u}_{1}\cdot\nabla)\\\
\rho_{2}\bm{l}_{2}(\bm{u}_{2}\cdot\nabla)\end{pmatrix}=-\begin{pmatrix}\theta_{1}\rho_{1}\bm{l}_{1}\\\
\theta_{2}\rho_{2}\bm{l}_{2}\end{pmatrix}(\bm{v}\cdot\nabla),$
$\displaystyle\begin{pmatrix}-\rho_{1}A_{1}&O\\\
O&-\rho_{2}A_{2}\end{pmatrix}\begin{pmatrix}\bm{q}_{1}&Q_{11}&Q_{12}\\\
\bm{q}_{2}&Q_{21}&Q_{22}\end{pmatrix}\begin{pmatrix}a&-\rho_{1}\bm{l}_{1}^{\rm
T}(\bm{u}_{1}\cdot\nabla)&\rho_{2}\bm{l}_{2}^{\rm T}(\bm{u}_{2}\cdot\nabla)\\\
\bm{0}&O&O\\\ \bm{0}&O&O\end{pmatrix}$
$\displaystyle\quad=q_{0}\begin{pmatrix}-a\rho_{1}\bm{l}_{1}&\rho_{1}^{2}\bm{l}_{1}\bm{l}_{1}^{\rm
T}(\bm{u}_{1}\cdot\nabla)&-\rho_{1}\rho_{2}\bm{l}_{1}\bm{l}_{2}^{\rm
T}(\bm{u}_{2}\cdot\nabla)\\\
a\rho_{2}\bm{l}_{2}&-\rho_{1}\rho_{2}\bm{l}_{2}\bm{l}_{1}^{\rm
T}(\bm{u}_{1}\cdot\nabla)&\rho_{2}^{2}\bm{l}_{2}\bm{l}_{2}^{\rm
T}(\bm{u}_{2}\cdot\nabla)\end{pmatrix},$

so that

(5.14)
$\displaystyle\begin{pmatrix}-\rho_{1}A_{1}\Delta\partial_{t}\dot{\bm{\phi}}_{1}-\theta_{1}\rho_{1}\bm{l}_{1}(\bm{v}\cdot\nabla)\partial_{t}\dot{\zeta}\\\
-\rho_{2}A_{2}\Delta\partial_{t}\dot{{\bm{\phi}}}_{2}-\theta_{2}\rho_{2}\bm{l}_{2}(\bm{v}\cdot\nabla)\partial_{t}\dot{\zeta}\end{pmatrix}$
$\displaystyle=q_{0}a\begin{pmatrix}\rho_{1}\bm{l}_{1}\\\
-\rho_{2}\bm{l}_{2}\end{pmatrix}\Delta\dot{\zeta}+q_{0}\begin{pmatrix}-\rho_{1}^{2}\bm{l}_{1}\bm{l}_{1}^{\rm
T}(\bm{u}_{1}\cdot\nabla)&\rho_{1}\rho_{2}\bm{l}_{1}\bm{l}_{2}^{\rm
T}(\bm{u}_{2}\cdot\nabla)\\\ \rho_{1}\rho_{2}\bm{l}_{2}\bm{l}_{1}^{\rm
T}(\bm{u}_{1}\cdot\nabla)&-\rho_{2}^{2}\bm{l}_{2}\bm{l}_{2}^{\rm
T}(\bm{u}_{2}\cdot\nabla)\end{pmatrix}\Delta\begin{pmatrix}\dot{\bm{\phi}}_{1}\\\
\dot{\bm{\phi}}_{2}\end{pmatrix}.$

On the other hand, taking the Euclidean inner product of the first and the
second equations in (5.8) with $-\rho_{1}\bm{q}_{1}$ and $\rho_{2}\bm{q}_{2}$,
respectively, we obtain

$\begin{cases}\theta_{2}(\partial_{t}\dot{\zeta}+\bm{u}_{1}\cdot\nabla\dot{\zeta})+q_{0}\rho_{1}\bm{l}_{1}\cdot\Delta\dot{\bm{\phi}}_{1}=0,\\\
\theta_{1}(\partial_{t}\dot{\zeta}+\bm{u}_{2}\cdot\nabla\dot{\zeta})-q_{0}\rho_{2}\bm{l}_{2}\cdot\Delta\dot{{\bm{\phi}}}_{2}=0,\end{cases}$

which are equivalent to

(5.15)
$\begin{cases}\partial_{t}\dot{\zeta}+\bm{u}\cdot\nabla\dot{\zeta}+q_{0}\Delta(\rho_{1}\bm{l}_{1}\cdot\dot{\bm{\phi}}_{1}-\rho_{2}\bm{l}_{2}\cdot\dot{\bm{\phi}}_{2})=0,\\\
\theta_{1}\theta_{2}\bm{v}\cdot\nabla\dot{\zeta}-q_{0}\Delta(\theta_{1}\rho_{1}\bm{l}_{1}\cdot\dot{\bm{\phi}}_{1}+\theta_{2}\rho_{2}\bm{l}_{2}\cdot\dot{\bm{\phi}}_{2})=0.\end{cases}$

It follows from the second equation in (5.15) that

$\displaystyle\theta_{1}\rho_{1}\bm{l}_{1}\cdot\partial_{t}\dot{\bm{\phi}}_{1}+\theta_{2}\rho_{2}\bm{l}_{2}\cdot\partial_{t}\dot{\bm{\phi}}_{2}$
$\displaystyle=q_{0}^{-1}\theta_{1}\theta_{2}(\bm{v}\cdot\nabla)\Delta^{-1}\partial_{t}\dot{\zeta}.$

Therefore, we obtain

(5.16) $\displaystyle
a\partial_{t}\dot{\zeta}+(\bm{v}\cdot\nabla)(\theta_{1}\rho_{1}\bm{l}_{1}\cdot\partial_{t}\dot{\bm{\phi}}_{1}+\theta_{2}\rho_{2}\bm{l}_{2}\cdot\partial_{t}\dot{\bm{\phi}}_{2})$
$\displaystyle=-a\bigl{(}(\bm{u}\cdot\nabla)\dot{\zeta}+q_{0}\Delta(\rho_{1}\bm{l}_{1}\cdot\dot{\bm{\phi}}_{1}-\rho_{2}\bm{l}_{2}\cdot\dot{\bm{\phi}}_{2})\bigr{)}$
$\displaystyle\quad-\theta_{1}\theta_{2}(\bm{v}\cdot\nabla)^{2}\bigl{(}q_{0}^{-1}(\bm{u}\cdot\nabla)\Delta^{-1}\dot{\zeta}+(\rho_{1}\bm{l}_{1}\cdot\dot{\bm{\phi}}_{1}-\rho_{2}\bm{l}_{2}\cdot\dot{\bm{\phi}}_{2})\bigr{)}.$

We proceed to symmetrize the second term in the right-hand side of (5.14).

$\displaystyle q_{0}\begin{pmatrix}-\rho_{1}^{2}\bm{l}_{1}\bm{l}_{1}^{\rm
T}(\bm{u}_{1}\cdot\nabla)&\rho_{1}\rho_{2}\bm{l}_{1}\bm{l}_{2}^{\rm
T}(\bm{u}_{2}\cdot\nabla)\\\ \rho_{1}\rho_{2}\bm{l}_{2}\bm{l}_{1}^{\rm
T}(\bm{u}_{1}\cdot\nabla)&-\rho_{2}^{2}\bm{l}_{2}\bm{l}_{2}^{\rm
T}(\bm{u}_{2}\cdot\nabla)\end{pmatrix}\Delta\begin{pmatrix}\dot{\bm{\phi}}_{1}\\\
\dot{\bm{\phi}}_{2}\end{pmatrix}$
$\displaystyle=q_{0}\begin{pmatrix}-\rho_{1}^{2}\bm{l}_{1}\bm{l}_{1}^{\rm
T}&\rho_{1}\rho_{2}\bm{l}_{1}\bm{l}_{2}^{\rm T}\\\
\rho_{1}\rho_{2}\bm{l}_{2}\bm{l}_{1}^{\rm
T}&-\rho_{2}^{2}\bm{l}_{2}\bm{l}_{2}^{\rm
T}\end{pmatrix}(\bm{u}\cdot\nabla)\Delta\begin{pmatrix}\dot{\bm{\phi}}_{1}\\\
\dot{\bm{\phi}}_{2}\end{pmatrix}+q_{0}\begin{pmatrix}\theta_{1}\rho_{1}^{2}\bm{l}_{1}\bm{l}_{1}^{\rm
T}&\theta_{2}\rho_{1}\rho_{2}\bm{l}_{1}\bm{l}_{2}^{\rm T}\\\
-\theta_{1}\rho_{1}\rho_{2}\bm{l}_{2}\bm{l}_{1}^{\rm
T}&-\theta_{2}\rho_{2}^{2}\bm{l}_{2}\bm{l}_{2}^{\rm
T}\end{pmatrix}(\bm{v}\cdot\nabla)\Delta\begin{pmatrix}\dot{\bm{\phi}}_{1}\\\
\dot{\bm{\phi}}_{2}\end{pmatrix},$

where

$\displaystyle
q_{0}\begin{pmatrix}\theta_{1}\rho_{1}^{2}\bm{l}_{1}\bm{l}_{1}^{\rm
T}&\theta_{2}\rho_{1}\rho_{2}\bm{l}_{1}\bm{l}_{2}^{\rm T}\\\
-\theta_{1}\rho_{1}\rho_{2}\bm{l}_{2}\bm{l}_{1}^{\rm
T}&-\theta_{2}\rho_{2}^{2}\bm{l}_{2}\bm{l}_{2}^{\rm
T}\end{pmatrix}\Delta\begin{pmatrix}\dot{\bm{\phi}}_{1}\\\
\dot{\bm{\phi}}_{2}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\rho_{1}\bm{l}_{1}\\\
-\rho_{2}\bm{l}_{2}\end{pmatrix}q_{0}\Delta(\theta_{1}\rho_{1}\bm{l}_{1}\cdot\dot{\bm{\phi}}_{1}+\theta_{2}\rho_{2}\bm{l}_{2}\cdot\dot{\bm{\phi}}_{2})$
$\displaystyle=\theta_{1}\theta_{2}\begin{pmatrix}\rho_{1}\bm{l}_{1}\\\
-\rho_{2}\bm{l}_{2}\end{pmatrix}(\bm{v}\cdot\nabla)\dot{\zeta}.$

In the above calculation, we used the second equation in (5.15). Therefore,

$\displaystyle\begin{pmatrix}-\rho_{1}A_{1}\Delta\partial_{t}\dot{\bm{\phi}}_{1}-\theta_{1}\rho_{1}\bm{l}_{1}(\bm{v}\cdot\nabla)\partial_{t}\dot{\zeta}\\\
-\rho_{2}A_{2}\Delta\partial_{t}\dot{\bm{\phi}}_{2}-\theta_{2}\rho_{2}\bm{l}_{2}(\bm{v}\cdot\nabla)\partial_{t}\dot{\zeta}\end{pmatrix}$
$\displaystyle=q_{0}a\begin{pmatrix}\rho_{1}\bm{l}_{1}\\\
-\rho_{2}\bm{l}_{2}\end{pmatrix}\Delta\dot{\zeta}+\theta_{1}\theta_{2}\begin{pmatrix}\rho_{1}\bm{l}_{1}\\\
-\rho_{2}\bm{l}_{2}\end{pmatrix}(\bm{v}\cdot\nabla)^{2}\dot{\zeta}$
$\displaystyle\quad+q_{0}\begin{pmatrix}-\rho_{1}^{2}\bm{l}_{1}\bm{l}_{1}^{\rm
T}&\rho_{1}\rho_{2}\bm{l}_{1}\bm{l}_{2}^{\rm T}\\\
\rho_{1}\rho_{2}\bm{l}_{2}\bm{l}_{1}^{\rm
T}&-\rho_{2}^{2}\bm{l}_{2}\bm{l}_{2}^{\rm
T}\end{pmatrix}(\bm{u}\cdot\nabla)\Delta\begin{pmatrix}\dot{\bm{\phi}}_{1}\\\
\dot{\bm{\phi}}_{2}\end{pmatrix}.$

Summarizing the above calculations, if we define the symmetrizer
$\mathscr{A}_{0}^{\rm mod}$ by

(5.17) $\mathscr{A}_{0}^{\rm
mod}=\begin{pmatrix}a&\theta_{1}\rho_{1}\bm{l}_{1}^{\rm
T}(\bm{v}\cdot\nabla)&\theta_{2}\rho_{2}\bm{l}_{2}^{\rm
T}(\bm{v}\cdot\nabla)\\\
-\theta_{1}\rho_{1}\bm{l}_{1}(\bm{v}\cdot\nabla)&-\rho_{1}A_{1}\Delta&O\\\
-\theta_{2}\rho_{2}\bm{l}_{2}(\bm{v}\cdot\nabla)&O&-\rho_{2}A_{2}\Delta\end{pmatrix},$

then we obtain

$\displaystyle\mathscr{A}_{0}^{\rm mod}\partial_{t}\dot{\bm{U}}$
$\displaystyle=\begin{pmatrix}a\partial_{t}\dot{\zeta}+(\bm{v}\cdot\nabla)(\theta_{1}\rho_{1}\bm{l}_{1}\cdot\partial_{t}\dot{\bm{\phi}}_{1}+\theta_{2}\rho_{2}\bm{l}_{2}\cdot\partial_{t}\dot{\bm{\phi}}_{2})\\\
-\theta_{1}\rho_{1}\bm{l}_{1}(\bm{v}\cdot\nabla)\partial_{t}\dot{\zeta}-\rho_{1}A_{1}\Delta\partial_{t}\dot{\bm{\phi}}_{1}\\\
-\theta_{2}\rho_{2}\bm{l}_{2}(\bm{v}\cdot\nabla)\partial_{t}\dot{\zeta}-\rho_{2}A_{2}\Delta\partial_{t}\dot{\bm{\phi}}_{2}\end{pmatrix}$
$\displaystyle=a\begin{pmatrix}-\bm{u}\cdot\nabla&-q_{0}\rho_{1}\bm{l}_{1}^{\rm
T}\Delta&q_{0}\rho_{2}\bm{l}_{2}^{\rm T}\Delta\\\
q_{0}\rho_{1}\bm{l}_{1}\Delta&O&O\\\
-q_{0}\rho_{2}\bm{l}_{2}\Delta&O&O\end{pmatrix}\dot{\bm{U}}$
$\displaystyle\quad+\begin{pmatrix}-q_{0}^{-1}\theta_{1}\theta_{2}(\bm{v}\cdot\nabla)^{2}(\bm{u}\cdot\nabla)\Delta^{-1}&-\theta_{1}\theta_{2}\rho_{1}\bm{l}_{1}^{\rm
T}(\bm{v}\cdot\nabla)^{2}&\theta_{1}\theta_{2}\rho_{2}\bm{l}_{2}^{\rm
T}(\bm{v}\cdot\nabla)^{2}\\\
\theta_{1}\theta_{2}\rho_{1}\bm{l}_{1}(\bm{v}\cdot\nabla)^{2}&-q_{0}\rho_{1}^{2}\bm{l}_{1}\bm{l}_{1}^{\rm
T}(\bm{u}\cdot\nabla)\Delta&q_{0}\rho_{1}\rho_{2}\bm{l}_{1}\bm{l}_{2}^{\rm
T}(\bm{u}\cdot\nabla)\Delta\\\
-\theta_{1}\theta_{2}\rho_{2}\bm{l}_{2}(\bm{v}\cdot\nabla)^{2}&q_{0}\rho_{1}\rho_{2}\bm{l}_{2}\bm{l}_{1}^{\rm
T}(\bm{u}\cdot\nabla)\Delta&-q_{0}\rho_{2}^{2}\bm{l}_{2}\bm{l}_{2}^{\rm
T}(\bm{u}\cdot\nabla)\Delta\end{pmatrix}\dot{\bm{U}}.$

Therefore, $\dot{\bm{U}}$ satisfies the symmetric system (5.11) with a skew-
symmetric operator $\mathscr{A}$ defined by

$\displaystyle\mathscr{A}$
$\displaystyle=a\begin{pmatrix}\bm{u}\cdot\nabla&q_{0}\rho_{1}\bm{l}_{1}^{\rm
T}\Delta&-q_{0}\rho_{2}\bm{l}_{2}^{\rm T}\Delta\\\
-q_{0}\rho_{1}\bm{l}_{1}\Delta&O&O\\\
q_{0}\rho_{2}\bm{l}_{2}\Delta&O&O\end{pmatrix}$
$\displaystyle\quad+\begin{pmatrix}q_{0}^{-1}\theta_{1}\theta_{2}(\bm{v}\cdot\nabla)^{2}(\bm{u}\cdot\nabla)\Delta^{-1}&\theta_{1}\theta_{2}\rho_{1}\bm{l}_{1}^{\rm
T}(\bm{v}\cdot\nabla)^{2}&-\theta_{1}\theta_{2}\rho_{2}\bm{l}_{2}^{\rm
T}(\bm{v}\cdot\nabla)^{2}\\\
-\theta_{1}\theta_{2}\rho_{1}\bm{l}_{1}(\bm{v}\cdot\nabla)^{2}&q_{0}\rho_{1}^{2}\bm{l}_{1}\bm{l}_{1}^{\rm
T}(\bm{u}\cdot\nabla)\Delta&-q_{0}\rho_{1}\rho_{2}\bm{l}_{1}\bm{l}_{2}^{\rm
T}(\bm{u}\cdot\nabla)\Delta\\\
\theta_{1}\theta_{2}\rho_{2}\bm{l}_{2}(\bm{v}\cdot\nabla)^{2}&-q_{0}\rho_{1}\rho_{2}\bm{l}_{2}\bm{l}_{1}^{\rm
T}(\bm{u}\cdot\nabla)\Delta&q_{0}\rho_{2}^{2}\bm{l}_{2}\bm{l}_{2}^{\rm
T}(\bm{u}\cdot\nabla)\Delta\end{pmatrix}.$

For the positive symmetric system (5.11), we can apply the standard theory for
partial differential equations to show its well-posedness of the initial value
problem. Moreover, these considerations help us to analyze the nonlinear
problem (2.14)–(2.16).

##  6 Analysis of related operators

We go back to consider the nonlinear problem, that is, the Kakinuma model
(2.14)–(2.16). We introduce second order differential operators
$L_{1,ij}=L_{1,ij}(H_{1})$ $(i,j=0,1,\ldots,N)$ and
$L_{2,ij}=L_{2,ij}(H_{2},b)$ $(i,j=0,1,\ldots,N^{*})$ by

(6.1) $\displaystyle L_{1,ij}\varphi_{1,j}$
$\displaystyle=-\nabla\cdot\left(\frac{1}{2(i+j)+1}H_{1}^{2(i+j)+1}\nabla\varphi_{1,j}\right)+\frac{4ij}{2(i+j)-1}H_{1}^{2(i+j)-1}\varphi_{1,j},$
(6.2) $\displaystyle L_{2,ij}\varphi_{2,j}$
$\displaystyle=-\nabla\cdot\left(\frac{1}{p_{i}+p_{j}+1}H_{2}^{p_{i}+p_{j}+1}\nabla\varphi_{2,j}-\frac{p_{j}}{p_{i}+p_{j}}H_{2}^{p_{i}+p_{j}}\varphi_{2,j}\nabla
b\right)$
$\displaystyle\quad\,-\frac{p_{i}}{p_{i}+p_{j}}H_{2}^{p_{i}+p_{j}}\nabla
b\cdot\nabla\varphi_{2,j}+\frac{p_{i}p_{j}}{p_{i}+p_{j}-1}H_{2}^{p_{i}+p_{j}-1}(1+|\nabla
b|^{2})\varphi_{2,j}.$

Then, we have $(L_{k,ij})^{*}=L_{k,ji}$ for $k=1,2$, where $(L_{k,ij})^{*}$ is
the adjoint operator of $L_{k,ij}$ in $L^{2}(\mathbf{R}^{n})$. We also use
$\bm{u}_{k}$ and $w_{k}$ for $k=1,2$ defined by (4.2) and (4.3), which
represent approximately the horizontal and the vertical components of the
velocity field on the interface from the water region $\Omega_{k}(t)$,
respectively. Then, the Kakinuma model (2.14)–(2.16) can be written simply as

$\begin{cases}\displaystyle
H_{1}^{2i}\partial_{t}\zeta+\sum_{j=0}^{N}L_{1,ij}(H_{1})\phi_{1,j}=0\quad\mbox{for}\quad
i=0,1,\ldots,N,\\\ \displaystyle-
H_{2}^{p_{i}}\partial_{t}\zeta+\sum_{j=0}^{N^{*}}L_{2,ij}(H_{2},b)\phi_{2,j}=0\quad\mbox{for}\quad
i=0,1,\ldots,N^{*},\\\
\displaystyle-\rho_{1}\left\\{\sum_{j=0}^{N}H_{1}^{2j}\partial_{t}\phi_{1,j}+g\zeta+\frac{1}{2}\bigl{(}|\bm{u}_{1}|^{2}+w_{1}^{2}\bigr{)}\right\\}\\\
\displaystyle\qquad+\rho_{2}\left\\{\sum_{j=0}^{N^{*}}H_{2}^{p_{j}}\partial_{t}\phi_{2,j}+g\zeta+\frac{1}{2}\bigl{(}|\bm{u}_{2}|^{2}+w_{2}^{2}\bigr{)}\right\\}=0.\end{cases}$

Moreover, introducing
$\bm{\phi}_{1}=(\phi_{1,0},\phi_{1,1},\ldots,\phi_{1,N})^{\rm T}$,
$\bm{\phi}_{2}=(\phi_{2,0},\phi_{2,1},\ldots,\phi_{2,N^{*}})^{\rm T}$, and

(6.3)
$\begin{cases}\bm{l}_{1}(H_{1})=(1,H_{1}^{2},H_{1}^{4},\ldots,H_{1}^{2N})^{\rm
T},\quad L_{1}(H_{1})=(L_{1,ij}(H_{1}))_{0\leq i,j\leq N},\\\
\bm{l}_{2}(H_{2})=(1,H_{2}^{p_{1}},H_{2}^{p_{2}},\ldots,H_{2}^{p_{N^{*}}})^{\rm
T},\quad L_{2}(H_{2},b)=(L_{2,ij}(H_{2},b))_{0\leq i,j\leq N^{*}},\end{cases}$

we can write the Kakinuma model (2.14)–(2.16) more simply as

(6.4)
$\begin{cases}\displaystyle\bm{l}_{1}(H_{1})\partial_{t}\zeta+L_{1}(H_{1})\bm{\phi}_{1}=\bm{0},\\\
\displaystyle-\bm{l}_{2}(H_{2})\partial_{t}\zeta+L_{2}(H_{2},b)\bm{\phi}_{2}=\bm{0},\\\
\displaystyle-\rho_{1}\left\\{{\bm{l}}_{1}(H_{1})\cdot\partial_{t}\bm{\phi}_{1}+g\zeta+\frac{1}{2}\bigl{(}|\bm{u}_{1}|^{2}+w_{1}^{2}\bigr{)}\right\\}\\\
\displaystyle\qquad+\rho_{2}\left\\{\bm{l}_{2}(H_{2})\cdot\partial_{t}\bm{\phi}_{2}+g\zeta+\frac{1}{2}\bigl{(}|\bm{u}_{2}|^{2}+w_{2}^{2}\bigr{)}\right\\}=0.\end{cases}$

By eliminating $\partial_{t}\zeta$ from the Kakinuma model, we obtain
$N+N^{*}+1$ scalar relations

$\begin{cases}\displaystyle\sum_{j=0}^{N}(L_{1,ij}(H_{1})\phi_{1,j}-H_{1}^{2i}L_{1,0j}(H_{1})\phi_{1,j})=0\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
\displaystyle\sum_{j=0}^{N^{*}}(L_{2,ij}(H_{2},b)\phi_{2,j}-H_{2}^{p_{i}}L_{2,0j}(H_{2},b)\phi_{2,j})=0\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\
\displaystyle\sum_{j=0}^{N}L_{1,0j}(H_{1})\phi_{1,j}+\sum_{j=0}^{N^{*}}L_{2,0j}(H_{2},b)\phi_{2,j}=0.\end{cases}$

These are compatibility conditions for the existence of the solution to the
Kakinuma model, and exactly the same as the compatibility conditions
(2.18)–(2.20). Introducing furthermore linear operators
$\mathcal{L}_{1,i}=\mathcal{L}_{1,i}(H_{1})$ $(i=0,1,\ldots,N)$ acting on
$\bm{\varphi}_{1}=(\varphi_{1,0},\ldots,\varphi_{1,N})^{\rm T}$ and
$\mathcal{L}_{2,i}=\mathcal{L}_{2,i}(H_{2},b)$ $(i=0,1,\ldots,N^{*})$ acting
on $\bm{\varphi}_{2}=(\varphi_{2,0},\ldots,\varphi_{2,N^{*}})^{\rm T}$ by

(6.5)
$\begin{cases}\displaystyle\mathcal{L}_{1,0}(H_{1})\bm{\varphi}_{1}=\sum_{j=0}^{N}L_{1,0j}(H_{1})\varphi_{1,j},\\\
\displaystyle\mathcal{L}_{1,i}(H_{1})\bm{\varphi}_{1}=\sum_{j=0}^{N}(L_{1,ij}(H_{1})\varphi_{1,j}-H_{1}^{2i}L_{1,0j}(H_{1})\varphi_{1,j})\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
\displaystyle\mathcal{L}_{2,0}(H_{2},b)\bm{\varphi}_{2}=\sum_{j=0}^{N^{*}}L_{2,0j}(H_{2},b)\varphi_{2,j},\\\
\displaystyle\mathcal{L}_{2,i}(H_{2},b)\bm{\varphi}_{2}=\sum_{j=0}^{N^{*}}(L_{2,ij}(H_{2},b)\varphi_{2,j}-H_{2}^{p_{i}}L_{2,0j}(H_{2},b)\varphi_{2,j})\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\end{cases}$

the compatibility conditions can be written simply as

(6.6)
$\begin{cases}\mathcal{L}_{1,i}(H_{1})\bm{\phi}_{1}=0\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
\mathcal{L}_{2,i}(H_{2},b)\bm{\phi}_{2}=0\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\
\mathcal{L}_{1,0}(H_{1})\bm{\phi}_{1}+\mathcal{L}_{2,0}(H_{2},b)\bm{\phi}_{2}=0.\end{cases}$

We proceed to derive evolution equations for $\bm{\phi}_{1}$ and
$\bm{\phi}_{2}$. To this end, we differentiate the above compatibility
conditions with respect to $t$ and use equations of the Kakinuma model to
eliminate $\partial_{t}\zeta$. Then, we obtain

(6.7)
$\begin{cases}\mathcal{L}_{1,i}(H_{1})\partial_{t}\bm{\phi}_{1}=F_{1,i}\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
\mathcal{L}_{2,i}(H_{2},b)\partial_{t}\bm{\phi}_{2}=F_{2,i}\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\
\mathcal{L}_{1,0}(H_{1})\partial_{t}\bm{\phi}_{1}+\mathcal{L}_{2,0}(H_{2},b)\partial_{t}\bm{\phi}_{2}=F_{3},\end{cases}$

where

$\begin{cases}F_{1,i}=-\frac{\partial\mathcal{L}_{1,i}}{\partial
H_{1}}(H_{1})[\mathcal{L}_{1,0}(H_{1})\bm{\phi}_{1}]\bm{\phi}_{1}\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\ F_{2,i}=-\frac{\partial\mathcal{L}_{2,i}}{\partial
H_{2}}(H_{2},b)[\mathcal{L}_{2,0}(H_{2},b)\bm{\phi}_{2}]\bm{\phi}_{2}\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\ F_{3}=-\frac{\partial\mathcal{L}_{1,0}}{\partial
H_{1}}(H_{1})[\mathcal{L}_{1,0}(H_{1})\bm{\phi}_{1}]\bm{\phi}_{1}-\frac{\partial\mathcal{L}_{2,0}}{\partial
H_{2}}(H_{2},b)[\mathcal{L}_{2,0}(H_{2},b)\bm{\phi}_{2}]\bm{\phi}_{2}.\end{cases}$

Here, we note that $F_{3}$ can be written in divergence form as

$F_{3}=\nabla\cdot\biggl{\\{}(\mathcal{L}_{1,0}(H_{1})\bm{\phi}_{1})\sum_{j=0}^{N}H_{1}^{2j}\nabla\phi_{1,j}+(\mathcal{L}_{2,0}(H_{2},b)\bm{\phi}_{2})\sum_{j=0}^{N^{*}}H_{2}^{p_{j}}\nabla\phi_{2,j}\biggr{\\}}.$

On the other hand, the last equation in the Kakinuma model can be written as

(6.8)
$-\rho_{1}{\bm{l}}_{1}(H_{1})\cdot\partial_{t}{\bm{\phi}}_{1}+\rho_{2}{\bm{l}}_{2}(H_{2})\cdot\partial_{t}{\bm{\phi}}_{2}=F_{4},$

where

$F_{4}=\rho_{1}\biggl{\\{}g\zeta+\frac{1}{2}\bigl{(}|{\bm{u}}_{1}|^{2}+w_{1}^{2}\bigr{)}\biggr{\\}}-\rho_{2}\biggl{\\{}g\zeta+\frac{1}{2}\bigl{(}|{\bm{u}}_{2}|^{2}+w_{2}^{2}\bigr{)}\biggr{\\}}.$

In view of these evolution equations (6.7)–(6.8) for $\bm{\phi}_{1}$ and
$\bm{\phi}_{2}$, we will consider the following equations for
$\bm{\varphi}_{1}$ and $\bm{\varphi}_{2}$.

(6.9)
$\begin{cases}\mathcal{L}_{1,i}(H_{1})\bm{\varphi}_{1}=f_{1,i}\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
\mathcal{L}_{2,i}(H_{2},b)\bm{\varphi}_{2}=f_{2,i}\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\
\mathcal{L}_{1,0}(H_{1})\bm{\varphi}_{1}+\mathcal{L}_{2,0}(H_{2},b)\bm{\varphi}_{2}=\nabla\cdot\bm{f}_{3},\\\
-\rho_{1}\bm{l}_{1}(H_{1})\cdot\bm{\varphi}_{1}+\rho_{2}\bm{l}_{2}(H_{2})\cdot\bm{\varphi}_{2}=f_{4}.\end{cases}$

In the following we will use the notation
$\bm{\varphi}_{1}^{\prime}=(\varphi_{1,1},\ldots,\varphi_{1,N})^{\rm T}$ and
$\bm{\varphi}_{2}^{\prime}=(\varphi_{2,1},\ldots,\varphi_{2,N^{*}})^{\rm T}$,
and we put $\bm{f}_{1}^{\prime}=(f_{1,1},\ldots,f_{1,N})^{\rm T}$ and
$\bm{f}_{2}^{\prime}=(f_{2,1},\ldots,f_{2,N^{*}})^{\rm T}$.

######  Lemma 6.1.

Let $c_{0}$ and $c_{1}$ be positive constants. There exists a positive
constant $C=C(c_{0},c_{1})$ depending only on $c_{0}$ and $c_{1}$ such that
for any $H_{1},H_{2},\nabla b\in L^{\infty}(\mathbf{R}^{n})$ satisfying
$H_{1}(x),H_{2}(x)\geq c_{0}$ and $|\nabla b(x)|\leq c_{1}$, any regular
solution $(\bm{\varphi}_{1},\bm{\varphi}_{2})$ to (6.9) satisfies

$\displaystyle\rho_{1}(\|\nabla\varphi_{1,0}\|_{L^{2}}^{2}+\|{\bm{\varphi}}_{1}^{\prime}\|_{H^{1}}^{2})+\rho_{2}(\|\nabla\varphi_{2,0}\|_{L^{2}}^{2}+\|{\bm{\varphi}}_{2}^{\prime}\|_{H^{1}}^{2})$
$\displaystyle\leq C\biggl{(}-\sum_{j=0}^{N}(\nabla
f_{4},\frac{1}{2j+1}H_{1}^{2j+1}\nabla\varphi_{1,j})_{L^{2}}$
$\displaystyle\qquad\quad+\rho_{1}({\bm{f}}_{1}^{\prime},{\bm{\varphi}}_{1}^{\prime})_{L^{2}}+\rho_{2}({\bm{f}}_{2}^{\prime},{\bm{\varphi}}_{2}^{\prime})_{L^{2}}+\rho_{2}(\nabla\cdot{\bm{f}}_{3},{\bm{l}}_{2}(H_{2})\cdot{\bm{\varphi}}_{2})_{L^{2}}\biggr{)}.$

###### Proof..

We introduce a dummy variable $\eta$ by

$\eta=-\mathcal{L}_{1,0}(H_{1})\bm{\varphi}_{1}.$

Then, we can rewrite the equations in (6.9) as

$\begin{cases}\eta\bm{l}_{1}(H_{1})+L_{1}(H_{1})\bm{\varphi}_{1}=\bm{f}_{1}=(0,f_{1,1},\ldots,f_{1,N})^{\rm
T},\\\
-\eta\bm{l}_{2}(H_{2})+L_{2}(H_{2},b)\bm{\varphi}_{2}=\bm{f}_{2}=(0,f_{2,1},\ldots,f_{2,N^{*}})^{\rm
T}+(\nabla\cdot\bm{f}_{3})\bm{l}_{2}(H_{2}),\\\
-\rho_{1}\bm{l}_{1}(H_{1})\cdot\bm{\varphi}_{1}+\rho_{2}\bm{l}_{2}(H_{2})\cdot\bm{\varphi}_{2}=f_{4},\end{cases}$

that is,

$\begin{pmatrix}0&-\rho_{1}{\bm{l}}_{1}(H_{1})^{\rm
T}&\rho_{2}{\bm{l}}_{2}(H_{2})^{\rm T}\\\
\rho_{1}{\bm{l}}_{1}(H_{1})&\rho_{1}L_{1}(H_{1})&O\\\
-\rho_{2}{\bm{l}}_{2}(H_{2})&O&\rho_{2}L_{2}(H_{2},b)\end{pmatrix}\begin{pmatrix}\eta\\\
{\bm{\varphi}}_{1}\\\ {\bm{\varphi}}_{2}\end{pmatrix}=\begin{pmatrix}f_{4}\\\
\rho_{1}{\bm{f}}_{1}\\\ \rho_{2}{\bm{f}}_{2}\end{pmatrix}.$

By taking the $L^{2}$-inner product of this equation with
$(\eta,\bm{\varphi}_{1},\bm{\varphi}_{2})^{\rm T}$, we see that

$\displaystyle\rho_{1}(L_{1}(H_{1})\bm{\varphi}_{1},\bm{\varphi}_{1})_{L^{2}}+\rho_{2}(L_{2}(H_{2},b)\bm{\varphi}_{2},\bm{\varphi}_{2})_{L^{2}}$
$\displaystyle=(f_{4},\eta)_{L^{2}}+\rho_{1}(\bm{f}_{1},\bm{\varphi}_{1})_{L^{2}}+\rho_{2}(\bm{f}_{2},\bm{\varphi}_{2})_{L^{2}}$
$\displaystyle=-\sum_{j=0}^{N}(\nabla
f_{4},\frac{1}{2j+1}H_{1}^{2j+1}\nabla\varphi_{1,j})_{L^{2}}$
$\displaystyle\quad\;+\rho_{1}(\bm{f}_{1}^{\prime},\bm{\varphi}_{1}^{\prime})_{L^{2}}+\rho_{2}(\bm{f}_{2}^{\prime},\bm{\varphi}_{2}^{\prime})_{L^{2}}+\rho_{2}(\nabla\cdot\bm{f}_{3},\bm{l}_{2}(H_{2})\cdot{\bm{\varphi}}_{2})_{L^{2}}.$

Here, by direct calculation we have

(6.10) $\displaystyle(L_{1}(H_{1})\bm{\varphi}_{1},\bm{\varphi}_{1})_{L^{2}}$
$\displaystyle=\sum_{i,j=0}^{N}(L_{1,ij}(H_{1})\varphi_{1,j},\varphi_{1,i})_{L^{2}}$
$\displaystyle=\int_{\mathbf{R}^{n}}\\!{\rm
d}\bm{x}\\!\int_{0}^{H_{1}}\left\\{\left|\sum_{i=0}^{N}(z^{2i}\nabla\varphi_{1,i})\right|^{2}+\left(\sum_{i=0}^{N}2iz^{2i-1}\varphi_{1,i}\right)^{2}\right\\}{\rm
d}z$ $\displaystyle\simeq\int_{\mathbf{R}^{n}}\\!{\rm
d}\bm{x}\\!\int_{0}^{H_{1}}\sum_{i=0}^{N}\left(z^{4i}|\nabla\varphi_{1,i}|^{2}+i^{2}z^{4i-2}\varphi_{1,i}^{2}\right){\rm
d}z$
$\displaystyle\simeq\int_{\mathbf{R}^{n}}\sum_{i=0}^{N}\left(H_{1}^{4i+1}|\nabla\varphi_{1,i}|^{2}+i^{2}H_{1}^{4i-1}\varphi_{1,i}^{2}\right){\rm
d}\bm{x},$

where we used the fact that $\\{z^{2i}\\}_{i=0,\dots,N}$ and
$\\{z^{2i-1}\\}_{i=1,\dots,N}$ are both linearly independent. We have also

$\displaystyle(L_{2}(H_{2},b)\bm{\varphi}_{2},\bm{\varphi}_{2})_{L^{2}}=\sum_{i,j=0}^{N^{*}}(L_{2,ij}(H_{2},b)\varphi_{2,j},\varphi_{2,i})_{L^{2}}$
$\displaystyle=\int_{\mathbf{R}^{n}}\\!{\rm
d}\bm{x}\\!\int_{0}^{H_{2}}\left\\{\left|\sum_{i=0}^{N^{*}}(z^{p_{i}}\nabla\varphi_{2,i}-p_{i}z^{p_{i}-1}\varphi_{2,i}\nabla
b)\right|^{2}+\left(\sum_{i=0}^{N^{*}}p_{i}z^{p_{i}-1}\varphi_{2,i}\right)^{2}\right\\}{\rm
d}z.$

If $\\{z^{p_{i}},z^{p_{i}-1}\\}_{i=0,\dots,N}$ are linearly independent, then
we have

(6.11)
$\displaystyle(L_{2}(H_{2},b)\bm{\varphi}_{2},\bm{\varphi}_{2})_{L^{2}}$
$\displaystyle\simeq\int_{\mathbf{R}^{n}}\\!{\rm
d}\bm{x}\\!\int_{0}^{H_{2}}\sum_{i=0}^{N^{*}}\left\\{\left(z^{2p_{i}}|\nabla\varphi_{2,i}|^{2}+p_{i}^{2}z^{2p_{i}-2}|\nabla
b|^{2}\varphi_{2,i}^{2}\right)+p_{i}^{2}z^{2p_{i}-2}\varphi_{2,i}^{2}\right\\}{\rm
d}z$
$\displaystyle\simeq\int_{\mathbf{R}^{n}}\sum_{i=0}^{N^{*}}\left\\{H_{2}^{2p_{i}+1}|\nabla\varphi_{2,i}|^{2}+p_{i}^{2}H_{2}^{2p_{i}-1}(1+|\nabla
b|^{2})\varphi_{2,i}^{2}\right\\}{\rm d}\bm{x}.$

Otherwise, for example, in the case $p_{i}=i$ $(i=0,\dots,N)$ we obtain

(6.12)
$\displaystyle(L_{2}(H_{2},b)\bm{\varphi}_{2},\bm{\varphi}_{2})_{L^{2}}$
$\displaystyle=\int_{\mathbf{R}^{n}}\\!{\rm
d}\bm{x}\\!\int_{0}^{H_{2}}\left\\{\left|\sum_{i=0}^{N^{*}-1}z^{i}(\nabla\varphi_{2,i}-(i+1)\varphi_{2,i+1}\nabla
b)+z^{N^{*}}\nabla\varphi_{2,N^{*}}\right|^{2}+\left(\sum_{i=0}^{N^{*}}p_{i}z^{p_{i}-1}\varphi_{2,i}\right)^{2}\right\\}{\rm
d}z$ $\displaystyle\simeq\int_{\mathbf{R}^{n}}\\!{\rm
d}\bm{x}\\!\int_{0}^{H_{2}}\left\\{\sum_{i=0}^{N^{*}-1}(z^{2i}|\nabla\varphi_{2,i}-(i+1)\varphi_{2,i+1}\nabla
b|^{2}+z^{2N^{*}}|\nabla\varphi_{2,N^{*}}|^{2}+\sum_{i=0}^{N^{*}}i^{2}z^{2(i-1)}\varphi_{2,i}^{2}\right\\}{\rm
d}z$
$\displaystyle\simeq\int_{\mathbf{R}^{n}}\left\\{\sum_{i=0}^{N^{*}-1}(H_{2}^{2i+1}|\nabla\varphi_{2,i}-(i+1)\varphi_{2,i+1}\nabla
b|^{2}+H_{2}^{2N^{*}+1}|\nabla\varphi_{2,N^{*}}|^{2}+\sum_{i=0}^{N^{*}}i^{2}H_{2}^{2i-1}\varphi_{2,i}^{2}\right\\}{\rm
d}\bm{x}.$

A similar estimate holds in other cases. These estimates give the desired one.
∎

Although this lemma gives an a priori bound of the solution to (6.9), the
equations in (6.9) do not have a good symmetry. In order to give an existence
theorem to (6.9) with robust elliptic estimates, it is better to rewrite them
in a symmetric form by introducing a good unknown variable. We introduce
scalar functions $\varphi_{1}$ and $\varphi_{2}$ by

(6.13)
$\varphi_{1}={\bm{l}}_{1}(H_{1})\cdot{\bm{\varphi}}_{1},\qquad\varphi_{2}={\bm{l}}_{2}(H_{2})\cdot{\bm{\varphi}}_{2}.$

We also introduce second order differential operators $P_{1,i}(H_{1})$
$(i=1,\ldots,N)$ and $Q_{1}(H_{1})$ acting on $\mathbf{R}^{N}$-valued
functions $\bm{\varphi}_{1}^{\prime}=(\varphi_{1,1},\ldots,\varphi_{1,N})^{\rm
T}$ and $P_{2,i}(H_{2},b)$ $(i=1,\ldots,N^{*})$ and $Q_{2}(H_{2})$ acting on
$\mathbf{R}^{N^{*}}$-valued functions
$\bm{\varphi}_{2}^{\prime}=(\varphi_{2,1},\ldots,\varphi_{2,N^{*}})^{\rm T}$
by

(6.14) $\begin{cases}\displaystyle
P_{1,i}(H_{1}){\bm{\varphi}}_{1}^{\prime}=\sum_{j=1}^{N}\bigl{\\{}\bigl{(}L_{1,ij}(H_{1})-H_{1}^{2i}L_{1,0j}(H_{1})\bigr{)}\varphi_{1,j}\\\
\makebox[68.88867pt]{}-\bigl{(}L_{1,i0}(H_{1})-H_{1}^{2i}L_{1,00}(H_{1})\bigr{)}\bigl{(}H_{1}^{2j}\varphi_{1,j}\bigr{)}\bigr{\\}},\\\
\displaystyle
Q_{1}(H_{1}){\bm{\varphi}}_{1}^{\prime}=\sum_{j=1}^{N}\bigl{\\{}L_{1,0j}(H_{1})\varphi_{1,j}-L_{1,00}(H_{1})\bigl{(}H_{1}^{2j}\varphi_{1,j}\bigr{)}\bigr{\\}},\\\
\end{cases}$

and

(6.15) $\begin{cases}\displaystyle
P_{2,i}(H_{2},b){\bm{\varphi}}_{2}^{\prime}=\sum_{j=1}^{N^{*}}\bigl{\\{}\bigl{(}L_{2,ij}(H_{2},b)-H_{2}^{p_{i}}L_{2,0j}(H_{2},b)\bigr{)}\varphi_{2,j}\\\
\makebox[77.49976pt]{}-\bigl{(}L_{2,i0}(H_{2},b)-H_{2}^{p_{i}}L_{2,00}(H_{2},b)\bigr{)}\bigl{(}H_{2}^{p_{j}}\varphi_{2,j}\bigr{)}\bigr{\\}},\\\
\displaystyle
Q_{2}(H_{2},b){\bm{\varphi}}_{2}^{\prime}=\sum_{j=1}^{N^{*}}\bigl{\\{}L_{2,0j}(H_{2},b)\varphi_{2,j}-L_{2,00}(H_{2},b)\bigl{(}H_{2}^{p_{j}}\varphi_{2,j}\bigr{)}\bigr{\\}},\end{cases}$

respectively, and put

$\begin{cases}P_{1}(H_{1})\bm{\varphi}_{1}^{\prime}=(P_{1,1}(H_{1})\bm{\varphi}_{1}^{\prime},\ldots,P_{1,N}(H_{1})\bm{\varphi}_{1}^{\prime})^{\rm
T},\\\
P_{2}(H_{2},b)\bm{\varphi}_{2}^{\prime}=(P_{2,1}(H_{2},b)\bm{\varphi}_{2}^{\prime},\ldots,P_{2,N^{*}}(H_{2},b)\bm{\varphi}_{2}^{\prime})^{\rm
T}.\end{cases}$

Then, we see easily that $P_{1}(H_{1})$ and $P_{2}(H_{2},b)$ are symmetric in
$L^{2}(\mathbf{R}^{n})$ and that

$\displaystyle\mathcal{L}_{1,i}(H_{1})\bm{\varphi}_{1}=\begin{cases}Q_{1}(H_{1})\bm{\varphi}_{1}^{\prime}+L_{1,00}(H_{1})(\bm{l}_{1}(H_{1})\cdot\bm{\varphi}_{1})&\mbox{for}\quad
i=0,\\\
P_{1,i}(H_{1})\bm{\varphi}_{1}^{\prime}+((Q_{1}(H_{1}))^{*}(\bm{l}_{1}(H_{1})\cdot\bm{\varphi}_{1}))_{i}&\mbox{for}\quad
i=1,\ldots,N,\end{cases}$
$\displaystyle\mathcal{L}_{2,i}(H_{2},b)\bm{\varphi}_{2}=\begin{cases}Q_{2}(H_{2},b)\bm{\varphi}_{2}^{\prime}+L_{2,00}(H_{2},b)(\bm{l}_{2}(H_{2})\cdot\bm{\varphi}_{2})&\mbox{for}\quad
i=0,\\\
P_{2,i}(H_{2},b){\bm{\varphi}}_{2}^{\prime}+((Q_{2}(H_{2},b))^{*}(\bm{l}_{2}(H_{2})\cdot\bm{\varphi}_{2}))_{i}&\mbox{for}\quad
i=1,\ldots,N^{*},\\\ \end{cases}$

where $Q^{*}$ denotes an adjoint operator of $Q$ in $L^{2}(\mathbf{R}^{n})$.
Therefore, we can rewrite (6.9) as

$\begin{cases}P_{1}(H_{1})\bm{\varphi}_{1}^{\prime}+(Q_{1}(H_{1}))^{*}\varphi_{1}=\bm{f}_{1},\\\
P_{2}(H_{2},b)\bm{\varphi}_{2}^{\prime}+(Q_{2}(H_{2},b))^{*}\varphi_{2}=\bm{f}_{2},\\\
Q_{1}(H_{1})\bm{\varphi}_{1}^{\prime}+L_{1,00}(H_{1})\varphi_{1}+Q_{2}(H_{2},b)\bm{\varphi}_{2}^{\prime}+L_{2,00}(H_{2},b)\varphi_{2}=\nabla\cdot\bm{f}_{3},\\\
-\rho_{1}\varphi_{1}+\rho_{2}\varphi_{2}=f_{4}.\end{cases}$

These equations for
$(\bm{\varphi}_{1}^{\prime},\varphi_{1},\bm{\varphi}_{2}^{\prime},\varphi_{2})$
do not have yet any good symmetry. But, it follows from the last equation that

$\rho_{2}\varphi_{2}=\rho_{1}\varphi_{1}+f_{4}.$

Using this we can remove $\varphi_{2}$ from the equations and obtain

$\begin{cases}\rho_{1}P_{1}(H_{1})\bm{\varphi}_{1}^{\prime}+\rho_{1}(Q_{1}(H_{1}))^{*}\varphi_{1}=\rho_{1}\bm{F}_{1},\\\
\rho_{2}P_{2}(H_{2},b)\bm{\varphi}_{2}^{\prime}+\rho_{1}(Q_{2}(H_{2},b))^{*}\varphi_{1}=\rho_{2}\bm{F}_{2},\\\
\rho_{1}Q_{1}(H_{1})\bm{\varphi}_{1}^{\prime}+\rho_{1}Q_{2}(H_{2},b)\bm{\varphi}_{2}^{\prime}+\rho_{1}\left(L_{1,00}(H_{1})+\frac{\rho_{1}}{\rho_{2}}L_{2,00}(H_{2},b)\right)\varphi_{1}=\rho_{1}\nabla\cdot\bm{F}_{3},\end{cases}$

where

(6.16)
$\bm{F}_{1}=\bm{f}_{1},\quad\bm{F}_{2}=\bm{f}_{2}-\frac{1}{\rho_{2}}(Q_{2}(H_{2},b))^{*}f_{4},\quad\bm{F}_{3}=\bm{f}_{3}+\frac{1}{\rho_{2}}H_{2}\nabla
f_{4}.$

These equations for
$(\bm{\varphi}_{1}^{\prime},\bm{\varphi}_{2}^{\prime},\varphi_{1})$ have a
good symmetry and can be written in the matrix form

(6.17) $\mathscr{P}(\zeta,b)\begin{pmatrix}\bm{\varphi}_{1}^{\prime}\\\
\bm{\varphi}_{2}^{\prime}\\\
\varphi_{1}\end{pmatrix}=\begin{pmatrix}\rho_{1}\bm{F}_{1}\\\
\rho_{2}\bm{F}_{2}\\\ \rho_{1}\nabla\cdot\bm{F}_{3}\end{pmatrix},$

where

(6.18)
$\mathscr{P}(\zeta,b)=\begin{pmatrix}\rho_{1}P_{1}(H_{1})&O&\rho_{1}(Q_{1}(H_{1}))^{*}\\\
O&\rho_{2}P_{2}(H_{2},b)&\rho_{1}(Q_{2}(H_{2},b))^{*}\\\
\rho_{1}Q_{1}(H_{1})&\rho_{1}Q_{2}(H_{2},b)&\rho_{1}\left(L_{1,00}(H_{1})+\frac{\rho_{1}}{\rho_{2}}L_{2,00}(H_{2},b)\right)\end{pmatrix},$

which is symmetric in $L^{2}(\mathbf{R}^{n})$. Moreover,
$\mathscr{P}(\zeta,b)$ is positive in $L^{2}(\mathbf{R}^{n})$ as shown in the
following lemma.

######  Lemma 6.2.

Let $c_{0},c_{1}$ be positive constants. There exists a positive constant
$C=C(c_{0},c_{1})$ depending only on $c_{0}$ and $c_{1}$ such that if
$\zeta,b\in W^{1,\infty}(\mathbf{R}^{n})$ satisfy $H_{1}(x),H_{2}(x)\geq
c_{0}$ and $H_{1}(x)+|\nabla H_{1}(x)|+|\nabla b(x)|\leq c_{1}$, then for any
$\widetilde{\bm{\varphi}}=(\bm{\varphi}_{1}^{\prime},\bm{\varphi}_{2}^{\prime},\varphi_{1})^{\rm
T}$ we have

$(\mathscr{P}(\zeta,b)\widetilde{\bm{\varphi}},\widetilde{\bm{\varphi}})_{L^{2}}\geq
C^{-1}\bigl{(}\rho_{1}\|\bm{\varphi}_{1}^{\prime}\|_{H^{1}}^{2}+\rho_{2}\|\bm{\varphi}_{2}^{\prime}\|_{H^{1}}^{2}+\rho_{1}\|\nabla\varphi_{1}\|_{L^{2}}^{2}\bigr{)}.$

###### Proof..

Given
$\widetilde{\bm{\varphi}}=(\bm{\varphi}_{1}^{\prime},\bm{\varphi}_{2}^{\prime},\varphi_{1})^{\rm
T}$, we define $\varphi_{1,0}$ and $\varphi_{2,0}$ by

$\varphi_{1,0}=\varphi_{1}-\sum_{j=1}^{N}H_{1}^{2j}\varphi_{1,j},\qquad\varphi_{2,0}=\frac{\rho_{1}}{\rho_{2}}\varphi_{1}-\sum_{j=1}^{N^{*}}H_{2}^{p_{j}}\varphi_{2,j}$

and put
$\bm{\varphi}_{1}=(\varphi_{1,0},\varphi_{1,1},\ldots,\varphi_{1,N})^{\rm T}$
and
$\bm{\varphi}_{2}=(\varphi_{2,0},\varphi_{2,1},\ldots,\varphi_{2,N^{*}})^{\rm
T}$. Then, we have
$\varphi_{1}=\bm{l}_{1}(H_{1})\cdot\bm{\varphi}_{1}=\frac{\rho_{2}}{\rho_{1}}\bm{l}_{2}(H_{2})\cdot\bm{\varphi}_{2}$,
so that

$\rho_{1}\bm{l}_{1}(H_{1})\cdot\bm{\varphi}_{1}-\rho_{2}\bm{l}_{2}(H_{2})\cdot\bm{\varphi}_{2}=0.$

We also define $\bm{F}_{1}=(F_{1,1},\ldots,F_{1,N})^{\rm T}$,
$\bm{F}_{2}=(F_{2,1},\ldots,F_{2,N^{*}})^{\rm T}$, and $F_{3}$ by

$\begin{pmatrix}\bm{F}_{1}\\\ \bm{F}_{2}\\\
F_{3}\end{pmatrix}=\mathscr{P}(\zeta,b)\begin{pmatrix}\bm{\varphi}_{1}^{\prime}\\\
\bm{\varphi}_{2}^{\prime}\\\ \varphi_{1}\end{pmatrix}.$

Then, we have

$\begin{cases}\rho_{1}\mathcal{L}_{1,i}(H_{1})\bm{\varphi}_{1}=F_{1,i}\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
\rho_{2}\mathcal{L}_{2,i}(H_{2},b)\bm{\varphi}_{2}=F_{2,i}\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\
\rho_{1}\mathcal{L}_{1,0}(H_{1})\bm{\varphi}_{1}+\rho_{2}\mathcal{L}_{2,0}(H_{2},b)\bm{\varphi}_{2}=F_{3}.\end{cases}$

Now, we introduce a dummy variable $\eta$ by

$\eta=-\mathcal{L}_{1,0}(H_{1})\bm{\varphi}_{1}.$

Then, it follows from the above equations that

$\begin{cases}-\rho_{1}\bm{l}_{1}(H_{1})\cdot\bm{\varphi}_{1}+\rho_{2}\bm{l}_{2}(H_{2})\cdot\bm{\varphi}_{2}=0,\\\
\rho_{1}(\eta\bm{l}_{1}(H_{1})+L_{1}(H_{1})\bm{\varphi}_{1})=\bm{f}_{1},\\\
\rho_{2}(-\eta\bm{l}_{2}(H_{2})+L_{2}(H_{2},b)\bm{\varphi}_{2})=\bm{f}_{2}+\frac{\rho_{2}}{\rho_{1}}\bm{l}_{2}(H_{2})F_{3},\\\
\end{cases}$

where $\bm{f}_{1}=(0,F_{1,1},\ldots,F_{1,N})^{\rm T}$ and
$\bm{f}_{2}=(0,F_{2,1},\ldots,F_{2,N^{*}})^{\rm T}$. These equations can be
written in the matrix form

$\begin{pmatrix}0&-\rho_{1}\bm{l}_{1}(H_{1})^{\rm
T}&\rho_{2}\bm{l}_{2}(H_{2})^{\rm T}\\\
\rho_{1}\bm{l}_{1}(H_{1})&\rho_{1}L_{1}(H_{1})&O\\\
-\rho_{2}\bm{l}_{2}(H_{2})&O&\rho_{2}L_{2}(H_{2},b)\end{pmatrix}\begin{pmatrix}\eta\\\
\bm{\varphi}_{1}\\\ \bm{\varphi}_{2}\end{pmatrix}=\begin{pmatrix}0\\\
\bm{f}_{1}\\\
\bm{f}_{2}+\frac{\rho_{2}}{\rho_{1}}\bm{l}_{2}(H_{2})F_{3}\end{pmatrix}.$

By taking the $L^{2}$-inner product of this equation with
$(\eta,\bm{\varphi}_{1},\bm{\varphi}_{2})^{\rm T}$ we see that

$\displaystyle\rho_{1}(L_{1}(H_{1})\bm{\varphi}_{1},\bm{\varphi}_{1})_{L^{2}}+\rho_{2}(L_{2}(H_{2},b)\bm{\varphi}_{2},\bm{\varphi}_{2})_{L^{2}}$
$\displaystyle=(\bm{f}_{1},\bm{\varphi}_{1})_{L^{2}}+(\bm{f}_{2},\bm{\varphi}_{2})_{L^{2}}+\frac{\rho_{2}}{\rho_{1}}(\bm{l}_{2}(H_{2})F_{3},\bm{\varphi}_{2})_{L^{2}}$
$\displaystyle=(\bm{F}_{1},\bm{\varphi}_{1}^{\prime})_{L^{2}}+(\bm{F}_{2},\bm{\varphi}_{2}^{\prime})_{L^{2}}+(F_{3},\varphi_{1})_{L^{2}}$
$\displaystyle=(\mathscr{P}(\zeta,b)\widetilde{\bm{\varphi}},\widetilde{\bm{\varphi}})_{L^{2}},$

which gives, by (6.10) and (6.11) or (6.12),

$\displaystyle(\mathscr{P}(\zeta,b)\widetilde{\bm{\varphi}},\widetilde{\bm{\varphi}})_{L^{2}}\simeq\rho_{1}(\|\bm{\varphi}_{1}^{\prime}\|_{H^{1}}^{2}+\|\nabla\varphi_{1,0}\|_{L^{2}}^{2})+\rho_{2}(\|\bm{\varphi}_{2}^{\prime}\|_{H^{1}}^{2}+\|\nabla\varphi_{2,0}\|_{L^{2}}^{2}).$

Since
$\|\nabla\varphi_{1}\|_{L^{2}}^{2}\lesssim\|\bm{\varphi}_{1}^{\prime}\|_{H^{1}}^{2}+\|\nabla\varphi_{1,0}\|_{L^{2}}^{2}$,
we obtain the desired estimate. ∎

By this lemma, the explicit expression (6.18) of the operator
$\mathscr{P}(\zeta,b)$, and the standard theory of elliptic partial
differential equations, we can obtain the following lemma.

######  Lemma 6.3.

Let $\rho_{1},\rho_{2},h_{1},h_{2},c_{0},M$ be positive constants and $m$ an
integer such that $m>\frac{n}{2}+1$. There exists a positive constant
$C=C(\rho_{1},\rho_{2},h_{1},h_{2},c_{0},m)$ such that if $\zeta$ and $b$
satisfy

$\begin{cases}\|\zeta\|_{H^{m}}+\|b\|_{W^{m,\infty}}\leq M,\\\
H_{1}(\bm{x})=h_{1}-\zeta(\bm{x})\geq c_{0},\quad
H_{2}(\bm{x})=h_{2}+\zeta(\bm{x})-b(\bm{x})\geq
c_{0}\quad\mbox{for}\quad\bm{x}\in\mathbf{R}^{n},\end{cases}$

then for any $\bm{F}_{1},\bm{F}_{2}\in H^{k-1}$ and $\bm{F}_{3}\in H^{k}$ with
$k\in\\{0,1,\ldots,m-1\\}$ there exists a solution
$(\bm{\varphi}_{1}^{\prime},\bm{\varphi}_{2}^{\prime},\varphi_{1})$ of (6.17)
satisfying

$\|(\bm{\varphi}_{1}^{\prime},\bm{\varphi}_{2}^{\prime})\|_{H^{k+1}}+\|\nabla\varphi_{1}\|_{H^{k}}\leq
C\left(\|(\bm{F}_{1},\bm{F}_{2})\|_{H^{k-1}}+\|\bm{F}_{3}\|_{H^{k}}\right).$

Moreover, the solution is unique up to an additive constant to $\varphi_{1}$.

We proceed to consider solvability to (6.9). Given
$\bm{f}_{1}^{\prime},\bm{f}_{2}^{\prime},\bm{f}_{3},f_{4}$, we define
$\bm{F}_{1},\bm{F}_{2},\bm{F}_{3}$ by (6.16), for which there exists a
solution $(\bm{\varphi}_{1}^{\prime},\bm{\varphi}_{2}^{\prime},\varphi_{1})$
to (6.17), define $\varphi_{1,0}$ and $\varphi_{2,0}$ by

$\varphi_{1,0}=\varphi_{1}-\sum_{j=1}^{N}H_{1}^{2j}\varphi_{1,j},\qquad\varphi_{2,0}=\frac{\rho_{1}}{\rho_{2}}\varphi_{1}-\sum_{j=1}^{N^{*}}H_{2}^{p_{j}}\varphi_{2,j}+\frac{1}{\rho_{2}}f_{4},$

and put
$\bm{\varphi}_{1}=(\varphi_{1,0},\varphi_{1,1},\ldots,\varphi_{1,N})^{\rm T}$
and
$\bm{\varphi}_{2}=(\varphi_{2,0},\varphi_{2,1},\ldots,\varphi_{2,N^{*}})^{\rm
T}$. Then, we see that $(\bm{\varphi}_{1},\bm{\varphi}_{2})$ is a solution to
(6.9). More precisely, we obtain the following lemma.

######  Lemma 6.4.

Under the hypothesis of Lemma 6.3, for any
$\bm{f}_{1}^{\prime}=(f_{1,1},\ldots,f_{1,N})^{\rm T}$,
$\bm{f}_{2}^{\prime}=(f_{2,1},\ldots,$ $f_{2,N^{*}})^{\rm T}$, $\bm{f}_{3}$,
and $f_{4}$ satisfying $\bm{f}_{1}^{\prime},\bm{f}_{2}^{\prime}\in H^{k-1}$
and $\bm{f}_{3},\nabla f_{4}\in H^{k}$ with $k\in\\{0,1,\ldots,m-1\\}$, there
exists a solution $(\bm{\varphi}_{1},\bm{\varphi}_{2})$ to (6.9) satisfying

$\|(\bm{\varphi}_{1}^{\prime},\bm{\varphi}_{2}^{\prime})\|_{H^{k+1}}+\|(\nabla\varphi_{1,0},\nabla\varphi_{2,0})\|_{H^{k}}\leq
C\bigl{(}\|(\bm{f}_{1}^{\prime},\bm{f}_{2}^{\prime})\|_{H^{k-1}}+\|(\bm{f}_{3},\nabla
f_{4})\|_{H^{k}}\bigr{)},$

where $C=C(\rho_{1},\rho_{2},h_{1},h_{2},c_{0},m)$. Moreover, the solution is
unique up to an additive constant of the form
$(\mathcal{C}\rho_{2},\mathcal{C}\rho_{1})$ to
$(\varphi_{1,0},\varphi_{2,0})$.

##  7 Construction of the solution

In this section, we will prove Theorem 2.1 one of the main theorems in this
paper. One possible strategy to construct the solution of the initial value
problem to the Kakinuma model (2.14)–(2.16) would consist in firstly
transforming the equations into a quasilinear positive symmetric system, that
is, a quasilinear version of the positive symmetric system (5.11), secondly
applying the method of parabolic regularization to construct the solution of
the transformed system, and finally to show that the solution to the
transformed system is in fact the solution of the Kakinuma model if we further
impose the compatibility conditions (2.18)–(2.18) on the initial data. Here,
in order to avoid the heavy computations that would be involved when following
this strategy, we find it more convenient to instead apply the method of
parabolic regularization to the Kakinuma model directly.

###  7.1 Parabolic regularization of the equations

We remind that the Kakinuma model (2.14)–(2.16) can be written compactly as
(6.4), that is,

(7.1)
$\begin{cases}\displaystyle\bm{l}_{1}(H_{1})\partial_{t}\zeta+L_{1}(H_{1})\bm{\phi}_{1}=\bm{0},\\\
\displaystyle-\bm{l}_{2}(H_{2})\partial_{t}\zeta+L_{2}(H_{2},b)\bm{\phi}_{2}=\bm{0},\\\
\displaystyle-\rho_{1}\bm{l}_{1}(H_{1})\cdot\partial_{t}\bm{\phi}_{1}+\rho_{2}\bm{l}_{2}(H_{2})\cdot\partial_{t}\bm{\phi}_{2}=F,\end{cases}$

where $\bm{\phi}_{1}=(\phi_{1,0},\phi_{1,1},\ldots,\phi_{1,N})^{\rm T}$,
$\bm{\phi}_{2}=(\phi_{2,0},\phi_{2,1},\ldots,\phi_{2,N^{*}})^{\rm T}$,
$\bm{l}_{k}$ and $L_{k}$ for $k=1,2$ are defined in (6.3), and

(7.2)
$F=\rho_{1}\left\\{g\zeta+\frac{1}{2}\left(|\bm{u}_{1}|^{2}+w_{1}^{2}\right)\right\\}-\rho_{2}\left\\{g\zeta+\frac{1}{2}\left(|\bm{u}_{2}|^{2}+w_{2}^{2}\right)\right\\}.$

Here $\bm{u}_{k}$ and $w_{k}$ for $k=1,2$ are defined by (4.2) and (4.3)
respectively. We regularize the Kakinuma model by adding artificial viscosity
terms as

(7.3)
$\begin{cases}\displaystyle\bm{l}_{1}(H_{1})(\partial_{t}\zeta-\varepsilon\Delta\zeta)+L_{1}(H_{1})\bm{\phi}_{1}=\bm{0},\\\
\displaystyle-\bm{l}_{2}(H_{2})(\partial_{t}\zeta-\varepsilon\Delta\zeta)+L_{2}(H_{2},b)\bm{\phi}_{2}=\bm{0},\\\
\displaystyle-\rho_{1}\bm{l}_{1}(H_{1})\cdot(\partial_{t}\bm{\phi}_{1}-\varepsilon\Delta\bm{\phi}_{1})+\rho_{2}{\bm{l}}_{2}(H_{2})\cdot(\partial_{t}\bm{\phi}_{2}-\varepsilon\Delta\bm{\phi}_{2})=F.\end{cases}$

We are going to show the existence of the solution to the initial value
problem for this regularized Kakinuma model under the initial conditions

(7.4)
$(\zeta,\bm{\phi}_{1},\bm{\phi}_{2})|_{t=0}=(\zeta_{(0)},\bm{\phi}_{1(0)},\bm{\phi}_{2(0)}).$

For this regularized Kakinuma model, the compatibility conditions for the
existence of the solution have the same form as the original Kakinuma model,
that is,

(7.5)
$\begin{cases}\mathcal{L}_{1,i}(H_{1})\bm{\phi}_{1}=0\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
\mathcal{L}_{2,i}(H_{2},b)\bm{\phi}_{2}=0\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\
\mathcal{L}_{1,0}(H_{1})\bm{\phi}_{1}+\mathcal{L}_{2,0}(H_{2},b)\bm{\phi}_{2}=0,\end{cases}$

where $\mathcal{L}_{1,i}(H_{1})$ for $i=0,1,\ldots,N$ and
$\mathcal{L}_{2,i}(H_{2},b)$ for $i=0,1,\ldots,N^{*}$ are defined in (6.5).
Here, we note the identities

$\begin{cases}[\partial_{t},\mathcal{L}_{1,i}(H_{1})]\bm{\phi}_{1}=f_{1,i}(\zeta,\bm{\phi}_{1})\partial_{t}\zeta\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
[\partial_{t},\mathcal{L}_{2,i}(H_{2},b)]\bm{\phi}_{2}=f_{2,i}(\zeta,\bm{\phi}_{2},b)\partial_{t}\zeta\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\
[\partial_{t},\mathcal{L}_{1,0}(H_{1})]\bm{\phi}_{1}+[\partial_{t},\mathcal{L}_{2,0}(H_{2},b)]\bm{\phi}_{2}=-\nabla\cdot(\bm{v}\partial_{t}\zeta),\end{cases}$

where $\bm{v}=\bm{u}_{2}-\bm{u}_{1}$ and

$\begin{cases}\displaystyle
f_{1,i}(\zeta,\bm{\phi}_{1})\displaystyle=-\sum_{j=0}^{N}\left\\{\frac{2i}{2j+1}H_{1}^{2(i+j)}\Delta\phi_{1,j}+4ijH_{1}^{2(i+j-1)}\phi_{1,j}\right\\},\\\
\displaystyle
f_{2,i}(\zeta,\bm{\phi}_{2},b)\displaystyle=\sum_{j=0}^{N^{*}}\left\\{\frac{p_{i}}{p_{j}+1}H_{2}^{p_{i}+p_{j}}\Delta\phi_{2,j}-\frac{p_{i}p_{j}}{p_{j}}H_{2}^{p_{i}+p_{j}-1}\nabla\cdot(\phi_{2,j}\nabla
b)\right.\\\
\displaystyle\phantom{f_{2,i}(\zeta,\bm{\phi}_{2},b)=\sum\Big{\\{}}-p_{i}H_{2}^{p_{i}+p_{j}-1}\nabla
b\cdot\nabla\phi_{2,j}+p_{i}p_{j}H_{2}^{p_{i}+p_{j}-2}(1+|\nabla
b|^{2})\phi_{2,j}\biggr{\\}},\end{cases}$

and

$\begin{cases}[\Delta,\mathcal{L}_{1,i}(H_{1})]\bm{\phi}_{1}=f_{1,i}(\zeta,\bm{\phi}_{1})\Delta\zeta+\widetilde{f}_{1,i}(\zeta,\bm{\phi}_{1})\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
[\Delta,\mathcal{L}_{2,i}(H_{2},b)]\bm{\phi}_{2}=f_{2,i}(\zeta,\bm{\phi}_{2},b)\Delta\zeta+\widetilde{f}_{2,i}(\zeta,\bm{\phi}_{2},b)\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\
[\Delta,\mathcal{L}_{1,0}(H_{1})]\bm{\phi}_{1}+[\Delta,\mathcal{L}_{2,0}(H_{2},b)]\bm{\phi}_{2}=-\nabla\cdot(\bm{v}\Delta\zeta)+f_{3}(\zeta,\bm{\phi}_{1},\bm{\phi}_{2},b),\end{cases}$

where

$\begin{cases}\displaystyle\widetilde{f}_{1,i}(\zeta,\bm{\phi}_{1})=\sum_{l=1}^{n}\left\\{[\partial_{l},\mathcal{L}_{1,i}(H_{1})]\partial_{l}\bm{\phi}_{1}+(\partial_{l}\zeta)\partial_{l}f_{1,i}(\zeta,\bm{\phi}_{1})\right\\},\\\
\displaystyle\widetilde{f}_{2,i}(\zeta,\bm{\phi}_{2},b)=\sum_{l=1}^{n}\biggl{\\{}[\partial_{l},\mathcal{L}_{2,i}(H_{2},b)]\partial_{l}\bm{\phi}_{2}+(\partial_{l}\zeta)f_{2,i}(\zeta,\bm{\phi}_{2},b)-\partial_{l}((\partial_{l}b)f_{2,i}(\zeta,\bm{\phi}_{2},b))\\\
\displaystyle\phantom{\widetilde{f}_{2,i}(\zeta,\bm{\phi}_{2},b)=\sum_{l=1}^{n}\biggl{\\{}}+\sum_{j=0}^{N^{*}}\partial_{l}\biggl{(}-\frac{p_{i}p_{j}}{(p_{i}+p_{j})p_{j}}H_{2}^{p_{i}+p_{j}}\nabla\cdot(\phi_{2,j}\nabla\partial_{l}b)-\frac{p_{i}}{p_{i}+p_{j}}H_{2}^{p_{i}+p_{j}}\nabla\partial_{l}b\cdot\nabla\phi_{2,j}\\\
\displaystyle\phantom{\widetilde{f}_{2,i}(\zeta,\bm{\phi}_{2},b)=\sum_{l=1}^{n}\biggl{\\{}+\sum_{j=0}^{N^{*}}\partial_{l}\biggl{(}}+\frac{p_{i}p_{j}}{p_{i}+p_{j}-1}H_{2}^{p_{i}+p_{j}-1}2(\nabla
b\cdot\nabla\partial_{l}b)\biggr{)}\biggr{\\}},\\\ \displaystyle
f_{3}(\zeta,\bm{\phi}_{1},\bm{\phi}_{2},b)=\sum_{l=1}^{n}\biggl{\\{}[\partial_{l},\mathcal{L}_{1,0}(H_{1})]\partial_{l}\bm{\phi}_{1}+[\partial_{l},\mathcal{L}_{2,0}(H_{2},b)]\partial_{l}\bm{\phi}_{2}\\\
\displaystyle\phantom{f_{3}(\zeta,\bm{\phi}_{1},\bm{\phi}_{2},b)=\sum_{l=1}^{n}\biggl{\\{}}+\nabla\cdot\biggl{(}-(\partial_{l}\zeta)(\partial_{l}\bm{v})+\partial_{l}\biggl{(}(\partial_{l}b)\bm{u}_{2}+\sum_{j=1}^{N^{*}}H_{2}^{p_{j}}\phi_{2,j}\nabla\partial_{l}b\biggr{)}\biggr{)}\biggr{\\}}.\end{cases}$

We also note that $f_{3}(\zeta,\bm{\phi}_{1},\bm{\phi}_{2},b)$ can be written
in a divergence form as

$f_{3}(\zeta,\bm{\phi}_{1},\bm{\phi}_{2},b)=\nabla\cdot\bm{f}_{3}(\zeta,\bm{\phi}_{1},\bm{\phi}_{2},b),$

where

$\displaystyle\bm{f}_{3}(\zeta,\bm{\phi}_{1},\bm{\phi}_{2},b)$
$\displaystyle=\sum_{l=1}^{n}\left\\{(\partial_{l}\zeta)\sum_{j=0}^{N}H_{1}^{2j}\nabla\partial_{l}\phi_{1,j}+\sum_{j=1}^{N^{*}}H_{2}^{p_{j}}(\partial_{l}\phi_{2,j})\nabla\partial_{l}b\right.$
$\displaystyle\phantom{=\sum_{l=1}^{n}\Bigl{\\{}}+(\partial_{l}b-\partial_{l}\zeta)\sum_{j=0}^{N^{*}}\bigl{(}H_{2}^{p_{j}}\nabla\partial_{l}\phi_{2,j}-p_{j}H_{2}^{p_{j}-1}(\partial_{l}\phi_{2,j})\nabla
b\bigr{)}$
$\displaystyle\phantom{=\sum_{l=1}^{n}\Bigl{\\{}}\left.-(\partial_{l}\zeta)(\partial_{l}\bm{v})+\partial_{l}\left((\partial_{l}b)\bm{u}_{2}+\sum_{j=1}^{N^{*}}H_{2}^{p_{j}}\phi_{2,j}\nabla\partial_{l}b\right)\right\\}.$

Therefore, applying the operator $\partial_{t}-\varepsilon\Delta$ to (7.5) we
obtain

(7.6)
$\begin{cases}\mathcal{L}_{1,i}(H_{1})(\partial_{t}\bm{\phi}_{1}-\varepsilon\Delta\bm{\phi}_{1})=-f_{1,i}(\zeta,\bm{\phi}_{1})(\partial_{t}\zeta-\varepsilon\Delta\zeta)+\varepsilon\widetilde{f}_{1,i}(\zeta,\bm{\phi}_{1})\\\
\makebox[180.00027pt]{}\mbox{for}\quad i=1,2,\ldots,N,\\\
\mathcal{L}_{2,i}(H_{2},b)(\partial_{t}\bm{\phi}_{2}-\varepsilon\Delta\bm{\phi}_{2})=-f_{2,i}(\zeta,\bm{\phi}_{2},b)(\partial_{t}\zeta-\varepsilon\Delta\zeta)+\varepsilon\widetilde{f}_{2,i}(\zeta,\bm{\phi}_{2},b)\\\
\makebox[180.00027pt]{}\mbox{for}\quad i=1,2,\ldots,N^{*},\\\
\mathcal{L}_{1,0}(H_{1})(\partial_{t}\bm{\phi}_{1}-\varepsilon\Delta\bm{\phi}_{1})+\mathcal{L}_{2,0}(H_{2},b)(\partial_{t}\bm{\phi}_{2}-\varepsilon\Delta\bm{\phi}_{2})\\\
\makebox[100.00015pt]{}=\nabla\cdot\bigl{(}\bm{v}(\partial_{t}\zeta-\varepsilon\Delta\zeta)+\varepsilon\bm{f}_{3}(\zeta,\bm{\phi}_{1},\bm{\phi}_{2},b)\bigr{)}.\end{cases}$

On the other hand, we have $N+N^{*}+2$ evolution equations for one scalar
function $\zeta$. To select an appropriate evolution equation for $\zeta$, we
will use the notation defined by (5.13). We note that they depend on the
unknown functions $H_{1}$ and $H_{2}$. Taking Euclidean inner products of the
first and the second equations in (7.3) with $\rho_{1}\bm{q}_{1}$ and
$\rho_{2}\bm{q}_{2}$, respectively, adding the resulting equations, and using
the relation
$-\rho_{1}\bm{l}_{1}\cdot\bm{q}_{1}+\rho_{2}\bm{l}_{2}\cdot\bm{q}_{2}=1$, we
obtain

(7.7) $\partial_{t}\zeta-\varepsilon\Delta\zeta=G_{0},$

where

$G_{0}=\rho_{1}\bm{q}_{1}\cdot
L_{1}(H_{1})\bm{\phi}_{1}+\rho_{2}\bm{q}_{2}\cdot
L_{2}(H_{2},b)\bm{\phi}_{2}.$

Plugging this into (7.6) and noting the last equation in (7.3), we have

(7.8)
$\begin{cases}\mathcal{L}_{1,i}(H_{1})(\partial_{t}\bm{\phi}_{1}-\varepsilon\Delta\bm{\phi}_{1})=-f_{1,i}(\zeta,\bm{\phi}_{1})G_{0}+\varepsilon\widetilde{f}_{1,i}(\zeta,\bm{\phi}_{1})\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
\mathcal{L}_{2,i}(H_{2},b)(\partial_{t}\bm{\phi}_{2}-\varepsilon\Delta\bm{\phi}_{2})=-f_{2,i}(\zeta,\bm{\phi}_{2},b)G_{0}+\varepsilon\widetilde{f}_{2,i}(\zeta,\bm{\phi}_{2},b)\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\
\mathcal{L}_{1,0}(H_{1})(\partial_{t}\bm{\phi}_{1}-\varepsilon\Delta\bm{\phi}_{1})+\mathcal{L}_{2,0}(H_{2},b)(\partial_{t}\bm{\phi}_{2}-\varepsilon\Delta\bm{\phi}_{2})\\\
\makebox[100.00015pt]{}=\nabla\cdot\bigl{(}\bm{v}G_{0}+\varepsilon\bm{f}_{3}(\zeta,\bm{\phi}_{1},\bm{\phi}_{2},b)\bigr{)},\\\
-\rho_{1}\bm{l}_{1}(H_{1})\cdot(\partial_{t}\bm{\phi}_{1}-\varepsilon\Delta\bm{\phi}_{1})+\rho_{2}{\bm{l}}_{2}(H_{2})\cdot(\partial_{t}\bm{\phi}_{2}-\varepsilon\Delta\bm{\phi}_{2})=F.\end{cases}$

Therefore, thanks to Lemma 6.4 we obtain

(7.9)
$\begin{cases}\partial_{t}\bm{\phi}_{1}-\varepsilon\Delta\bm{\phi}_{1}=\bm{G}_{1},\\\
\partial_{t}\bm{\phi}_{2}-\varepsilon\Delta\bm{\phi}_{2}=\bm{G}_{2},\end{cases}$

where $\bm{G}_{1}=(G_{1,0},G_{1,1},\ldots,G_{1,N})^{\rm T}$ and
$\bm{G}_{2}=(G_{2,0},G_{2,1},\ldots,G_{2,N^{*}})^{\rm T}$ are defined as a
solution to the following equations

(7.10)
$\begin{cases}\mathcal{L}_{1,i}(H_{1})\bm{G}_{1}=-f_{1,i}(\zeta,\bm{\phi}_{1})G_{0}+\varepsilon\widetilde{f}_{1,i}(\zeta,\bm{\phi}_{1})\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
\mathcal{L}_{2,i}(H_{2},b)\bm{G}_{2}=-f_{2,i}(\zeta,\bm{\phi}_{2},b)G_{0}+\varepsilon\widetilde{f}_{2,i}(\zeta,\bm{\phi}_{2},b)\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\
\mathcal{L}_{1,0}(H_{1})\bm{G}_{1}+\mathcal{L}_{2,0}(H_{2},b)\bm{G}_{2}=\nabla\cdot\left(\bm{v}G_{0}+\varepsilon\bm{f}_{3}(\zeta,\bm{\phi}_{1},\bm{\phi}_{2},b)\right),\\\
-\rho_{1}\bm{l}_{1}(H_{1})\cdot\bm{G}_{1}+\rho_{2}\bm{l}_{2}(H_{2})\cdot\bm{G}_{2}=F.\end{cases}$

Precisely speaking, $(\bm{G}_{1},\bm{G}_{2})$ are defined uniquely up to an
additive constant of the form $(\mathcal{C}\rho_{2},\mathcal{C}\rho_{1})$ to
$(G_{1,0},G_{2,0})$. However, this indeterminacy does not cause any
difficulties in the following arguments.

######  Remark 7.1.

The equations in (7.9) are valid even in the case $\varepsilon=0$, that is,
any regular solutions to the Kakinuma model (2.14)–(2.15) satisfy (7.9) with
$\varepsilon=0$. Particularly, $\partial_{t}\bm{\phi}_{k}(\bm{x},0)$ for
$k=1,2$ can be expressed in term of the initial data
$(\zeta_{(0)},\bm{\phi}_{1(0)},\bm{\phi}_{2(0)})$ and the bottom topography
$b$.

###  7.2 Existence of the solution to the regularized problem

######  Lemma 7.2.

Let $g,\rho_{1},\rho_{2},h_{1},h_{2},c_{0}$ be positive constants and $m$ an
integer such that $m>\frac{n}{2}+1$. For any initial data
$(\zeta_{(0)},\bm{\phi}_{1(0)},\bm{\phi}_{2(0)})$ and bottom topography $b$
satisfying

$\begin{cases}\zeta_{(0)},\nabla\phi_{1,0(0)},\nabla\phi_{2,0(0)}\in
H^{m},\quad\bm{\phi}_{1(0)}^{\prime},\bm{\phi}_{2(0)}^{\prime}\in
H^{m+1},\quad b\in W^{m+2,\infty},\\\ h_{1}-\zeta_{(0)}(\bm{x})\geq
c_{0},\quad h_{2}+\zeta_{(0)}(\bm{x})-b(\bm{x})\geq
c_{0}\quad\mbox{for}\quad\bm{x}\in\mathbf{R}^{n},\end{cases}$

and for any $\varepsilon>0$ there exists a maximal existence time
$T_{\varepsilon}\in(0,+\infty]$ such that the initial value problem (7.7),
(7.9), and (7.4) has a unique solution
$(\zeta^{\varepsilon},{\bm{\phi}}_{1}^{\varepsilon},{\bm{\phi}}_{2}^{\varepsilon})$
satisfying

$\zeta^{\varepsilon},\nabla\phi_{1,0}^{\varepsilon},\nabla\phi_{2,0}^{\varepsilon}\in
C([0,T_{\varepsilon});H^{m}),\qquad{\bm{\phi}}_{1}^{\varepsilon\prime},{\bm{\phi}}_{2}^{\varepsilon\prime}\in
C([0,T_{\varepsilon});H^{m+1}).$

###### Proof..

We evaluate the right-hand sides of the equations, that is, the terms $G_{0}$,
$\bm{G}_{1}$, and $\bm{G}_{2}$. To this end, suppose that
$(\zeta,\bm{\phi}_{1},\bm{\phi}_{2})$ and $b$ satisfy

(7.11)
$\begin{cases}\|(\zeta,\nabla\phi_{1,0},\nabla\phi_{2,0})\|_{H^{m}}+\|(\bm{\phi}_{1}^{\prime},\bm{\phi}_{2}^{\prime})\|_{H^{m+1}}+\|b\|_{W^{m+2,\infty}}\leq
M,\\\ h_{1}-\zeta(\bm{x})\geq c_{1},\quad h_{2}+\zeta(\bm{x})-b(\bm{x})\geq
c_{1}\quad\mbox{for}\quad\bm{x}\in\mathbf{R}^{n}.\end{cases}$

Then, we see that

$\|G_{0}\|_{H^{m-1}}+\|(\bm{f}_{1}^{\prime},\bm{f}_{2}^{\prime},\bm{f}_{3})\|_{H^{m-1}}+\|(\widetilde{\bm{f}}_{1}^{\prime},\widetilde{\bm{f}}_{2}^{\prime})\|_{H^{m-2}}+\|F\|_{H^{m}}\leq
C(M,c_{1}),$

where
$\bm{f}_{1}^{\prime}=(f_{1,1}(\zeta,\bm{\phi}_{1}),\ldots,f_{1,N}(\zeta,\bm{\phi}_{1}))$
and so on. Therefore, by Lemma 6.4 we have

$\|(\nabla G_{1,0},\nabla
G_{2,0})\|_{H^{m-1}}+\|({\bm{G}}_{1}^{\prime},{\bm{G}}_{2}^{\prime})\|_{H^{m}}\leq
C(M,c_{1},\varepsilon),$

where we notice for further use that $C(M,c_{1},\varepsilon)$ is bounded
uniformly with respect to $\varepsilon\in(0,1]$. We obtain the desired result
by the standard theory on the heat equation. ∎

######  Lemma 7.3.

Suppose that the initial data
$(\zeta_{(0)},\bm{\phi}_{1(0)},\bm{\phi}_{2(0)})$ and the bottom topography
$b$ satisfy the hypotheses in Lemma 7.2 and the compatibility conditions
(7.5). Then, the solution
$(\zeta^{\varepsilon},\bm{\phi}_{1}^{\varepsilon},\bm{\phi}_{2}^{\varepsilon})$
constructed in Lemma 7.2 satisfies the regularized Kakinuma model (7.3).

###### Proof..

By the construction of the solution, we easily see that it satisfies (7.8) and
in particular the last equation in (7.3). Therefore, it is sufficient to show
that it satisfies also the first two equations in (7.3). By (7.7) and (7.8),
we have

$\begin{cases}(\partial_{t}-\varepsilon\Delta)(\mathcal{L}_{1,i}(H_{1})\bm{\phi}_{1})=0\quad\mbox{for}\quad
i=1,2,\ldots,N,\\\
(\partial_{t}-\varepsilon\Delta)(\mathcal{L}_{2,i}(H_{2},b)\bm{\phi}_{2})=0\quad\mbox{for}\quad
i=1,2,\ldots,N^{*},\\\
(\partial_{t}-\varepsilon\Delta)\left(\mathcal{L}_{1,0}(H_{1})\bm{\phi}_{1}+\mathcal{L}_{2,0}(H_{2},b)\bm{\phi}_{2}\right)=0,\end{cases}$

so that by the uniqueness of the solution to the initial value problem of the
heat equation, if the initial data satisfy the compatibility conditions (7.5),
then the solution satisfies also (7.5) for all $t\in[0,T_{\varepsilon})$.
Particularly, we obtain

$\begin{cases}-\bm{l}_{1}(H_{1})(\mathcal{L}_{1,0}(H_{1})\bm{\phi}_{1})+L_{1}(H_{1})\bm{\phi}_{1}=\bm{0},\\\
-\bm{l}_{2}(H_{2})(\mathcal{L}_{2,0}(H_{2},b)\bm{\phi}_{2})+L_{2}(H_{2},b)\bm{\phi}_{2}=\bm{0},\end{cases}$

so that by the last equation in the compatibility conditions (7.5) we have

(7.12)
$\begin{cases}-\bm{l}_{1}(H_{1})(\mathcal{L}_{1,0}(H_{1})\bm{\phi}_{1})+L_{1}(H_{1})\bm{\phi}_{1}=\bm{0},\\\
\bm{l}_{2}(H_{2})(\mathcal{L}_{1,0}(H_{1})\bm{\phi}_{1})+L_{2}(H_{2},b)\bm{\phi}_{2}=\bm{0}.\end{cases}$

Taking Euclidean inner products of the first and the second equations with
$\rho_{1}\bm{q}_{1}$ and $\rho_{2}\bm{q}_{2}$, respectively, adding the
resulting equations, and using the relation
$-\rho_{1}\bm{l}_{1}\cdot\bm{q}_{1}+\rho_{2}\bm{l}_{2}\cdot\bm{q}_{2}=1$, we
obtain

$\mathcal{L}_{1,0}(H_{1})\bm{\phi}_{1}+\rho_{1}\bm{q}_{1}\cdot
L_{1}(H_{1})\bm{\phi}_{1}+\rho_{2}\bm{q}_{2}\cdot
L_{2}(H_{2},b)\bm{\phi}_{2}=0,$

which together with (7.7) implies

$\mathcal{L}_{1,0}(H_{1})\bm{\phi}_{1}=-(\partial_{t}\zeta-\varepsilon\Delta\zeta).$

Plugging this into (7.12), we see that the solution satisfies the first two
equations in (7.3). ∎

###  7.3 Uniform bound of the solution to the regularized problem

We proceed to derive estimates on solutions
$(\zeta^{\varepsilon},\bm{\phi}_{1}^{\varepsilon},\bm{\phi}_{2}^{\varepsilon})$
to the regularized Kakinuma model (7.3) uniform with respect to the
regularized parameter $\varepsilon\in(0,1]$ and for a time interval
independent of $\varepsilon$. To this end, we make use of a good symmetric
structure of the Kakinuma model based on the analysis of Section 5.1. In order
to simplify the notation we write $(\zeta,\bm{\phi}_{1},\bm{\phi}_{2})$ in
place of
$(\zeta^{\varepsilon},\bm{\phi}_{1}^{\varepsilon},\bm{\phi}_{2}^{\varepsilon})$.

In view of (6.1) and (6.2) we decompose $L_{1}(H_{1})\bm{\phi}_{1}$ and
$L_{2}(H_{2},b)\bm{\phi}_{2}$ into their principal parts and the remainder
parts as

(7.13) $\displaystyle
L_{1}(H_{1})\bm{\phi}_{1}=-A_{1}(H_{1})\Delta\bm{\phi}_{1}+\bm{l}_{1}(H_{1})(\bm{u}_{1}\cdot\nabla\zeta)+L_{1}^{\rm
low}(H_{1})\bm{\phi}_{1},$ (7.14) $\displaystyle
L_{2}(H_{2},b)\bm{\phi}_{2}=-A_{2}(H_{2})\Delta\bm{\phi}_{2}-\bm{l}_{2}(H_{2})(\bm{u}_{2}\cdot\nabla\zeta)+L_{2}^{\rm
low}(H_{2},b)\bm{\phi}_{2},$

where the matrices $A_{1}(H_{1})$ and $A_{2}(H_{2})$ are given by (5.10),