A Friedlander type estimate for Stokes operators

C. Denis and A.F.M. ter Elst

  * Abstract. Let $\Omega\subset\mathds{R}^{d}$ be a bounded open connected set with Lipschitz boundary. Let $A^{N}$ and $A^{D}$ be the Neumann Stokes operator and Dirichlet Stokes operator on $\Omega$, respectively. Further let $\lambda_{1}^{N}\leq\lambda_{2}^{N}\leq\ldots$ and $\lambda_{1}^{D}\leq\lambda_{2}^{D}\leq\ldots$ be the eigenvalues of $A^{N}$ and $A^{D}$ repeated with multiplicity, respectively. Then

$\lambda_{n+1}^{N}<\lambda_{n}^{D}$

for all $n\in\mathds{N}$.

††footnotetext:  Mathematics Subject Classification. 47A75, 35P15. Keywords.
Eigenvalues, counting function, Stokes operator.

##  1 Introduction

If $\lambda_{1}^{D}\leq\lambda_{2}^{D}\leq\ldots$ are the eigenvalues of the
Dirichlet Laplacian and $\lambda_{1}^{N}\leq\lambda_{2}^{N}\leq\ldots$ the
eigenvalues of the Neumann Laplacian on a domain $\Omega$ with Lipschitz
boundary, then

$\lambda_{n+1}^{N}<\lambda_{n}^{D}$ (1)

for all $n\in\mathds{N}$. This result has a long history, starting with Payne
[Pay] and Pólya [Pól]. In 1991, Friedlander [Fri] proved that
$\lambda_{n+1}^{N}\leq\lambda_{n}^{D}$ for a $C^{1}$ domain, using the
Dirichlet-to-Neumann operator. The method was modified and extended by Filonov
[Fil] to prove a strict inequality on more general domains. A different and
independent proof was obtained by Arendt–Mazzeo [AM2]. On the other hand,
Mazzeo [Maz] gave a counter example for manifolds. Under certain geometric
conditions on three dimensional domains Hansson [Han] proved (1) for Dirichlet
and Neumann eigenvalues of the Heisenberg Laplacian. This was extended to
general Heisenberg groups and fairly general domains by Frank–Laptev [FL].
Other results are for powers $(-\Delta)^{m}$ of the Laplacian by Provenzano
[Pro] and a comparison between Robin and Dirichlet eigenvalues by
Gesztesy–Mitrea [GM] or for certain operators with mixed boundary conditions
by Lotoreichik–Rohleder [LR]. In [BRS] Behrndt–Rohleder–Stadler considered
external domains and added a potential to the Laplacian.

The aim of this paper is to prove a similar result for the Stokes operator.
Let $\Omega\subset\mathds{R}^{d}$ be an open bounded connected set. One can
define two versions of the Stokes operator: a Dirichlet version and a Neumann
version. The Dirichlet Stokes operator $A^{D}$ is defined in the closure
$L_{2,\sigma}(\Omega)$ of solenoidal test functions $\\{u\in
C_{c}^{\infty}(\Omega,\mathds{C}^{d}):\mathop{\rm div}u=0\\}$ in
$L_{2}(\Omega,\mathds{C}^{d})$. If $u,f\in L_{2,\sigma}(\Omega)$, then $u\in
D(A^{D})$ and $A^{D}u=f$ if and only if $u\in
H^{1}_{0}(\Omega,\mathds{C}^{d})$, $\mathop{\rm div}u=0$ and there exists a
$\pi\in L_{2}(\Omega)$ such that $f=-\Delta u+\nabla\pi$ in
$H^{-1}(\Omega,\mathds{C}^{d})$. Then $A^{D}$ is a positive self-adjoint
operator in $L_{2,\sigma}(\Omega)$ with discrete spectrum. Let
$\lambda_{1}^{D}\leq\lambda_{2}^{D}\leq\ldots$ be the eigenvalues of $A^{D}$
repeated with multiplicity.

The Neumann Stokes operator $A^{N}$ is defined in the Hilbert space $\\{u\in
L_{2}(\Omega,\mathds{C}^{d}):\mathop{\rm div}u=0\\}$. In Proposition 3.5 we
shall give a precise definition of the self-adjoint operator $A^{N}$. Loosely
speaking, $u\in D(A^{N})$ and $A^{N}u=f$ if $u\in
H^{1}(\Omega,\mathds{C}^{d})$, $\mathop{\rm div}u=0$ and there exists a
$\pi\in L_{2}(\Omega)$ such that $f=-\Delta u+\nabla\pi$ in
$H^{-1}(\Omega,\mathds{C}^{d})$ and $\partial_{\nu}u=({\mathop{\rm
Tr\,}}\pi)\nu$, where $\partial_{\nu}$ is the outward normal derivative. Also
$A^{N}$ is a positive self-adjoint operator with discrete spectrum. Let
$\lambda_{1}^{N}\leq\lambda_{2}^{N}\leq\ldots$ be the eigenvalues of $A^{D}$
repeated with multiplicity. The main result of this paper is that

$\lambda_{n+1}^{N}<\lambda_{n}^{D}$ (2)

for all $n\in\mathds{N}$.

Although the Dirichlet and Neumann Laplacians are self-adjoint operators in
the same Hilbert space $L_{2}(\Omega)$, this is not the case for the Stokes
Dirichlet operator $A^{D}$ and Stokes Neumann operator $A^{N}$. The Stokes
Dirichlet operator $A^{D}$ is a self-adjoint operator in the orthogonal
complement in $L_{2}(\Omega,\mathds{C}^{d})$ of the space

$\\{\nabla\pi:\pi\in H^{1}(\Omega)\\},$

whilst the Stokes Neumann operator $A^{N}$ is a self-adjoint operator in the
orthogonal complement in $L_{2}(\Omega,\mathds{C}^{d})$ of the space

$\\{\nabla\pi:\pi\in H^{1}_{0}(\Omega)\\}.$

Nevertheless, we relate the eigenvalues in (2).

As in Friedlander’s original paper, the main ingredient for our proof is the
Dirichlet-to-Neumann operator, adapted to the Stokes problem. We use
sesquilinear form methods as in [AM1] and [AM2]. In Section 2 we prove an
abstract theorem, Theorem 2.2, on eigenvalue comparison of Dirichlet-type and
Neumann-type operators associated with forms, which relies (and expands) on a
combination of the papers [AM1], [AM2], [AEKS] and the BSc(Honours) thesis
[Gor]. As an immediate application we obtain in Theorem 2.8 the classical
result for the Laplacian on a Lipschitz domain. Section 3 begins with the
description of the Dirichlet and Neumann versions of the Stokes operator on a
Lipschitz domain, as well as of the related Dirichlet-to-Neumann operator. We
then apply our abstract eigenvalue theorem, Theorem 2.2, to obtain a Stokes
version of the Friedlander eigenvalue theorem.

##  2 An abstract Friedlander comparison theorem

We recall the basic notions of form methods as developed in [AE] and [AEKS],
which extend the classical form methods of Kato [Kat] and Lions [Lio]. Since
we only use self-adjoint graphs in this paper, we restrict our attention to
symmetric forms.

Let $V$, $K$ be Hilbert spaces and $\gothic{a}\colon V\times V\to\mathds{C}$
be a sesquilinear form. The form $\gothic{a}$ is called continuous if there is
an $M>0$ such that $|\gothic{a}(u,v)|\leq M\,\|u\|_{V}\,\|v\|_{V}$ for all
$u,v\in V$. The form $\gothic{a}$ is called positive if $\gothic{a}(u,u)\geq
0$ for all $u\in V$ and it is called symmetric if
$\gothic{a}(v,u)=\overline{\gothic{a}(u,v)}$ for all $u,v\in V$. Let $j\colon
V\to K$ be a continuous linear map. In Lions [Lio], the map $j$ is an
inclusion map. We say that $\gothic{a}$ is $j$-elliptic if there are
$\omega,\mu>0$ such that

$\mathop{\rm
Re}\gothic{a}(u,u)+\omega\,\|j(u)\|_{K}^{2}\geq\mu\,\|u\|_{V}^{2}$

for all $u\in V$. The graph associated with $(\gothic{a},j)$ is the subspace

$\\{(j(u),f)\in K\times K:u\in V\mbox{ and }\gothic{a}(u,v)=(f,j(v))_{K}\mbox{
for all }v\in V\\}$

in $K\times K$.

There are a number of sufficient conditions to conclude that the graph
associated with $(\gothic{a},j)$ is a self-adjoint graph if $\gothic{a}$ is
continuous and symmetric. For a summary of the terminology in self-adjoint
graphs we refer to [AEKS] Section 3 and [BE] Section 2. In this paper we use
the following two sufficient conditions.

######  Proposition 2.1.

Let $V$, $K$ be Hilbert spaces and $\gothic{a}\colon V\times V\to\mathds{C}$
be a continuous symmetric sesquilinear form. Let $j\colon V\to K$ be a
continuous linear map. Suppose at least one of the following conditions is
valid.

  * (a)

The form $\gothic{a}$ is $j$-elliptic.

  * (b)

There exists a Hilbert space $H$ and a compact linear map $i\colon V\to H$
such that $\gothic{a}$ is $i$-elliptic.

Then the graph $A$ associated with $(\gothic{a},j)$ is a self-adjoint graph
which is lower bounded.

Moreover, if $j$ is compact, then $A$ has compact resolvent and in particular
a discrete spectrum.

###### Proof.

‘(a)’. This follows from [AE] Theorem 2.1 applied to the Hilbert space
$\overline{j(V)}$.

‘(b)’. See [AEKS] Theorems 4.5 and 4.15.

For the last statement, see [AEKS] Proposition 4.8. ∎

Now we are able to formulate the abstract eigenvalue theorem. We denote by
$A^{\circ}$ the single valued part of a self-adjoint graph $A$. For an
immediate example we refer the reader to the beginning of the proof of Theorem
2.8.

######  Theorem 2.2.

Let $V$, $H$ and $K$ be Hilbert spaces. Suppose that $V$ is embedded in $H$
and that the inclusion map $i\colon V\to H$ is compact. Let $j\colon V\to K$
be a compact linear map. Let $\gothic{a}\colon V\times V\to\mathds{C}$ be a
positive symmetric continuous $i$-elliptic sesquilinear form. Let $V_{D}=\ker
j$. Let $A^{N}$ be the self-adjoint operator in $\overline{V}\subset H$
associated with $\gothic{a}$ and let $A^{D}$ be the self-adjoint operator in
$\overline{V_{D}}\subset H$ associated with $\gothic{a}|_{V_{D}\times V_{D}}$,
where the closures are in $H$. Further, for all $\lambda\in\mathds{R}$ define
$\gothic{b}_{\lambda}\colon V\times V\to\mathds{C}$ by

$\gothic{b}_{\lambda}(u,v)=\gothic{a}(u,v)-\lambda\,(u,v)_{H}.$

Let ${\cal N}_{\lambda}$ be the self-adjoint graph associated with
$(\gothic{b}_{\lambda},j)$. Suppose

  * (I)

the operator $A^{N}$ has no eigenvector in $V_{D}$ and

  * (II)

$\displaystyle\dim\mathop{\rm span}\\{\varphi\in D({\cal N}_{\lambda}):({\cal
N}_{\lambda}^{\circ}\varphi,\varphi)_{K}=0\\}=\infty$ for all
$\lambda\in(0,\infty)$.

Let $\lambda_{1}^{N}\leq\lambda_{2}^{N}\leq\ldots$ and
$\lambda_{1}^{D}\leq\lambda_{2}^{D}\leq\ldots$ be the eigenvalues of $A^{N}$
and $A^{D}$ repeated with multiplicity, respectively. Then

$\lambda_{n+1}^{N}<\lambda_{n}^{D}$

for all $n\in\mathds{N}$.

###### Proof.

Since $\gothic{a}$ is $i$-elliptic there are $\omega,\delta>0$ such that

$\gothic{a}(u,u)+\omega\,\|u\|_{H}^{2}\geq\delta\,\|u\|_{V}^{2}$ (3)

for all $u\in V$. Let $\widehat{A}^{D}$ be the self-adjoint graph in $H$
associated with $(\gothic{a}|_{V_{D}\times V_{D}},i|_{V_{D}})$. Then
$A^{D}=(\widehat{A}^{D})^{\circ}$, the singular part of $\widehat{A}^{D}$ and
$\widehat{A}^{D}$ is positive. Moreover, $\widehat{A}^{D}$ has compact
resolvent since $i$ is compact.

For all $\mu\in\mathds{R}$ define the form $\gothic{a}_{\mu}\colon V\times
V\to\mathds{C}$ by

$\gothic{a}_{\mu}(u,v)=\gothic{a}(u,v)-\mu\,(j(u),j(v))_{K}.$

Then $\gothic{a}_{\mu}$ is $i$-elliptic. Let $A_{\mu}$ be the self-adjoint
graph associated with $(\gothic{a}_{\mu},i)$. Then $A_{\mu}$ has compact
resolvent by Proposition 2.1. Let
$\lambda_{1}(\mu)\leq\lambda_{2}(\mu)\leq\ldots$ be the eigenvalues of
$A_{\mu}$, repeated with multiplicity. Then
$\lambda_{n}\colon\mathds{R}\to\mathds{R}$ is a function for all
$n\in\mathds{N}$. Since $\mu\mapsto\gothic{a}_{\mu}(u)$ is a decreasing
function for all $u\in V$, it follows from the mini-max theorem that
$\lambda_{n}$ is a decreasing function for all $n\in\mathds{N}$.

######  Lemma 2.3.

$\lim_{\mu\to-\infty}(I+A_{\mu})^{-1}=(I+\widehat{A}^{D})^{-1}$ in ${\cal
L}(H)$.

###### Proof.

Let $\mu_{1},\mu_{2},\ldots\in(-\infty,0)$ and suppose that
$\lim_{n\to\infty}\mu_{n}=-\infty$. We shall prove that
$\lim_{n\to\infty}\|(I+A_{\mu_{n}})^{-1}-(I+\widehat{A}^{D})^{-1}\|_{{\cal
L}(H)}=0$. Suppose not. Then there are $\varepsilon>0$ and
$f_{1},f_{2},\ldots\in H$ such that

$\|(I+A_{\mu_{n}})^{-1}f_{n}-(I+\widehat{A}^{D})^{-1}f_{n}\|_{H}>\varepsilon\,\|f_{n}\|_{H}$
(4)

for all $n\in\mathds{N}$. Without loss of generality we may assume that
$\|f_{n}\|_{H}=1$ for all $n\in\mathds{N}$. Passing to a subsequence if
necessary, there exists an $f\in H$ such that $\lim f_{n}=f$ weakly in $H$.
Then $\lim(I+\widehat{A}^{D})^{-1}f_{n}=(I+\widehat{A}^{D})^{-1}f$ in $H$
since $(I+\widehat{A}^{D})^{-1}$ is compact. For all $n\in\mathds{N}$ set
$u_{n}=(I+A_{\mu_{n}})^{-1}f_{n}$.

Let $n\in\mathds{N}$. Then

$\gothic{a}(u_{n},v)-\mu_{n}\,(j(u_{n}),j(v))_{K}+(u_{n},v)_{H}=(f_{n},v)_{H}$
(5)

for all $v\in V$. Choosing $v=u_{n}$ gives

$\gothic{a}(u_{n},u_{n})+|\mu_{n}|\,\|j(u_{n})\|_{K}^{2}+\|u_{n}\|_{H}^{2}=\mathop{\rm
Re}(f_{n},u_{n})_{H}\leq\|f_{n}\|_{H}\,\|u_{n}\|_{H}=\|u_{n}\|_{H}.$

So $\|u_{n}\|_{H}\leq 1$ and then also $\gothic{a}(u_{n},u_{n})\leq 1$ and
$|\mu_{n}|\,\|j(u_{n})\|_{K}^{2}\leq 1$. Hence $\lim j(u_{n})=0$ in $K$.
Moreover, $\delta\,\|u_{n}\|_{V}^{2}\leq(1+\omega)$ for all $n\in\mathds{N}$
by (3). Therefore the sequence $(u_{n})_{n\in\mathds{N}}$ is bounded in $V$.
Passing to a subsequence if necessary, there exists a $u\in V$ such that $\lim
u_{n}=u$ weakly in $V$. Then $\lim u_{n}=u$ strongly in $H$ since the
embedding $i$ is compact. Similarly $j(u)=\lim j(u_{n})$ in $K$. But $\lim
j(u_{n})=0$ in $K$. So $j(u)=0$ and $u\in V_{D}$.

Let $v\in V_{D}$. Taking the limit $n\to\infty$ in (5) gives
$\gothic{a}(u,v)+(u,v)_{H}=(f,v)_{H}$. So $u\in D(\widehat{A}^{D})$ and
$u=(I+\widehat{A}^{D})^{-1}f$. We proved that $\lim u_{n}=u$ in $H$. Hence

$\lim_{n\to\infty}(I+A_{\mu_{n}})^{-1}f_{n}-(I+\widehat{A}^{D})^{-1}f_{n}=(I+\widehat{A}^{D})^{-1}f-(I+\widehat{A}^{D})^{-1}f=0$

in $H$. This contradicts (4). ∎

######  Proposition 2.4.

If $n\in\mathds{N}$, then
$\lim_{\mu\to-\infty}\lambda_{n}(\mu)=\lambda_{n}^{D}$.

###### Proof.

Let $n\in\mathds{N}$. Let $\varepsilon>0$. By Lemma 2.3 there exists an
$M\in\mathds{N}$ such that

$((I+\widehat{A}^{D})^{-1}f,f)_{H}-\varepsilon\,\|f\|_{H}^{2}\leq((I+A_{\mu})^{-1}f,f)_{H}\leq((I+\widehat{A}^{D})^{-1}f,f)_{H}+\varepsilon\,\|f\|_{H}^{2}$

for all $f\in H$ and $\mu\in(-\infty,-M]$. Note that the max-min theorem gives

$(1+\lambda_{n}(\mu))^{-1}=\max_{\scriptstyle W\subset H\atop\scriptstyle\dim
W=n}\min_{\scriptstyle f\in
W\atop\scriptstyle\|f\|_{H}=1}((I+A_{\mu})^{-1}f,f)_{H}$

and a similar expression is valid for $\widehat{A}^{D}$. So
$(1+\lambda_{n}^{D})^{-1}-\varepsilon\leq(1+\lambda_{n}(\mu))^{-1}\leq(1+\lambda_{n}^{D})^{-1}+\varepsilon$
for all $\mu\in(-\infty,-M]$ and the proposition follows. ∎

The next lemma is a version of the Birman–Schwinger principle.

######  Lemma 2.5.

Let $\lambda,\mu\in\mathds{R}$. Then $\lambda\in\sigma(A_{\mu})$ if and only
if $\mu\in\sigma({\cal N}_{\lambda})$.

###### Proof.

‘$\Rightarrow$’. Suppose $\lambda\in\sigma(A_{\mu})$. Then there exists a
$u\in D(A_{\mu})$ such that $u\neq 0$ and $A_{\mu}u=\lambda\,u$. Then

$\gothic{a}(u,v)-\mu\,(j(u),j(v))_{K}=(\lambda\,u,v)_{H}$ (6)

for all $v\in V$. Hence
$\gothic{b}_{\lambda}(u,v)=\gothic{a}(u,v)-\lambda\,(u,v)_{H}=(\mu\,j(u),j(v))_{K}$
for all $v\in V$. Therefore $j(u)\in D({\cal N}_{\lambda})$ and
$(j(u),\mu\,j(u))\in{\cal N}_{\lambda}$. It remains to show that $j(u)\neq 0$.
Suppose that $j(u)=0$. Then it follows from (6) that
$\gothic{a}(u,v)=(\lambda\,u,v)_{H}$ for all $v\in V$. So $u\in D(A^{N})$ and
$A^{N}u=\lambda\,u$. Hence $u$ is an eigenvector for $A^{N}$ and also $u\in
V_{D}$. By Assumption (I) this is not possible. So $j(u)\neq 0$ and
$\mu\in\sigma({\cal N}_{\lambda})$.

‘$\Leftarrow$’. The proof is similar, even easier. ∎

######  Lemma 2.6.

Let $n\in\mathds{N}$. Then $\lambda_{n}$ is strictly decreasing.

###### Proof.

We already know that $\lambda_{n}$ is decreasing. Let
$\mu_{1},\mu_{2}\in\mathds{R}$ with $\mu_{1}<\mu_{2}$. Suppose that
$\lambda_{n}(\mu_{1})=\lambda_{n}(\mu_{2})$. Write
$\lambda=\lambda_{n}(\mu_{1})$. Let $\mu\in[\mu_{1},\mu_{2}]$. Since
$\lambda_{n}$ is decreasing, it follows that $\lambda_{n}(\mu)=\lambda$. So
$\lambda\in\sigma(A_{\mu})$. Hence $\mu\in\sigma({\cal N}_{\lambda})$ by Lemma
2.5. So $[\mu_{1},\mu_{2}]\subset\sigma({\cal N}_{\lambda})$. This is a
contradiction since ${\cal N}_{\lambda}$ has a discrete spectrum by
Proposition 2.1. ∎

######  Lemma 2.7.

Let $\lambda>0$. Then $\sigma({\cal
N}_{\lambda})\cap(-\infty,0)\neq\emptyset$.

###### Proof.

Let ${\cal N}_{\lambda}^{\circ}$ be the single valued part of ${\cal
N}_{\lambda}$. Suppose $\sigma({\cal N}_{\lambda})\cap(-\infty,0)=\emptyset$.
Then ${\cal N}_{\lambda}^{\circ}$ is a positive self-adjoint operator.

Let $\varphi\in D({\cal N}_{\lambda})$ and suppose that $({\cal
N}_{\lambda}^{\circ}\varphi,\varphi)_{K}=0$. Let $\psi\in D({\cal
N}_{\lambda}^{\circ})$. Then

$|({\cal N}_{\lambda}^{\circ}\varphi,\psi)_{K}|\leq\|({\cal
N}_{\lambda}^{\circ})^{1/2}\varphi\|_{K}\,\|({\cal
N}_{\lambda}^{\circ})^{1/2}\psi\|_{K}=({\cal
N}_{\lambda}^{\circ}\varphi,\varphi)_{K}^{1/2}\,\|({\cal
N}_{\lambda}^{\circ})^{1/2}\psi\|_{K}=0.$

Hence ${\cal N}_{\lambda}^{\circ}\varphi=0$ and $\varphi\in\ker({\cal
N}_{\lambda}^{\circ})$. By Assumption (II)

$\dim\mathop{\rm span}\\{\varphi\in D({\cal N}_{\lambda}):({\cal
N}_{\lambda}^{\circ}\varphi,\varphi)_{K}=0\\}=\infty.$

So $\dim\ker({\cal N}_{\lambda}^{\circ})=\infty$. But the operator ${\cal
N}_{\lambda}^{\circ}$ has compact resolvent by Proposition 2.1. This is a
contradiction. ∎

Now we complete the proof of Theorem 2.2. Let $n\in\mathds{N}$. Choose
$\lambda=\lambda_{n}^{D}$. Since $\lambda_{n}$ is strictly decreasing by Lemma
2.6, it follows from Proposition 2.4 that

$0\leq\lambda_{n}^{N}=\lambda_{n}(0)<\lambda_{n}^{D}=\lambda.$

Let $\mu=\min\sigma({\cal N}_{\lambda})$. Then $\mu<0$ by Lemma 2.7. Moreover,
$\lambda\in\sigma(A_{\mu})$ by Lemma 2.5. Hence there exists a
$j\in\mathds{N}$ such that $\lambda_{j}(\mu)=\lambda$. Using again Lemma 2.6
and Proposition 2.4 it follows that
$\lambda_{n}(\mu)<\lambda_{n}^{D}=\lambda=\lambda_{j}(\mu)$. Therefore $j\geq
n+1$. Then

$\lambda_{n+1}^{N}=\lambda_{n+1}(0)<\lambda_{n+1}(\mu)\leq\lambda_{j}(\mu)=\lambda=\lambda_{n}^{D},$

where the strict inequality follows from Lemma 2.6. ∎

The first main example is the classical version of Friedlander’s theorem.

######  Theorem 2.8.

Let $\Omega\subset\mathds{R}^{d}$ be a bounded open set with Lipschitz
boundary and $d\geq 2$. Let $\Delta^{N}$ and $\Delta^{D}$ be the Neumann and
Dirichlet Laplacians on $\Omega$, respectively. Further let
$\lambda_{1}^{N}\leq\lambda_{2}^{N}\leq\ldots$ and
$\lambda_{1}^{D}\leq\lambda_{2}^{D}\leq\ldots$ be the eigenvalues of
$-\Delta^{N}$ and $-\Delta^{D}$ repeated with multiplicity, respectively. Then

$\lambda_{n+1}^{N}<\lambda_{n}^{D}$

for all $n\in\mathds{N}$.

###### Proof.

Choose $V=H^{1}(\Omega)$, $H=L_{2}(\Omega)$ and $K=L_{2}(\partial\Omega)$. Let
$i\colon V\to H$ be the inclusion map and $j={\mathop{\rm Tr\,}}\colon V\to
K$. Then $i$ and $j$ are compact. Moreover, $\ker j=H^{1}_{0}(\Omega)$. Define
$\gothic{a}\colon V\times V\to\mathds{C}$ by

$\gothic{a}(u,v)=\int_{\Omega}\nabla u\cdot\overline{\nabla v}.$

Then $\gothic{a}$ is a positive symmetric continuous $i$-elliptic sesquilinear
form. Further, $-\Delta^{N}$ is the self-adjoint operator associated with
$\gothic{a}$ and $-\Delta^{D}$ is the self-adjoint operator associated with
$\gothic{a}|_{H^{1}_{0}(\Omega)\times H^{1}_{0}(\Omega)}=\gothic{a}|_{\ker
j\times\ker j}$.

It remains to verify Conditions (I) and (II) of Theorem 2.2.

‘(I)’. Let $u\in D(A^{N})\cap V_{D}$ and $\lambda\in\mathds{R}$. Suppose that
$A^{N}u=\lambda\,u$. Then $u\in H^{1}_{0}(\Omega)$. Let $\tilde{u}\in
H^{1}(\mathds{R}^{d})$ be the extension of $u$ by zero. Then

$\int_{\mathds{R}^{d}}\nabla\tilde{u}\cdot\overline{\nabla
v}=\int_{\Omega}\nabla
u\cdot\overline{\nabla(v|_{\Omega})}=(\lambda\,u,v|_{\Omega})_{L_{2}(\Omega)}=(\lambda\,\tilde{u},v)_{L_{2}(\mathds{R}^{d})}$

for all $v\in H^{1}(\mathds{R}^{d})$. So
$-\widetilde{\Delta}^{N}\tilde{u}=\lambda\,\tilde{u}$, where
$\widetilde{\Delta}^{N}$ is the Laplacian on $\mathds{R}^{d}$. Since
$\tilde{u}$ vanishes on a non-empty open set, it follows from the unique
continuation property that $\tilde{u}=0$. (See for example [RS] Theorem
XIII.57.) Then also $u=0$. So $A^{N}$ has no eigenvector in $V_{D}$.

‘(II)’. Let $\lambda\in(0,\infty)$. Let $\omega\in\mathds{R}^{d}$ and suppose
that $|\omega|^{2}=\lambda$. Define $u\in H^{2}(\Omega)$ by
$u(x)=e^{i\omega\cdot x}$. Then $-\Delta u=\lambda u$ as distribution and

$\displaystyle\gothic{b}_{\lambda}(u,v)$ $\displaystyle=$
$\displaystyle\int_{\Omega}\nabla u\cdot\overline{\nabla
v}-\lambda\int_{\Omega}u\,\overline{v}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\nabla u\cdot\overline{\nabla
v}+\int_{\Omega}(\Delta u)\,\overline{v}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\mathop{\rm div}(\overline{v}\,\nabla u)$
$\displaystyle=$ $\displaystyle\int_{\partial\Omega}\nu\cdot{\mathop{\rm
Tr\,}}(\overline{v}\,\nabla
u)=\int_{\partial\Omega}i(\omega\cdot\nu)\,({\mathop{\rm
Tr\,}}u)\,\overline{{\mathop{\rm Tr\,}}v}=(i(\omega\cdot\nu)\,{\mathop{\rm
Tr\,}}u,j(v))_{K}$

for all $v\in C_{b}^{\infty}(\Omega)$, where we used the divergence theorem
(see [Alt] Theorem A8.8). So by density (see for example [Maz] Theorem
1.1.6/2) one deduces that

$\gothic{b}_{\lambda}(u,v)=(i(\omega\cdot\nu)\,{\mathop{\rm Tr\,}}u,j(v))_{K}$

for all $v\in H^{1}(\Omega)$. Hence ${\mathop{\rm Tr\,}}u\in D({\cal
N}_{\lambda})$ and $({\mathop{\rm Tr\,}}u,i(\omega\cdot\nu)\,{\mathop{\rm
Tr\,}}u)\in{\cal N}_{\lambda}$. Moreover,

$({\cal
N}_{\lambda}^{\circ}j(u),j(u))_{L_{2}(\Omega,\mathds{C}^{d})}=\gothic{b}_{\lambda}(u,u)=\int_{\partial\Omega}i(\omega\cdot\nu)\,|{\mathop{\rm
Tr\,}}u|^{2}=i\int_{\partial\Omega}(\omega\cdot\nu)=i\int_{\Omega}\mathop{\rm
div}\omega=0.$

Therefore

${\mathop{\rm Tr\,}}u\in\\{\varphi\in D({\cal N}_{\lambda}):({\cal
N}_{\lambda}^{\circ}\varphi,\varphi)_{L_{2}(\partial\Omega)}=0\\}.$

It is then obvious that

$\dim\mathop{\rm span}\\{\varphi\in D({\cal N}_{\lambda}):({\cal
N}_{\lambda}^{\circ}\varphi,\varphi)_{L_{2}(\partial\Omega)}=0\\}=\infty.$

This completes the proof of the theorem. ∎

##  3 The Stokes operators

In this section we describe the Dirichlet and the Neumann version of the
Stokes operator on a Lipschitz domain with the aid of sesquilinear forms. We
then apply Theorem 2.2 (abstract eigenvalues comparison theorem) to obtain a
Stokes version of Friedlander’s eigenvalues theorem. For independent interest
we present a description of the Stokes type Dirichlet-to-Neumann operator.

Let $\Omega\subset\mathds{R}^{d}$ be a bounded open connected set with
Lipschitz boundary. We assume that $d\geq 2$. We denote by $\nu$ the outward
normalised normal. Let $H=L_{2}(\Omega,\mathds{C}^{d})$,

$V=\\{u\in H^{1}(\Omega,\mathds{C}^{d}):\mathop{\rm div}u=0\\}$

and $K=L_{2}(\partial\Omega,\mathds{C}^{d})$. Here and below $\mathop{\rm
div}$ is the divergence operator acting on distributions. Then the inclusion
$i\colon V\to H$ is compact. Let $j={\mathop{\rm Tr\,}}\colon V\to K$ be the
trace operator. Then also $j$ is compact. Define the form $\gothic{a}\colon
V\times V\to\mathds{C}$ by

$\gothic{a}(u,v)=\int_{\Omega}\nabla u\cdot\overline{\nabla v}.$

Then $\gothic{a}$ is continuous and $i$-elliptic. Let $A^{N}$ be the self-
adjoint operator in $\overline{V}\subset H$ associated with $\gothic{a}$ and
let $A^{D}$ be the self-adjoint operator in $\overline{V_{D}}\subset H$
associated with $\gothic{a}|_{V_{D}\times V_{D}}$, where the closures are in
$H$ and $V_{D}=\ker j$. We first characterise $A^{D}$ and $A^{N}$.

Define

$C^{\infty}_{c,\sigma}(\Omega)=\\{u\in
C_{c}^{\infty}(\Omega,\mathds{C}^{d}):\mathop{\rm div}u=0\\},$

the space of solenoidal test functions, and let $L_{2,\sigma}(\Omega)$ be the
closure of $C^{\infty}_{c,\sigma}(\Omega)$ in $L_{2}(\Omega,\mathds{C}^{d})$.

######  Proposition 3.1.

  * (a)

The closure of the space $V_{D}$ in $L_{2}(\Omega,\mathds{C}^{d})$ is
$L_{2,\sigma}(\Omega)$.

  * (b)

$V_{D}^{\perp}=L_{2,\sigma}(\Omega)^{\perp}=\\{\nabla\pi:\pi\in
H^{1}(\Omega)\\}$, where the orthogonal complement is in
$L_{2}(\Omega,\mathds{C}^{d})$.

  * (c)

Let $u,f\in L_{2,\sigma}(\Omega)$. Then $u\in D(A^{D})$ and $A^{D}u=f$ if and
only if $u\in H^{1}_{0}(\Omega,\mathds{C}^{d})$, $\mathop{\rm div}u=0$ and
there exists a $\pi\in L_{2}(\Omega)$ such that $f=-\Delta u+\nabla\pi$ in
$H^{-1}(\Omega,\mathds{C}^{d})$.

###### Proof.

‘(a)’. This follows from [Soh] Lemma II.2.2.3.

‘(b)’. The first equality is obvious. The second one follows from [Tem]
Theorem 1.1.4.

‘(c)’. Clearly $V_{D}=\\{w\in H^{1}_{0}(\Omega,\mathds{C}^{d}):\mathop{\rm
div}w=0\\}$. ‘$\Rightarrow$’. Since $A^{D}$ is the operator associated with
$\gothic{a}|_{V_{D}\times V_{D}}$ it follows that $u\in V_{D}$ and
$\gothic{a}(u,v)=(f,v)_{L_{2}(\Omega,\mathds{C}^{d})}$ for all $v\in V_{D}$.
So $u\in H^{1}_{0}(\Omega,\mathds{C}^{d})$ and $\mathop{\rm div}u=0$. Consider
$(f+\Delta u)\in H^{-1}(\Omega,\mathds{C}^{d})$. If $v\in
H^{1}_{0}(\Omega,\mathds{C}^{d})$ and $\mathop{\rm div}v=0$, then

$(f+\Delta u)(v)=(f,v)_{L_{2}(\Omega,\mathds{C}^{d})}+\langle\Delta
u,v\rangle_{H^{-1}(\Omega,\mathds{C}^{d})\times
H^{1}_{0}(\Omega,\mathds{C}^{d})}=(f,v)_{L_{2}(\Omega,\mathds{C}^{d})}-\gothic{a}(u,v)=0.$

Since $\Omega$ has a Lipschitz boundary, it follows from [Soh] Lemma II.2.1.1b
that there exists a $\pi\in L_{2}(\Omega)$ such that $f+\Delta u=\nabla\pi$ in
$H^{-1}(\Omega,\mathds{C}^{d})$.

‘$\Leftarrow$’. Obviously $u\in V_{D}$. If $v\in V_{D}$, then $v\in
H^{1}_{0}(\Omega,\mathds{C}^{d})$ and $\mathop{\rm div}v=0$. So

$\displaystyle\gothic{a}(u,v)$ $\displaystyle=$ $\displaystyle-\langle\Delta
u,v\rangle_{H^{-1}(\Omega,\mathds{C}^{d})\times
H^{1}_{0}(\Omega,\mathds{C}^{d})}=\langle
f-\nabla\pi,v\rangle_{H^{-1}(\Omega,\mathds{C}^{d})\times
H^{1}_{0}(\Omega,\mathds{C}^{d})}$ $\displaystyle=$
$\displaystyle(f,v)_{L_{2}(\Omega,\mathds{C}^{d})}+(\pi,\mathop{\rm
div}v)_{L_{2}(\Omega)}=(f,v)_{L_{2}(\Omega,\mathds{C}^{d})}$

and the proposition follows. ∎

We next consider the operator $A^{N}$. In order to obtain a description as in
Proposition 3.1(c) for $A^{D}$, we first need a kind of normal derivative. If
$\Omega$ is smooth, $u\in H^{2}(\Omega,\mathds{C}^{d})$, $\pi\in
H^{1}(\Omega)$ and $f\in L_{2}(\Omega,\mathds{C}^{d})$ with $-\Delta
u+\nabla\pi=f$ in $H^{-1}(\Omega,\mathds{C}^{d})$, then

$\displaystyle\langle\partial_{\nu}u-({\mathop{\rm
Tr\,}}\pi)\,\nu,{\mathop{\rm Tr\,}}\Phi\rangle_{H^{-1/2}\times H^{1/2}}$
$\displaystyle=$ $\displaystyle\int_{\Omega}\nabla
u\cdot\overline{\nabla\Phi}+\langle\Delta u,\Phi\rangle_{H^{-1}\times
H^{1}}-\langle\nabla\pi,\Phi\rangle_{H^{-1}\times
H^{1}}-\int_{\Omega}\pi\,\overline{\mathop{\rm div}\Phi}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\nabla
u\cdot\overline{\nabla\Phi}-\int_{\Omega}\pi\,\overline{\mathop{\rm
div}\Phi}-\int_{\Omega}f\cdot\overline{\Phi}$

for all $\Phi\in H^{1}(\Omega,\mathds{C}^{d})$. We next define a weak version
of this normal derivative.

######  Lemma 3.2.

Let $u\in H^{1}(\Omega,\mathds{C}^{d})$, $\pi\in L_{2}(\Omega)$ and $f\in
L_{2}(\Omega,\mathds{C}^{d})$. Suppose that $-\Delta u+\nabla\pi=f$ in
$H^{-1}(\Omega,\mathds{C}^{d})$. Then there exists a unique $F\in
H^{-1/2}(\partial\Omega,\mathds{C}^{d})$ such that

$\langle F,\varphi\rangle_{H^{-1/2}(\partial\Omega,\mathds{C}^{d})\times
H^{1/2}(\partial\Omega,\mathds{C}^{d})}=\int_{\Omega}\nabla
u\cdot\overline{\nabla\Phi}-\int_{\Omega}\pi\,\overline{\mathop{\rm
div}\Phi}-\int_{\Omega}f\cdot\overline{\Phi}$

for all $\varphi\in H^{1/2}(\partial\Omega,\mathds{C}^{d})$ and $\Phi\in
H^{1}(\Omega,\mathds{C}^{d})$ with ${\mathop{\rm Tr\,}}\Phi=\varphi$.

###### Proof.

Let $\Phi\in H^{1}(\Omega,\mathds{C}^{d})$ and suppose that ${\mathop{\rm
Tr\,}}\Phi=0$. Then

$\displaystyle\int_{\Omega}\nabla
u\cdot\overline{\nabla\Phi}-\int_{\Omega}\pi\,\overline{\mathop{\rm
div}\Phi}-\int_{\Omega}f\cdot\overline{\Phi}$ $\displaystyle=$
$\displaystyle\langle-\Delta u,\Phi\rangle_{H^{-1}\times
H^{1}}+\langle\nabla\pi,\Phi\rangle_{H^{-1}\times
H^{1}}-\int_{\Omega}f\cdot\overline{\Phi}$ $\displaystyle=$
$\displaystyle\langle-\Delta
u+\nabla\pi-f,\Phi\rangle_{H^{-1}(\Omega,\mathds{C}^{d})\times
H^{1}(\Omega,\mathds{C}^{d})}=0.$

Hence there exists a unique $F\colon
H^{1/2}(\partial\Omega,\mathds{C}^{d})\to\mathds{C}$ such that

$F(\varphi)=\int_{\Omega}\nabla
u\cdot\overline{\nabla\Phi}-\int_{\Omega}\pi\,\overline{\mathop{\rm
div}\Phi}-\int_{\Omega}f\cdot\overline{\Phi}$

for all $\varphi\in H^{1/2}(\partial\Omega,\mathds{C}^{d})$ and $\Phi\in
H^{1}(\Omega,\mathds{C}^{d})$ with ${\mathop{\rm Tr\,}}\Phi=\varphi$.
Obviously $F$ is anti-linear and

$|F({\mathop{\rm
Tr\,}}\Phi)|\leq(\|u\|_{H^{1}(\Omega,\mathds{C}^{d})}+\|\pi\|_{L_{2}(\Omega)}+\|f\|_{L_{2}(\Omega,\mathds{C}^{d})})\|\Phi\|_{H^{1}(\Omega,\mathds{C}^{d})}$

for all $\Phi\in H^{1}(\Omega,\mathds{C}^{d})$. Since ${\mathop{\rm
Tr\,}}\colon H^{1}(\Omega,\mathds{C}^{d})\to
H^{1/2}(\partial\Omega,\mathds{C}^{d})$ has a continuous right-inverse, the
lemma follows. ∎

We denote the element $F\in H^{-1/2}(\partial\Omega,\mathds{C}^{d})$ in Lemma
3.2 by $\partial_{\nu}(u,\pi)$. So

$\langle\partial_{\nu}(u,\pi),{\mathop{\rm
Tr\,}}\Phi\rangle_{H^{-1/2}(\partial\Omega,\mathds{C}^{d})\times
H^{1/2}(\partial\Omega,\mathds{C}^{d})}=\int_{\Omega}\nabla
u\cdot\overline{\nabla\Phi}-\int_{\Omega}\pi\,\overline{\mathop{\rm
div}\Phi}-\int_{\Omega}f\cdot\overline{\Phi}$

for all $u\in H^{1}(\Omega,\mathds{C}^{d})$, $\pi\in L_{2}(\Omega)$, $f\in
L_{2}(\Omega,\mathds{C}^{d})$ and $\Phi\in H^{1}(\Omega,\mathds{C}^{d})$ with
$-\Delta u+\nabla\pi=f$ in $H^{-1}(\Omega,\mathds{C}^{d})$.

Let $u\in H^{1}(\Omega,\mathds{C}^{d})$, $\pi\in L_{2}(\Omega)$, $f\in
L_{2}(\Omega,\mathds{C}^{d})$ and $c\in\mathds{C}$ with $-\Delta
u+\nabla\pi=f$ in $H^{-1}(\Omega,\mathds{C}^{d})$. Then $-\Delta
u+\nabla(\pi+c\,\mathds{1}_{\Omega})=f$ in $H^{-1}(\Omega,\mathds{C}^{d})$ and

$\partial_{\nu}(u,\pi+c\,\mathds{1}_{\Omega})=\partial_{\nu}(u,\pi)-c\,\nu$
(7)

by a straightforward calculation.

We also need to characterise the range ${\mathop{\rm Tr\,}}V$ of the form
domain under the trace map. Note that ${\mathop{\rm Tr\,}}\colon
H^{1}(\Omega,\mathds{C}^{d})\to H^{1/2}(\partial\Omega,\mathds{C}^{d})$. If
$u\in V$, then the divergence theorem gives

$0=\int_{\Omega}\mathop{\rm div}u=\int_{\partial\Omega}\nu\cdot{\mathop{\rm
Tr\,}}u.$

Define

$L_{2,0}(\partial\Omega,\mathds{C}^{d})=\\{\varphi\in
L_{2}(\partial\Omega,\mathds{C}^{d}):\int_{\partial\Omega}\varphi\cdot\nu=0\\}.$

Then ${\mathop{\rm Tr\,}}(V)\subset H^{1/2}(\partial\Omega,\mathds{C}^{d})\cap
L_{2,0}(\partial\Omega,\mathds{C}^{d})$. We next show that one actually has an
equality.

######  Lemma 3.3.

${\mathop{\rm Tr\,}}(V)=H^{1/2}(\partial\Omega,\mathds{C}^{d})\cap
L_{2,0}(\partial\Omega,\mathds{C}^{d})$.

###### Proof.

We only have to show ‘$\supset$’. Let $\varphi\in
H^{1/2}(\partial\Omega,\mathds{C}^{d})\cap
L_{2,0}(\partial\Omega,\mathds{C}^{d})$. There exists a $u\in
H^{1}(\Omega,\mathds{C}^{d})$ such that ${\mathop{\rm Tr\,}}u=\varphi$. Then
$\int_{\Omega}\mathop{\rm div}u=\int_{\partial\Omega}\nu\cdot{\mathop{\rm
Tr\,}}u=\int_{\partial\Omega}\nu\cdot\varphi=0$. Hence by [Soh] Lemma
II.2.1.1a there exists a $w\in H^{1}_{0}(\Omega,\mathds{C}^{d})$ such that
$\mathop{\rm div}w=\mathop{\rm div}u$. Set $v=u-w$. Then $v\in
H^{1}(\Omega,\mathds{C}^{d})$ and $\mathop{\rm div}v=0$. So $v\in V$.
Moreover, ${\mathop{\rm Tr\,}}v={\mathop{\rm Tr\,}}u=\varphi$. ∎

Although we do not need the next density lemma, we state it for completeness.

######  Lemma 3.4.

The range ${\mathop{\rm Tr\,}}(V)$ is dense in
$L_{2,0}(\partial\Omega,\mathds{C}^{d})$.

###### Proof.

Let $\varphi\in L_{2,0}(\partial\Omega,\mathds{C}^{d})$. Since
$H^{1/2}(\partial\Omega,\mathds{C}^{d})$ is dense in
$L_{2}(\partial\Omega,\mathds{C}^{d})$, there is a sequence
$(\varphi_{n})_{n\in\mathds{N}}$ in $H^{1/2}(\partial\Omega,\mathds{C}^{d})$
such that $\lim\varphi_{n}=\varphi$ in $L_{2}(\partial\Omega,\mathds{C}^{d})$.
Then
$\lim\int_{\partial\Omega}\varphi_{n}\cdot\nu=\int_{\partial\Omega}\varphi\cdot\nu=0$.
Note that $\nu\in L_{\infty}(\partial\Omega,\mathds{C}^{d})\subset
L_{2}(\partial\Omega,\mathds{C}^{d})$. Hence as above there is a sequence
$(\nu_{n})_{n\in\mathds{N}}$ in $H^{1/2}(\partial\Omega,\mathds{C}^{d})$ such
that $\lim\nu_{n}=\nu$ in $L_{2}(\partial\Omega,\mathds{C}^{d})$. Then
$\lim\int_{\partial\Omega}\nu_{n}\cdot\nu=\sigma(\partial\Omega)\neq 0$, so we
may assume that $\int_{\partial\Omega}\nu_{n}\cdot\nu\neq 0$ for all
$n\in\mathds{N}$. For all $n\in\mathds{N}$ define
$\tilde{\varphi}_{n}=\varphi_{n}-c_{n}\,\nu_{n}$, where
$c_{n}=\Big{(}\int_{\partial\Omega}\nu_{n}\cdot\nu\Big{)}^{-1}\int_{\partial\Omega}\varphi_{n}\cdot\nu$.
Then $\tilde{\varphi}_{n}\in H^{1/2}(\partial\Omega,\mathds{C}^{d})\cap
L_{2,0}(\partial\Omega,\mathds{C}^{d})={\mathop{\rm Tr\,}}(V)$ by Lemma 3.3.
Moreover, $\lim\tilde{\varphi}_{n}=\varphi$ in
$L_{2}(\partial\Omega,\mathds{C}^{d})$. ∎

Now we are able to characterise the operator $A^{N}$.

######  Proposition 3.5.

  * (a)

The closure of the space $V$ in $L_{2}(\Omega,\mathds{C}^{d})$ is $\\{u\in
L_{2}(\Omega,\mathds{C}^{d}):\mathop{\rm div}u=0\\}$.

  * (b)

$V^{\perp}=\\{\nabla\pi:\pi\in H^{1}_{0}(\Omega)\\}$, where the orthogonal
complement is in $L_{2}(\Omega,\mathds{C}^{d})$.

  * (c)

Let $u,f\in\overline{V}^{(H)}$. Then $u\in D(A^{N})$ and $A^{N}u=f$ if and
only if $u\in H^{1}(\Omega,\mathds{C}^{d})$, $\mathop{\rm div}u=0$ and there
exists a $\pi\in L_{2}(\Omega)$ such that $f=-\Delta u+\nabla\pi$ in
$H^{-1}(\Omega,\mathds{C}^{d})$ and $\partial_{\nu}(u,\pi)=0$.

###### Proof.

‘(a)’. See [MMW] Lemma 2.1.

‘(b)’. By definition of $V$ we obtain that $V^{\perp}$ is the closure of
$\\{\nabla\pi:\pi\in C_{c}^{\infty}(\Omega)\\}$ in
$L_{2}(\Omega,\mathds{C}^{d})$. Obviously the latter is the closure of the
space $\\{\nabla\pi:\pi\in H^{1}_{0}(\Omega)\\}$ in
$L_{2}(\Omega,\mathds{C}^{d})$. By [Maz] Corollary 1.1.11 and the Dirichlet
Poincaré inequality the space $\\{\nabla\pi:\pi\in H^{1}_{0}(\Omega)\\}$ in
closed in $L_{2}(\Omega,\mathds{C}^{d})$.

‘(c)’. ‘$\Rightarrow$’. By definition $u\in V$ and
$\gothic{a}(u,v)=(f,v)_{L_{2}(\Omega,\mathds{C}^{d})}$ for all $v\in V$. So
$u\in H^{1}(\Omega,\mathds{C}^{d})$ and $\mathop{\rm div}v=0$. Consider
$(f+\Delta u)\in H^{-1}(\Omega,\mathds{C}^{d})$. If $v\in
H^{1}_{0}(\Omega,\mathds{C}^{d})$ and $\mathop{\rm div}v=0$, then

$(f+\Delta u)(v)=(f,v)_{L_{2}(\Omega,\mathds{C}^{d})}+\langle\Delta
u,v\rangle_{H^{-1}(\Omega,\mathds{C}^{d})\times
H^{1}_{0}(\Omega,\mathds{C}^{d})}=(f,v)_{L_{2}(\Omega,\mathds{C}^{d})}-\gothic{a}(u,v)=0.$

Since $\Omega$ has a Lipschitz boundary, it follows from [Soh] Lemma II.2.1.1b
that there exists a $\pi\in L_{2}(\Omega)$ such that $f+\Delta u=\nabla\pi$ in
$H^{-1}(\Omega,\mathds{C}^{d})$.

If $v\in V$, then

$\displaystyle\langle\partial_{\nu}(u,\pi),{\mathop{\rm
Tr\,}}v\rangle_{H^{-1/2}(\partial\Omega,\mathds{C}^{d})\times
H^{1/2}(\partial\Omega,\mathds{C}^{d})}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\nabla u\cdot\overline{\nabla
v}-\int_{\Omega}\pi\,\overline{\mathop{\rm
div}v}-\int_{\Omega}f\cdot\overline{v}$ $\displaystyle=$
$\displaystyle\gothic{a}(u,v)-(f,v)_{L_{2}(\Omega,\mathds{C}^{d})}=0.$

So by Lemma 3.3 one deduces that
$\langle\partial_{\nu}(u,\pi),\varphi\rangle_{H^{-1/2}(\partial\Omega,\mathds{C}^{d})\times
H^{1/2}(\partial\Omega,\mathds{C}^{d})}=0$ for all $\varphi\in
H^{1/2}(\partial\Omega,\mathds{C}^{d})\cap
L_{2,0}(\partial\Omega,\mathds{C}^{d})$.

Fix $\varphi_{0}\in H^{1/2}(\partial\Omega,\mathds{C}^{d})$ such that
$(\varphi_{0},\nu)_{L_{2}(\partial\Omega,\mathds{C}^{d})}\neq 0$. Let
$c\in\mathds{C}$ be such that

$\langle\partial_{\nu}(u,\pi)-c\,\nu,\varphi_{0}\rangle_{H^{-1/2}(\partial\Omega,\mathds{C}^{d})\times
H^{1/2}(\partial\Omega,\mathds{C}^{d})}=0.$

Then it follows from (7) that

$\langle\partial_{\nu}(u,\pi+c\,\mathds{1}_{\Omega}),\varphi\rangle_{H^{-1/2}(\partial\Omega,\mathds{C}^{d})\times
H^{1/2}(\partial\Omega,\mathds{C}^{d})}=0$

for all
$\varphi\in\\{\varphi_{0}\\}\cup\Big{(}H^{1/2}(\partial\Omega,\mathds{C}^{d})\cap
L_{2,0}(\partial\Omega,\mathds{C}^{d})\Big{)}$ and then by linearity for all

$\varphi\in\mathop{\rm
span}\bigg{(}\\{\varphi_{0}\\}\cup\Big{(}H^{1/2}(\partial\Omega,\mathds{C}^{d})\cap
L_{2,0}(\partial\Omega,\mathds{C}^{d})\Big{)}\bigg{)}=H^{1/2}(\partial\Omega,\mathds{C}^{d}).$

So $\partial_{\nu}(u,\pi+c\,\mathds{1}_{\Omega})=0$. Replacing $\pi$ by
$\pi+c\,\mathds{1}_{\Omega}$ completes the proof of the implication
‘$\Rightarrow$’.

‘$\Leftarrow$’. Let $u\in H^{1}(\Omega,\mathds{C}^{d})$ and
$f\in\overline{V}^{(H)}$. Suppose that $\mathop{\rm div}u=0$ and there exists
a $\pi\in L_{2}(\Omega)$ such that $f=-\Delta u+\nabla\pi$ and
$\partial_{\nu}(u,\pi)=0$. Then $u\in V$ and

$\int_{\Omega}\nabla u\cdot\overline{\nabla
v}-\int_{\Omega}\pi\,\overline{\mathop{\rm
div}v}-\int_{\Omega}f\cdot\overline{v}=\langle\partial_{\nu}(u,\pi),{\mathop{\rm
Tr\,}}v\rangle_{H^{-1/2}(\partial\Omega,\mathds{C}^{d})\times
H^{1/2}(\partial\Omega,\mathds{C}^{d})}=0$

for all $v\in H^{1}(\Omega,\mathds{C}^{d})$. Hence if $v\in V$, then
$\mathop{\rm div}v=0$ and

$\gothic{a}(u,v)=\int_{\Omega}\nabla u\cdot\overline{\nabla
v}=\int_{\Omega}f\cdot\overline{v}=(f,v)_{L_{2}(\Omega,\mathds{C}^{d})}.$

So $u\in\mathop{\rm dom}(A^{N})$ and $A^{N}u=f$. ∎

We next turn to a Stokes version of the Dirichlet-to-Neumann operator. For all
$\lambda\in\mathds{R}$ define $\gothic{b}_{\lambda}\colon V\times
V\to\mathds{C}$ by

$\gothic{b}_{\lambda}(u,v)=\gothic{a}(u,v)-\lambda\,(u,v)_{H}.$

Let ${\cal N}_{\lambda}$ be the self-adjoint graph in
$L_{2}(\partial\Omega,\mathds{C}^{d})$ associated with
$(\gothic{b}_{\lambda},j)$. We call ${\cal N}_{\lambda}$ the Stokes Dirichlet-
to-Neumann graph. The singular part ${\cal N}_{\lambda}^{\circ}$ is a self-
adjoint operator in the Hilbert space $\overline{D({\cal N}_{\lambda})}$,
where the closure is in $L_{2}(\partial\Omega,\mathds{C}^{d})$. Note that
$\overline{D({\cal N}_{\lambda})}\subset\overline{j(V)}=\overline{{\mathop{\rm
Tr\,}}V}=L_{2,0}(\partial\Omega,\mathds{C}^{d})$ by Lemma 3.4. We next
characterise ${\cal N}_{\lambda}^{\circ}$ in the case where $\lambda$ is not
an eigenvalue of $A^{D}$.

######  Proposition 3.6.

Let $\lambda\in\mathds{R}\setminus\sigma(A^{D})$. Let $\varphi,\psi\in
L_{2,0}(\partial\Omega,\mathds{C}^{d})$. Then the following are equivalent.

  * (i)

$\varphi\in D({\cal N}_{\lambda}^{\circ})$ and ${\cal
N}_{\lambda}^{\circ}\varphi=\psi$.

  * (ii)

There exist $u\in H^{1}(\Omega,\mathds{C}^{d})$ and $\pi\in L_{2}(\Omega)$
such that

    * •

$\varphi={\mathop{\rm Tr\,}}u$,

    * •

$\mathop{\rm div}u=0$,

    * •

$-\Delta u+\nabla\pi=\lambda\,u$ in $H^{-1}(\Omega,\mathds{C}^{d})$, and

    * •

$\psi=\partial_{\nu}(u,\pi)$.

###### Proof.

‘(i)$\Rightarrow$(ii)’. By definition there exists a $u\in V$ such that
${\mathop{\rm Tr\,}}u=\varphi$ and
$\gothic{b}_{\lambda}(u,v)=(\psi,{\mathop{\rm
Tr\,}}v)_{L_{2}(\partial\Omega,\mathds{C}^{d})}$ for all $v\in V$. Consider
$(-\Delta u-\lambda\,u)\in H^{-1}(\Omega,\mathds{C}^{d})$. If $v\in
H^{1}_{0}(\Omega,\mathds{C}^{d})$ and $\mathop{\rm div}v=0$, then

$(-\Delta u-\lambda\,u)(v)=\gothic{b}_{\lambda}(u,v)=(\psi,{\mathop{\rm
Tr\,}}v)_{L_{2}(\partial\Omega,\mathds{C}^{d})}=0.$

Since $\Omega$ has a Lipschitz boundary, it follows from [Soh] Lemma II.2.1.1b
that there exists a $\pi\in L_{2}(\Omega)$ such that $-\Delta
u-\lambda\,u=-\nabla\pi$ in $H^{-1}(\Omega,\mathds{C}^{d})$.

Let $v\in V$. Then

$\displaystyle(\psi,{\mathop{\rm
Tr\,}}v)_{L_{2}(\partial\Omega,\mathds{C}^{d})}$ $\displaystyle=$
$\displaystyle\gothic{b}_{\lambda}(u,v)=\int_{\Omega}\nabla
u\cdot\overline{\nabla v}-\lambda\int_{\Omega}u\cdot\overline{v}$
$\displaystyle=$ $\displaystyle\int_{\Omega}\nabla u\cdot\overline{\nabla
v}-\lambda\int_{\Omega}u\cdot\overline{v}-\int_{\Omega}\pi\,\overline{\mathop{\rm
div}v}=\langle\partial_{\nu}(u,\pi),{\mathop{\rm
Tr\,}}v\rangle_{H^{-1/2}\times H^{1/2}}.$

So $\langle\partial_{\nu}(u,\pi)-\psi,\tau\rangle_{H^{-1/2}\times H^{1/2}}=0$
for all $\tau\in{\mathop{\rm
Tr\,}}(V)=H^{1/2}(\partial\Omega,\mathds{C}^{d})\cap
L_{2,0}(\partial\Omega,\mathds{C}^{d})$ by Lemma 3.3. Now it follows as at the
end of the proof of the implication ‘$\Rightarrow$’ in the proof of
Proposition 3.5(c) that there exists a $c\in\mathds{C}$ such that
$\partial_{\nu}(u,\pi+c\,\mathds{1}_{\Omega})=\psi$.

‘(ii)$\Rightarrow$(i)’. Let $u$ and $\pi$ be as in (ii). Note that $u\in V$.
Let $v\in V$. Then

$\displaystyle\gothic{b}_{\lambda}(u,v)$ $\displaystyle=$
$\displaystyle\int_{\Omega}\nabla u\cdot\overline{\nabla
v}-\lambda\int_{\Omega}u\cdot\overline{v}-\int_{\Omega}\pi\,\overline{\mathop{\rm
div}v}$ $\displaystyle=$
$\displaystyle\langle\partial_{\nu}(u,\pi),{\mathop{\rm
Tr\,}}v\rangle_{H^{-1/2}\times H^{1/2}}=\langle\psi,{\mathop{\rm
Tr\,}}v\rangle_{H^{-1/2}\times H^{1/2}}=(\psi,{\mathop{\rm
Tr\,}}v)_{L_{2}(\partial\Omega,\mathds{C}^{d})}$

and (i) follows. ∎

Now we are able to state and prove the main theorem of this paper. For an
explicit description of the Neumann Stokes operator and Dirichlet Stokes
operator we refer to Propositions 3.5(c) and 3.1(c).

######  Theorem 3.7.

Let $\Omega\subset\mathds{R}^{d}$ be a bounded open connected set with
Lipschitz boundary. Let $A^{N}$ and $A^{D}$ be the Neumann Stokes operator and
Dirichlet Stokes operator on $\Omega$, respectively. Further let
$\lambda_{1}^{N}\leq\lambda_{2}^{N}\leq\ldots$ and
$\lambda_{1}^{D}\leq\lambda_{2}^{D}\leq\ldots$ be the eigenvalues of $A^{N}$
and $A^{D}$ repeated with multiplicity, respectively. Then

$\lambda_{n+1}^{N}<\lambda_{n}^{D}$

for all $n\in\mathds{N}$.

###### Proof.

We have to verify Conditions (I) and (II) of Theorem 2.2.

‘(I)’. Let $u\in D(A^{N})\cap\ker j$ and suppose that $u$ is an eigenvector
for $A^{N}$ with eigenvalue $\lambda\in\mathds{R}$. Then $u\in
H^{1}_{0}(\Omega,\mathds{C}^{d})$ and $\mathop{\rm div}u=0$ in $\Omega$. Let
$\tilde{u}\in H^{1}(\mathds{R}^{d},\mathds{C}^{d})$ be the extension by zero
of $u$. Then $\nabla\tilde{u}=\widetilde{\nabla u}$, the extension by zero of
$\nabla u$. So $\mathop{\rm div}\tilde{u}=0$ in $\mathds{R}^{d}$. Let $v\in
H^{1}(\mathds{R}^{d},\mathds{C}^{d})$ and suppose that $\mathop{\rm div}v=0$
in $\mathds{R}^{d}$. Then

$\int_{\mathds{R}^{d}}\nabla\tilde{u}\cdot\overline{\nabla
v}=\int_{\Omega}\nabla
u\cdot\overline{\nabla(v|_{\Omega})}=\gothic{a}(u,v|_{\Omega})=(\lambda\,u,v|_{\Omega})_{L_{2}(\Omega,\mathds{C}^{d})}=(\lambda\,\tilde{u},v)_{L_{2}(\mathds{R}^{d},\mathds{C}^{d})}.$

So $\tilde{u}\in D(A_{\mathds{R}^{d}})$, the Stokes operator on
$\mathds{R}^{d}$, and $A_{\mathds{R}^{d}}\tilde{u}=\lambda\,\tilde{u}$. Then
$-\Delta\tilde{u}=\lambda\,\tilde{u}$ by [Soh] Lemma III.2.3.2. But
$\tilde{u}$ vanishes on a non-empty open set. So $\tilde{u}=0$ by the unique
continuation property, see for example [RS] Theorem XIII.57. Hence $u=0$ and
$u$ is not an eigenvector.

‘(II)’. Let $\lambda>0$. Let $\omega\in\mathds{R}^{d}$ and suppose that
$|\omega|^{2}=\lambda$. Let $b\in\mathds{R}^{d}$ be such that $|b|=1$ and
$b\cdot\omega=0$. Define $\tau\in H^{2}(\Omega)$ by $\tau(x)=e^{i\omega\cdot
x}$. Then $-\Delta\tau=|\omega|^{2}\,\tau=\lambda\,\tau$. Let $u=\tau\,b\in
H^{1}(\Omega,\mathds{C}^{d})$. Then $\mathop{\rm
div}u=i\tau\,(b\cdot\omega)=0$, so $u\in V$. Let $v\in V$. Then

$\displaystyle\gothic{b}_{\lambda}(u,v)$ $\displaystyle=$
$\displaystyle\int_{\Omega}\nabla u\cdot\overline{\nabla
v}-\lambda\int_{\Omega}u\cdot\overline{v}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\nabla\tau\cdot\overline{\nabla(b\cdot
v)}-\lambda\int_{\Omega}\tau\,\overline{(b\cdot v)}$ $\displaystyle=$
$\displaystyle-\int_{\Omega}(\Delta\tau)\,\overline{(b\cdot
v)}+\int_{\partial\Omega}(\partial_{\nu}\tau)\,\overline{{\mathop{\rm
Tr\,}}(b\cdot v)}-\lambda\int_{\Omega}\tau\,\overline{(b\cdot v)}$
$\displaystyle=$
$\displaystyle\int_{\partial\Omega}(\partial_{\nu}\tau)\,b\cdot\overline{{\mathop{\rm
Tr\,}}v}=\int_{\partial\Omega}i\,(\omega\cdot\nu)\,({\mathop{\rm
Tr\,}}\tau)\,b\cdot\overline{{\mathop{\rm Tr\,}}v}.$

Therefore $j(u)\in D({\cal N}_{\lambda})$ and
$(j(u),i\,(\omega\cdot\nu)\,({\mathop{\rm Tr\,}}\tau)\,b)\in{\cal
N}_{\lambda}$. Moreover,

$\displaystyle({\cal
N}_{\lambda}^{\circ}j(u),j(u))_{L_{2}(\Omega,\mathds{C}^{d})}$
$\displaystyle=$ $\displaystyle\gothic{b}_{\lambda}(u,u)$ $\displaystyle=$
$\displaystyle i\int_{\partial\Omega}(\omega\cdot\nu)\,|b|^{2}\,|{\mathop{\rm
Tr\,}}\tau|^{2}=i\,|b|^{2}\int_{\partial\Omega}(\omega\cdot\nu)=i\,|b|^{2}\int_{\Omega}\mathop{\rm
div}\omega=0.$

Consequently

$\dim\mathop{\rm span}\\{\varphi\in D({\cal N}_{\lambda}):({\cal
N}_{\lambda}^{\circ}\varphi,\varphi)_{L_{2}(\Omega,\mathds{C}^{d})}=0\\}=\infty$

as required. ∎

We finish with a small extension of the previous theorem.

######  Example 3.8.

Let $\Omega\subset\mathds{R}^{d}$ be a bounded open connected set with
Lipschitz boundary. Let $H=L_{2}(\Omega,\mathds{C}^{d})$,

$V=\\{u\in H^{1}(\Omega,\mathds{C}^{d}):\mathop{\rm div}u=0\\}$

and $K=L_{2}(\partial\Omega,\mathds{C}^{d})$. Then the inclusion $i\colon V\to
H$ is compact. Let $j={\mathop{\rm Tr\,}}\colon V\to K$ be the trace operator.
Then $j$ is compact. Fix $\alpha\in(-1,1]$. Define the form $\gothic{a}\colon
V\times V\to\mathds{C}$ by

$\gothic{a}(u,v)=\int_{\Omega}\nabla u\cdot\overline{\nabla
v}+\alpha\int_{\Omega}(\nabla u)^{\rm T}\cdot\overline{\nabla v},$

where the superscript ${\rm T}$ denotes the transpose. So Theorem 3.7
corresponds to the case $\alpha=0$. Obviously $\gothic{a}$ is continuous and
if $|\alpha|<1$, then $\gothic{a}$ is $i$-elliptic. If $\alpha=1$, then

$\mathop{\rm Re}\gothic{a}(u)=\frac{1}{2}\int_{\Omega}|\nabla u+(\nabla
u)^{\rm T}|^{2}$

for all $u\in V$. Since $\Omega$ has a Lipschitz boundary, Korn’s second
inequality states that there exists a $C>0$ such that

$\|u\|_{H^{1}(\Omega,\mathds{C}^{d})}^{2}\leq C\Big{(}\int_{\Omega}|\nabla
u+(\nabla u)^{\rm T}|^{2}+\|u\|_{L_{2}(\Omega,\mathds{C}^{d})}^{2}\Big{)}$

for all $u\in H^{1}(\Omega,\mathds{C}^{d})$. For an easy proof, see [Nit]
Section 3. Hence the form $\gothic{a}$ is also $i$-elliptic if $\alpha=1$. The
associated operator $A^{N}$ is the Neumann Stokes operator studied in [MMW]
and $A^{D}$ is the (Dirichlet) Stokes operator. We next show that the
conditions in Theorem 2.2 are valid.

‘(I)’. Let $u\in D(A^{N})\cap\ker j$ and suppose that $u$ is an eigenvector
for $A^{N}$ with eigenvalue $\lambda\in\mathds{R}$. Then as before $u\in
H^{1}_{0}(\Omega,\mathds{C}^{d})$ and $\mathop{\rm div}u=0$ in $\Omega$. Let
$\tilde{u}\in H^{1}(\mathds{R}^{d},\mathds{C}^{d})$ be the extension by zero
of $u$. Again $\nabla\tilde{u}=\widetilde{\nabla u}$, the extension by zero of
$\nabla u$, and $\mathop{\rm div}\tilde{u}=0$ in $\mathds{R}^{d}$. Let $v\in
H^{1}(\mathds{R}^{d},\mathds{C}^{d})$ and suppose that $\mathop{\rm div}v=0$
in $\mathds{R}^{d}$. Then

$\displaystyle\int_{\mathds{R}^{d}}\Big{(}\nabla\tilde{u}+\alpha\,(\nabla\tilde{u})^{\rm
T}\Big{)}\cdot\overline{\nabla v}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\Big{(}\nabla u+\alpha\,(\nabla u)^{\rm
T}\Big{)}\cdot\overline{\nabla(v|_{\Omega})}=\gothic{a}(u,v|_{\Omega})$
$\displaystyle=$
$\displaystyle(\lambda\,u,v|_{\Omega})_{L_{2}(\Omega,\mathds{C}^{d})}=(\lambda\,\tilde{u},v)_{L_{2}(\mathds{R}^{d},\mathds{C}^{d})}.$

Since $\mathop{\rm div}\tilde{u}=0$, one deduces that

$\int_{\mathds{R}^{d}}(\nabla\tilde{u})^{\rm T}\cdot\overline{\nabla v}=0.$

Therefore $\tilde{u}\in D(A_{\mathds{R}^{d}})$, the Stokes operator on
$\mathds{R}^{d}$, and $A_{\mathds{R}^{d}}\tilde{u}=\lambda\,\tilde{u}$. Now
one can argue as in the proof of Theorem 3.7 that $u=0$ and hence $u$ is not
an eigenvector.

‘(II)’. Let $\lambda>0$. Let $\omega\in\mathds{R}^{d}$ and suppose that
$|\omega|^{2}=\lambda$. As before let $b\in\mathds{R}^{d}$ and $\tau\in
H^{2}(\Omega)$ be such that $|b|=1$, $b\cdot\omega=0$ and
$\tau(x)=e^{i\omega\cdot x}$. Then
$\Delta\tau=-|\omega|^{2}\,\tau=-\lambda\,\tau$. Let $u=\tau\,b$. Then
$\mathop{\rm div}u=i\tau\,(b\cdot\omega)=0$, so $u\in V$. Let $v\in V$. Then

$\displaystyle\gothic{b}_{\lambda}(u,v)$ $\displaystyle=$
$\displaystyle\int_{\Omega}\nabla u\cdot\overline{\nabla
v}+\alpha\int_{\Omega}(\nabla u)^{\rm T}\cdot\overline{\nabla
v}-\lambda\int_{\Omega}u\cdot\overline{v}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\nabla\tau\cdot\overline{\nabla(b\cdot
v)}+\alpha\sum_{k,l=1}^{d}\int_{\Omega}i\,\omega_{k}\,\tau\,b_{l}\,\overline{\partial_{l}v_{k}}-\lambda\int_{\Omega}\tau\,\overline{(b\cdot
v)}$ $\displaystyle=$
$\displaystyle-\int_{\Omega}(\Delta\tau)\,\overline{(b\cdot
v)}+\int_{\partial\Omega}(\partial_{\nu}\tau)\,\overline{{\mathop{\rm
Tr\,}}(b\cdot v)}$ $\displaystyle\hskip
28.45274pt{}+\alpha\sum_{k,l=1}^{d}\int_{\Omega}\omega_{k}\,\omega_{l}\,\tau\,b_{l}\,\overline{v_{k}}+i\,\alpha\sum_{k,l=1}^{d}\int_{\partial\Omega}\omega_{k}\,({\mathop{\rm
Tr\,}}\tau)\,\nu_{l}\,b_{l}\,\overline{{\mathop{\rm
Tr\,}}v_{k}}-\lambda\int_{\Omega}\tau\,\overline{(b\cdot v)}$ $\displaystyle=$
$\displaystyle\int_{\partial\Omega}(\partial_{\nu}\tau)\,b\cdot\overline{{\mathop{\rm
Tr\,}}v}+i\,\alpha\int_{\partial\Omega}({\mathop{\rm Tr\,}}\tau)\,(\nu\cdot
b)\,\omega\cdot\overline{{\mathop{\rm Tr\,}}v},$

where we used in the last step that $b\cdot\omega=0$. Consequently $j(u)\in
D({\cal N}_{\lambda})$ and
$(j(u),(\partial_{\nu}\tau)\,b+i\alpha\,({\mathop{\rm Tr\,}}\tau)\,(\nu\cdot
b)\,\omega)\in{\cal N}_{\lambda}$. Moreover,

$({\cal
N}_{\lambda}^{\circ}j(u),j(u))_{L_{2}(\Omega,\mathds{C}^{d})}=\gothic{b}_{\lambda}(u,u)\\\
=i\int_{\partial\Omega}(\omega\cdot\nu)\,|b|^{2}\,|{\mathop{\rm
Tr\,}}\tau|^{2},$

where we used once more that $b\cdot\omega=0$. Then argue as in the proof of
Theorem 3.7.

## Acknowledgement

The authors wish to thank Sylvie Monniaux for helpful discussions. The authors
are most grateful for the hospitality and fruitful stay of the first-named
author at the University of Auckland and the second-named author at the Aix-
Marseille Université. This work is partly supported by the Aix-Marseille
Université, an NZ-EU IRSES counterpart fund and the Marsden Fund Council from
Government funding, administered by the Royal Society of New Zealand and the
EU Marie Curie IRSES program, project ‘AOS’, No. 318910.

## References

  * [Alt] Alt, H. W., Linear functional analysis.  Universitext. Springer-Verlag, Berlin etc., 2016.
  * [AE] Arendt, W. and Elst, A. F. M. ter, Sectorial forms and degenerate differential operators.  J. Operator Theory 67 (2012), 33–72.
  * [AEKS] Arendt, W., Elst, A. F. M. ter, Kennedy, J. B. and Sauter, M., The Dirichlet-to-Neumann operator via hidden compactness.  J. Funct. Anal. 266 (2014), 1757–1786.
  * [AM1] Arendt, W. and Mazzeo, R., Spectral properties of the Dirichlet-to-Neumann operator on Lipschitz domains.  Ulmer Seminare 12 (2007), 23–37.
  * [AM2] , Friedlander’s eigenvalue inequalities and the Dirichlet-to-Neumann semigroup.  Commun. Pure Appl. Anal. 11 (2012), 2201–2212.
  * [BE] Batty, C. J. K. and Elst, A. F. M. ter, On series of sectorial forms.  J. Evol. Equ. 14 (2014), 29–47.
  * [BRS] Behrndt, J., Rohleder, J. and Stadler, S., Eigenvalue inequalities for Schrödinger operators on unbounded Lipschitz domains.  J. Spectr. Theory 8 (2018), 493–508.
  * [Fil] Filonov, N., On an inequality between Dirichlet and Neumann eigenvalues for the Laplace operator.  St. Petersburg Math. J. 16 (2005), 413–416.
  * [FL] Frank, R. L. and Laptev, A., Inequalities between Dirichlet and Neumann eigenvalues on the Heisenberg group.  Int. Math. Res. Not. 2010 (2010), 2889–2902.
  * [Fri] Friedlander, L., Some inequalities between Dirichlet and Neumann eigenvalues.  Arch. Rational Mech. Anal. 116 (1991), 153–160.
  * [GM] Gesztesy, F. and Mitrea, M., Nonlocal Robin Laplacians and some remarks on a paper by Filonov on eigenvalue inequalities.  J. Differential Equations 247 (2009), 2871–2896.
  * [Gor] Gordon, G. D., Friedlander’s eigenvalue inequalities, 2014.  BSc(Honours) thesis, Department of Mathematics, University of Auckland, unpublished.
  * [Han] Hansson, A. M., An inequality between Dirichlet and Neumann eigenvalues of the Heisenberg Laplacian.  Comm. Partial Differential Equations 33 (2008), 2157–2163.
  * [Kat] Kato, T., Perturbation theory for linear operators.  Second edition, Grundlehren der mathematischen Wissenschaften 132. Springer-Verlag, Berlin etc., 1980.
  * [Lio] Lions, J. L., Lectures on elliptic partial differential equations.  Lectures on Mathematics 10. Tata Institute of Fundamental Research, Bombay, 1957.
  * [LR] Lotoreichik, V. and Rohleder, J., Eigenvalue inequalities for the Laplacian with mixed boundary conditions.  J. Differential Equations 263 (2017), 491–508.
  * [Maz] Maz’ya, V. G., Sobolev spaces with applications to elliptic partial differential equations.  Second edition, Grundlehren der mathematischen Wissenschaften 342. Springer-Verlag, Berlin etc., 2011.
  * [Maz] Mazzeo, R., Remarks on a paper of Friedlander concerning inequalities between Neumann and Dirichlet eigenvalues.  Internat. Math. Res. Notices 4 (1991), 41–48.
  * [MMW] Mitrea, M., Monniaux, S. and Wright, M., The Stokes operator with Neumann boundary conditions in Lipschitz domains.  J. Math. Sci. 176 (2011), 409–457.
  * [Nit] Nitsche, J. A., On Korn’s second inequality.  RAIRO Anal. Numér. 15 (1981), 237–248.
  * [Pay] Payne, L. E., Inequalities for eigenvalues of membranes and plates.  J. Rational Mech. Anal. 4 (1955), 517–529.
  * [Pól] Pólya, G., Remarks on the foregoing paper.  J. Math. Physics 31 (1952), 55–57.
  * [Pro] Provenzano, L., Inequalities between Dirichlet and Neumann eigenvalues of the polyharmonic operators.  Proc. Amer. Math. Soc. 147 (2019), 4813–4921.
  * [RS] Reed, M. and Simon, B., Methods of modern mathematical physics IV. Analysis of operators.  Academic Press, New York etc., 1978.
  * [Soh] Sohr, H., The Navier–Stokes equation. An elementary functional analytic approach.  Birkhäuser, Basel, 2001.
  * [Tem] Temam, R., Navier–Stokes equations. Theory and numerical analysis.  Studies in mathematics and its applications 2. North-Holland Publishing Company, Amsterdam, 1977.

C. Denis, Centre de Mathématiques et Informatique (CMI) Technopôle Château-
Gombert, 39, rue F. Joliot Curie, 13453 Marseille Cedex 13, France
E-mail address<EMAIL_ADDRESS>
A.F.M. ter Elst, Department of Mathematics, University of Auckland, Private
bag 92019, Auckland 1142, New Zealand
E-mail address<EMAIL_ADDRESS>


# Bayesian Committee Machine Potential for Isothermal-Isobaric Molecular
Dynamics Simulations

Soohaeng Yoo Willow  Department of Energy Science, Sungkyunkwan University,
Seobu-ro 2066, Suwon, 16419, Korea  Chang Woo Myung<EMAIL_ADDRESS>Department of Energy Science, Sungkyunkwan University, Seobu-ro 2066, Suwon,
16419, Korea

###### Abstract

In recent advancements in material simulations, the utilization of Sparse
Gaussian Process Regression (SGPR)-based machine learning potentials (MLPs)
has proven to be highly successful in diverse applications such as catalysis,
batteries, and solar cells. In the context of isothermal and isobaric
molecular dynamics simulations, achieving precise pressure estimates is
crucial for an accurate understanding of the system behavior under constant
pressure conditions. In this study, we introduce a novel kernel function
designed for estimating the virial term, a critical component for pressure
calculations in materials simulations. Our study reveals that the inclusion of
a virial prediction in the kernel function leads to significantly improved
accuracy in calculating the stress of a system. We present a kernel-based ML
potential that can be estimated via the Bayesian Committee Machine without the
need for additional training. This improvement allows us to calculate the
melting temperature of ice through Isobaric-Isenthalpy $NpH$ molecular
dynamics simulations of ice-liquid coexisting phases, employing the BCM
potential.

##  I Introduction

In the last two decades, the application of machine learning potentials (MLPs)
in the field of materials simulations has proven to be highly successful,
particularly for the accurate prediction of potential energy surfaces (PESs)
of the reference ab initio data,[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] such
as density functional theory (DFT),[13, 14, 15] many-body perturbation
(MP),[16, 17, 18] and coupled-cluster (CC) theories. Moreover, MLPs uphold an
affordable computational cost, allowing one to investigate many interesting
physical phenomena with currently accessible supercomputers. This breakthrough
has paved the way for large-scale and long-time simulations in various
scientific domains, including materials science, chemistry, and biology.
Several distinct machine learning methods have been employed for the
construction of MLPs. These include the atom-centered artificial neural
network method (ANN),[19, 20, 21, 22] Gaussian approximation potentials
(GAPs),[3, 23] graph neural network potentials (GNNs),[5, 6] gradient-domain
machine learning (GDML),[24] moment tensor potentials (MTP),[25, 26, 27] and
sparse Gaussian process regression (SGPR).[7, 8, 28, 9, 10, 11] The
versatility of these methodologies has led to their widespread adoption and
successful application across different scientific disciplines, promising
valuable insights and outcomes.

The Gaussian Process (GP) based MLP boasts in the active learning, uncertainty
prediction, and low data requirement.[29, 30, 31, 32] However, despite all
these advantages, the scalability of GP-based ML potential came under scrutiny
in particular for achieving universal potential due to its poor scaling. MLPs
based on GP face computational challenges when dealing with large datasets.
The training complexity increases dramatically to $\mathcal{O}(n^{3})$ due to
the necessity of inverting a kernel matrix, and prediction complexity grows
quadratically with the dataset size $n$. As a result, GPs are restricted to
handling datasets smaller than $10^{4}$. This led to the development of the
SGPR technique,[7, 8, 28, 9, 10] which uses a subset of data for efficiency.
However, SGPR still struggles with large numbers $m$ of inducing points,
showing computational demands of $\mathcal{O}(nm^{2})$. To mitigate this
limitation, the focus shifted to aggregation models such as the Bayesian
Committee Machine (BCM)[33] to address computational challenges. This involves
combining predictions from $p$ distinct sub-models, each trained on different
data sections, with a computational scaling of $\mathcal{O}(nm^{2}/p^{2})$.
This approach is deemed effective in enhancing the scalability of kernel-based
ab initio MLPs.

Molecular dynamics (MD) simulations serve as a powerful tool to complement
experimental findings by offering a detailed atomic-level understanding in
materials. By comparing the results of MD simulations with their corresponding
experimental observations, such as phase diagrams and radial distribution
functions, one can effectively test the accuracy and reliability of the PES.
For accurately performing isothermal-isobaric MD simulations, it is crucial to
correctly estimate the stress. In this context, improving the accuracy of
stress estimations of MLPs can enhance the reliability and effectiveness of MD
simulations. While other types of MLPs such as neural network MLPs [34]
incorporated the estimation of viral stress, a SGPR-based model has not yet
included the estimation of virial stress.

Zhang and coworkers in 2021[35] employed MD simulations with a SCAN-DFT[36]
based Deep potential model[37, 34] to predict ice-water phase diagram,
achieving a generally satisfactory agreement with experimental results.
Subsequently, Bore and Paesani in 2023[38] utilized the deep potential
model[34] based on MB-pol[39, 40, 41] to estimate the ice-water phase diagram
and showed improved agreement with experimental results. Surprisingly,
however, they observed that the phase diagram predicted by the MLP did not
mirror the one obtained using its reference potential, called MB-pol. To
enhance our comprehension of the disparity in predicting phase diagrams
between MLP and the reference model, we performed isobaric-isenthalpy ($NpH$)
MD simulations of ice-liquid coexisting phases to estimate the melting points
of the SGPR-based MLP.

##  II Theory

###  II.1 Virial Kernel prediction in Sparse Gaussian Process Regression
potential

The MLP energy is expressed in terms of a local chemical environment (LCE)
$\rho_{i}$ as

$E=\sum_{i=1}^{N}\mathcal{E}(\rho_{i}),$ (1)

where $\mathcal{E}$ is a fictional latent function of local energy. This
indicates that the MLP energy as well as force acting on atom $i$ depends on
the relative coordinates of surrounding atoms. The LCE of atom $i$ is defined
by

$\rho_{i}=\\{\vec{r}_{ij};j\neq i\quad|\vec{r}_{ij}|<r_{c}\\},$ (2)

where $\vec{r}_{ij}=\vec{r}_{j}-\vec{r}_{i}$, $r_{c}$ is a cutoff radius for
neighborhood relations.

The core component of Bayesian regression is building a covariance Kernel
$\mathcal{K}(\rho_{i},\rho_{j})$ that encodes the similarity between
$\rho_{i}$ and $\rho_{j}$. The Kernel is invariant with respect to the
symmetry operations that leaves the potential energy invariant, e.g.
translations and and rotations. Here, the smooth overlap of atomic positions
(SOAP)[23] is used for the descriptor of a similarity kernel between LCEs and
is defined as

$\mathcal{K}(\rho_{i},\rho_{j})=\int d\bm{R}~{}|\int
d\vec{r}\xi_{i}(\vec{r})\xi_{j}(\bm{R}\vec{r})|^{2},$ (3)

where $\bm{R}$ is the three-dimensional rotation operator, $\xi_{i}(\vec{r})$
is the atomic density in the neighborhood of $i$,

$\xi_{i}(\vec{r})=\sum_{j\in\rho_{i}}e^{\alpha|\vec{r}-\vec{r}_{ij}|^{2}},$
(4)

and $\alpha$ is a hyperparameter.

To improve the scalability, one can exploit the LCE similarities in the data
$X$. It is very likely that many LCE pairs in $X$ are similar:
$\mathcal{K}(\rho_{i},\rho_{j})\approx 1$. Let $z=\\{\chi_{j}\\}_{j=1}^{m}$
denote a reduced set of LCEs which are significantly distinct from each other
and are sufficient statistics for $X=\\{R_{1},\cdots,R_{n}\\}$. $z$ is called
the inducing set of LCEs. With the inducing set of LCEs, the latent energy
function of atom $i$ is expressed as

$\mathcal{E}(\rho_{i})=\sum_{j=1}^{m}w_{j}\mathcal{K}(\rho_{i},\chi_{j})$ (5)

where $\bm{w}=[w_{1},\cdots,w_{n}]$ are the weights for the descriptors in
$z$. The potential energy for $R$ becomes

$E(R)=\sum_{i}^{N}\mathcal{E}(\rho_{i})=\sum_{i,j}w_{j}\mathcal{K}(\rho_{i},\chi_{j}).$
(6)

The forces $f_{i}^{\mu}$ are obtained from

$f_{i}^{\mu}=\sum_{j=1}^{m}w_{j}\left[-\sum_{k=1}\frac{\partial\mathcal{K}(\rho_{k},\chi_{j})}{\partial
r_{i}^{\mu}}\right]=\sum_{j=1}^{m}w_{j}\mathcal{K}_{i}^{\mu}(R,\chi_{j}),$ (7)

where $\mathcal{K}_{i}^{\mu}(R,\chi_{j})$ is the covariance kernel for the
force on atom $i$ and $\mu=\\{x,y,z\\}$.

The previous SGPR potential did not include the stress tensor in the training
process.[7] Due to this, the SGPR potential often fails to accurately predict
the stress tensor of a given system. To improve this shortcoming, our strategy
is to incorporate the stress data into the training process of SGPR. The
instantaneous stress tensor $\mathcal{P}_{\alpha\beta}$ is defined as

$\mathcal{P}_{\alpha\beta}=k_{B}\mathcal{T}_{\alpha\beta}/V+\mathcal{W}_{\alpha\beta}/V,$
(8)

where $\mathcal{T}_{\alpha\beta}$ is the instantaneous temperature tensor. To
predict the stress tensor, we build the internal virial kernel function
$\mathcal{W}_{\alpha\beta}$ as below

$\displaystyle\mathcal{W}_{\alpha\beta}$ $\displaystyle=$
$\displaystyle\sum_{i}f_{i}^{\alpha}r_{i}^{\beta}$ (9) $\displaystyle=$
$\displaystyle\sum_{j}w_{j}\left[-\sum_{i,k}\frac{\partial\mathcal{K}(\rho_{k},\chi_{j})}{\partial
r_{i}^{\alpha}}r_{i}^{\beta}\right]$

The weight vector $\bm{w}$ is obtained by linear fitting of $E(R)$,
$f_{k}^{\mu}$, and $\mathcal{W}_{\alpha\beta}$ to the target energies, forces,
virial stresses $\bm{y}_{n}$ for a given dataset $X=\\{R_{1},\cdots,R_{n}\\}$.
In the Bayesian linear regression method, the mathematical expression that we
wish to solve becomes

$\bm{K}_{nm}\bm{w}_{m}=\bm{y}_{n},$ (10)

where $\bm{K}_{nm}$ is the covariance matrix between data X and inducing set
$z$ of LCEs. In the SGPR algorithm,[7] $\bm{w}_{m}$ can be the solution of the
following linear system:

$\begin{bmatrix}K_{nm}\\\ \sigma
L^{T}\end{bmatrix}\bm{w}_{m}=\begin{bmatrix}\bm{y}_{n}\\\
\bm{0}\end{bmatrix},$ (11)

where $L$ is the Cholesky factor of $K_{mm}$ ($K_{mm}=LL^{T}$, where $L$ is
lower triangular), $\bm{0}$ is a columnar vector of zeros with length $m$, and
$\sigma$ is the noise hyperparameter.

Figure 1: Schematic illustration of Bayesian committee machine potential. The
entire dataset is divided into the distinct subsets based on classification.
Within each subset, an SGPR-based MLP is trained. Subsequently, the universal
MLP for the entire dataset is estimated through the BCM. The weight factor
$\frac{\beta_{\alpha}}{\sigma_{\alpha}^{2}}$ ensures that the predictions from
confident models are more outspoken than those from less confident models.

###  II.2 Bayesian Committee Machine potential

The SGPR algorithm encounters difficulties in obtaining the solution
$\bm{w}_{m}$ for Eq. (11), which involves the inversion of a kernel matrix,
with large size of training datasets $X$ and inducing set $z$ of LCEs. In
order to mitigate these challenge, we employed the BCM[33] as a novel approach
to combine SGPR-MLPs which were separately trained on different subsets as
shown in Figure 1. Taking the example of ice polymorphs, instead of utilizing
all ice phases as a collective training dataset, we assign each ice phase to
an individual training dataset. Subsequently, we trained separate models for
each subset to create machine learning potential models tailored to specific
phases. We consolidated these individually trained model into an universal MLP
via BCM. This approach not only streamlines the training process but also
addresses the computational challenges associated with the inversion of large
kernel matrices, paving the way for an efficient and accurate universal MLPs.

With classification of the training dataset and inducing set into the distinct
subsets: $\\{n_{1},\cdots,n_{p}\\}$ and $\\{m_{1},\cdots,m_{p}\\}$, in the
BCM,[33] the predicted universal potential energy of the system becomes

$\displaystyle
E\approx\hat{S}\sum_{\alpha=1}^{p}\frac{\beta_{\alpha}}{\sigma_{\alpha}^{2}}E_{\alpha},$
(12)

where the potential energy $E_{\alpha}$ is calculated using $\alpha$-th subset
($R_{n\in\alpha}$ and $z_{\alpha}=\\{\chi_{j}\\}_{j\in\alpha})$ and the weight
vector $\bm{w}^{\alpha}$ of $\alpha$-th sub SGPR model. This energy can be
represented as:

$\displaystyle E_{\alpha}(R)$ $\displaystyle\approx$
$\displaystyle\sum_{i,(j\in m_{\alpha})}w_{j}^{\alpha}K(\rho_{i},\chi_{j}).$
(13)

The maximum value of the covariance loss for the $\alpha$-th subset, denoted
as $\sigma_{\alpha}^{2}$, is used to weight the $\alpha$-th committee
prediction. The covariance loss $s(\rho_{i})$ is defined as

$\sigma_{\alpha}^{2}=\mathrm{max}(1-K_{\rho_{i},m}K_{mm}^{-1}K_{\rho_{i},m}^{T}).$
(14)

Furthermore, we weight the prediction with $\beta_{\alpha}$ following the
concepts of robust BCM where the individual committee prediction is weighted
by the differential entropy, $\beta_{\alpha}=-\log(\sigma_{\alpha}^{2})$. The
final normalization factor $\hat{S}$ is thus given:

$\displaystyle\hat{S}$ $\displaystyle=$
$\displaystyle\left[\sum_{\alpha}^{p}\frac{\beta_{\alpha}}{\sigma_{\alpha}^{2}}\right]^{-1}.$
(15)

The significant benefit of employing the BCM with SGPR lies in its ability to
work with a reduced $(n/p)$-dimensional training dataset and a
$(m/p)$-dimensional inducing set, instead of the original $n$-dimensional
training dataset and $m$-dimensional inducing set, where
$n=\sum_{\alpha}n_{\alpha}$ and $m=\sum_{\alpha}m_{\alpha}$. Consequently, the
computational cost in BCM-based SGPR scales as from $\mathcal{O}(nm^{2})$ to
$p\;\mathcal{O}(nm^{2}/p^{3})$. Another advantage of the BCM is that a
particular subset, exhibiting greater similarity with a test configuration
$x^{*}$, contributes more significantly to energy, forces, and stress.

Compared to “product of experts” (POE),[42] in which the predicted universal
potential energy reads $E\approx\frac{1}{P}\sum_{\alpha}^{P}E_{\alpha}$, the
weight factor $\beta_{\alpha}/\sigma_{\alpha}^{2}$ in the BCM ensures that the
predictions from confident models showing higher similarity are given more
emphasis that those from less confident models.

Figure 2: Comparison of SGPR-MLP’s predicted pressure with reference values.
The testing datasets were gathered from molecular dynamics simulations of
systems where ice III and liquid coexist, using MB-pol water at pressures of
$P$ = 0.4 GPa and $P$ = 0.6 GPa in isobaric isothermal conditions.

##  III Results

###  III.1 Pressure prediction using virial Kernel method

In Section II.1, the weight vector $\bm{w}_{m}$ relies on the covariance
matrix $K_{nm}$ and corresponding target values for energy, forces, virial
stress denoted as $\bm{y}_{n}$. The weight vector’s values vary based on
whether the covariance matrix and corresponding target values include the
virial stress term. In this subsection, our objective is to compare
instantaneous pressure calculated using the weight vector obtained under
different conditions. This comparative analysis seeks to provide insights into
the impact of incorporating or omitting the virial stress term in covariance
matrix and target values on estimation of instantaneous pressure. To achieve
this, we explore two schemes. The first scheme (Scheme 1) involves
constructing $\bm{w}_{m}$ solely with $K_{nm}$ and $\bm{y}_{n}$ including only
both energies and forces. The second scheme (Scheme 2) incorporates virial
term into $K_{nm}$ and $\bm{y}_{n}$ in addition to energies and forces to
obtain $\bm{w}_{m}$.

To obtain the training and testing dataset, we performed isobaric isothermal
molecular dynamics simulations of the system where ice III and liquid coexist,
using MB-pol water at pressures of $P=0.4$ and $0.6$ GPa and temperature of
$T=255$ K. Once we trained the weight vector $\bm{w}_{m}$ based on the
training dataset, we estimated the pressures with the testing dataset (see
Figure 2) using Eqs. (8) and (9). Figure 2 shows that only Scheme 2 could
reproduce target instantaneous pressures accurately.

Figure 3: Ice-liquid melting points were calculated through $NpH$ MD
simulations of coexisting ice-liquid systems. The melting points obtained
using the MB-pol model are denoted by circle points, and those obtained with
the universal SGPR-MLP are indicated by star points. The melting points
calculated by Bore and Paesani[38] are represented by open square points.
Solid black lines represent the experimental melting lines.

###  III.2 Estimation of the coexistence line of ice-liquid

$NpH$ molecular dynamics simulations were employed to compute the melting
temperature ($T_{m}$) of ice within the coexisting ice-liquid phase.[43, 44,
45] $NpH$ MD simulations enable the spontaneous adjustment of temperature to
satisfy the condition
$\mu_{\mathrm{ice}}(P,T)_{T=T_{m}}=\mu_{\mathrm{liq}}(P,T)_{T=T_{m}}$.
Specifically, when the temperature of ice-liquid coexisting system is higher
than $T_{m}$, the chemical potential of the ice phase is higher than that of
the liquid phase:
$\mu_{\mathrm{ice}}(P,T)_{T>T_{m}}>\mu_{\mathrm{liq}}(P,T)_{T>T_{m}}$.
Consequently, the ice phase undergoes a process of melting with the absorption
of heat, leading to a reduction in the kinetic temperature of the coexisting
system towards the melting temperature $T_{m}$. One notable advantage of
employing $NpH$ MD simulations over $NVE$ MD simulation is the absence of the
stress-anisotropy problem, as the three principle components of the stress
tensor can be adjusted to the specified external pressure in the former.

We aim to estimate the melting temperatures of ice II, III, and V under low
pressure conditions (0.2 $\sim$ 0.6 GPa). To achieve this, on-the-fly active
learning using isobaric isothermal ($NpT$) MD simulations were performed to
create a SGPR MLP for each phases (ice II, ice III, and ice V) where ice and
liquid coexist. The MB-pol water model[41] was chosen for the reference model
in this study. All MD simulations were carried out using ASE package. We
employed a time step of 0.5 fs for the integration of the equations of motion
that were propagate according to combined Nose-Hoover and Parrinello-Rahman
dynamics.[46, 47, 48]

Three individually trained SGPR MLP models were merged into one universal SGPR
MLP model for ice-liquid coexisting systems using the BCM. $NpH$ MD
simulations with this universal model were performed to determine the melting
temperatures of ice II, III, and V. As shown in Figure S1, the $NpH$ MD
simulations maintained the conservation of enthalpy ($H$) while the total
energy as the sum of the potential energy and kinetic energy exhibited
fluctuations. Figure 3 shows that melting points of ice II, III, and V, as
obtained from isobaric isenthalpy ($NpH$) MD simulations with the MB-pol water
model, largely coincide with those calculated from enhanced-coexistence
simulations.[38] The universal SGPR-MLP gives melting points that exhibit a
slight shift towards higher values compared to those obtained using the
reference potential (MB-pol). This disparity arise due to slight noise in
energy and forces during training of SGPR-MLP. In general, aside from noting a
shift in the phase diagram, it is expected that the well-trained SGPR-MLP can
faithfully reproduce physical and thermodynamic properties of materials using
the reference potential models.

Deep neural network ML models[35, 38] based on DFT and MB-pol predict that ice
III is metastable when compared to ice II. Our SGPR ML model based on MB-pol
exhibits the same trend, indicating that ice III is indeed metastable in
comparison to ice II. The discrepancy between these ML models and experimental
observations may arise from the omission of entropy resulting from proton
disorder in ice III.

##  IV Discussion

The focus of our investigation lies in the critical aspect of achieving
precise pressure estimates, particularly in the context of isothermal and
isobaric molecular dynamics simulations, as accurate pressure information is
paramount for a thorough understanding of system behavior. A pivotal
contribution of our work involves the introduction of a novel kernel function
explicitly tailored for estimating the virial term, a key factor in pressure
calculations during molecular dynamics simulations. Our findings underscore
the significance of incorporating a virial term correction within the kernel
function, leading to a substantial improvement in the accuracy of system
pressure calculations. This refinement proves instrumental in enhancing the
overall reliability and precision of the MLP, particularly when simulating
material systems under isothermal and isobaric conditions.

Furthermore, the refined MLP with the virial term correction enables us to
extend its applicability to estimating the melting temperature of ice
polymorphs. We achieve this by employing isobaric-isenthalpy ($NpH$) molecular
dynamics simulations that capture the coexisting phases of ice and liquid. The
universal SGPR-MLP, now augmented with the virial term correction,
demonstrates its efficacy in providing accurate predictions of the melting
temperature for ice. This not only highlights the versatility of our proposed
methodology but also underscores its potential in advancing the understanding
of material behavior under specific thermodynamic conditions.

In the end, the study highlights the efficacy of the combining model (BCM)
generated through kernel-based methods, such as SGPR, for enhanced generality
and transferability. By employing the BCM, a universal SGPR-MLP can be built
by leveraging previously trained models from a database. In comparison to the
conventional “product of experts” (POE) method, the universal SGPR-MLP
obtained via the BCM demonstrates superior predictive accuracy for potential
energy, forces, and pressures. This indicates that the BCM offers a more
robust and reliable framework for broader applications in the field of
machine-learning potentials.

##  V Acknowledgements

SYW and CWM acknowledge the support from the National Research Foundation of
Korea (NRF) grant funded by the Korea government (MSIT) (No.
RS-2023-00222245). CWM acknowledges the support from the National Research
Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No.
NRF-2022R1C1C1010605). The authors are grateful for: the computational support
from the Korea Institute of Science and Technology Information (KISTI) for the
Nurion cluster (KSC-2023-CRE-0059,KSC-2023-CRE-0355,KSC-2023-CRE-0454).

## References

  * Unke  _et al._ [2021] O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko, and K.-R. Müller, Machine Learning Force Fields, Chem. Rev. 121, 10142 (2021).
  * Deringer  _et al._ [2021] V. L. Deringer, A. P. Bartók, N. Bernstein, D. M. Wilkins, M. Ceriotti, and G. Csányi, Gaussian Process Regression for Materials and Molecules, Chem. Rev. 121, 10073 (2021).
  * Bartók  _et al._ [2010] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons, Phys. Rev. Lett. 104, 136403 (2010).
  * Chen  _et al._ [2023] M. S. Chen, J. Lee, H.-Z. Ye, T. C. Berkelbach, D. R. Reichman, and T. E. Markland, Data-Efficient Machine Learning Potentials from Transfer Learning of Periodic Correlated Electronic Structure Methods: Liquid Water at AFQMC, CCSD, and CCSD(T) Accuracy, J. Chem. Theory Comput. 19, 4510 (2023).
  * Gilmer  _et al._ [2017] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl, Neural Message Passing for Quantum Chemistry (2017), arxiv:1704.01212 [cs] .
  * Batzner  _et al._ [2022] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials, Nat Commun 13, 2453 (2022).
  * Hajibabaei  _et al._ [2021a] A. Hajibabaei, C. W. Myung, and K. S. Kim, Sparse Gaussian process potentials: Application to lithium diffusivity in superionic conducting solid electrolytes, Phys. Rev. B 103, 214102 (2021a).
  * Hajibabaei and Kim [2021] A. Hajibabaei and K. S. Kim, Universal Machine Learning Interatomic Potentials: Surveying Solid Electrolytes, J. Phys. Chem. Lett. 12, 8115 (2021).
  * Hajibabaei  _et al._ [2021b] A. Hajibabaei, M. Ha, S. Pourasad, J. Kim, and K. S. Kim, Machine Learning of First-Principles Force-Fields for Alkane and Polyene Hydrocarbons, J. Phys. Chem. A 125, 9414 (2021b).
  * Myung  _et al._ [2022] C. W. Myung, A. Hajibabaei, J.-H. Cha, M. Ha, J. Kim, and K. S. Kim, Challenges, Opportunities, and Prospects in Metal Halide Perovskites from Theoretical and Machine Learning Perspectives, Adv.Energy Mater. 12, 2202279 (2022).
  * Vandermause  _et al._ [2022] J. Vandermause, Y. Xie, J. S. Lim, C. J. Owen, and B. Kozinsky, Active learning of reactive Bayesian force fields applied to heterogeneous catalysis dynamics of H/Pt, Nat Commun 13, 5183 (2022).
  * Vandermause  _et al._ [2020] J. Vandermause, S. B. Torrisi, S. Batzner, Y. Xie, L. Sun, A. M. Kolpak, and B. Kozinsky, On-the-fly active learning of interpretable Bayesian force fields for atomistic rare events, npj Comput Mater 6, 20 (2020).
  * Laasonen  _et al._ [1993] K. Laasonen, M. Sprik, M. Parrinello, and R. Car, ”Ab initio” liquid water, J. Chem. Phys. 99, 9080 (1993).
  * Richards  _et al._ [2016] W. D. Richards, T. Tsujimura, L. J. Miara, Y. Wang, J. C. Kim, S. P. Ong, I. Uechi, N. Suzuki, and G. Ceder, Design and synthesis of the superionic conductor Na10SnP2S12, Nat Commun 7, 11009 (2016).
  * Boero  _et al._ [1998] M. Boero, M. Parrinello, and K. Terakura, First Principles Molecular Dynamics Study of Ziegler-Natta Heterogeneous Catalysis, J. Am. Chem. Soc. 120, 2746 (1998).
  * Willow  _et al._ [2015] S. Y. Willow, M. A. Salim, K. S. Kim, and S. Hirata, Ab initio molecular dynamics of liquid water using embedded-fragment second-order many-body perturbation theory towards its accurate property prediction, Sci Rep 5, 14358 (2015).
  * Del Ben  _et al._ [2013] M. Del Ben, M. Schönherr, J. Hutter, and J. VandeVondele, Bulk Liquid Water at Ambient Temperature and Pressure from MP2 Theory, J. Phys. Chem. Lett. 4, 3753 (2013).
  * Del Ben  _et al._ [2015] M. Del Ben, J. Hutter, and J. VandeVondele, Forces and stress in second order Møller-Plesset perturbation theory for condensed phase systems within the resolution-of-identity Gaussian and plane waves approach, J. Chem. Phys. 143, 102803 (2015).
  * Behler and Parrinello [2007] J. Behler and M. Parrinello, Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces, Phys. Rev. Lett. 98, 146401 (2007).
  * Behler [2011] J. Behler, Atom-centered symmetry functions for constructing high-dimensional neural network potentials, J. Chem. Phys. 134, 074106 (2011).
  * Pun  _et al._ [2019] G. P. P. Pun, R. Batra, R. Ramprasad, and Y. Mishin, Physically informed artificial neural networks for atomistic modeling of materials, Nat Commun 10, 2339 (2019).
  * Eckhoff and Behler [2021] M. Eckhoff and J. Behler, High-dimensional neural network potentials for magnetic systems using spin-dependent atom-centered symmetry functions, Npj Comput. Mater. 7, 170 (2021).
  * Bartók  _et al._ [2013] A. P. Bartók, R. Kondor, and G. Csányi, On representing chemical environments, Phys. Rev. B 87, 184115 (2013).
  * Chmiela  _et al._ [2017] S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt, and K.-R. Müller, Machine learning of accurate energy-conserving molecular force fields, Sci. Adv. 3, e1603015 (2017).
  * Shapeev [2016] A. V. Shapeev, Moment Tensor Potentials: A Class of Systematically Improvable Interatomic Potentials, Multiscale Model. Simul. 14, 1153 (2016).
  * Podryabinkin and Shapeev [2017] E. V. Podryabinkin and A. V. Shapeev, Active learning of linearly parametrized interatomic potentials, Comput. Mater. Sci. 140, 171 (2017).
  * Novikov  _et al._ [2021] I. S. Novikov, K. Gubaev, E. V. Podryabinkin, and A. V. Shapeev, The MLIP package: Moment tensor potentials with MPI and active learning, Mach. Learn.: Sci. Technol. 2, 025002 (2021).
  * Quiñonero-Candela and Rasmussen [2005] J. Quiñonero-Candela and C. E. Rasmussen, A Unifying View of Sparse Approximate Gaussian Process Regression, J. Mach. Learn. Res. 6, 1939 (2005).
  * Rasmussen and Williams [2005] C. E. Rasmussen and C. K. I. Williams,  _Gaussian Processes for Machine Learning_ (The MIT Press, 2005).
  * Shahriari  _et al._ [2016] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. De Freitas, Taking the Human Out of the Loop: A Review of Bayesian Optimization, Proc. IEEE 104, 148 (2016).
  * Lawrence [2005] N. Lawrence, Probabilistic Non-linear Principal Component Analysis with Gaussian Process Latent Variable Models, J. Mach. Learn. Res. 6, 1783 (2005).
  * Álvarez and Lawrence [2011] M. A. Álvarez and N. D. Lawrence, Computationally Efficient Convolved Multiple Output Gaussian Processes, J. Mach. Learn. Res. 12, 1459 (2011).
  * Tresp [2000] V. Tresp, A Bayesian Committee Machine, Neural Computation 12, 2719 (2000).
  * Wang  _et al._ [2018] H. Wang, L. Zhang, J. Han, and W. E, DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics, Comput. Phys. Commun. 228, 178 (2018).
  * Zhang  _et al._ [2021] L. Zhang, H. Wang, R. Car, and W. E, Phase Diagram of a Deep Potential Water Model, Phys. Rev. Lett. 126, 236001 (2021).
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly Constrained and Appropriately Normed Semilocal Density Functional, Phys. Rev. Lett. 115, 036402 (2015).
  * Zhang  _et al._ [2019] L. Zhang, D.-Y. Lin, H. Wang, R. Car, and W. E, Active learning of uniformly accurate interatomic potentials for materials simulation, Phys. Rev. Mater. 3, 023804 (2019).
  * Bore and Paesani [2023] S. L. Bore and F. Paesani, Realistic phase diagram of water from ”first principles” data-driven quantum simulations, Nat Commun 14, 3349 (2023).
  * Babin  _et al._ [2013] V. Babin, C. Leforestier, and F. Paesani, Development of a “First Principles” Water Potential with Flexible Monomers: Dimer Potential Energy Surface, VRT Spectrum, and Second Virial Coefficient, J. Chem. Theory Comput. 9, 5395 (2013).
  * Babin  _et al._ [2014] V. Babin, G. R. Medders, and F. Paesani, Development of a “First Principles” Water Potential with Flexible Monomers. II: Trimer Potential Energy Surface, Third Virial Coefficient, and Small Clusters, J. Chem. Theory Comput. 10, 1599 (2014).
  * Riera  _et al._ [2023] M. Riera, C. Knight, E. F. Bull-Vulpe, X. Zhu, H. Agnew, D. G. A. Smith, A. C. Simmonett, and F. Paesani, MBX: A many-body energy and force calculator for data-driven many-body simulations, J. Chem. Phys. 159, 054802 (2023).
  * Hinton [2002] G. E. Hinton, Training Products of Experts by Minimizing Contrastive Divergence, Neural Comput. 14, 1771 (2002).
  * Yoo  _et al._ [2004] S. Yoo, X. C. Zeng, and J. R. Morris, The melting lines of model silicon calculated from coexisting solid–liquid phases, The Journal of Chemical Physics 120, 1654 (2004).
  * Yoo  _et al._ [2009] S. Yoo, X. C. Zeng, and S. S. Xantheas, On the phase diagram of water with density functional theory potentials: The melting temperature of ice Ih with the Perdew–Burke–Ernzerhof and Becke–Lee–Yang–Parr functionals, The Journal of Chemical Physics 130, 221102 (2009).
  * Yoo and Xantheas [2011] S. Yoo and S. S. Xantheas, Communication: The effect of dispersion corrections on the melting temperature of liquid water, The Journal of Chemical Physics 134, 121105 (2011).
  * Melchionna  _et al._ [1993] S. Melchionna, G. Ciccotti, and B. L. Holian, Hoover NPT dynamics for systems varying in shape and size, Mol. Phys. 78, 533 (1993).
  * Melchionna [2000] S. Melchionna, Constrained systems and statistical distribution, Phys. Rev. E 61, 6165 (2000).
  * Holian  _et al._ [1990] B. L. Holian, A. J. De Groot, W. G. Hoover, and C. G. Hoover, Time-reversible equilibrium and nonequilibrium isothermal-isobaric simulations with centered-difference Stoermer algorithms, Phys. Rev. A 41, 4552 (1990).