# Weyl asymptotics for Poincaré-Steklov eigenvalues in a domain with Lipschitz
boundary

Grigori Rozenblum  Chalmers Univ. of Technology; The Euler International
Mathematical Institute and St.Petersburg State University
<EMAIL_ADDRESS>
###### Abstract.

We justify the Weyl asymptotic formula for the eigenvalues of the Poincaré-
Steklov spectral problem for a domain bounded by a Lipschitz surface.

The work of G.R. was performed at the Saint Petersburg Leonhard Euler
International Mathematical Institute and supported by the Ministry of Science
and Higher Education (Agreement No. 075–15–2022–287).

In the memory of Mikhail Birman and Mikhail Solomyak, my teachers.

##  1\. Introduction

###  1.1. The Poincaré-Steklov problem.

Let $\Omega\subset{\mathbb{R}}^{d+1}$ be a connected bounded open set. We
suppose that the boundary $\Sigma=\partial\Omega$ is connected. The classical
_Poincaré-Steklov (P-S) eigenvalue problem_ consists in the study of the
spectral properties for the Laplacian in $\Omega$ with the spectral parameter
entering in the boundary condition,

$-\Delta
u(x)=0,\,x\in\Omega;\partial_{\nu}u(x)-\lambda^{-1}u(x)=0,\,x\in\Sigma,$ (1.1)

where $\partial_{\nu}$ is the derivative in the direction of the exterior
normal $\nu(x)$ at $x\in\Sigma.$ Further on, we assume that the boundary
$\Sigma$ is Lipschitz and we equip $\Sigma$ with the surface measure induced
by the Lebesgue measure in ${\mathbb{R}}^{d+1},$ which coincides with the
$d$-dimensional Hausdorff measure.

This problem, often attributed to and named after V.A. Steklov, [40], was, in
fact, first considered by H. Poincaré in 1896, see [32], in relation to the
analysis of tidal waves. Another early application, dealing with liquid waves
in an open container was studied by D. Hilbert in [22].

The P-S spectral problem and its numerous generalizations found a lot of
applications in physics and technology (see, e.g., [7] and the quite recent
review [17], for a, far from complete, bibliography on this topic.) Studies
dealing with this problem continue up to now. The Steklov eigenvalue problem
appears in quite a few physical fields, such as fluid mechanics,
electromagnetism, elasticity, etc. It has applications for the study of
various kinds of wave phenomena (see, e.g. [42]), as well as in the seismology
and tomography. Mathematically, this problem keeps being in the center of
interest in spectral geometry and approximations. MathSciNet shows more than
600 publications where various facets of Steklov type spectral problems are
dealt with.

One of traditional topics in this field is the study of the behavior of
eigenvalues of the P-S problem. In the earliest paper by L. Sandgren on the
asymptotics of these eigenvalues (see [39]), inspired by Lars
G$\overset{{}_{{}_{\circ}}}{\mathrm{a}}$rding and $\AA$ke Plejel, with co-
operation of Lars Hörmander, the author considered a moderately smooth domain
$\Omega$ in a Riemannian manifold. The boundary of the domain was supposed be
of the class $C^{2}$, but the coefficients of the main operator were set to
belong to $C^{1}.$ Under these conditions, the asymptotic formula of H. Weyl
type was obtained. Further on, the conditions imposed on the domain and
coefficients were being gradually relaxed. Finally, in 2006, M.S. Agranovich
[2] established the Weyl type formula for the case of a Lipschitz boundary,
being, however, ’almost smooth’ in the sense that this boundary should be at
least $C^{1}$ outside a closed set of zero surface measure; the symmetric
second order elliptic operator ${\mathcal{L}}$ replacing the Laplacian in
(1.1) is supposed to be of divergence form with continuous coefficients. Even
earlier, probably, since [6], M.S. Agranovich started popularizing the problem
of extending the eigenvalue asymptotics results to domains with Lipschitz
boundary, without additional assumptions. This question was morally supported
by a flow of impressive results on various properties of elliptic boundary
problems in Lipschitz domains by A.P.Calderón, B.Dahlberg, C.Kenig, S.Hoffman,
M.Taylor, M.Mitrea, E.Fabes, G.Verchota, M.S.Agranovich himself, and many
others. It turned out that the setting of Lipschitz boundary problems and the
typical results are usually more or less the same as for a more smooth
boundary, while the methods used in analysis need to be essentially new – and
often quite complicated. This latter circumstance is, partially, caused by the
fact that certain important results in potential theory break down when
passing from $C^{1}$ to Lipschitz surfaces.

Meanwhile, the study of the P-S problem continued. We address the Readers to
[8], [9], [10], [18], and references therein, where various aspects of this
problem were investigated. Important results were obtained even for boundaries
considerably less regular than the Lipschitz ones, see, e.g., [43]. Quite
extensive became studies in the spectral geometry relating the geometric
properties of the boundary to spectral properties. Especially rich were the
results in the two-dimensional case. Here, it was known since long ago, see
[35], that in the infinitely smooth case, the Steklov eigenvalues are, faster
than any negative power of their sequential number, close to those for the
disk with the same perimeter. Quite a lot of further results in the two-
dimensional case were obtained since then. A description and huge bibliography
can be found in [20], [17], and [29]. We just mention here that in a series of
papers, starting with [30], for the case of a piecewise smooth Lipschitz
boundary, the authors succeeded in describing how the corners influence the
deviation of the eigenvalues from their behavior in the smooth case.
Meanwhile, in the Lipschitz case, many important properties of the P-S (or the
Neumann-to-Dirichlet, N-to-D, $N\\!\\!\\!D$ operator) were established, see,
e.g., [8], [9], [10], [19]; a number of important Steklov-type problems
arising in hydrodynamics were considered recently in [42].

Finally, after all these years, the essential progress in the eigenvalue
asymptotics for a Lipschitz boundary was made in [26], where, in the two-
dimensional case, the Weyl eigenvalue asymptotics for the problem (1.1) was
proved. Moreover, in [26], the boundary of the domain may be even somewhat
more rough than Lipschitz. The method in [26] is based upon some deep results
in the theory of conformal mappings, therefore, seemingly, it cannot be
extended to higher dimensions. Thus, the eigenvalue asymptotics problem for
Lipschitz boundary remained unresolved. In the recent fundamental review [17]
and in [19], this problem is listed among unsolved and challenging.

The author of the present paper, together with Grigory Tashchiyan, spent quite
a lot of time attacking this problem, being initially inspired by M.S.
Agranovich. Our hope was based upon our result in [37], where, for integral
operators of the type of single layer potential on a Lipschitz surface, a Weyl
type eigenvalue asymptotic formula was proved. Since the P-S operator can be
expressed in a simple way via the single layer and double layer potentials
(see, e.g. [4], [3], and in a very general setting, [43]), it would be
sufficient, for example, to know that the double layer potential is a compact
operator. Unfortunately, for a Lipschitz surface, one should not expect this
property of the double layer potential to hold; in fact, the opposite is
known. Our efforts to circumvent this obstacle led to some partial results
(see a mentioning in [42]); we, however, did not consider these partial
results deserving being published (a brief description of our efforts here, is
given in Appendix, in hope that someone might be interested and more lucky –
if successful, such proof will, probably, be more esthetic than the present
one).

A different approach is used in _this_ paper, not based directly upon the
results in [37] but rather on the _idea_ in [37], a quite natural one, of
approximating the Lipschitz surface by smooth surfaces and tracing how the
corresponding compact Neumann-to-Dirichlet operators converge. In the study of
solvability of elliptic boundary problems in Lipschitz domains, such
approximation was, probably, first used by G.Verchota in [45]. In fact, a
somewhat different realization of this approximation idea is implemented here.
By a change of variables, the P-S problem in a Lipschitz domain for an
elliptic operator in divergence form with coefficients continuous at the
boundary, is relocated to the problem for another elliptic operator, this time
in a smooth domain. Unfortunately, under such transformation, the coefficients
of the operator _cease to be continuous_ even at the boundary, but stay only
bounded. Functions in $L_{\infty}$ cannot be approximated by smooth functions
in $L_{\infty}$ norm. We, however, construct an approximation of the
coefficients in a certain $L_{p}$-related norm, $p<\infty,$ and after a series
of further transformations, we succeed in establishing a sufficiently strong
convergence of operators describing the P-S spectrum, which enables us to
perform the passage to the limit in eigenvalue asymptotic formulas.

The reasoning in the paper is, to a large extent, based upon the ideas the
author absorbed many years ago while being a student of M.Sh. Birman and M.Z.
Solomyak; the paper is dedicated to their memory. The author thanks Prof.
Tatiana Suslina for useful discussions and Prof. Jean Lagacé who acquainted
him with an early version of the paper [26], which encouraged our efforts.
Permanent discussions with Prof. Grigory Tashchiyan for at least 15 years
contributed a lot to a better understanding of the problem and stimulated the
author, who is deeply grateful to Grigory for his involvement.

###  1.2. Setting

We are going to study the distribution of eigenvalues and singular numbers of
various kinds of compact operators. By $n_{\pm}(\lambda,{\mathbf{K}})$ we
denote the distribution function of positive (negative) eigenvalues of the
_self-adjoint_ operator ${\mathbf{K}}$; for an arbitrary compact operator
${\mathbf{K}}$, $n(\lambda,{\mathbf{K}})$ denotes the counting function of
singular numbers of ${\mathbf{K}}.$ In the case when the operator is described
by a quadratic form $\mathfrak{B}[u]$ in the Hilbert space with norm defined
by means of the quadratic form $\mathfrak{A}[u]$, we replace ${\mathbf{K}}$ by
the ratio of these quadratic forms in the notation $n(\lambda,.)$,
$n_{\pm}(\lambda,.)$. For $\theta>0$, we denote by
${\mathbf{n}}^{\sup}_{\pm}(\theta,{\mathbf{K}})$ the quantity
$\limsup_{\lambda\to 0}\lambda^{\theta}n_{\pm}(\lambda,{\mathbf{K}}).$ Without
the subscript $\pm$, these notations concern the distribution of singular
numbers of these operators. Without the superscript $\sup$, these symbols
denote the limits $\lim_{\lambda\to
0}\lambda^{\theta}n_{\pm}(\lambda,{\mathbf{K}}),$ provided these limits exist.
Finally, if some operator or the ratio of quadratic form is defined in the
displayed formula with tag (X.Y), the above notation is used with this tag
replacing the notation of the operator, e.g.,
${\mathbf{n}}^{\sup}(\theta,(X.Y))$, etc.

We consider a bounded domain $\Omega\subset{\mathbb{R}}^{d+1}$ with connected
Lipschitz boundary $\Sigma$. (The connectedness condition can be easily
removed, at the cost of certain notational complications.) This means that
$\Sigma=\partial\Omega$ can be covered by a finite collection of co-ordinate
neighborhoods ${\mathcal{U}}_{\iota}$ such that the portion of $\Sigma$ in
${\mathcal{U}}_{\iota}$ can be, in a properly rotated Euclidean co-ordinate
system, $(x^{\prime},x_{d+1})\in{\mathcal{U}}_{\iota}$, represented by the
equation
$x_{d+1}=\psi^{(\iota)}(x^{\prime}),\,x^{\prime}\in{\mathcal{V}}\subset{\mathbb{R}}^{d}$
with Lipschitz function $\psi^{(\iota)}$ (sometimes such surfaces are called
strongly Lipschitz in the literature, see, e.g., [23].) It is convenient to
use a global co-ordinate system in a neighborhood of the boundary, see Sect.3.

In [39], [2], and some other papers, a more general setting is considered,
namely the spectral problem (1.1), with the Laplacian in the equation replaced
by some formally self-adjoint second order elliptic differential operator
${\mathcal{L}}$ in divergence form, with the normal derivative
$\partial_{\nu}$ replaced by the derivative along the conormal associated with
${\mathcal{L}}$ and with a weight function present on the right-hand side of
the eigenvalue equation. For us, it is also convenient to consider such more
general, the _weighted_ $N\\!\\!\\!D$, Neumann-to-Dirichlet, problem. Namely,
for a uniformly positive definite (elliptic) real matrix-function
$\mathfrak{a}(x)=(a_{j,k}(x),\,j,k=1,\dots,d+1)$, $\mathfrak{a}\in
L_{\infty}(\Omega),$ we denote by ${\mathcal{L}}={\mathcal{L}}_{\mathfrak{a}}$
the formal differential operator

${\mathcal{L}}\equiv{\mathcal{L}}_{\mathfrak{a}}=-\sum_{j,k}\partial_{j}a_{jk}(x)\partial_{k}.$
(1.2)

The weighted $N\\!\\!\\!D$ problem for ${\mathcal{L}}$ is the eigenvalue
problem

${\mathcal{L}}u(x)=0,\,x\in\Omega;\partial_{\mathfrak{a}}u(x)=\lambda^{-1}\rho(x)u(x),\,x\in\Sigma;\,\partial_{\mathfrak{a}}=\sum_{j,k}a_{j,k}\nu_{k}\partial_{x_{j}},$
(1.3)

with a real function $\rho(x).$ The conormal differential operator
$\partial_{\mathfrak{a}}$ has symbol
$\imath\langle\mathfrak{a}(x)\xi,\nu_{x}\rangle.$

A direct setting of the problem (1.3) requires an explicit description of the
space of functions $u$ where the equation is considered. Under very weak
conditions imposed here on $\Omega$ and $\mathfrak{a},$ this task is very
hard. Instead of this, we follow [39], [2] in considering the $N\\!\\!\\!D$
problem in the variational form, where all spaces under consideration admit
explicit description, see Section 2. It is noted in [2] that such variational
approach follows the way of reasoning of H.Poincaré, V.A.Steklov and other
researchers of that time, for whom it is the variational setting of the
eigenvalue problem, based upon considering energy balance, was the primary
point of the eigenvalue analysis, while the differential Euler-Lagrange
equation appeared as a secondary object.

###  1.3. Asymptotic formulas. The main result

The asymptotic eigenvalue formula for the problem (1.2) can be written in the
following way (see, e.g., [2]). In a fixed orthogonal co-ordinate system, we
associate with the matrix $\mathfrak{a}(x)$ its sesquilinear form
$\boldsymbol{\mathfrak{a}}_{x}$:
$\boldsymbol{\mathfrak{a}}_{x}(\xi)=\sum_{j,k}a_{jk}(x)\xi_{j}\xi_{k};$ let
$\boldsymbol{\mathfrak{a}}_{x}(\xi,\eta)$ be the corresponding bilinear form.
For $x\in\Sigma,$ $\xi^{\prime}\in\mathrm{T}^{*}_{x}\Sigma$ and the normal
vector $\nu$ at $x,$ we denote by $\beta(x;\xi^{\prime})$ the positive square
root of
$\boldsymbol{\mathfrak{a}}_{x}(\nu)\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime})-\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime},\nu)^{2}$,

$\beta(x,\xi^{\prime})=(\boldsymbol{\mathfrak{a}}_{x}(\nu)\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime})-\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime},\nu)^{2})^{\frac{1}{2}};$
(1.4)

this is a function positively homogeneous in $\xi^{\prime}$ of order $1$. Then
the density $\alpha_{\pm}(x)$ is defined as

$\alpha_{\pm}(x)={\rm
vol}\,\\{\xi^{\prime}\in\mathrm{T}^{*}_{x}\Sigma:\beta(x;\xi^{\prime})<\rho_{\pm}(x)\\},$
(1.5)

and the asymptotic formulas to be proved are

${\mathbf{n}}_{\pm}(d,\ref{LSteklov})=(2\pi)^{-d}\int_{\Sigma}\alpha_{\pm}(x)d\mu_{\Sigma}(x),$
(1.6)

where $\mu_{\Sigma}$ is the surface measure on $\Sigma$ generated by the
embedding of $\Sigma$ into ${\mathbb{R}}^{d+1}.$

For a _continuous_ matrix-function $\mathfrak{a},$ such definition of the
eigenvalue problem (1.3) and the understanding of the formula (1.6) do not
cause confusion, even if the matrix $\mathfrak{a}$ is continuous only at the
boundary $\Sigma.$ In our case, the matrix $\mathfrak{a}(x)$ may be
discontinuous. However, under the conditions imposed below, the limit values
of $a_{jk}$ at the boundary are well defined almost everywhere in $\Sigma,$
and therefore these formulas (1.6) make sense.

In [2], the asymptotic formula (1.6) was justified under the following
conditions, see Theorem on P.242.

######  Theorem 1.1.

Let $\Omega\subset{\mathbb{R}}^{d+1}$ be a bounded domain; suppose that
$\Sigma=\partial\Omega$ is a surface of class $C^{1}$ and $\rho\in
L_{\infty}(\Sigma).$ Suppose finally that the coefficients $a_{jk}$ are
continuous in $\overline{\Omega}$. Then the asymptotics (1.6) holds. Moreover,
the same conclusion is correct if $\Sigma$ is Lipschitz but belongs to $C^{1}$
in a neighborhood of the support of $\rho,$ outside a closed set of zero
surface measure. The function $\rho$ may belong to $L_{d}(\Sigma)$ for $d>1$
and any $L_{q}(\Sigma)$ with $q>1,$ for $d=1.$

The _elementary_ formulation of the main result of our present paper is the
following.

######  Theorem 1.2.

Let $\Omega\subset{\mathbb{R}}^{d+1}$ be a bounded domain with Lipschitz
boundary $\Sigma.$ Suppose that the operator ${\mathcal{L}}_{\mathfrak{a}}$ is
uniformly elliptic in $\Omega$, the coefficients $a_{jk}$ are continuous in
$\overline{\Omega}$. Suppose that $\rho\in L_{d}(\Sigma)$ for $d>1$ and $\rho$
belongs to the Orlicz class $L\log L(\Sigma)$ for $d=1.$ Then the asymptotic
formula (1.6) is valid.

The continuity condition in Theorem 1.2 can be considerably relaxed. We say
that a function $a(x)\in L_{\infty}(\Omega)$ is _continuous at the boundary_
if there exists a continuous function $a_{b}\in C(\Sigma),$
$\Sigma=\partial\Omega$ such that for any $x^{\prime}\in\Sigma$,

$\lim_{\delta\to
0}\|a(.)-a_{b}(x^{\prime})\|_{L_{\infty}(B(x^{\prime},\delta)\cap\Omega)}=0,$
(1.7)

where $B(x^{\prime},\delta)$ is the ball in ${\mathbb{R}}^{d+1}$ with radius
$\delta$, centered at $x^{\prime}.$ Since the functions in
$L_{\infty}(\Omega)$ are defined up to a set of measure $0$, we may identify
in the above definition the continuous function $a_{b}$ with the restriction
of $a$ to $\Sigma$.

Having this definition, the formulation of the more general result that we are
going to prove here, is as follows.

######  Theorem 1.3.

Let $\Omega\subset{\mathbb{R}}^{d+1}$ be a bounded domain with Lipschitz
boundary $\Sigma.$ Suppose that the operator ${\mathcal{L}}_{\mathfrak{a}}$ is
uniformly elliptic in $\Omega$, the coefficients $a_{jk}$ are bounded and
measurable and they are continuous at the boundary in the sense of (1.7). Then
the asymptotic formula (1.6) is valid.

######  Remark.

The conditions imposed on the coefficients of the operator ${\mathcal{L}}$ can
be further relaxed, allowing them to be unbounded in a certain sense and
admitting certain degeneration of ellipticity. We do not pursue this line
here, in order to keep the elementary level of the presentation.

We describe here briefly the structure of the proof. After certain
transformations, the weighted $N\\!\\!\\!D$ problem in the domain with
Lipschitz boundary is relocated to a similar problem for a uniformly elliptic
operator ${\mathcal{L}}$ but with bounded coefficients, possibly,
discontinuous at the boundary, in a domain with smooth boundary. This operator
is approximated by operators ${\tilde{\mathcal{L}}}$ with smooth coefficients
which converge to the coefficients of ${\mathcal{L}}$ in a certain $L_{p}$ \-
based norm, with some $p<\infty$. Then it is proved that the compact operators
describing the P-S spectrum for ${\tilde{\mathcal{L}}}$ converge so strongly
that it becomes possible to pass to the limit in the formulas for coefficients
${\mathbf{n}}(d,.)$ in the Steklov eigenvalue asymptotics for
${\tilde{\mathcal{L}}}.$

##  2\. The variational representation

The setting of the eigenvalue problem under very weak regularity conditions
imposed on the coefficients and the domain, is variational. As usual, for more
regular data, the variational setting is equivalent to the classical, ’strong’
one.

###  2.1. The weighted $N\\!\\!\\!D$ operator

Since the earliest paper on the topic by Lenart Sandgren [39] in 1955, the
most commonly used method of the study of Steklov eigenvalues is the
variational one. Its main advantage is its robustness with respect to various
perturbations of the problem. We too will need a version of the variational
setting of the Steklov eigenvalue problem (see, however, Appendix where we
discuss a different approach).

First of all, since we are going to use the instruments of the perturbation
theory for compact operators, we will need to pass from the unbounded
Dirichlet-to-Neumann operators $D\\!\\!\\!N$ (which are often considered) to
the N-to-D ones, $N\\!\\!\\!D,$ since the latter operators are compact (this
corresponds to the placement of the spectral parameter in (1.3)). A minor
inconvenience here is that the Neumann problem for the Laplacian (as well as
for other second order elliptic operators without zero order terms) has an
eigenfunction with zero eigenvalue, so the Neumann operator is not invertible.
There are several standard ways to deal with this circumstance. Say, often,
the Neumann operator is restricted to the codimension 1 subspace of functions
orthogonal to constants. In our considerations, it is more convenient to avoid
hitting the zero eigenvalue by means of adding a nontrivial nonnegative
function ${\mathbf{v}}_{0}(x)$ to the operator ${\mathcal{L}}$, namely to
consider the operator ${\mathcal{L}}+{\mathbf{v}}_{0}(x)$. Such perturbation
is relatively compact with respect to self-adjoint realizations of
${\mathcal{L}}$, and, by the usual perturbation reasoning, see, e.g., Lemma
1.3 in [13], it does not influence the asymptotics of the spectrum. One might
have set ${\mathbf{v}}_{0}\equiv 1,$ but this does not simplify the matter.
For economy of symbols, we use the notation ${\mathcal{L}}$ for the operator
with ${\mathbf{v}}_{0}$ already absorbed, further on. At a proper moment, we
will use the freedom in choosing ${\mathbf{v}}_{0}.$

We denote by ${\mathbf{T}}\equiv{\mathbf{T}}_{\mathfrak{a}}$ the operator
${\mathcal{L}}$ with _Neumann_ boundary conditions on $\Sigma$. Namely, it
corresponds to the quadratic form

${\mathbf{a}}_{0}[u]=\int_{\Omega}(\boldsymbol{\mathfrak{a}}[u](x)+{\mathbf{v}}_{0}(x)|u(x)|^{2})dx,\,\boldsymbol{\mathfrak{a}}[u](x):=\sum_{j,k}a_{j,k}(x)\partial_{j}u(x)\partial_{k}\bar{u}(x),$
(2.1)

with domain $u\in H^{1}(\Omega);$ this form is equivalent to the square of the
standard norm $\|u\|_{H^{1}(\Omega)}$ in the Sobolev space $H^{1}(\Omega).$
Being considered in $L_{2}(\Omega)$, this quadratic form defines the positive
self-adjoint operator ${\mathbf{T}}$ satisfying

${\mathbf{a}}_{0}[u]=({\mathbf{T}}^{\frac{1}{2}}u,{\mathbf{T}}^{\frac{1}{2}}u),\,u\in
H^{1}(\Omega)={\mathcal{D}}({\mathbf{T}}^{\frac{1}{2}}).$

(This fact, which we will use systematically, that the domain of a closed
positive quadratic form coincides with the domain of the square root of the
operator defined by this quadratic form, is a standard property in the
abstract spectral theory, see, e.g., [27], Theorem 2.23 in Sect. VI. For
quadratic forms which are only sectorial (the case we do not need), this
equality constitutes the famous _Kato’s square root problem_ which has been
the topic of deep extensive studies since 1980-s.)

Further on, we will sometimes use the notion of the spectrum of the ratio of
quadratic forms in the Hilbert space. Actually, the spectrum of the ratio
$\frac{\mathfrak{B}[f]}{\mathfrak{A}[f]}$ with $f\in\mathfrak{X}$ is a
shorthand for the spectrum of the operator defined by the quadratic form
$\mathfrak{B}[f]$ in the Hilbert space $\mathfrak{X}$ equipped with the norm
defined by the quadratic form $\mathfrak{A}[f]$.

For a real function $\rho_{0}(x),x\in\Sigma$ we consider the quadratic form

$\boldsymbol{\boldsymbol{\rho}}_{0}[f]=\int_{\Sigma}\rho_{0}(x)|f(x)|^{2}d\mu_{\Sigma}(x).$
(2.2)

Let $\rho_{0}$ belong to $L_{d}(\Sigma).$ By Kondrashov’s trace theorem, for
$d>1,$ the trace mapping $u\mapsto u_{|_{\Sigma}},$ defined initially on
continuous functions $u\in H^{1}(\Omega)$, extends by continuity to the
bounded mapping $\gamma:u\mapsto u_{|_{\Sigma}},$ $\gamma:H^{1}(\Omega)\to
L_{q}(\Sigma),$ $q=\frac{2d}{d-1}.$ Therefore, by the Hölder inequality, if
$\rho_{0}\in L_{d}(\Sigma)$, with $d^{-1}+2q^{-1}=1$, the quadratic form (2.2)
satisfies

$\boldsymbol{\boldsymbol{\rho}}_{0}[\gamma u]\leq
C\|\rho_{0}\|_{L_{d}(\Sigma)}{\mathbf{a}}_{0}[u],\,u\in H^{1}(\Omega).$ (2.3)

For $d=1,$ a somewhat more complicated reasoning (see, e.g., [1]) gives

$\boldsymbol{\boldsymbol{\rho}}_{0}[\gamma u]\leq C\|\rho_{0}\|_{L\log
L(\Sigma)}{\mathbf{a}}_{0}[u],\,u\in H^{1}(\Omega),$ (2.4)

where ${L\log L(\Sigma)}$ is the Orlicz space consisting of functions $u$ on
$\Sigma$ for which $|u|(1+|\log|u||)\in L_{1}(\Sigma),$ with the Luxemburg
norm.

The following variational representation is known since [39], see e.g., [2],
[43], and sources cited there: the Steklov spectrum is described by the ratio
of quadratic forms $\frac{\boldsymbol{\boldsymbol{\rho}}_{0}[\gamma
u]}{{\mathbf{a}}_{0}[u]}$ considered on the space of $L_{2}-$solutions of the
equation ${\mathcal{L}}u=0.$ Since $\boldsymbol{\boldsymbol{\rho}}_{0}[\gamma
u]=0$ for $u\in\overset{\circ}{H}{}^{1}(\Omega),$ one can, as this is usually
done, drop the condition ${\mathcal{L}}u=0$ in the variational setting and
consider the spectrum of the ratio

$\frac{\boldsymbol{\boldsymbol{\rho}}_{0}[\gamma
u]}{{\mathbf{a}}_{0}[u]},\,u\in H^{1}(\Omega).$ (2.5)

###  2.2. A freedom in the choice of the weight function $\rho_{0}$.

Estimates (2.3), (2.4) mean that the quadratic form (2.2) is bounded in
$H^{1}(\Omega)$, the equality in (2.2) can be therefore extended by continuity
to all traces of functions in $H^{1}(\Omega)$ and therefore the form
$\boldsymbol{\boldsymbol{\rho}}_{0}$ defines a bounded operator
${\mathbf{S}}={\mathbf{S}}[{\rho_{0}}]$ in $H^{1}(\Omega)$, with norm
estimated by $\|\rho_{0}\|_{L_{d}(\Sigma)},$ resp., $\|\rho_{0}\|_{L\log
L(\Sigma)}.$ If the function $\rho_{0}$ is real-valued, this operator is self-
adjoint. In fact, this operator is compact and its eigenvalues satisfy a power
order estimate.

######  Theorem 2.1.

Let $\Omega\subset{\mathbb{R}}^{d+1},$ $d>1,$ be a bounded domain with
Lipschitz boundary $\Sigma$ and let $\rho_{0}\in L_{d}(\Sigma)$. Then for the
singular numbers of ${\mathbf{S}}={\mathbf{S}}{[\rho_{0}]}$ the estimate holds

$n(\lambda,{\mathbf{S}}{[\rho_{0}]})\leq
C(\Omega,\mathfrak{a})\lambda^{-d}\|\rho_{0}\|_{L_{d}(\Sigma)}^{d}.$ (2.6)

If the function $\rho_{0}$ is real-valued, estimates of the form (2.6) hold,
separately, for positive and negative eigenvalues,

$n_{\pm}(\lambda,{\mathbf{S}}{[\rho_{0}]})\leq
C(\Omega,\mathfrak{a})\lambda^{-d}\|(\rho_{0})_{\pm}\|_{L_{d}(\Sigma)}^{d}.$
(2.7)

For $d=1$ the above results hold, with $L_{d}(\Sigma)$ replaced by $L\log
L(\Sigma).$

###### Proof.

For $d>1,$ we use the estimate for the eigenvalues of Birman-Schwinger type
operators with singular measures, obtained recently in [38]. We consider some
bounded open set $\Omega^{\prime}$ containing $\Omega$ strictly inside and
apply Theorem 3.3 in [38], for the particular choice of
${\mathbf{N}}=d+1,l=1,s=d,d>1,$ $\mu$ being the surface measure on the surface
$\Sigma$. This theorem states the following. Let
$\Omega^{\prime}\subset{\mathbb{R}}^{{\mathbf{N}}}$ be a bounded open set,
$\mu$ be a compactly supported measure in $\Omega^{\prime}$ satisfying, for
$X\in\Sigma,$ $r\leq\operatorname{diam\,}(\Sigma)$, $\Sigma=\hbox{{\rm
supp}}\,\mu$,

$\mu(B(X,r))\leq{\mathcal{A}}r^{s},\,s>{\mathbf{N}}-2(=d-1),\rho\in
L_{\theta,\mu},\,\theta=\frac{s}{s+2-{\mathbf{N}}},$ (2.8)

where $B(X,r)$ is the ball with radius $r$ centered at $X.$ Then, for the
particular case $s=d,$ $\theta=d,$ for the operator
${\mathbf{S}}={\mathbf{S}}[\rho_{0},\Omega^{\prime}]$ in
$\overset{\circ}{H}{}^{1}(\Omega^{\prime})$ defined by the quadratic form
$\boldsymbol{\boldsymbol{\rho}}_{0}$ in (2.2), the estimates (2.6), (2.7)
hold.

We use this result in the following way. First, note that for a compact
Lipschitz surface $\Sigma$ of dimension $d$ with $\mu$ being the
$d$-dimensional Hausdorff measure on $\Sigma$, the condition (2.8) is
satisfied with $s=d$. Next, by the Calderón-Stein extension theorem, (see,
e.g., Theorem 5.24 in [1]), there exists a bounded extension operator
${\mathcal{E}}:H^{1}(\Omega)\to\overset{\circ}{H}{}^{1}(\Omega^{\prime}),$

$\|{\mathcal{E}}v\|_{\overset{\circ}{H}{}^{1}(\Omega^{\prime})}\leq
C(\Omega,\Omega^{\prime})\|v\|_{H^{1}(\Omega)}.$

Therefore, if on some subspace $\sc\mbox{L}\hskip 1.0pt\subset H^{1}(\Omega)$
of dimension $n,$ we have

$\frac{\boldsymbol{\boldsymbol{\rho}}_{0}[\gamma
v]}{\|v\|^{2}_{H^{1}(\Omega)}}\geq\lambda,\,v\in\sc\mbox{L}\hskip 1.0pt,$

it follows that

$\frac{\boldsymbol{\boldsymbol{\rho}}_{0}[\gamma{\mathcal{E}}v]}{\|{\mathcal{E}}v\|^{2}_{\overset{\circ}{H}{}^{1}(\Omega^{\prime})}}=\frac{\boldsymbol{\boldsymbol{\rho}}_{0}[\gamma
v]}{\|{\mathcal{E}}v\|^{2}_{\overset{\circ}{H}{}^{1}(\Omega^{\prime})}}\geq
C(\Omega,\Omega^{\prime})^{-2}\lambda.$

By the variational principle, this implies the inequality for the counting
functions of eigenvalues of operators,

$n_{+}(\lambda,{\mathbf{S}}[\rho_{0},\Omega])\leq
n_{+}(C(\Omega,\Omega^{\prime})^{-2}\lambda,{\mathbf{S}}[\rho_{0},\Omega^{\prime}]),$

which gives (2.7) with ’+’ sign. Other estimates follow in a similar way.
For $d=1,$ the reasoning is the same, just we use the Orlicz estimate in [36]
instead of the eigenvalue estimate in [38]. ∎

######  Remark.

In fact, some stronger results which, however, are not needed for our present
topic, hold. Actually, the extension theorem for the class $H^{1}(\Omega)$,
which we used, is valid under somewhat weaker restrictions than the Lipschitz
property. Namely, in [25], the class of $(\epsilon-\delta)$ \- domains has
been introduced (see the definition in [25], p.73) called also _locally
uniform domains_ in the literature. This class contains all Lipschitz domains,
but some other, less regular ones, as well (it is mentioned in [25] that an
$(\epsilon-\delta)$-domain in ${\mathbb{R}}^{d+1}$ may have a fractal boundary
of any Haussdorff dimension in $[d,d+1)$, particular examples being the Von
Koch snowflakes. On the other hand, an exterior cusp prevents the domain from
being locally uniform). For an $(\epsilon-\delta)$ domain, a bounded extension
operator ${\mathcal{E}}$ for Sobolev spaces exists. Therefore, if the
$d$-dimensional Hausdorff measure on $\partial\Omega$ satisfies (2.8) with
$s=d$ and $\Omega$ is an $(\epsilon-\delta)$-domain then for the operator
${\mathbf{S}}[\rho_{0},\Omega]$ with a function $\rho_{0}\in L_{d}(\mu)$,
which can be considered as a natural generalization of the $N\\!\\!\\!D$
operator, the estimate (2.6), (2.7) is valid.

The operator ${\mathbf{S}}\equiv{\mathbf{S}}{[\rho_{0}]}$ has the same
eigenvalues as the operator of the P-S problem in the smooth case. In the
course of our reasoning, we will introduce some more operators with the same
spectrum.

The eigenvalue estimates in Theorem 2.1, with coefficient in front of the
power term depending on the integral norm of the functional parameter, enable
one to establish the property of the passage to limit in asymptotic
eigenvalues formulas. This way of reasoning, starting in the classical studies
by M.Sh. Birman and M.Z. Solomyak in 1960-s -1970-s (see, e.g., [13], Lemma
1.5), is now the standard tool in establishing the eigenvalue asymptotics for
singular problems by means of approximating by more regular ones. We reproduce
here this extraordinary lemma, which we use at least on three occasions in the
course of the paper.

######  Lemma 2.2.

Let ${\mathbf{K}}$ be a compact operator, $\theta>0,$ and for sufficiently
small $\epsilon$ one can split ${\mathbf{K}}$ into the sum,
${\mathbf{K}}={\mathbf{K}}_{\epsilon}+{\mathbf{K}}_{\epsilon}^{\prime}$ so
that for ${\mathbf{K}}_{\epsilon}$ the singular numbers asymptotic formula
holds,

${\mathbf{n}}(\theta,{\mathbf{K}}_{\epsilon})=A_{\epsilon},$ (2.9)

and the operators ${\mathbf{K}}_{\epsilon}^{\prime}$ are asymptotically small
in the sense

${\mathbf{n}}^{\sup}(\theta,{\mathbf{K}}_{\epsilon}^{\prime})\to 0,$ (2.10)

as $\epsilon\to 0$. Then the limit $\lim_{\epsilon\to 0}A_{\epsilon}=A$ exists
and the asymptotics holds,

${\mathbf{n}}(\theta,{\mathbf{K}})=A.$ (2.11)

An analogous statement is valid for the positive/negative eigenvalues of a
_self-adjoint_ operator ${\mathbf{K}}$. One should replace in the above
formulation ${\mathbf{n}}$ by ${\mathbf{n}}_{\pm}$ in (2.9) and (2.11) (but
not in (2.10)).

Among most recent examples of applying this way of reasoning, one can cite
[38], Sect.6 and [26]. We formulate here the particular statement concerning
the eigenvalue asymptotics for the operator ${\mathbf{S}}[{\rho_{0}}];$ it
follows from estimates (2.6), (2.7)

######  Corollary 2.3.

Suppose that for the weights $\rho_{0}$ in some set ${\mathcal{Y}}\subset
L_{d}(\Sigma)$, $d>1,$ the asymptotic formula (1.6) is established. Then this
formula is valid for all $\rho_{0}\in\overline{{\mathcal{Y}}},$ the closure of
${\mathcal{Y}}$ in the norm of $L_{d}(\Sigma).$ The same statement is valid
for the dimension $d=1,$ with $L_{d}(\Sigma)$ replaced by the Orlicz space
$L\log L(\Sigma).$

In particular, this means that we can restrict ourselves to considering as
${\mathcal{Y}}$, the set of continuous, and, later, smooth functions
$\rho_{0}$, which are dense in the corresponding space with integral norm.

We follow now [39], [41], [46], [3], [37] and other papers where the
variational method was used for the study of spectral problems containing a
weight, including the P-S type spectral problems. Namely, once order sharp
eigenvalue _estimates_ involving an integral norm of the weight function are
obtained, it is possible to restrict the further study of eigenvalue
_asymptotics_ to considering _sign-definite_ weight functions only. We do not
need to repeat this standard reasoning in detail, but just describe briefly
the common scheme. Starting with a non-sign-definite weight function
$\rho_{0}$, we first approximate it in a proper integral norm by a function
$\rho$ which has ’separated’ positive and negative parts, ${\rm
dist}\,(\hbox{{\rm supp}}\,\rho_{+},\hbox{{\rm supp}}\,\rho_{-})>0.$ After
this, the study of the distribution of the positive, resp., negative spectrum
of the problem is performed by cutting-away the part of the weight with wrong
sign.

Consequently, we restrict ourselves to non-negative weight functions
$\rho_{0}$ further on.

###  2.3. The asymptotic coefficient

The next preparatory step consists in obtaining a convenient expression for
the coefficient in the asymptotic formula for the eigenvalues of the P-S
problem so that it is possible to pass to the limit in this expression as the
coefficients of the operator converge in a proper sense.

The asymptotic coefficient, known to be correct in the smooth case and aimed
for in the Lipschitz case, is given by (1.6), (1.4), (1.5).

For our further needs, it will be more convenient to use a somewhat different
representation of the coefficient in (1.6). Namely, for a given $x\in\Sigma,$
we can represent $\boldsymbol{\mathfrak{a}}_{x}(\xi)$,
$\boldsymbol{\mathfrak{a}}_{x}(\xi,\eta)$ as matrix products

$\boldsymbol{\mathfrak{a}}_{x}(\xi)=\xi^{*}\mathfrak{a}(x)\xi,\,\,\boldsymbol{\mathfrak{a}}_{x}(\xi,\eta)=\xi^{*}\mathfrak{a}(x)\eta=\eta^{*}\mathfrak{a}(x)\xi,$

(where $\xi^{*}$ denotes the row-vector, matrix-adjoint to the column-vector
$\xi$.) In these notations, we have

$\displaystyle\beta(x,\xi^{\prime})^{2}=\nu^{*}\mathfrak{a}(x)\nu{\xi^{\prime}}^{*}\mathfrak{a}(x)\xi^{\prime}-{\xi^{\prime}}^{*}\mathfrak{a}(x)\nu\nu^{*}\mathfrak{a}(x)\xi^{\prime}$
(2.12)
$\displaystyle=\nu^{*}\mathfrak{a}(x)\nu\langle\mathfrak{a}(x)\xi^{\prime},\xi^{\prime}\rangle-\langle[(\nu^{*}\mathfrak{a}(x))^{*}(\nu^{*}\mathfrak{a}(x))]\xi^{\prime},\xi^{\prime}\rangle:=\langle\Theta(x)\xi^{\prime},\xi^{\prime}\rangle$

where
$\Theta(x)=(\nu^{*}\mathfrak{a}(x)\nu)\mathfrak{a}(x)-\mathfrak{a}(x)\nu\nu^{*}\mathfrak{a}(x)$
(note that the last matrix here has rank 1, together with $\nu\nu^{*}$.). Now
let $\Pi_{x}$ be the embedding of $\mathrm{T}^{*}_{x}\Sigma$ into
${\mathbb{R}}^{d+1},$ and $\Pi_{x}^{*}$ be the adjoint operator, the
projection of ${\mathbb{R}}^{d+1}$ to $\mathrm{T}^{*}_{x}\Sigma.$ Then (2.12)
can be written as

$\beta(x,\xi^{\prime},\nu)^{2}=\langle{\Pi_{x}^{*}\Theta_{x}\Pi_{x}}\xi^{\prime},\xi^{\prime}\rangle,$
(2.13)

therefore,

$\alpha_{\pm}(x)={\rm
vol}\,\\{\xi^{\prime}\in{\mathbb{R}}^{d}:0<\langle\Theta^{\prime}_{x})\xi^{\prime},\xi^{\prime}\rangle^{\frac{1}{2}}<(\rho_{0})_{\pm}(x)\\}=\omega_{d}(\rho_{0})_{\pm}(x)^{d}\det(\Theta^{\prime}_{x})^{-1/2},$
(2.14)

with $d\times d$ matrix $\Theta^{\prime}_{x}=\Pi_{x}^{*}\Theta_{x}\Pi_{x}$ The
meaning of the matrix $\Theta^{\prime}_{x}$ can be understood in the following
way. We represent the $(d+1)\times(d+1)$ matrix $\Theta$ in an orthogonal
frame where one of basis vectors is the normal vector $\nu_{x}.$ Next, we
strike out in this matrix the row and the column corresponding to $\nu_{x}$.
What remains is just our matrix $\Theta^{\prime}_{x}$. Using (2.14), we can
establish certain estimates for the quantity
${\mathbf{n}}_{\pm}(d,\ref{LSteklov}),$ uniform in a class of elliptic
operators.

It follows from the Cauchy-Schwartz inequality that

$\boldsymbol{\mathfrak{a}}_{x}(\nu)\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime})-\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime},\nu)^{2}=\langle\mathfrak{a}(x)\nu,\nu\rangle\langle\mathfrak{a}(x)\xi^{\prime},\xi^{\prime}\rangle-\langle\mathfrak{a}(x)\xi^{\prime},\nu\rangle^{2}\geq
0.$ (2.15)

For a positive matrix $\mathfrak{a},$ the exact equality in the Cauchy-
Schwartz inequality (2.15) is possible only if the vectors $\xi^{\prime}$ and
$\nu$ are parallel. The latter may never happen, therefore the inequality in
(2.15) is strict; by homogeneity we obtain

$\boldsymbol{\mathfrak{a}}_{x}(\nu)\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime})-\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime},\nu)^{2}\geq{\mathbf{k}}(\xi^{\prime}/|\xi^{\prime}|)\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime}),$
(2.16)

with some ${\mathbf{k}}(\xi^{\prime}/|\xi^{\prime}|)>0.$

The sharp constant ${\mathbf{k}}(\xi^{\prime}/|\xi^{\prime}|)$ in (2.16)
depends continuously on $\xi^{\prime}\neq 0;$ by the compactness of the unit
sphere, we have

$\beta_{x}(\xi^{\prime})^{2}=\boldsymbol{\mathfrak{a}}_{x}(\nu)\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime})-\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime},\nu)^{2}\geq{\mathbf{k}}(\xi^{\prime}/|\xi^{\prime}|)\boldsymbol{\mathfrak{a}}_{x}(\xi^{\prime})\geq
C\mathfrak{k}(x)|\xi^{\prime}|^{2},$ (2.17)

where $\mathfrak{k}(x)$ is the ellipticity constant of the matrix
$\mathfrak{a}(x)$ at the point $x\in\Sigma.$ We suppose that the matrix
$\mathfrak{a}(x)$ is uniformly elliptic,
$\mathfrak{a}(x)\geq{\mathbf{k}}_{0}>0.$ Then (2.17) implies

$\alpha_{\pm}(x)\leq C{\mathbf{k}}_{0}^{-d}\omega_{d}\rho_{\pm}(x)^{d}.$
(2.18)

In the equivalent representation, the estimate for the matrix
$\Theta^{\prime}_{x}$ also follows, $\det(\Theta^{\prime}_{x})\geq
C{\mathbf{k}}_{0}^{-d}.$

The last estimates enable us to establish the following convergence property.

######  Lemma 2.4.

Let $\mathfrak{a}_{s}(x)$, $s\in[0,1)$ be a family of Hermitian matrix
functions on $\Sigma$, such that they satisfy ellipticity estimates
$\mathfrak{a}_{s}(x)\geq{\mathbf{k}}_{0}$ uniformly in $x,s,$ and are also
uniformly bounded, $|\mathfrak{a}_{s}(x)|\leq{\mathbf{k}}_{1}$. Suppose that
the matrices $\mathfrak{a}_{s}$ converges as $s\to 0$ to $\mathfrak{a}_{0}$ in
$L_{q}(\Sigma),$ $q<\infty.$ Denote by $\varrho_{s}(x)$ the function

$\varrho_{s}(x)=\det(\Theta_{s}^{\prime}(x))^{-\frac{1}{2}},$

where $\Theta_{s}^{\prime}$ is the matrix $\Theta^{\prime}$ in (2.12),
corresponding to the coefficient matrix $\mathfrak{a}_{s}$. Then

$\varrho_{s}(x)\to\varrho_{0}(x)$

in $L_{q}(\Sigma).$

###### Proof.

The statement follows easily from the fact that $\varrho_{s}(x)^{2}$ is a
rational function of the entries of the matrix $\mathfrak{a}_{s}(x)$ and these
matrices are separated from zero. ∎

The following result shows that the expression for the coefficient in the
asymptotics of P-S eigenvalues endures the convergence of the coefficient
matrix in the integral metric.

######  Lemma 2.5.

Let $\mathfrak{a}_{s}(x)$ be a family of matrices, as in Lemma 2.4 and
$\rho_{0}$ be a bounded function on $\Sigma.$ Then
$\int\varrho_{s,\pm}(x)\rho_{0}(x)d\Sigma\to\int\varrho_{\pm}(x)\rho_{0}(x)d\Sigma.$

The statement follows by the passage to limit, using Lemma 2.4.

###  2.4. D-to-N and N-to-D operators as pseudodifferential ones

It is well known that for the Laplacian, the $D\\!\\!\\!N$ operator equals, up
to lower order terms, $(-\Delta_{\Sigma})^{\frac{1}{2}}$, for a smooth
boundary; an interesting discussion can be found in [19]. For general elliptic
operators, in [4], [5], M.S.Agranovich presented an explanation of the
coefficient $\beta(x,\xi^{\prime})$ for the case of a smooth (or ’almost
smooth’) surface $\Sigma.$ We discuss now the latter formula. Suppose first
that the coefficients matrix $\mathfrak{a}(x)$ does not depend on
$x\in\Omega.$ Consider the fundamental solution $R(x-y)$ for the operator
${\mathcal{L}},$ so that ${\mathcal{L}}_{x}R(x-y)=\delta(x-y).$ With this
fundamental solution the classical potential operators are associated, namely,
the single layer potential operator S ,

$f\mapsto\sc\mbox{S}\hskip
1.0pt_{f}(x)=\int_{\Sigma}R(x-y)f(y)d\mu(y),\,x\in\Sigma,$

and the ’direct value’ of the conormal derivative of the single layer
potential (often called the Neumann-Poincaré operator)

$\sc\mbox{D}\hskip 1.0pt^{\prime}:\,f\mapsto\sc\mbox{D}\hskip
1.0pt^{\prime}_{f}(x)=\int_{\Sigma}\partial_{\mathfrak{a}(x)}R(x-y)f(y)d\mu(y),\,x\in\Sigma,$

the latter integral understood in the principal value sense. The adjoint to
$\sc\mbox{D}\hskip 1.0pt^{\prime}_{f}$ is the integral operator
$\sc\mbox{D}\hskip 1.0pt,$ called the direct value of the double layer
potential,

$\sc\mbox{D}\hskip 1.0pt:\,f\mapsto\sc\mbox{D}\hskip
1.0pt_{f}(x)=\int_{\Sigma}\partial_{\mathfrak{a}(y)}R(x-y)f(y)d\mu(y),\,x\in\Sigma.$

Then the $N\\!\\!\\!D$ operator is, up to lower order terms, the composition

$N\\!\\!\\!D=(\frac{1}{2}+\sc\mbox{D}\hskip
1.0pt^{\prime})^{-1}\sc\mbox{S}\hskip 1.0pt=\sc\mbox{S}\hskip
1.0pt(\frac{1}{2}+\sc\mbox{D}\hskip 1.0pt)^{-1}.$ (2.19)

For a smooth boundary $\Sigma$, the operators S and D are order $-1$
pseudodifferential operators on $\Sigma.$ Since S is a restriction of the
order $-2$ pseudodifferential operator to the boundary, its symbol is
calculated according to the usual rules of pseudodifferential calculus,

$\mathfrak{s}(x,\xi^{\prime})=\frac{1}{2\pi}\int_{-\infty}^{\infty}[\mathfrak{a}(\xi^{\prime},\xi_{d+1})]^{-1}d\xi_{d+1},\,x\in\Sigma,$
(2.20)

in co-ordinates where $x_{d+1}$ is directed along the normal vector $\nu_{x}$
to the tangent plane at the point $x\in\Sigma$ and the corresponding co-
ordinates $\xi^{\prime},\xi_{d+1},$ with $\xi^{\prime}$ in the cotangent plane
and $\xi_{d+1}$ directed along $\nu_{x}$, see, e.g., [4], (4.28). Calculations
by $\eqref{D2N rest}$ give the expression for
$\beta(x,\xi^{\prime})=\mathfrak{s}(x,\xi^{\prime})^{-1}$ in (1.4) as the
principal symbol of the operator $D\\!\\!\\!N=N\\!\\!\\!D^{-1}.$ Similar
considerations cover the case of variable smooth coefficients, where the
fundamental solution takes the form $R(x,y),$
${\mathcal{L}}_{x}R(x,y)=\delta(x-y)$ and in (2.20) one should change
$\mathfrak{a}(\xi^{\prime},\xi_{d+1})^{-1}$ to
$\mathfrak{a}(x^{\prime},0,\xi^{\prime},\xi_{d+1})^{-1},$ see [3], Sect.3.4.
For a non-smooth, Lipschitz surface, the function $\beta(x,\xi^{\prime})$
still can be calculated by the above formula, but it is not a symbol of a nice
pseudodifferential operator any more.

The principal symbol of the double layer operator D can be also calculated,
see, e.g., [44], but the we do not need an explicit expression here.

##  3\. Some geometry considerations

In this Section we introduce a global co-ordinate system in an interior collar
neighborhood of the boundary $\Sigma=\partial\Omega.$ After the passage to
this co-ordinate system and a special co-ordinates change, the Poincaré-
Steklov problem in a Lipschitz domain, for the equation with coefficients,
_continuous at the boundary_ , reduces to an equivalent problem in a domain
with smooth boundary, but for an elliptic operator with _discontinuous_
bounded measurable coefficients, fortunately, still with a certain very weak
continuity property.

###  3.1. A smooth surface approximating the boundary

By the Rademacher theorem, a Lipschitz function has derivatives almost
everywhere and these derivatives are bounded. Therefore, the Lipschitz surface
$\Sigma\subset{\mathbb{R}}^{d+1}$ has a tangent plane almost everywhere. This
fact enables one to describe explicitly the surface measure on $\Sigma$
generated by the Lebesgue measure on ${\mathbb{R}}^{d+1}.$ Namely, on the
local Lipschitz graph
$x_{d+1}=\psi(x^{\prime}),\,x^{\prime}\in\sc\mbox{V}\hskip
1.0pt\subset{\mathbb{R}}^{d},$ the surface measure is given by

$\mu_{\Sigma}(E)=\int_{\pi
E}\sqrt{1+|\nabla\psi(x^{\prime})|^{2}}dx^{\prime},$ (3.1)

where $\pi E$ is the projection of the set $E\subset\Sigma$ to
${\mathbb{R}}^{d}.$ It is important to note that for a Lipschitz function
$\psi,$ the integrand in (3.1) is a bounded function, an algebraic expression
of the partial derivatives of $\psi.$

A detailed study of geometric properties of Lipschitz domains is presented in
[23]. In particular, for a domain $\Omega\subset{\mathbb{R}}^{d+1}$ with
Lipschitz boundary, the distributional gradient ${\rm grad}\,(\chi_{\Omega})$
of the characteristic function of $\Omega$ is a vector measure which is
absolutely continuous with respect to the measure $\mu_{\Sigma}$ in (3.1),

${\rm grad}\,(\chi_{\Omega})=-{\nu}(x)\mu_{\Sigma}$ (3.2)

with the vector-valued density ${\nu}$, $|{\nu}(x)|=1$ coinciding with the
normal at $x$ almost everywhere on $\Sigma.$ (Such vector measure can be well
defined for a class of domains even somewhat less regular as well, see [23].)

Starting from [39], and then in [2], [41], [3], etc., a localization was used,
reducing the spectral analysis for operators related with the P-S problem in
$\Omega$ to the ones in Lipschitz cylinders. In our study, it is more
convenient to define the surface $\Sigma$ and perform related constructions
globally.

To do this, we will need the classical result on the existence of a smooth
vector field transversal to the Lipschitz surface. Namely, there exists a
$C^{\infty}$-vector field $\boldsymbol{\Gamma}(x),$ $x\in{\mathbb{R}}^{d+1}$
such that at the points $x\in\Sigma=\partial\Omega,$ those points where the
tangent plane exists, $|\boldsymbol{\Gamma}(x)|=1$ and
$\boldsymbol{\Gamma}(x)$ forms an acute angle with the normal $\nu(x)$ to
$\Sigma$ at $x,$ not exceeding some constant which is strictly smaller than
$\pi/2.$ The construction of such field can be found, e.g., in [23],
Proposition 2.3 (in fact a similar construction was used as long ago as in
1983 by A.P. Calderón, see [16]).

For further reasoning, we need one more geometric construction. We suspect
that it might have been well known since long ago, but were not able to locate
a reference in the literature. Therefore, we describe it here.

######  Lemma 3.1.

Let $\Omega\subset{\mathbb{R}}^{d+1}$ be a bounded domain with Lipschitz
boundary $\Sigma.$ Let $\boldsymbol{\Gamma}$ be the transversal smooth vector
field, as above. In a neighborhood of $\Sigma$ denote by $\ell_{x}$ the
integral curve of this vector field passing through $x$. Then, for a certain
interior collar neighborhood U of $\Sigma$ there exists a smooth surface
$\boldsymbol{\boldsymbol{\Sigma}}\subset\sc\mbox{U}\hskip 1.0pt$, which has
exactly one intersection point with every integral curve $\ell$ of the field
$\boldsymbol{\Gamma}$ in $\sc\mbox{U}\hskip 1.0pt,$ and this intersection is
transversal.

###### Proof.

In some two-sided collar neighborhood $\widetilde{\sc\mbox{U}\hskip 1.0pt}$ of
$\Sigma,$ the integral curves $\ell_{x^{\prime}},$ $x^{\prime}\in\Sigma,$ of
the vector field $\boldsymbol{\Gamma}$ form a one-dimensional foliation. We
define the function $\mathfrak{f}(x),$ $x\in\widetilde{\sc\mbox{U}\hskip
1.0pt},$ as the distance of the point $x$ to the intersection point
$x^{\prime}\in\Sigma$ along the integral curve $\ell_{x^{\prime}}$ passing
through $x$ (with proper sign, positive inside $\Sigma.$) Then
$x=(x^{\prime},t)$, $t=\mathfrak{f}(x)$ form a new co-ordinate system in
$\widetilde{\sc\mbox{U}\hskip 1.0pt}.$ In these co-ordinates the function
$\mathfrak{f}$ is continuous in $\widetilde{\sc\mbox{U}\hskip 1.0pt}$,
differentiable almost everywhere and has derivatives of any order in $t$
variable, in other words, along $\ell_{x^{\prime}}$. Moreover, $\mathfrak{f}$
is strictly decaying in $t$ variable, i.e., along the curves
$\ell_{x^{\prime}}$.

We take a smooth non-negative cut-off function
$\chi_{\epsilon}(x),\,\operatorname{diam\,}\hbox{{\rm
supp}}\,\chi_{\epsilon}<\epsilon,$ $\int\chi_{\epsilon}(x)dx=1,$ with
sufficiently small $\epsilon$, and consider the convolution, the smooth
function $\mathfrak{f}_{\epsilon}=\mathfrak{f}*\chi_{\epsilon}.$ Sufficiently
close to $\Sigma$, this function is smooth, growing in $t$ variable, and it
has the same sign as $\mathfrak{f}$ outside the $\epsilon$-neighborhood of
$\Sigma$. We consider the level surface $\boldsymbol{\boldsymbol{\Sigma}}$ for
$\mathfrak{f}_{\epsilon},$ $\mathfrak{f}_{\epsilon}(x)=\delta$ for some
$\delta>0.$ By construction, being the level surface of a smooth function
without stationary points, this surface is smooth, has only one intersection
point with each of trajectories $\ell_{x^{\prime}}$ and is uniformly
nontangent with these trajectories. Due to the monotonicity of $\mathfrak{f},$
the level $\delta$ can be chosen in such way that the surface
$\boldsymbol{\boldsymbol{\Sigma}}$ lies strictly inside $\Omega;$ we may also
suppose that ${\rm
dist}\,(\Sigma,\boldsymbol{\boldsymbol{\Sigma}})>{\mathbf{r}}_{0}>0.$ ∎

###  3.2. Change of variables

We introduce now new co-ordinates in $\Omega$ near $\Sigma$. Namely, for the
smooth surface $\boldsymbol{\boldsymbol{\Sigma}}$, just constructed, for a
point $x\in\widetilde{\sc\mbox{U}\hskip 1.0pt},$ we take now
$x=(x^{\prime},x_{d+1}),$ where
$x^{\prime}\in\boldsymbol{\boldsymbol{\Sigma}}$ is the intersection point of
the curve $\ell_{x}$ with $\boldsymbol{\boldsymbol{\Sigma}}$ and $x_{d+1}$ is
the distance from the point $x\in\sc\mbox{U}\hskip 1.0pt$ to the surface
$\boldsymbol{\boldsymbol{\Sigma}}$ along the integral curve $\ell_{x};$ we set
$x_{d+1}$ to be positive outside $\boldsymbol{\boldsymbol{\Sigma}}$. According
to the properties of the surface $\boldsymbol{\boldsymbol{\Sigma}},$ this is a
smooth co-ordinates change from the initial Euclidean co-ordinates.

In these co-ordinates, the initial surface $\Sigma$ is described _globally_ by
the equation $x_{d+1}=\boldsymbol{\psi}(x^{\prime})$,
$x^{\prime}\in\boldsymbol{\boldsymbol{\Sigma}},$ with some Lipschitz function
$\boldsymbol{\psi}(x^{\prime})>{\mathbf{r}}_{0}>0$ for all
$x^{\prime}\in\boldsymbol{\boldsymbol{\Sigma}}$ (recall that the surface
$\boldsymbol{\boldsymbol{\Sigma}}$ lies strictly inside $\overline{\Omega}$).
We denote by $\sc\mbox{C}\hskip 1.0pt_{0}$ the truncated cylinder with base
$\boldsymbol{\boldsymbol{\Sigma}}$ and top $\Sigma$:

$\sc\mbox{C}\hskip
1.0pt_{0}=\\{(x^{\prime},x_{d+1}):x^{\prime}\in\boldsymbol{\boldsymbol{\Sigma}},x_{d+1}\in[0,\boldsymbol{\psi}(x^{\prime})]\\}\subset{\Omega},$

which is diffeomorphic to a collar neighborhood of $\Sigma$ in $\Omega,$ and
denote by $\Omega_{0}$ the complement of this cylinder in $\Omega$,

$\Omega_{0}=\Omega\setminus\sc\mbox{C}\hskip 1.0pt_{0}$

In this setting the domain $\Omega$ can be considered as a smooth manifold
with Lipschitz boundary $\Sigma$,

$\Omega=\Omega_{0}\cup\sc\mbox{C}\hskip 1.0pt_{0}$

having a truncated cylindrical exit $\sc\mbox{C}\hskip
1.0pt_{0}=\sc\mbox{C}\hskip 1.0pt_{0}[\boldsymbol{\psi}]$ with base section
$\boldsymbol{\boldsymbol{\Sigma}}$ and boundary $\Sigma$.

Our next aim now is to perform a change of variables $\boldsymbol{\Phi}$ of
$\Omega$ to the standard domain
$\boldsymbol{\Omega}=\Omega_{0}\cup\sc\mbox{C}\hskip 1.0pt$ where C is the
_straight cylinder_ , $\sc\mbox{C}\hskip
1.0pt=[0,1]\times\boldsymbol{\boldsymbol{\Sigma}}.$ This change of variables
is the identity mapping, $\boldsymbol{\Phi}(x)=x$ on $\Omega_{0}$ and an
adjoining part of the cylinder $\sc\mbox{C}\hskip
1.0pt_{0}[\boldsymbol{\psi}]$, $x_{d+1}<t_{0}$, while near the Lipschitz
boundary, it is glued together, by means of proper cut-offs, with the
_stretching_ along the $x_{d+1}$ variable, i.e., along the integral curves of
$\boldsymbol{\Gamma}$ namely, with the mapping

$\boldsymbol{\Phi}(x)\equiv\boldsymbol{\Phi}(x^{\prime},x_{d+1})=(y^{\prime},y_{d+1})\equiv\left(x^{\prime},\frac{x_{d+1}}{\boldsymbol{\psi}(x^{\prime})}\right),\,x^{\prime}\in\boldsymbol{\boldsymbol{\Sigma}},$
(3.3)

so, at the top, $\Sigma$, we have $y_{d+1}=1$. Thus, in new co-ordinates, the
Lipschitz boundary $\Sigma$ of the domain $\Omega$ is transformed to the
smooth boundary of the domain $\boldsymbol{\Omega}$:
$y^{\prime}\in\boldsymbol{\boldsymbol{\Sigma}},y_{d+1}=1.$

We denote by $x=\boldsymbol{\Psi}(y)$ the mapping inverse to
$\boldsymbol{\Phi}$ (it, obviously, exists and it is the inverse stretching
along the $y_{d+1}$ variable near
$\boldsymbol{\boldsymbol{\Sigma}}\times\\{1\\}$). The mapping
$\boldsymbol{\Phi}$ is is a bilipschitz mapping of $\Omega$ to
$\boldsymbol{\Omega}$. Its restriction to $\Sigma$ is a bilipschitz mapping to
$\boldsymbol{\boldsymbol{\Sigma}}$. Such changes of variables leave invariant
the Sobolev spaces $H^{1}$ and for the weak derivatives the usual chain rule
holds, see, e.g., [48], Theorem 2.2.2.

The observation, crucial for our further reasoning is the following.

######  Lemma 3.2.

  1. (1)

The components of the mappings $\boldsymbol{\Phi},$ $\boldsymbol{\Psi}$ are
continuous, their first order derivatives in $x^{\prime}$, resp., $y^{\prime}$
variables and all their derivatives in ${x_{d+1}},$ resp., $y_{d+1}$ variables
are bounded functions, smooth in $x_{d+1},$ resp., $y_{d+1}$ variable; they
are continuous as functions of these variables with values in
$L_{\infty}(\boldsymbol{\boldsymbol{\Sigma}})$, resp. in $L_{\infty}(\Sigma).$

  2. (2)

The Jacobian matrices of the mappings $\boldsymbol{\Psi}$ and
$\boldsymbol{\Phi}$ are bounded; their entries are continuous in $y_{d+1},$
resp., $x_{d+1},$ variables as functions with values in
$L_{\infty}(\boldsymbol{\boldsymbol{\Sigma}})$, resp. in $L_{\infty}(\Sigma)$.

All these properties follow automatically from the definition of the mapping
$\boldsymbol{\Psi}$ and the chain rule for derivatives, taking into account
that the first order derivatives of a Lipschitz function are bounded.

##  4\. Transformations and approximation of the P-S operator.

In the course of the proof of our main result we will often transform
operators under study while preserving their spectra, so that the initial and
the approximating operator can be comfortably compared. Our considerations
take place in the domain $\Omega$ with cylindrical exit, as presented in the
previous section. Having the quadratic forms
$\boldsymbol{\boldsymbol{\rho}}_{0}$, ${\mathbf{a}}_{0}$, we find out what
happens with them under the transformation $\boldsymbol{\Psi}.$

###  4.1. Transformation of quadratic forms

We study here how the quadratic forms defining the spectrum of the Poincare-
Steklov problem transform under the change of variables, as described in Sect.
3; we set here $v(y)=u(x),$ $y=\boldsymbol{\Phi}(x)$,
$x=\boldsymbol{\Psi}(y)$, $x\in\Omega,$ $y\in\boldsymbol{\Omega}.$

First, we consider the quadratic form

$\boldsymbol{\boldsymbol{\rho}}_{0}[u]=\boldsymbol{\boldsymbol{\rho}}[v]\equiv\int_{\boldsymbol{\boldsymbol{\Sigma}}}\rho_{0}(\boldsymbol{\Psi}(y))|v(y^{\prime},1)|^{2}|\det
J(y^{\prime},1)|d\mu_{\boldsymbol{\boldsymbol{\Sigma}}}(y^{\prime}),$

where $J(y^{\prime},1)$ is the Jacobian matrix of the mapping
$\boldsymbol{\Psi}$ and $\mu_{\boldsymbol{\boldsymbol{\Sigma}}}$ is now the
natural measure on $\boldsymbol{\boldsymbol{\Sigma}},$ a smooth surface in
${\mathbb{R}}^{d+1}$. This quadratic form can be written as

$\boldsymbol{\boldsymbol{\rho}}[v]=\int_{\boldsymbol{\boldsymbol{\Sigma}}}{\rho}(y^{\prime})|v(y^{\prime},1)|^{2}d\mu_{\boldsymbol{\boldsymbol{\Sigma}}}(y^{\prime}),$

with density
${\rho}(y^{\prime})=\rho_{0}(\boldsymbol{\Psi}(y^{\prime},1))|\det
J(y^{\prime},1)|$.

The quadratic form ${\mathbf{a}}_{0}[u],$ see (2.1), is transformed to

${\mathbf{a}}_{0}[u]\equiv{{\mathbf{a}}}[v]=\int_{\boldsymbol{\Omega}}\sum_{j,k=1}^{d+1}\check{a}_{j,k}(y)\partial_{j}v(y)\overline{\partial_{k}{v(y)}}dy+\int_{\boldsymbol{\Omega}}\check{{\mathbf{v}}}(y)|v(y)|^{2}dy,$
(4.1)

$v\in H^{1}(\boldsymbol{\Omega}).$ The transformed coefficients
$\check{a}_{j,k}(y),\check{{\mathbf{v}}}(y)$, in the quadratic forms (4.1),
can be calculated, using (3.3), but we do not need here their particular
explicit expression at this point. What is important is the kind of dependence
of these coefficients on the variables $y_{k}$.

Recall that the coefficients $a_{j,k}(x)$ are bounded and are continuous at
$\Sigma$ in the sense of (1.7). This property persists after the continuous
change of variables $x=\boldsymbol{\Psi}(y)$, so the functions
$a_{j,k}(\boldsymbol{\Psi}(y))$ are continuous at the new, smooth, boundary
$\boldsymbol{\boldsymbol{\Sigma}}$ in the sense of (1.7). Next, to obtain
$\check{a}_{j,k}(y)$, we need to multiply $a_{j,k}(\boldsymbol{\Psi}(y))$ by
some algebraic combinations of derivatives of the mapping $\boldsymbol{\Psi},$
namely, the derivatives of $\boldsymbol{\Psi}$, appearing in the chain rule
when passing from $x$-derivatives to $y$-derivatives and also by the Jacobian
of the mapping $\boldsymbol{\Psi}$ arising in the passage from $dx$ in the
integral to $dy$. The multiplication by these derivatives which are bounded
but, generally, not continuous, destroys the continuity of the boundary values
of $a_{j,k}(\boldsymbol{\Psi}(y))$ for $y\in\boldsymbol{\boldsymbol{\Sigma}}$.
However, these derivatives are still continuous in $y_{d+1}$ variable as
functions with values in $L_{\infty}(\boldsymbol{\boldsymbol{\Sigma}})$. We
introduce the following definition.

######  Definition 4.1.

Let $q(y)=q(y^{\prime},y_{d+1})$ be a bounded measurable function in the
variables $y^{\prime}\in\boldsymbol{\boldsymbol{\Sigma}},$
$y_{d+1}\in(1-\delta^{\prime},1].$ Let ${{\mathbf{q}}}(y^{\prime})\in
L_{\infty}(\boldsymbol{\boldsymbol{\Sigma}})$ be fixed. Denote by
${{\mathbf{q}}^{\diamond}}(y)\equiv{{\mathbf{q}}^{\diamond}}(y^{\prime},y_{d+1})$
the continuation of ${{\mathbf{q}}}(y^{\prime})$ as a function not depending
on $y_{d+1},$
${{\mathbf{q}}^{\diamond}}(y^{\prime},y_{d+1}):={{\mathbf{q}}}(y^{\prime}),$
$y^{\prime}\in\boldsymbol{\boldsymbol{\Sigma}}.$ If for any
$y^{\prime}\in\boldsymbol{\boldsymbol{\Sigma}},$

$\lim_{\delta\to
0}\|q-{\mathbf{q}}^{\diamond}\|_{L_{\infty}(B((y^{\prime},1),\delta)}=0,$
(4.2)

we say that $q(y)$ has $L_{\infty}(\boldsymbol{\boldsymbol{\Sigma}})$ limit
${{\mathbf{q}}}(y^{\prime})$ at $\boldsymbol{\boldsymbol{\Sigma}}$ in the
variable $y_{d+1}\to 1$. If this is the case, we redefine the initial function
$q(y)$ on $\boldsymbol{\boldsymbol{\Sigma}},$ by setting
$q(y^{\prime},1)={\mathbf{q}}(y^{\prime})$, thus changing $q$ on a set of
$d+1$-dimensional measure $0$. After this, we say that $q$ is
$L_{\infty}(\boldsymbol{\boldsymbol{\Sigma}})$-continuous at
$\boldsymbol{\boldsymbol{\Sigma}}.$

One can easily see that the product of a function
$L_{\infty}(\boldsymbol{\boldsymbol{\Sigma}})$-continuous at
$\boldsymbol{\boldsymbol{\Sigma}}$ and a function continuous at
$\boldsymbol{\boldsymbol{\Sigma}}$ is again a function
$L_{\infty}(\boldsymbol{\boldsymbol{\Sigma}})$-continuous at
$\boldsymbol{\boldsymbol{\Sigma}}.$ Therefore, by Lemma 3.2, the coefficients
$\check{a}_{j,k}$ of the quadratic form ${\mathbf{a}}[v]$ in (4.1) are
$L_{\infty}(\boldsymbol{\boldsymbol{\Sigma}})$-continuous at
$\boldsymbol{\boldsymbol{\Sigma}}.$ This ’very week’ continuity property
proves to be sufficient for our approximation construction.

We recall that the ratio of quadratic forms
$\boldsymbol{\boldsymbol{\rho}}_{0}/{\mathbf{a}}_{0}$ defines an operator
whose spectrum coincides with the one of the P-S operator. The domain of the
quadratic forms is $H^{1}(\boldsymbol{\Omega})$. Therefore, our
transformations can be summarized in the following way.

######  Proposition 4.2.

The spectrum of the variational Poincaré-Steklov problem (1.3) coincides with
the spectrum of the ratio

$\frac{{\boldsymbol{\boldsymbol{\rho}}}[v]}{{{\mathbf{a}}}[v]},\,v\in
H^{1}(\boldsymbol{\Omega}),$ (4.3)

in a domain $\boldsymbol{\Omega}$ with smooth boundary
$\boldsymbol{\boldsymbol{\Sigma}}$, and where the quadratic form
${\mathbf{a}}[v]$ is elliptic and has bounded coefficients, which are
$L_{\infty}(\boldsymbol{\boldsymbol{\Sigma}})$ \- continuous in $y_{d+1}$
variable at the boundary
$\boldsymbol{\boldsymbol{\Sigma}}=\partial\boldsymbol{\Omega}.$

From now on, we may forget about the initial Lipschitz problem and study, from
scratch, the eigenvalue distribution for the ratio (4.3). To simplify
notations, we drop further on the ’check’ over the symbols and, again, denote
the main domain by $\Omega$ and its, now smooth, boundary by $\Sigma$
($=\Sigma\times\\{1\\}$ in our $y$ co-ordinates). Recall that we have paid for
the smoothness of the boundary by weakening the continuity property of the
coefficient matrix at the boundary, from continuity to $L_{\infty}(\Sigma)$
continuity. In the notations of the quadratic forms
$\boldsymbol{\boldsymbol{\rho}}_{0},$ ${\mathbf{a}}_{0}$ we drop the subscript
’0’.

###  4.2. An operator representation of the eigenvalue problem.

We recall that we consider the case of a non-negative function $\rho.$

The quadratic form ${\mathbf{a}}[v]$ with domain $H^{1}(\Omega)$ defines the
self-adjoint positive ’Neumann’ operator
${\mathbf{T}}={\mathbf{T}}_{\mathfrak{a}}$ in $L_{2}(\Omega).$ This elliptic
operator acts, formally, as

${\mathbf{T}}v\equiv{\mathbf{T}}_{\mathfrak{a}}v=-\sum_{j,k}\partial_{j}a_{j,k}\partial_{k}v+{\mathbf{v}}v,$

in the sense of distributions, but its domain is rather hard to describe due
to the lack of continuity of the coefficients. Fortunately, neither the domain
nor the exact action are needed in our reasoning. Instead, since
${\mathbf{T}}$ is a positive self-adjoint elliptic operator with bounded
coefficients, the domain of the square root ${\mathbf{T}}^{\frac{1}{2}}$ can
be described explicitly: it coincides with $H^{1}(\Omega)$ and
${\mathbf{T}}^{\frac{1}{2}}$ is an isomorphism from $H^{1}(\Omega)$ to
$L_{2}(\Omega)$.

It was already stated that we are free in the choice of the function
${\mathbf{v}}$. Now we fix this choice, supposing that ${\mathbf{v}}\in
C_{0}^{\infty}(\Omega)$ (after the above co-ordinates change).

Since
${\mathbf{a}}[v]=({\mathbf{T}}^{\frac{1}{2}}v,{\mathbf{T}}^{\frac{1}{2}}v),$
we can set $f={\mathbf{T}}^{\frac{1}{2}}v\in L_{2}(\Omega)$ and therefore the
ratio (4.3) transforms to

$\frac{\boldsymbol{\boldsymbol{\rho}}[{\mathbf{T}}^{-\frac{1}{2}}f]}{\|f\|_{L_{2}(\Omega)}^{2}}=\frac{\int_{\Sigma}|(\gamma{\mathbf{T}}^{-\frac{1}{2}}f)(y^{\prime})|^{2}\rho(y^{\prime})d\mu_{\Sigma}(y^{\prime})}{\|f\|_{L_{2}(\Omega)}^{2}},\,f\in
L_{2}(\Omega).$ (4.4)

Here and further on, we denote by $\gamma:f\mapsto\gamma f$ the trace operator
acting as the restriction to $\Sigma$ of a function $f$ defined on the domain
$\Omega$ with cylindric exit C having the boundary lid $\Sigma$ identified
with $\Sigma\times\\{1\\}$. We recall that the boundary $\Sigma$ is smooth.
Therefore, this restriction is known to act as a bounded operator from
$H^{s}(\Omega)$ to $H^{s-1/2}(\Sigma)$ for all _positive_ $s$ _except the
half-integer ones_ ($s\in{\mathbb{N}}+\frac{1}{2}).$ In particular, this
implies that the operator $\gamma{\mathbf{T}}^{-\frac{1}{2}}$ is bounded as
acting from $L_{2}(\Omega)$ to $L_{2}(\Sigma)$ (and even to
$H^{\frac{1}{2}}(\Sigma)$). Following Sect.2, we restrict ourselves to
$\rho\geq 0$ and $\rho\in L_{\infty},$ therefore, the operator
$\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}:L_{2}(\Omega)\to
L_{2}(\Sigma)$ is bounded as well.

As a result, the numerator in (4.4) can be written as

$\displaystyle{\boldsymbol{\boldsymbol{\rho}}[{\mathbf{T}}^{-\frac{1}{2}}f]}=(\rho(y^{\prime})^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}f,\rho(y^{\prime})^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}f)_{L_{2}(\Sigma)}=$
$\displaystyle\left(\left[\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}\right]^{*}[\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}]f,f\right)_{L_{2}(\Omega)}.$

In this way, we have reduced our spectral problem to the study of the spectrum
of the self-adjoint operator

${\mathbf{G}}_{\mathfrak{a}}=\left[\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}\right]^{*}\left[\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}\right],$

acting in $L_{2}(\Omega)$. Since the nonzero discrete spectrum of the product
does not change under the cyclic permutation of the factors, the nonzero
eigenvalues of ${\mathbf{G}}_{\mathfrak{a}}$ are the same as the ones of the
self-adjoint operator

${\mathbf{H}}_{\mathfrak{a}}\equiv\left[\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}\right]\left[\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}\right]^{*},$
(4.5)

acting in $L_{2}(\Sigma).$

It is more convenient to have a somewhat different representation of the
operator ${\mathbf{H}}_{\mathfrak{a}}.$

######  Lemma 4.3.

The following equality holds

${\mathbf{H}}_{\mathfrak{a}}=\left[\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}\right]\left[\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}\right]^{*}=\rho^{\frac{1}{2}}\gamma(\gamma{\mathbf{T}}^{-1})^{*}\rho^{\frac{1}{2}}=\rho^{\frac{1}{2}}\gamma(\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-1})^{*}.$
(4.6)

###### Proof.

We start with the equality

$[\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}]{\mathbf{T}}^{-\frac{1}{2}}=\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-1}:L_{2}(\Omega)\to
L_{2}(\Sigma),$ (4.7)

where we have the product of two bounded operators on the left. Next we take
adjoint of both parts in (4.7),

${\mathbf{T}}^{-\frac{1}{2}}[\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-\frac{1}{2}}]^{*}=(\rho^{\frac{1}{2}}\gamma{\mathbf{T}}^{-1})^{*}:L_{2}(\Sigma)\to
L_{2}(\Omega).$ (4.8)

This is a bounded operator acting from $L_{2}(\Sigma)$ to $H^{1}(\Omega)$,
again a product of two bounded operators. Therefore, the restriction operator
$\gamma$ is well defined on the range of the operator in (4.8) and we can
apply $\gamma$ on the left, and then multiply by the bounded function
$\rho^{\frac{1}{2}},$ which gives us (4.6). ∎

###  4.3. Smooth approximation

Our study of P-S eigenvalues will be based upon a special smooth approximation
of the coefficient matrix $\mathfrak{a}.$

######  Lemma 4.4.

Let $\mathfrak{a}\in L_{\infty}(\Omega)$ be a Hermitian matrix-function,
continuous at $\Sigma$ as a function of $y_{d+1}$ variable with values in
$L_{\infty}(\Sigma)$ in the sense of Definition 4.2. Then for any $\epsilon>0$
and $p<\infty,$ there exist a $C^{\infty}(\overline{\Omega})$-smooth matrix-
function $\tilde{\mathfrak{a}}(y),$ $y\in\bar{\Omega},$ such that the
following conditions are satisfied

  1. (1)

$\|\mathfrak{a}(.,1)-\tilde{\mathfrak{a}}(.,1)\|_{L_{p}(\Sigma)}<\epsilon$;

  2. (2)

the function $\mathfrak{a}(y)-\tilde{\mathfrak{a}}(y)$ is continuous at
$\Sigma$ as a function of $y_{d+1}$ variable with values in
$L_{\infty}(\Sigma);$

  3. (3)

$\|\tilde{\mathfrak{a}}\|_{L_{\infty}(\Omega)}\leq
c\|\mathfrak{a}\|_{L_{\infty}(\Omega)}$;

  4. (4)

if the matrix $\mathfrak{a}$ is uniformly positive in $\Omega,$ then,
additionally, $\|\tilde{\mathfrak{a}}^{-1}\|_{L_{\infty}(\Omega)}\leq
c\|\mathfrak{a}^{-1}\|_{L_{\infty}(\Omega)}.$

with the constant $c$ not depending on $\mathfrak{a}.$

The condition (3) means that the $L_{\infty}$ norm for the approximating
matrix can be arranged to be not depending on $\epsilon$, being controlled, by
the $L_{\infty}$ norm of the initial matrix $\mathfrak{a};$ by condition (4),
in the elliptic case, the same can be arranged for the ellipticity constant of
the approximating matrix. Note that the approximation quality in Lemma 4.4 is
required only at the boundary $\Sigma$.

###### Proof.

We discuss the ellipticity case, $\mathfrak{a}(y)\geq C>0,$ first. Consider
the constant matrix $\mathfrak{a}_{1}(y)=C_{1}\mathbf{I}_{(d+1)\times(d+1)}$,
where $\|\mathfrak{a}^{-1}\|_{L_{\infty}(\Omega)}^{-1}\leq
C_{1}\leq\|\mathfrak{a}\|_{L_{\infty}(\Omega)}.$ Next, to construct
$\tilde{\mathfrak{a}}(y)$ near the boundary, we take $\delta>0$ such that for
any ball $B(y^{\prime},\delta),\,y^{\prime}\in\Sigma$ the inequality
$\|\mathfrak{a}(y)-\mathfrak{a}^{\diamond}(y)\|_{L_{\infty}(B(y^{\prime},\delta))}<\epsilon$
holds (where ${\mathfrak{a}}^{\diamond}(y)$ is the continuation of
$\mathfrak{a}(y^{\prime},1)$ from $\Sigma$ to the $\delta$-neighborhood of
$\Sigma,$
${\mathfrak{a}}^{\diamond}(y^{\prime},y_{d+1})=\mathfrak{a}(y^{\prime},1).$)
Next we approximate the matrix $\mathfrak{a}(.,1)$ by a smooth matrix
$\tilde{\mathfrak{a}}_{0}\in L_{p}(\Sigma)\cap L_{\infty}(\Sigma)$ such that
$\|\tilde{\mathfrak{a}}_{0}(.)-\mathfrak{a}(.,1)\|_{L_{p}(\Sigma)}<\epsilon,$
thus the condition 1) in Lemma 4.4 is satisfied. Also we require that

$\|\tilde{\mathfrak{a}}_{0}\|_{L_{\infty}(\Sigma)}\leq
C\|\mathfrak{a}\|_{L_{\infty}(\Omega)},\,\|\tilde{\mathfrak{a}}_{0}^{-1}\|_{L_{\infty}(\Sigma)}\leq
C^{\prime}\|\mathfrak{a}^{-1}\|_{L_{\infty}(\Omega)};$ (4.9)

all these conditions can be satisfied, for example, by setting
$\tilde{\mathfrak{a}}_{0}(.)={\mathbf{V}}(\tau)\mathfrak{a}(.,1)$, where
${\mathbf{V}}(\tau)$ is the heat semigroup on $\Sigma,$ since
${\mathbf{V}}(\tau)\to\boldsymbol{1}$ strongly in $L_{p}(\Sigma),$ $p<\infty$
as $\tau\to+0$ and ${\mathbf{V}}(\tau)$ is bounded in $L_{\infty},$ $\tau\geq
0.$ We consider now the matrix ${\tilde{\mathfrak{a}}_{0}}^{\diamond}(y),$ the
continuation of $\tilde{\mathfrak{a}}_{0}(y^{\prime}),$ namely
${\tilde{\mathfrak{a}}_{0}}^{\diamond}(y^{\prime},y_{d+1})=\tilde{\mathfrak{a}}_{0}(y^{\prime}).$
This smooth matrix, not depending on $y_{d+1}$ is, of course, continuous in
the $y_{d+1}$ variable as a function with values in $L_{\infty}(\Sigma)$.
Therefore, the difference,
$\mathfrak{a}(y)-{\tilde{\mathfrak{a}}_{0}}^{\diamond}(y)$ is also continuous
at $\Sigma$ as a function of $y_{d+1}$ with values in $L_{\infty}(\Sigma)$.
This grants the condition 2) in Lemma. Conditions 3) and 4) are satisfied due
to (4.9). Finally, we glue together the matrix
${\tilde{\mathfrak{a}}_{0}}^{\diamond}$ defined in a neighborhood of $\Sigma$
and the matrix $\mathfrak{a}_{1}(y)$ defined outside such neighborhood by
means of proper cut-off functions, and this matrix satisfies all conditions of
Lemma.

In the case when the matrix $\mathfrak{a}$ is not supposed to be uniformly
positive definite, the simplified reasoning succeeds: namely we just skip all
estimates concerning $\mathfrak{a}^{-1}$, therefore, we may simply set
$\mathfrak{a}_{1}(y)=0$. ∎

We will use Lemma 4.4 twice; first time when approximating the elliptic
coefficients matrix $\mathfrak{a}(y),$ where the control of the ellipticity
constant is important, and the second time for the approximation of the matrix
$\mathfrak{b}=\mathfrak{a}-\tilde{\mathfrak{a}}$ when proving the required
operator convergence of ${\mathbf{H}}_{\tilde{\mathfrak{a}}}$ to
${\mathbf{H}}_{\mathfrak{a}}$.

###  4.4. The approximating operator

Consider, for a given $\epsilon>0$, the approximating coefficient matrix
$\tilde{\mathfrak{a}}$ constructed in Section 4.3 for the elliptic coefficient
matrix $\mathfrak{a}$ and the corresponding elliptic differential operator
${\tilde{\mathcal{L}}}={\mathcal{L}}_{\tilde{\mathfrak{a}}},$ _with the same
zero order term_ ${\mathbf{v}}\in C^{\infty}(\overline{\Omega}).$ With the
coefficient matrix $\tilde{\mathfrak{a}}$ we associate the Neumann operator
$\tilde{{\mathbf{T}}}={\mathbf{T}}_{\tilde{\mathfrak{a}}},$ by means of the
quadratic form $\boldsymbol{\tilde{\mathfrak{a}}}[v],$ $v\in H^{1}(\Omega)$.

Since the boundary of $\Omega$ is smooth, as well as the coefficients matrix
$\tilde{\mathfrak{a}}$, some more can be said about the properties of the
operator (4.5) for this matrix. The corresponding compact operator, which we
denote by ${{\mathbf{H}}}_{\tilde{\mathfrak{a}}},$ has, by Lemma 4.3, the form

${{\mathbf{H}}}_{\tilde{\mathfrak{a}}}=\rho^{\frac{1}{2}}\gamma(\gamma{\mathbf{T}}_{\tilde{\mathfrak{a}}}^{-1})^{*}\rho^{\frac{1}{2}}.$
(4.10)

First of all, for the eigenvalues of this operator, which coincide with
eigenvalues of the corresponding $N\\!\\!\\!D$ operator, the asymptotic
formula holds, of the form (1.6), (1.4), (1.5), with the matrix
$\tilde{\mathfrak{a}}$ replacing $\mathfrak{a},$ due to the existing results,
see Theorem 1.1. Next, since the coefficients of
${\mathbf{T}}_{\tilde{\mathfrak{a}}}$ are smooth, the domain and action of the
operator ${\mathbf{T}}_{\tilde{\mathfrak{a}}}$ can be described explicitly.
Namely, by the standard elliptic regularity results, the domain
${\mathcal{D}}(\tilde{{\mathbf{T}}})$ of $\tilde{{\mathbf{T}}}$ is the
subspace in the Sobolev space $H^{2}(\Omega)$, consisting of functions $v$
satisfying the classical Neumann boundary conditions,
$\sum_{k}\tilde{a}_{j,k}\nu_{k}\partial_{j}v|_{\partial\Omega}=0,$ where
$\nu_{k}$ are the components of the unit normal vector to
$\Sigma=\partial\Omega.$

######  Lemma 4.5.

For a smooth matrix $\tilde{\mathfrak{a}},$ for any $f\in L_{2}(\Sigma)$, the
function
$v=(\rho^{\frac{1}{2}}\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}\rho^{\frac{1}{2}}f$
satisfies in $\Omega$ the elliptic equation

$\tilde{{\mathcal{L}}}v\equiv\tilde{{\mathcal{L}}}(\rho^{\frac{1}{2}}\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}\rho^{\frac{1}{2}}f=0\,\,\mathrm{in}\,\,\Omega.$
(4.11)

###### Proof.

Let $g\in C^{\infty}_{0}(\Omega)$ be a smooth function with compact support in
$\Omega.$ Consider the sesquilinear form
${\mathbf{I}}[f,g]=(\tilde{{\mathcal{L}}}(\rho^{\frac{1}{2}}\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}\rho^{\frac{1}{2}}f,g)_{L_{2}(\Omega)}.$
The function $g$ belongs to the domain of the self-adjoint operator
$\tilde{{\mathbf{T}}},$ therefore,

${\mathbf{I}}[f,g]=((\rho^{\frac{1}{2}}\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}\rho^{\frac{1}{2}}f,\tilde{{\mathbf{T}}}{}g)_{L_{2}(\Omega)}=(\rho^{\frac{1}{2}}f,\rho^{\frac{1}{2}}\gamma\tilde{{\mathbf{T}}}{}^{-1}\tilde{{\mathbf{T}}g})_{L_{2}(\Sigma)}=(f,\rho\gamma
g)_{L_{2}(\Sigma)}.$ (4.12)

Since $g=0$ on $\Sigma,$ the last expression in (4.12) is zero. Therefore, the
function
$\tilde{{\mathcal{L}}}(\rho^{\frac{1}{2}}\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}\rho^{\frac{1}{2}}f$
is orthogonal to $C^{\infty}_{0}(\Omega),$ and this implies that
$\tilde{{\mathcal{L}}}(\rho^{\frac{1}{2}}\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}\rho^{\frac{1}{2}}f$
is zero in $\Omega.$ ∎

###  4.5. ${\mathbf{H}}_{\mathfrak{a}}$ as a pseudodifferential operator

Here we establish an important property of the approximating operator
${\mathbf{H}}_{\tilde{\mathfrak{a}}},$ needed further for evaluating the rate
of approximation of spectra.

######  Proposition 4.6.

For a smooth approximating matrix $\tilde{\mathfrak{a}}$ and a smooth weight
$\rho$, the operator ${\mathbf{H}}_{\tilde{\mathfrak{a}}}$ is an order $-1$
pseudodifferential operator on the boundary $\Sigma.$

###### Proof.

When finding, in the smooth case, an expression for the operator
${\mathbf{H}}_{\tilde{\mathfrak{a}}},$ we study its component
$(\rho^{\frac{1}{2}}\gamma\tilde{{\mathbf{T}}}^{-1})^{*}$ first. Recall that
$\tilde{{\mathbf{T}}}{}^{-1}$ is the inverse of the realization of the
operator ${\tilde{\mathcal{L}}}$ with Neumann boundary conditions,
$\gamma_{1}u=0$, where $\gamma_{1}$ is the conormal derivative. Continue the
coefficients $\tilde{\mathfrak{a}}$ outside $\Omega$ in a smooth way. We
denote by ${\mathbf{R}}$ the fundamental solution for ${\tilde{\mathcal{L}}}$;
so that ${\tilde{\mathcal{L}}}{\mathbf{R}}-1,$
${\mathbf{R}}{\tilde{\mathcal{L}}}-1$ are infinitely smoothing operators;
${\mathbf{R}}$ is the order $-2$ pseudodifferential operator with symbol
${\mathbf{r}}(y,\eta)$ having the principal term
${\mathbf{r}}_{-2}(y,\eta)=\boldsymbol{\tilde{\mathfrak{a}}}_{y}(\eta)^{-1},$
where
$\boldsymbol{\tilde{\mathfrak{a}}}_{y}(\eta)=\sum_{j,k}\tilde{a}_{j,k}(y)\eta_{j}\eta_{k}$
is the principal symbol of the operator ${\tilde{\mathcal{L}}}$. Let $d>1$
first. Then the operator ${\mathbf{R}}$ acts on functions on $\Omega$ as an
integral operator with kernel $R(y,z-y)$ having leading singularity of order
$1-d$ at the diagonal $y=z$. Namely, the principal term, $R_{1-d}(y,y-z),$ the
leading singularity in $R(y,z-y)$ is the Fourier transform of the principal
symbol ${\mathbf{r}}_{-2}(y,\eta)$,

$R_{1-d}(y,y-z)={\mathcal{F}}_{\eta\to y-z}{\mathbf{r}}_{-2}(y,\eta).$

For the dimension $d+1=2$, the kernel $R(y,z-y)$ has logarithmic singularity
at the diagonal.

Let $f$ be a function in $L_{2}(\Omega)$. In order to find a representation
for $u=\tilde{{\mathbf{T}}}{}^{-1}f,$ we set
$u_{0}=\int_{\Omega}R(y,y-z)f(z)dz$. This function satisfies the equation
${\tilde{\mathcal{L}}}u_{0}=f$ in $\Omega$ but not the Neumann boundary
condition $\gamma_{1}u=0$, where $\gamma_{1}$ is the conormal derivative at
the boundary, corresponding to the elliptic operator ${\tilde{\mathcal{L}}}$.
We construct the correction $u_{1},$ such that ${\tilde{\mathcal{L}}}u_{1}=0,$
$\gamma_{1}u_{1}=-\gamma_{1}u_{0}$. To do this, we consider the _Green
function_ ${\mathcal{G}}_{1}(y,z)$ for the Neumann problem, in other words,
the integral kernel of the solution operator $\boldsymbol{G}_{1}$ of the
Neumann problem for ${\tilde{\mathcal{L}}}$ . This means that for a smooth
function $h$ on $\Sigma$,

${\tilde{\mathcal{L}}}\boldsymbol{G}_{1}h:={\tilde{\mathcal{L}}}_{y}\int_{\Sigma}{\mathcal{G}}_{1}(y,z)h(z)dz=0,\,y\in\Omega,$
(4.13)

and

$\lim_{y\to
y_{0}\in\Sigma}\gamma_{1,y}\int_{\Sigma}{\mathcal{G}}_{1}(y,z)h(z)dz=h(y_{0}).$
(4.14)

The operator $\boldsymbol{G}_{1}:C^{\infty}(\Sigma)\to C^{\infty}(\Omega)$ is
a Poisson operator in the Boutet-de-Monvel algebra. By (4.14), it satisfies
$\gamma_{1}\boldsymbol{G}_{1}:C^{\infty}(\Sigma)\to
C^{\infty}(\Sigma)=\boldsymbol{1},$ i.e., it gives the identity operator.
Therefore the solution of the boundary problem ${\tilde{\mathcal{L}}}u=f$ in
$\Omega,$ $\gamma_{1}u=0$ on $\Sigma$ can be expressed as

$\displaystyle u(y)=u_{0}(y)+u_{1}(y)=$
$\displaystyle\int_{\Omega}R(y,y-z)f(z)dz-\int_{\Sigma}{\mathcal{G}}_{1}(y,z)\gamma_{1}\left(\int_{\Omega}R(y,y-z)f(z)dz\right)$
$\displaystyle=({\mathbf{R}}f)(y)-(\boldsymbol{G}_{1}\gamma_{1}{\mathbf{R}}f)(y).$

As a result, the operator $\tilde{{\mathbf{T}}}{}^{-1}$ can be represented as

$\tilde{{\mathbf{T}}}{}^{-1}={\mathbf{R}}-\boldsymbol{G}_{1}\gamma_{1}{\mathbf{R}}=(1-\boldsymbol{G}_{1}\gamma_{1}){\mathbf{R}};$

such representation, in a somewhat different setting, can be found, e.g., in
[21], Theorem 9.20. There, after the standard straightening of the boundary,
the operator ${\mathbf{R}}$ is treated as a truncated pseudodifferential
operator,

${\mathbf{R}}=\chi_{\Omega}\mathfrak{R}{\mathbf{e}}_{\Omega},$ (4.15)

where ${\mathbf{e}}_{\Omega}$ is the operator of extension by zero of a
function in $\Omega$ to the whole ${\mathbb{R}}^{d+1}$, $\mathfrak{R}$ is a
pseudodifferential operator in ${\mathbb{R}}^{d+1}$ with symbol
${\mathbf{r}}(y,\eta)$ and $\chi_{\Omega}$ is the restriction of functions in
${\mathbb{R}}^{d+1}$ to $\Omega.$

We also need an expression for the operator $\tilde{{\mathbf{T}}}_{0}^{-1}$,
the resolvent of the _Dirichlet_ problem for the equation
${\tilde{\mathcal{L}}}u=f$ in $\Omega.$ In a similar way,

$\tilde{{\mathbf{T}}_{0}}{}^{-1}={\mathbf{R}}-\boldsymbol{G}_{0}\gamma{\mathbf{R}}=(1-\boldsymbol{G}_{0}\gamma){\mathbf{R}},$
(4.16)

where $\boldsymbol{G}_{0}$ is the Poisson operator (the Green function)
solving the nonhomogeneous Dirichlet boundary problem for
${\tilde{\mathcal{L}}}$ in $\Omega$: ${\tilde{\mathcal{L}}}u=0,\,\gamma u=h;$
this means that $\gamma\boldsymbol{G}_{0}=\boldsymbol{1}.$

We pass now to $\gamma\tilde{{\mathbf{T}}}{}^{-1}.$ Since
$\gamma{\tilde{{\mathbf{T}}}_{0}{}^{-1}}=0,$ we can write

$\displaystyle\gamma\tilde{{\mathbf{T}}}{}^{-1}=\gamma(\tilde{{\mathbf{T}}}{}^{-1}-\tilde{{\mathbf{T}}}_{0}{}^{-1})=\gamma(\boldsymbol{G}_{0}\gamma{\mathbf{R}}-\boldsymbol{G}_{1}\gamma_{1}{\mathbf{R}})$
(4.17)
$\displaystyle=\gamma\boldsymbol{G}_{0}\gamma{\mathbf{R}}-\gamma\boldsymbol{G}_{1}\gamma_{1}{\mathbf{R}}=(\gamma\boldsymbol{G}_{0}\gamma-\gamma\boldsymbol{G}_{1}\gamma_{1}){\mathbf{R}}.$

We recall now that $\gamma\boldsymbol{G}_{0}=\boldsymbol{1}.$ Further on, the
function $u_{1}=\boldsymbol{G}_{1}\gamma_{1}{\mathbf{R}}f$ satisfies the
equation ${\tilde{\mathcal{L}}}u_{1}=0$ in $\Omega$, therefore $\gamma u_{1}$
and $\gamma_{1}u_{1}$ are connected by the Neumann-to Dirichlet operator
$N\\!\\!\\!D$, $\gamma u_{1}=N\\!\\!\\!D\gamma_{1}u_{1}$. We set all this into
(4.17) and obtain

$\gamma\tilde{{\mathbf{T}}}^{-1}=(\gamma-N\\!\\!\\!D\gamma_{1}){\mathbf{R}}.$
(4.18)

Consequently,

$(\gamma\tilde{{\mathbf{T}}}^{-1})^{*}=(\gamma{\mathbf{R}})^{*}-(N\\!\\!\\!D\gamma_{1}{\mathbf{R}})^{*}.$

Finally, we have

$\gamma(\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}=\gamma(\gamma{\mathbf{R}})^{*}-\gamma(N\\!\\!\\!D\gamma_{1}{\mathbf{R}})^{*}.$
(4.19)

Both terms on the right in (4.19) belong to the Boutet-de-Monvel algebra, see
e.g., [21], where the composition rules are described in detail. In this
setting, ${\mathbf{R}}$ is a truncated pseudodifferential operator of order
$-2$, $\gamma$ is a trace operator of zero order, further,
$(\gamma{\mathbf{R}})^{*}$ is a Poisson operator, and as a result,
$\gamma(\gamma{\mathbf{R}})^{*}$ is a pseudodifferential operator of order
$-1.$ Similarly, $\gamma_{1}$ is a trace operator of order $1$, the Neumann-to
Dirichlet operator $N\\!\\!\\!D$ is an order $-1$ pseudodifferential operator
on the boundary, therefore, the second term in (4.19) is a pseudodifferential
operator of order $-1$ as well. Its principal symbol can be expressed in a
standard way, algebraically, via the principal symbols of ${\mathbf{R}}$ and
$N\\!\\!\\!D$, as we will see in Sect.7. Finally, the multiplication by the
smooth function $\rho^{\frac{1}{2}}$ produces, again, a pseudodifferential
operator. ∎

##  5\. Operator perturbations

Our approach to establishing the asymptotic formula for Steklov eigenvalues,
similarly to [37], is based upon the operator approximation. This
approximation idea was successfully used by M.S. Agranovich in [2] on the base
of the variational setting of the problem. In fact, for more smooth
boundaries, the ones of class $C^{1,1}$, a similar perturbational approach was
even used in [39], the earliest paper on the Steklov eigenvalue asymptotics.

Namely, if the coefficients matrix $\mathfrak{a}(y)$ is continuous, it can be
approximated in $C(\Omega)$, both from above and from below, by smooth
elliptic matrices $\tilde{\mathfrak{a}}_{\pm}(y)$, and the required closeness
of spectra of operators ${\mathbf{T}}$ and $\tilde{{\mathbf{T}}}$ follows in
[39] and in [2] from rather simple monotonicity estimates. The weakening of
conditions imposed on the boundary in _this_ paper, actually, the passage from
the _continuous_ derivatives of the function $\psi$ defining the boundary to
the function having only _bounded_ derivatives, seemingly, a minor one, is, in
fact, rather essential, the resulting coefficients $a_{j,k}$ are not
continuous any more and they cannot be approximated in the $C$ metric by
smooth ones. Instead, we use an approximation in a weaker, $L_{p}$ sense, as
in Lemma 4.4, which turns out to be sufficient.

###  5.1. Approximation of the operator

In this section, we find the expression for the difference of operators
describing the eigenvalues of the N-to-D operators. Let
$\mathfrak{a},\tilde{\mathfrak{a}}$ be the matrices of coefficients of the
operators ${\mathcal{L}},\tilde{{\mathcal{L}}}$, described in Sect.4, so that
$\mathfrak{a},\tilde{\mathfrak{a}}{}^{-1}$ belong to $L_{\infty}(\Omega)$,
$\tilde{\mathfrak{a}},\tilde{\mathfrak{a}}^{-1}\in
C^{\infty}(\overline{\Omega})$ and
$\tilde{\mathfrak{a}}|_{\Sigma}-\mathfrak{a}|_{\Sigma}$ is small in the sense
of Lemma 4.4.

Consider the operators ${\mathbf{T}},\tilde{{\mathbf{T}}}$, the Neumann
operators for ${\mathcal{L}},{\tilde{\mathcal{L}}}$.

######  Lemma 5.1.

Under the above conditions, the following factorization is valid:

${\mathbf{T}}^{-1}-\tilde{{\mathbf{T}}}{}^{-1}={\mathbf{X}}^{*}{\mathbf{Y}},$
(5.1)

where ${\mathbf{X}},{\mathbf{Y}}$ are bounded operators acting from
$L_{2}(\Omega)$ to $L_{2}(\Omega)\otimes{\mathbb{C}}^{d+1}$.

${\mathbf{X}}=\nabla{\mathbf{T}}^{-1},{\mathbf{Y}}=\mathfrak{b}(y)\nabla\tilde{{\mathbf{T}}}{}^{-1},$
(5.2)

with $\mathfrak{b}=(\tilde{\mathfrak{a}}-\mathfrak{a}).$

###### Proof.

Essentially, our Lemma 5.1 is an analogy of Lemma 8.1 in [13]. For the sake of
completeness, we present the detailed proof of (5.1) in our setting.

The equality (5.1) is equivalent to

$({\mathbf{T}}^{-1}f,g)_{L_{2}(\Omega)}-(\tilde{{\mathbf{T}}}{}^{-1}f,g)_{L_{2}(\Omega)}=\int_{\Omega}\langle(\tilde{\mathfrak{a}}-\mathfrak{a})\nabla\tilde{{\mathbf{T}}}{}^{-1}f,\nabla{\mathbf{T}}^{-1}g\rangle
dy,\,f,g\in L_{2}(\Omega),$ (5.3)

where the angle brackets denote the scalar product in ${\mathbb{C}}^{d+1}.$

In our conditions, the domains of the operators
${\mathbf{T}}^{\frac{1}{2}},\tilde{{\mathbf{T}}}{}^{\frac{1}{2}}$ coincide,
both are equal to the Sobolev space $H^{1}(\Omega).$ It follows, in
particular, that ${\mathcal{D}}(\tilde{{\mathbf{T}}})\subset
H^{2}(\Omega)\subset{\mathcal{D}}({\mathbf{T}}^{\frac{1}{2}}).$ Let $u,v$ be
arbitrary functions in
${\mathcal{D}}({\mathbf{T}}^{\frac{1}{2}})={\mathcal{D}}(\tilde{{\mathbf{T}}}{}^{\frac{1}{2}})=H^{1}(\Omega).$
Consider the equality

$\displaystyle(\tilde{{\mathbf{T}}}{}^{\frac{1}{2}}u,\tilde{{\mathbf{T}}}{}^{\frac{1}{2}}v)_{L_{2}(\Omega)}-({\mathbf{T}}^{\frac{1}{2}}u,{\mathbf{T}}^{\frac{1}{2}}v)_{L_{2}(\Omega)}=\tilde{{\mathbf{a}}}[u,v]-{\mathbf{a}}[u,v]=$
(5.4) $\displaystyle\int_{\Omega}\langle\tilde{\mathfrak{a}}\nabla u,\nabla
v\rangle dy-\int_{\Omega}\langle\mathfrak{a}\nabla u,\nabla v\rangle dy$
$\displaystyle=\int_{\Omega}\langle(\tilde{\mathfrak{a}}-\mathfrak{a})\nabla
u,\nabla v\rangle dy.$

We set here $u=\tilde{{\mathbf{T}}}{}^{-1}f,$ $v={\mathbf{T}}^{-1}g,$ where
$f,g$ are arbitrary elements in $L_{2}(\Omega)$. These functions $u,v$ belong
to $H^{1}(\Omega)$, therefore (5.4) is satisfied. Such substitution leads to
(5.3):

$(\tilde{{\mathbf{T}}}^{\frac{1}{2}}u,\tilde{{\mathbf{T}}}^{\frac{1}{2}}v)_{L_{2}(\Omega)}=(\tilde{{\mathbf{T}}}^{\frac{1}{2}}\tilde{{\mathbf{T}}}^{-1}f,\tilde{{\mathbf{T}}}^{\frac{1}{2}}{{\mathbf{T}}}^{-1}g)_{L_{2}(\Omega)}=(f,{{\mathbf{T}}}^{-1}g)_{L_{2}(\Omega)}=({{\mathbf{T}}}^{-1}f,g)_{L_{2}(\Omega)},$

and, similarly, for
$({\mathbf{T}}^{\frac{1}{2}}u,{\mathbf{T}}^{\frac{1}{2}}v),$

$\displaystyle({{\mathbf{T}}}^{\frac{1}{2}}u,\tilde{{\mathbf{T}}}^{\frac{1}{2}}v)_{L_{2}(\Omega)}=({\mathbf{T}}^{\frac{1}{2}}\tilde{{\mathbf{T}}}^{-1}f,{\mathbf{T}}^{\frac{1}{2}}{{\mathbf{T}}}^{-1}g)_{L_{2}(\Omega)}=$
$\displaystyle({\mathbf{T}}^{\frac{1}{2}}\tilde{{\mathbf{T}}}^{-1}f,{\mathbf{T}}^{-\frac{1}{2}}g)_{L_{2}(\Omega)}=(\tilde{{\mathbf{T}}}^{-1}f,g)_{L_{2}(\Omega)}.$

∎

Using (5.2), we arrive at the representation of the difference of the
operators ${\mathbf{H}}_{a}$ and ${\mathbf{H}}_{\tilde{\mathfrak{a}}}$, see
(4.5), (4.10):

${\mathbf{H}}_{\mathfrak{a}}-{\mathbf{H}}_{\tilde{\mathfrak{a}}}\equiv\rho^{\frac{1}{2}}\gamma(\gamma({\mathbf{T}}^{-1}-\tilde{{\mathbf{T}}}{}^{-1}))^{*}\rho^{\frac{1}{2}}=\rho^{\frac{1}{2}}\gamma[\gamma{\mathbf{X}}^{*}{\mathbf{Y}}]^{*}\rho^{\frac{1}{2}}.$
(5.5)

We need a somewhat different representation for the operator in (5.5).

######  Lemma 5.2.

The following equality is valid

${\mathbf{H}}_{\mathfrak{a}}-{\mathbf{H}}_{\tilde{\mathfrak{a}}}=\rho^{\frac{1}{2}}(\gamma{\mathbf{Y}}^{*})(\gamma{\mathbf{X}}^{*})^{*}\rho^{\frac{1}{2}}.$
(5.6)

###### Proof.

The operator
$\gamma{\mathbf{X}}^{*}{\mathbf{Y}}=(\gamma{\mathbf{X}}^{*})({\mathbf{Y}})$ is
a product of two bounded operators, therefore,
$(\gamma{\mathbf{X}}^{*}{\mathbf{Y}})^{*}={\mathbf{Y}}^{*}(\gamma{\mathbf{X}}^{*})^{*}.$
We apply the operator $\gamma$ from the left, multiply by the bounded function
$\rho^{\frac{1}{2}},$ and obtain (5.6). ∎

We will use expression (5.6) to evaluate the singular numbers of the
difference ${\mathbf{H}}_{\mathfrak{a}}-{\mathbf{H}}_{\tilde{\mathfrak{a}}}.$
Since the function $\rho$ is bounded, the multiplication by
$\rho^{\frac{1}{2}}$ preserves spectral estimates, therefore, it suffices to
drop this weight in further estimates.

###  5.2. Spectral estimates for the operator $\gamma{\mathbf{X}}^{*}$

We represent the operator ${\mathbf{X}}=\nabla{\mathbf{T}}^{-1}$ in the
following way as a product of two bounded operators,

${\mathbf{X}}=(\nabla{\mathbf{T}}^{-\frac{1}{2}}){\mathbf{T}}^{-\frac{1}{2}},{\mathbf{X}}^{*}={\mathbf{T}}^{-\frac{1}{2}}(\nabla{\mathbf{T}}^{-\frac{1}{2}})^{*}.$

This gives

$\gamma{\mathbf{X}}^{*}=(\gamma{\mathbf{T}}^{-\frac{1}{2}})(\nabla{\mathbf{T}}^{-\frac{1}{2}})^{*}.$

The operator $(\nabla{\mathbf{T}}^{-\frac{1}{2}})$ is bounded in
$L_{2}(\Omega)$ since ${\mathbf{T}}^{-\frac{1}{2}}:L_{2}(\Omega)\to
H^{1}(\Omega)$ and $\nabla$ is bounded as acting from $H^{1}(\Omega)$ to
$L_{2}(\Omega).$ The norm of $\nabla{\mathbf{T}}^{-\frac{1}{2}}$ is controlled
by the ellipticity constant of ${\mathcal{L}}$ and does not depend on the
approximation $\tilde{\mathfrak{a}}$. Therefore, $s$-numbers of the operator
$\gamma{\mathbf{X}}^{*}$ are majorated by $s$-numbers of
$\gamma{\mathbf{T}}^{-\frac{1}{2}},$

${\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{X}}^{*})\leq
C{\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{T}}^{-\frac{1}{2}}).$

Next we have

$n(\lambda,\gamma{\mathbf{T}}^{-\frac{1}{2}})=n(\lambda^{2},(\gamma{\mathbf{T}}^{-\frac{1}{2}})^{*}(\gamma{\mathbf{T}}^{-\frac{1}{2}})).$

Therefore,

${\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{T}}^{-\frac{1}{2}})={\mathbf{n}}^{\sup}(d,(\gamma{\mathbf{T}}^{-\frac{1}{2}})^{*}(\gamma{\mathbf{T}}^{-\frac{1}{2}}))$
(5.7)

The expression on the right in (5.7) characterizes the counting function for
the eigenvalues of the ratio

$\frac{\|\gamma{\mathbf{T}}^{-\frac{1}{2}}f\|^{2}_{L_{2}(\Sigma)}}{\|f\|^{2}_{L_{2}(\Omega)}},\,f\in
L_{2}(\Omega).$ (5.8)

We set $f={\mathbf{T}}^{\frac{1}{2}}u$, $u\in H^{1}(\Omega)$ in (5.8) and
obtain the ratio

$\frac{\|\gamma
u\|^{2}_{L_{2}(\Sigma)}}{\|{\mathbf{T}}^{\frac{1}{2}}u\|_{L_{2}(\Omega)}}=\frac{\int_{\Sigma}|u(y^{\prime},1)|^{2}d\mu_{\Sigma}(y^{\prime})}{{\mathbf{a}}[u]}\leq\frac{C^{-1}\int_{\Sigma}|u(y^{\prime},1)|^{2}d\mu_{\Sigma}(y^{\prime})}{\|u\|^{2}_{H^{1}(\Omega)}},\,u\in
H^{1}(\Omega),$ (5.9)

since ${\mathbf{a}}[u]\geq C\|u\|^{2}_{H^{1}(\Omega)}$, with constant
determined by ellipticity constant of ${\mathcal{L}},$ the ratio (5.9) is
majorated by the spectral problem considered in Theorem 2.1. From this
theorem, taking into account (5.7), we obtain

${\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{T}}^{-\frac{1}{2}})\leq C,$ (5.10)

with constant $C$ depending only on the $L_{\infty}$ norm of the matrix
$\mathfrak{a}^{-1},$ therefore, returning to the original problem, only on the
ellipticity constant of the initial operator and the Lipshitz norm of the
function $\psi.$

##  6\. Spectral estimates for the operator $\gamma{\mathbf{Y}}^{*}$. Rough
estimates

In the study of the second term in the factorization (5.5), (5.6), namely,
$\gamma{\mathbf{Y}}[\mathfrak{b}]^{*}$, where we set
${\mathbf{Y}}[\mathfrak{b}]=\mathfrak{b}\nabla\tilde{{\mathbf{T}}}{}^{-1},$ we
will need to obtain an asymptotic singular numbers estimate for
${\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{Y}}^{*})=\limsup_{\lambda\to
0}{\lambda^{2d}n(\lambda,\gamma{\mathbf{Y}}^{*})}.$ This will be done in two
steps. First, in this section, for an _arbitrary_ matrix function
$\mathfrak{b}$, $L_{\infty}(\Sigma)$-continuous at $\Sigma\times\\{1\\},$ we
find an estimate of ${\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{Y}}^{*})$ in terms
of certain integral norm of $\mathfrak{b}(.,1).$ The constant in this estimate
will, unfortunately, depend, in an uncontrollable manner, on the approximating
matrix $\tilde{\mathfrak{a}}$. However, this estimate enables us to restrict
ourselves to smooth matrices $\mathfrak{b}$, using Lemma 2.2. Then, in the
next section, for a, _now smooth,_ matrix $\mathfrak{b}$, we can use the
pseudodifferential calculus to establish the _asymptotic estimate_ for
singular numbers of $\gamma{\mathbf{Y}}^{*},$ containing the integral norm of
$\mathfrak{b}$ but depending now only on the ellipticity constant and some
algebraic norm of the approximating matrix $\tilde{\mathfrak{a}}$ and its
inverse, or, what is the same, on these bounds for the initial matrix
$\mathfrak{a}$. Finally, we collect our estimates and establish the asymptotic
singular numbers bound for $\mathfrak{b}=\mathfrak{a}-\tilde{\mathfrak{a}}.$

The important point here that we have already got rid of the non-smoothness of
the coefficients of the operator, the latter stayed behind in the operator
${\mathbf{X}},$ therefore our reasoning uses essentially the smoothness of the
coefficients $\tilde{\mathfrak{a}}(y).$

As a preparation, we note that we can suppose that $\mathfrak{b}=0$ outside an
(arbitrarily small) neighborhood of the boundary $\Sigma$ of the domain
$\Omega.$ This observation has been used many times in papers on our topic,
including [39], [2], [41]. A simple explanation is that if ${\rm
dist}\,(\hbox{{\rm supp}}\,\mathfrak{b},\partial\Omega)>0$ then the integral
kernel of $\tilde{{\mathbf{T}}}{}^{-1}$ is a smooth function on the set
$\hbox{{\rm supp}}\,\mathfrak{b}\times\Sigma$ (which is separated from the
diagonal) and the eigenvalues of the corresponding operator decay faster than
any power,

${\mathbf{n}}^{\sup}(\theta,\gamma(\mathfrak{b}\nabla\tilde{{\mathbf{T}}}{}^{-1})^{*})=0,\,\mbox{for}\,\mbox{any}\,\theta>0,\,\hbox{{\rm
supp}}\,{\mathfrak{b}}\cap\Sigma=\varnothing.$ (6.1)

Therefore, we may suppose that $\mathfrak{b}$ is supported in a conveniently
small neighborhood of $\Sigma$ where it possesses the properties discussed in
Section 3.

######  Proposition 6.1.

Let $\tilde{\mathfrak{a}}(y),y\in\Omega,$ be a smooth elliptic matrix function
and $\mathfrak{b}(y)\in L_{\infty}(\Omega).$ Suppose that $\mathfrak{b}$ is
$L_{\infty}(\Sigma)$\- continuous at $\Sigma\times\\{1\\}$ and $\mathfrak{b}$
is zero outside some neighborhood of $\Sigma\times\\{1\\}$. Then

${\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{Y}}[\mathfrak{b}]^{*})\leq
C(\tilde{\mathfrak{a}})\|\mathfrak{b}(.,1)\|_{2d+2}^{\frac{d}{2}}\|\mathfrak{b}(.,1)\|_{\infty}^{\frac{d}{2}}.$
(6.2)

###### Proof.

We consider the operator ${\mathbf{Z}}[\mathfrak{b}]=\gamma{\mathbf{Y}}^{*}$,
where
${\mathbf{Y}}={\mathbf{Y}}[\mathfrak{b}]=\mathfrak{b}\nabla\tilde{{\mathbf{T}}}^{-1}$.
Since we are studying the singular numbers of
${\mathbf{Z}}[\mathfrak{b}]:L_{2}(\Omega)\to L_{2}(\Sigma)$, we can consider
the adjoint operator ${\mathbf{Z}}[\mathfrak{b}]^{*}:L_{2}(\Sigma)\to
L_{2}(\Omega)$ instead. This operator admits a convenient representation,
namely,

${\mathbf{Z}}[\mathfrak{b}]^{*}=(\gamma{\mathbf{Y}}^{*})^{*}=\left[\gamma(\mathfrak{b}\nabla\tilde{{\mathbf{T}}}^{-1})^{*}\right]^{*}=(\mathfrak{b}\nabla)(\gamma\tilde{{\mathbf{T}}}^{-1})^{*}.$
(6.3)

In fact, starting with the identity

$\mathfrak{b}\nabla\tilde{{\mathbf{T}}}^{-1}=(\mathfrak{b}\nabla\tilde{{\mathbf{T}}}^{-\frac{1}{2}})\tilde{{\mathbf{T}}}^{-\frac{1}{2}},$

we obtain

$(\mathfrak{b}\nabla\tilde{{\mathbf{T}}}^{-1})^{*}=\tilde{{\mathbf{T}}}^{-\frac{1}{2}}(\mathfrak{b}\nabla\tilde{{\mathbf{T}}}^{-\frac{1}{2}})^{*};$

further, we have

$\gamma(\mathfrak{b}\nabla\tilde{{\mathbf{T}}}^{-1})^{*}=\left(\gamma\tilde{{\mathbf{T}}}^{-\frac{1}{2}}\right)\left(\mathfrak{b}\nabla\tilde{{\mathbf{T}}}^{-\frac{1}{2}}\right)^{*}.$

And now we take adjoints,

$\displaystyle{\mathbf{Z}}(\mathfrak{b})^{*}\equiv\left[\gamma(\mathfrak{b}\nabla\tilde{{\mathbf{T}}}^{-1})^{*}\right]^{*}=$
$\displaystyle\left(\mathfrak{b}\nabla\tilde{{\mathbf{T}}}^{-\frac{1}{2}}\right)^{**}\left(\gamma\tilde{{\mathbf{T}}}^{-\frac{1}{2}}\right)^{*}=\mathfrak{b}\nabla\tilde{{\mathbf{T}}}^{-\frac{1}{2}}(\gamma\tilde{{\mathbf{T}}}^{-\frac{1}{2}})^{*}=\mathfrak{b}\nabla(\gamma\tilde{{\mathbf{T}}}^{-1})^{*},$

since $\mathfrak{b}$ is a bounded function. We consider now the function
$w(y)=(\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}g\in L_{2}(\Omega)$ for $g\in
L_{2}(\Sigma).$ By Lemma 4.5, this function satisfies the second order
elliptic equation ${\tilde{\mathcal{L}}}w=0$ in $\Omega.$ Its restriction to
$\Sigma$ is $\gamma w=\gamma(\gamma\tilde{{\mathbf{T}}}^{-1})^{*}g$, and,
since by Proposition 4.6 the operator
$\gamma(\gamma\tilde{{\mathbf{T}}}^{-1})^{*}$ is an order $-1$
pseudodifferential operator, the inequality holds

$\|\gamma w\|_{H^{1}(\Sigma)}\leq C\|g\|_{L_{2}(\Sigma)}.$ (6.4)

By the elliptic regularity property, for the solution $w$ of the second order
elliptic equation ${\tilde{\mathcal{L}}}w=0$ in $\Omega$ with the Dirichlet
boundary condition $\gamma w$, the Sobolev spaces estimate holds:

$\|w\|_{H^{s+\frac{1}{2}}(\Omega)}\asymp\|\gamma w\|_{H^{s}(\Sigma)},$ (6.5)

now for _all_ $s\in{\mathbb{R}}^{1}.$ Note here, that the constants concealed
in the ${}^{\prime}\asymp^{\prime}$ symbol in (6.5) depend on the particular
value of $s$ and, what is important, on the operator ${\tilde{\mathcal{L}}}$,
and may deteriorate when ${\tilde{\mathcal{L}}}$ is changing, while the
derivatives of the coefficients of ${\tilde{\mathcal{L}}}$ grow in the process
of approximation, even with the ellipticity constant preserved.

Now we pass to estimating the singular numbers of the operator
$(\gamma{\mathbf{Y}}^{*})^{*}$, whose _squares_ , by (6.3), are described by
the ratio

$\frac{\int_{\Omega}|\mathfrak{b}(y)\nabla(\gamma\tilde{{\mathbf{T}}}^{-1})^{*}g|^{2}dy}{\|g\|^{2}_{L_{2}(\Sigma)}}=\frac{\int_{\Omega}|\mathfrak{b}(y)\nabla
w(y)|^{2}dy}{\|g\|^{2}_{L_{2}(\Sigma)}},$ (6.6)

where, recall, $w=(\gamma\tilde{{\mathbf{T}}}^{-1})^{*}g$, $g\in
L_{2}(\Sigma).$

We use the relations (6.5) (with $s=1$) and (6.4). As a result, the
eigenvalues of the ratio (6.6) are majorized by the eigenvalues of the ratio

$\frac{\int_{\Omega}|\mathfrak{b}(y)\nabla
w(y)|^{2}dy}{\|w\|^{2}_{H^{3/2}(\Omega)}},\,{\tilde{\mathcal{L}}}w=0\,\mbox{{in}}\,\Omega.$
(6.7)

On the next step, while evaluating the eigenvalues of the ratio (6.7), we use
the following weighted estimate for solutions of elliptic equations, see [41],
Lemma 3.3.

######  Lemma 6.2.

Let a function $w\in H^{s}(\Omega),\,s>0,$ be a solution of an elliptic
equation ${\tilde{\mathcal{L}}}w=0$ with smooth coefficients in a bounded
domain $\Omega$ and let $r(y)$ be the distance from the point $y\in\Omega$ to
the boundary of $\Omega.$ Then, for any number $\kappa\geq 0$ such that
$s+\kappa$ is an integer,

$\int_{\Omega}r(y)^{2\kappa}|\nabla_{s+\kappa}w(y)|^{2}dy\leq
C\|w\|^{2}_{H^{s}(\Omega)},$ (6.8)

for a constant $C=C(s,\kappa,\Omega)$ not depending on $w.$

The constant in (6.8) depends, of course, on the operator
$\tilde{{\mathcal{L}}}$ as well, but this dependence is not mentioned in [41]
– although this fact is of no importance at the moment, as long as the
operator is fixed.

We also have the obvious inequality,

$\int_{\Omega}r(y)^{2\kappa}|\nabla w(y)|^{2}dy\leq
C\|w\|^{2}_{H^{s}(\Omega)},$

for $s\geq 1.$

We consider the case $d>1$ first, thus excluding temporarily the case of a
two-dimensional domain $\Omega.$ We apply Lemma 6.2 for the particular values
$s=\frac{3}{2},$ $\kappa=\frac{1}{2},$ $s+\kappa=2,$ and thus reduce our task
to estimating the eigenvalues of the ratio

$\frac{\int_{\Omega}|\mathfrak{b}(y)\nabla
w(y)|^{2}dy}{\int_{\Omega}r(y)(|\nabla_{2}w|^{2}+|\nabla
w|^{2}dy},{\tilde{\mathcal{L}}}w=0;$ (6.9)

where, in our case, $r(y)=1-y_{d+1}$ If we drop the condition
${\tilde{\mathcal{L}}}w=0$ in (6.9), the eigenvalue counting function may only
increase and the result will give us an upper eigenvalue estimate for the
ratio (6.9).

The denominator in (6.9) is the quadratic form of a degenerate elliptic
operator. Such kind of spectral problems was considered in the series of
papers of M.Z. Solomyak and I.L. Vulis in 1970-s, see, especially, [46], where
order sharp eigenvalue estimates and formulas for asymptotics have been
proved. We need only estimates; we might have cited the paper [41], where the
eigenvalue estimates for problems of the type (6.9) were established in an
even more general setting, namely, for a domain with _piecewise smooth_
boundary, see there Lemma 4.1. However, for our, more simple, case, since the
boundary is smooth, we refer to a more easily accessible and more elementary
paper [46]. In order to do this, we make some more transformations of our
spectral problem. Namely, we majorize the matrix $\mathfrak{b}(y)$ by its
matrix norm $b(y)=|\mathfrak{b}(y)|=(\sum|b_{j,k}(y)|^{2})^{\frac{1}{2}}.$
After this, we replace the gradient of $w$, $\nabla
w=(\partial_{1}w,\dots,\partial_{d+1}w)$ by an arbitrary vector function with
$d+1$ components ${\mathbf{w}}_{\iota}$. This widens the set of functions
where the variational ratio is considered, therefore, we arrive at the ratio

$\frac{\sum_{\iota}\int_{\Omega}b(y)^{2}|{\mathbf{w}}_{\iota}(y)|^{2}dy}{\sum_{\iota}\int_{\Omega}r(y)(|\nabla{\mathbf{w}}_{\iota}(y)|^{2}+|{\mathbf{w}}_{\iota}|^{2})dy}.$
(6.10)

In this way, the spectral problem (6.9) splits into the direct sum of $d+1$
identical scalar spectral problems,

$\frac{\int_{\Omega}b(y)^{2}|{\mathbf{w}}(y)|^{2}dy}{\int_{\Omega}r(y)(|\nabla{\mathbf{w}}(y)|^{2}+|{\mathbf{w}}|^{2})dy},$
(6.11)

exactly of the form, considered in [46]. The result which we are going to use
here is the combination of Lemma 5.1 and Lemma 5.4 in [46].

We cite and further discuss these results, tailored for our special case.
Recall that the boundary $\Sigma$ corresponds to $y_{d+1}=1,$
$r(y)=1-y_{d+1},$ and the corresponding change is made in the formulation.

Set $m=d+1.$ Let ${\mathbf{b}}[{\mathbf{w}}]$ be the quadratic form in the
cylinder $\sc\mbox{C}\hskip 1.0pt\subset\Omega\subset{\mathbb{R}}^{m}$:

${\mathbf{b}}[{\mathbf{w}}]=\int_{\sc\mbox{C}\hskip
1.0pt}r(y)^{\beta}b(y)^{2}{\mathbf{w}}(y)|^{2}dy,$ (6.12)

and

${\mathbf{a}}^{(\alpha)}[{\mathbf{w}}]=\int_{\Omega}r(y)^{\alpha}(|\nabla{\mathbf{w}}|^{2}+|{\mathbf{w}}|^{2})dy.$
(6.13)

######  Lemma 6.3.

(Lemma 5.1, a), in [46].) Let the coefficient $b(y^{\prime},y_{m}),$
$1-\delta\leq y_{m}\leq 1,$ not depend on $y_{m}$, namely,
$b(y^{\prime},y_{m})=\check{b}(y^{\prime}),\,y^{\prime}\in\Sigma,\check{b}\in
L_{\infty}(\Sigma).$ Suppose that ${\mathbf{k}}:=2-\alpha+\beta>\frac{2}{m}$
(this case is called the _strong degeneration_ case in [46]). Set
$\theta=\frac{m-1}{{\mathbf{k}}}(=\frac{d}{{\mathbf{k}}})$. Then for the
spectrum of the ratio

$\frac{{\mathbf{b}}[{\mathbf{w}}]}{{\mathbf{a}}^{(\alpha)}[{\mathbf{w}}]},{\mathbf{w}}\in
H^{1}(\Omega),$ (6.14)

the estimate holds

${\mathbf{n}}^{\sup}(\theta,\ref{VS3})\leq
C\|\check{b}\|_{L_{2m}(\Sigma)}^{\frac{m-1}{2}}\|\check{b}\|_{L_{\infty}(\Sigma)}^{\theta-\frac{m-1}{2}}.$
(6.15)

The crucial importance of this estimate is that it involves the integral norm
of $\check{b}.$

The proof of this Lemma in [46] is based upon the separation of variables in
the cylinder, thorough bookkeeping of the eigenvalues of the separated problem
and finally using the result on the asymptotics of eigenvalues of an elliptic
boundary problem with singular weight.

Note that in [46], the weight function in the form ${\mathbf{b}}$ is denoted
by $b(y)$, while it is $b(y)^{2}$ in (6.12); the corresponding change is made
in the formulation of Lemma 6.3.

The second lemma relaxes the condition
$b(y^{\prime},y_{m})=\check{b}(y^{\prime})$ used in Lemma 6.3. We formulate it
in our terms:

######  Lemma 6.4.

(Lemma 5.4 in [46].) Let all conditions of Lemma 6.3, except
$b=\check{b}(y^{\prime})$, be fulfilled, the latter being replaced by

$b(y^{\prime},y_{m})\,\mbox{is}\,L_{\infty}(\Sigma)\,\mbox{-continuous}\,\mbox{at}\,\Sigma\times\\{1\\}.$
(6.16)

Then the asymptotic estimate (6.15) holds, with $\check{b}=b(.,1)$.

In fact, a stronger statement is formulated in Lemma 5.4 in [46], namely, that
the asymptotic bounds for ${\mathbf{n}}_{\pm}$ for the functions $b$ and
$\check{b}$ coincide, but we need only the upper estimate (6.15).

In [46], only a short sketch of the proof of Lemma 5.4 is given, with
reference to ’standard tools of the variational method’ (which might, in fact,
have been considered standard at that glory period of the variational method,
but are, probably, not that standard nowadays). In more detail, and in much
more generality, the reasoning, explaining the passage from
$y_{m}$-independent coefficient $b$ to the one satisfying (6.16), is given in
[41], however, formally, the condition $b\in C(\Omega)$ was imposed there. In
fact, only the condition (6.16) was actually used in [41]. A short but,
hopefully, sufficient explanation is the following.

By the Ky Fan inequality, it suffices to prove that if (6.16) holds and
$b(.,1)=0$ then ${\mathbf{n}}^{\sup}(\theta,\eqref{VS3})=0$. Fix some
$\boldsymbol{\epsilon}>0$ and find $l>0$ such that
$|b(y)|<\boldsymbol{\epsilon}$ for $y_{m}\geq 1-l.$ The quadratic form
${{\mathbf{b}}[{\mathbf{w}}]}$ splits into the sum of two forms,

$\displaystyle{\mathbf{b}}[{\mathbf{w}}]=\int_{y_{m}\geq
1-l}b(y)^{2}|{\mathbf{w}}(y)|^{2}dy+\int_{y_{m}<1-l}b(y)^{2}|{\mathbf{w}}(y)|^{2}dy\leq$
(6.17) $\displaystyle
C\left(\boldsymbol{\epsilon}^{2}\int_{\Omega}|{\mathbf{w}}(y)|^{2}dy+l^{-\delta}\int_{\Omega}|y_{m}-1|^{\delta}|{\mathbf{w}}(y)|^{2}dy\right)\equiv{\mathbf{b}}_{\boldsymbol{\epsilon}}[{\mathbf{w}}]+{\mathbf{b}}_{l}[{\mathbf{w}}],\,\delta>0.$

For the first term in the second line in (6.17), Lemma 6.3 with
$\check{b}=\boldsymbol{\epsilon}$ applies by monotonicity, which gives the
estimate
${\mathbf{n}}_{\pm}^{\sup}(\theta,{\mathbf{b}}_{\boldsymbol{\epsilon}}/{\mathbf{a}})\leq
C\boldsymbol{\epsilon}^{2\theta}.$ For the eigenvalues of the operator
described by the second term, ${\mathbf{b}}_{l}[{\mathbf{w}}],$ we apply Lemma
6.3 with some $\beta=\delta>0,$ this means, with a different order in the
weight in the numerator (6.12). As a result, Lemma 6.3 gives for this term a
faster eigenvalues decay: (6.15) takes the form
${\mathbf{n}}_{\pm}^{\sup}(\theta^{\prime},{\mathbf{b}}_{\delta,L}/{\mathbf{a}})<\infty,$
$\theta^{\prime}=\frac{m-1}{1+\delta}<\theta,$ and therefore,
${\mathbf{n}}(\theta,{\mathbf{b}}_{\delta,l}/{\mathbf{a}})=0.$ After this, the
required equality ${\mathbf{n}}_{\pm}^{\sup}(\theta,\eqref{VS3})=0$ follows
due to the arbitrariness of $\boldsymbol{\epsilon}.$

We apply Lemma 6.4 for $m=d+1$, $\alpha=1$, $\beta=0,$
${\mathbf{k}}=2-\alpha+\beta=1$, $\theta=d$, to the ratio (6.11) and this
gives us the desired estimate for the singular numbers of
$\gamma{\mathbf{Y}}^{*}$

The above reasoning breaks down in the two-dimensional case, $m=d+1=2,$ since
here ${\mathbf{k}}=1$ and we have the equality ${\mathbf{k}}=\frac{2}{m}$
instead of the required inequality ${\mathbf{k}}>\frac{2}{m}.$ In this case,
called in [46] ’the intermediate degeneration’, the above scheme produces a
non-sharp order in the singular numbers estimate. This kind of complication
was handled in [41], Sect.4, in the following way which we adapt to our
situation. The idea is in using (6.8) for a different, larger, value of
$\kappa$, so that the resulting spectral problem becomes the one with strong
degeneration. The order of eigenvalue estimates obtained in this way does not
depend on the chosen value of $\kappa$. (Note that the dimension of the
enveloping space, denoted by $m$ in [46], is denoted by $m+1$ in [41],
therefore we change notations correspondingly when citing the latter paper.)

For $m=2,d=1,$ we choose the number $\kappa$ in (6.8) to be not $\frac{1}{2}$
but $\frac{3}{2},$ so, $s+\kappa=3$, and the weighted inequality (6.8) takes
the form

$\int_{\Omega}r(y)^{3}|\nabla_{3}w|^{2}dy\leq
C\|w\|^{2}_{H^{\frac{3}{2}}(\Omega)},\,{\tilde{\mathcal{L}}}w=0.$ (6.18)

Following the reasoning above, we arrive at estimating the singular numbers of
the ratio

$\frac{\int_{\Omega}|\mathfrak{b}(y)|^{2}|\nabla
w|^{2}dy}{\int_{\Omega}r(y)^{3}(|\nabla_{3}w(y)|^{2}+|w(y)|^{2})dy}.$ (6.19)

instead of (6.11). With this set of parameters, the problem (6.19) in
dimension $m=2$ is of a strong degeneration type, The result in [47] gives in
this case the estimate, see also [41], Lemma 4.1 or [42].

${\mathbf{n}}^{\sup}(1,\ref{W13})\leq
C(\mathfrak{a})\|\mathfrak{b}(.,1)\|_{L_{4}(\Sigma)}^{1/3}\|\mathfrak{b}(.,1)\|_{L_{\infty}(\Sigma)}^{2/3}.$
(6.20)

Formally, the proof of Lemma 4.1 in [41] requires $\mathfrak{b}$ to be
continuous, however the passage to the discontinuous $\mathfrak{b}$ which is
$L_{\infty}(\Sigma)$\- continuous at $\Sigma$ as function of $y_{d+1}$
variable is made identically with the above case of the dimension $d>1.$ A
more simple treatment of this case, for a smooth boundary, can be found in
[47].

Finally, we recall that the eigenvalues of the ratio (6.6) are squares of the
singular numbers of the operator ${\mathbf{Z}}[\mathfrak{b}],$ therefore, the
spectral estimate of order $d$ for (6.6) produces the singular numbers
estimate of order $2d$ for $\gamma{\mathbf{Y}}^{*}$ ∎

##  7\. Sharp estimates for ${\mathbf{Z}}=\gamma{\mathbf{Y}}^{*}$

We recall that the constant $C(\tilde{\mathfrak{a}})$ in (6.2) depends in a
non-controllable way on the approximating matrix $\tilde{\mathfrak{a}}$
(although, in an analogous situation, in [41] it was found that the constants
in this kind of eigenvalue estimates depend only on the bounds for some finite
collection of derivatives of $\tilde{\mathfrak{a}}.$) Now we, for the case of
a _smooth_ matrix $\mathfrak{b}(y),$ establish an estimate for
${\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{Y}}^{*})$ in terms of the norm of
$\mathfrak{b}$ in an integral metric and the bounds for the principal symbol
of the operator ${\tilde{\mathcal{L}}}.$

In this section we establish asymptotic bounds for singular numbers of the
operator ${\mathbf{Z}}$ with smooth matrix $\mathfrak{b}\in
C^{\infty}(\overline{\Omega})$, which depend only on the ellipticity constant
of the matrix $\tilde{\mathfrak{a}}$ (in other words, on the ellipticity
constant of the operator ${\mathcal{L}}$), its norm and certain integral norm
of the matrix $\mathfrak{b}(.,1)$.

###  7.1. The structure of ${\mathbf{Z}}^{*}{\mathbf{Z}}$

The squares of singular numbers of ${\mathbf{Z}}$ are eigenvalues of the self-
adjoint operator ${\mathbf{Z}}^{*}{\mathbf{Z}}$ acting in $L_{2}(\Sigma)$. We
discuss its structure in more detail now. We have

$\displaystyle{\mathbf{Z}}^{*}{\mathbf{Z}}=\left[(\mathfrak{b}(y)\nabla)(\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}\right]^{*}\left[\mathfrak{b}(y)\nabla(\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}\right]$
(7.1)
$\displaystyle=\gamma\tilde{{\mathbf{T}}}{}^{-1}[(\mathfrak{b}(y)\nabla)^{*}(\mathfrak{b}(y)\nabla)](\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}.$

Since all coefficient functions entering in ${\mathbf{Z}}$ are smooth,
operators composing (7.1) belong to the Boutet-de-Monvel algebra of
pseudodifferential operators (see, e.g., [21]). We consider separate terms
more closely.

The operator $\gamma\tilde{{\mathbf{T}}}{}^{-1}$ and its adjoint were
described in detail in Sect. 4.4 and 4.5. Since ${\mathbf{R}}$ is an integral
operator with Hermitian kernel $R(y,z)$ the operator in (4.18) is an integral
operator acting from $\Omega$ to $\Sigma$ with kernel

$R^{(1)}(y,z)=\gamma_{y}R(y,z)-N\\!\\!\\!D_{y}\gamma_{1,y}R(y,z),\,y\in\Sigma,\,z\in\Omega,$

where $N\\!\\!\\!D_{y}$, $\gamma_{y},$ $\gamma_{1,y}$ denote the N-to-D
operator $N\\!\\!\\!D,$ the trace $\gamma$ and the conormal derivative
$\gamma_{1}$ acting upon the $y$ variable at $y\in\Sigma,$
$\gamma_{1,y}=\gamma_{y}\partial_{\nu_{\tilde{\mathfrak{a}}}(y)}.$ Recall that
$N\\!\\!\\!D$ is an order $-1$ pseudodifferential operator at the boundary
with symbol $\beta(y^{\prime},\eta^{\prime}),$ see (1.4).

The adjoint operator $(\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}$ is, therefore,
an integral operator acting from $\Sigma$ to $\Omega$ with the adjoint kernel,

$R^{(1)*}(y,z)=R^{(1)}(z,y)=\gamma_{z}R(y,z)-N\\!\\!\\!D_{z}\gamma_{1,z}R(y,z),\,y\in\Omega,z\in\Sigma.$

Now we can describe how the operator
${\mathbf{Z}}=\mathfrak{b}(y)\nabla(\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}$
acts. By (7.1),
${\mathbf{Z}}=\mathfrak{b}(y)\nabla(\gamma\tilde{{\mathbf{T}}}{}^{-1})^{*}$ is
an integral operator acting from $\Sigma$ to $\Omega$ with integral kernel

${\mathcal{Z}}(y,z)=\mathfrak{b}(y)\nabla_{y}\gamma_{z}R(y,z)-\mathfrak{b}(y)\nabla_{y}N\\!\\!\\!D_{z}\gamma_{1,z}R(y,z),\,y\in\Omega,z\in\Sigma.$
(7.2)

Similarly, ${\mathbf{Z}}^{*}$ is an integral operator acting from $\Omega$ to
$\Sigma$ with the integral kernel

${\mathcal{Z}}^{*}(y,z)={\mathcal{Z}}(z,y).$ (7.3)

###  7.2. Composition

Now we collect the description of the entries in the operator
${\mathbf{W}}={\mathbf{Z}}^{*}{\mathbf{Z}}.$ By (7.2),(7.3),
${\mathbf{Z}}^{*}{\mathbf{Z}}$ is the composition

${\mathbf{Z}}^{*}{\mathbf{Z}}=((\mathfrak{b}\nabla)((\gamma-N\\!\\!\\!D\gamma_{1}){\mathbf{R}})^{*})^{*}(\mathfrak{b}\nabla)((\gamma-N\\!\\!\\!D\gamma_{1}){\mathbf{R}})^{*}$
(7.4)

is the integral operator with kernel

${\mathcal{W}}(y,z)=\int_{\Omega}{\mathcal{Z}}(y,\varsigma){\mathcal{Z}}(\varsigma,z)d\varsigma,$
(7.5)

with kernel ${\mathcal{Z}}$ given by (7.2), acting on $\Sigma.$

As elements in the Boutet-de-Montvel algebra,
$(\gamma-N\\!\\!\\!D\gamma_{1}){\mathbf{R}}$ is a trace operator acting from
$\Omega$ to $\Sigma,$ its adjoint
$((\gamma-N\\!\\!\\!D\gamma_{1}){\mathbf{R}})^{*}$ is a Poisson operator
acting from $\Sigma$ to $\Omega,$ and the whole composition, ${\mathbf{W}},$
has integral kernel of the form

$\displaystyle{\mathcal{W}}(y,z)=\int_{\Omega}{\mathcal{Z}}(y,\sigma){\mathcal{Z}}(\sigma,z)d\sigma=$
(7.6)
$\displaystyle\int_{\Omega}(\gamma_{y}-N\\!\\!\\!D_{y}(\gamma_{1})_{y})(\mathfrak{b}(\sigma)\nabla_{\sigma})^{*}R(y,\sigma)\mathfrak{b}(\sigma)\nabla_{\sigma}(\gamma_{z}-N\\!\\!\\!D_{z}(\gamma_{1})_{z})R(\sigma,z)d\sigma.$

By the Boutet de Monvel calculus, ${\mathbf{W}}$ is a pseudodifferential
operator of order $-1$ on $\Sigma.$

We need to find the dependence of the principal symbol of this operator on the
matrices $\mathfrak{b}(y)$ and $\mathfrak{a}(y),$ $y\in\Sigma.$ In this
evaluation, we may ignore terms of lower order, in particular, those appearing
when we commute the factors in the product in (7.6).

We start by considering the term
$(\mathfrak{b}(\sigma)\nabla_{\sigma})^{*}(\mathfrak{b}(\sigma)\nabla_{\sigma}).$
This is a pseudodifferential operator of order $2$ in $\Omega$ with principal
symbol satisfying

${\mathcal{X}}(\sigma,\varsigma)=\langle\mathfrak{b}(\sigma)\varsigma,\mathfrak{b}(\sigma)\varsigma\rangle=|\mathfrak{b}(\sigma)\varsigma|^{2}\leq|\mathfrak{b}(\sigma)|^{2}|\varsigma|^{2}$
(7.7)

Further on, the truncated self-adjoint pseudodifferential operator
${\mathbf{M}}_{\mathfrak{b}}=$
$(\mathfrak{b}\nabla{\mathbf{R}})^{*}(\mathfrak{b}\nabla{\mathbf{R}})$ of
order $-2$ has principal symbol

${\mathcal{M}}_{\mathfrak{b}}(\sigma,\varsigma)=\mathfrak{r}_{-2}(\sigma,\varsigma)^{2}|\mathfrak{b}(\sigma)\varsigma|^{2},$
(7.8)

where, recall,
$\mathfrak{r}_{-2}(\sigma,\varsigma)=\langle\tilde{\mathfrak{a}}(\sigma)\varsigma,\varsigma\rangle^{-1}$
is the principal symbol of the fundamental solution ${\mathbf{R}}$, so
${\mathbf{M}}_{\mathfrak{b}}$ is an operator of order $-2$. Next, we need to
make the restriction of ${\mathbf{M}}_{\mathfrak{b}}$ to the boundary,
$y\in\Sigma.$

The trace operators $\gamma-N\\!\\!\\!D\gamma_{1}$ acting from both sides upon
${\mathbf{M}}\equiv{\mathbf{M}}_{\mathfrak{b}}$ produce an order $-1$
pseudodifferential operator on the boundary with principal symbol, again,
containing the matrix $\mathfrak{b}^{*}\mathfrak{b}$ and a homogeneous symbol
depending algebraically on the matrix $\tilde{\mathfrak{a}}$ and its inverse
$\tilde{\mathfrak{a}}^{-1}.$ Its principal symbol can be calculated following
the Boutet de Monvel calculus rules, which involve only algebraical operations
with principal symbols. We however do not need to calculate this symbol or
write down its explicit expression. For our needs it is sufficient to describe
its properties, especially, on what data of our operator it depends. We
demonstrate it on one of the terms,
${\mathbf{M}}_{\mathfrak{b}}^{(\gamma)}:=\gamma(\gamma{\mathbf{M}}_{\mathfrak{b}})^{*}$,
which is the restriction of ${\mathbf{M}}_{\mathfrak{b}}$ to the boundary. Its
symbol can be evaluated as

$\displaystyle{\mathcal{M}}^{\gamma}(y,\eta^{\prime})=\frac{1}{2\pi}\int_{-\infty}^{\infty}{\mathcal{M}}_{\mathfrak{b}}(y,\eta^{\prime},\eta_{d+1})d\eta_{d+1}$
(7.9)
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}|\mathfrak{b}(y)\eta|^{2}(\tilde{\mathfrak{a}}(y)\eta,\eta)^{-4}d\eta_{d+1}\leq$
$\displaystyle|\mathfrak{b}(y)|^{2}|\tilde{\mathfrak{a}}(y)|^{-2}\frac{1}{2\pi}\int_{-\infty}^{\infty}(|\eta^{\prime}|^{2}+|\eta_{d+1}|^{2})^{-1}d\eta_{d+1}\leq
C|\tilde{\mathfrak{a}}(y)^{-2}||\mathfrak{b}(y)|^{2}|\eta^{\prime}|^{-1},\,y\in\Sigma,$

the expression (7.9) is calculated in local co-ordinates at the point
$y\in\Sigma$, where $y_{d+1}$ axis is directed along the normal to $\Sigma.$

In a similar way, other terms in the symbol of ${\mathbf{Z}}^{*}{\mathbf{Z}}$
can be estimated, using the representation (7.4), following the composition
rules in the Boutet de Monvel algebra, see, e.g., [21], Sect.9.5, 10.4. Here
we use the expression for the symbol of the conormal derivative,
$\imath\langle\tilde{\mathfrak{a}}(y)\eta,\nu\rangle$ and the symbol
$\beta(y,\eta^{\prime})$ of the $N\\!\\!\\!D$ operator, see (1.4). All these
terms are majorated by $|\mathfrak{b}(y)|^{2},$ with bounded dependence on
$|\tilde{\mathfrak{a}}(y)^{-1}|$ and $|\tilde{\mathfrak{a}}(y)|,$
$y\in\Sigma.$ Singular Green operators arising in the process of composition
of truncated pseudodifferential operators give no contribution to the
principal symbol of the composition

As a result we have

######  Proposition 7.1.

The principal symbol $\mathfrak{W}(y,\eta^{\prime})$, $(y,\eta^{\prime})\in
T^{*}\Sigma$ of the order $-1$ pseudodifferential operator
${\mathbf{W}}={\mathbf{Z}}^{*}{\mathbf{Z}}$ admits the estimate

$|\mathfrak{W}(y,\eta^{\prime})|\leq
C|\mathfrak{b}(y)|^{2}A(|\tilde{\mathfrak{a}}(y)^{-1}|,|\tilde{\mathfrak{a}}(y)|)|\eta^{\prime}|^{-1},$
(7.10)

with a function $A(s_{1},s_{2})$ bounded on bounded intervals separated from
zero.

Now we apply the asymptotic estimate for singular numbers of negative order
pseudodifferential operators. This formula was obtained by M.Birman and
M.Solomyak in [14], Theorem 1, for operators in a domain in the Euclidean
space and then carried over to manifolds in [15]. The proof in [14], II, a
very technical one, having been published in the, now quite obscure, Russian
journal, was almost unaccessible to Western researchers, although the result
was rather widely cited. Fortunately, quite recently, R.Ponge in [33], Sect.6,
proposed a rather soft proof of a special case of the main results in [14],
concerning operators with smooth symbols, which fits our needs. We arrive at
the following asymptotic estimate for the singular numbers of the operator
${\mathbf{Z}}.$

######  Proposition 7.2.

Let $\tilde{{\mathcal{L}}}$ be an elliptic operator with leading coefficients
matrix $\tilde{\mathfrak{a}}(y)\in C^{\infty}(\overline{\Omega}).$ Suppose
that for all $y\in\overline{\Omega},$

$|\tilde{\mathfrak{a}}(y)^{-1}|,|\tilde{\mathfrak{a}}(y)|\leq{\mathbf{C}}_{\mathfrak{a}}.$
(7.11)

Let $\mathfrak{b}(y)$ be a smooth symmetric matrix in
$C^{\infty}(\overline{\Omega}).$ Then for the singular numbers of the operator
${\mathbf{Z}}=\mathfrak{b}\nabla(\gamma\tilde{T}^{-1})^{*}$ the estimate holds

$\displaystyle{\mathbf{n}}^{\sup}(2d,{\mathbf{Z}})={\mathbf{n}}^{\sup}(d,{\mathbf{Z}}^{*}{\mathbf{Z}})\leq
A({\mathbf{C}}_{\mathfrak{a}})\int_{\Sigma}|\mathfrak{b}(y^{\prime})|^{2d}d\mu_{\Sigma}(y^{\prime})\leq$
(7.12) $\displaystyle
CA({\mathbf{C}}_{\mathfrak{a}})\|\mathfrak{b}\|_{L_{2d+2}(\Sigma)}^{d}\|\mathfrak{b}\|_{L_{\infty}(\Sigma)}^{d},$

with constant $A({\mathbf{C}}_{\mathfrak{a}})$ depending only on the
$L_{\infty}$ norms of $\tilde{\mathfrak{a}}^{-1}$ and $\tilde{\mathfrak{a}}$.

###### Proof.

As we just found, the operator ${\mathbf{Z}}^{*}{\mathbf{Z}}$ is an order $-1$
pseudodifferential operator on the boundary $\Sigma,$ with leading symbol
$\mathfrak{W}(y^{\prime},\eta^{\prime}),$
$(y^{\prime},\eta^{\prime})\in\mathrm{T}^{*}\Sigma$ majorated by
${\mathbf{C}}_{\mathfrak{a}}|\mathfrak{b}(y)|^{2}|\eta^{\prime}|^{-1}.$ By the
formula (23) in [14], for the eigenvalue counting function
$n(\lambda,{\mathbf{Z}}^{*}{\mathbf{Z}})$ the asymptotics holds (in our
notations)

${\mathbf{n}}(2d,{\mathbf{Z}})={\mathbf{n}}(d,{\mathbf{Z}}^{*}{\mathbf{Z}})=d^{-1}(2\pi)^{-d}\int_{\Sigma}\int_{\mathrm{S}^{*}\Sigma}\mathfrak{W}(y,\eta^{\prime})^{d}d\mu_{\Sigma}\omega(\eta^{\prime}),$
(7.13)

and $\omega$ is the standard volume form on the unit sphere $S^{d-1}.$ Thus,
our estimate for the principal symbol $\mathfrak{W}$ of the operator
${\mathbf{Z}}^{*}{\mathbf{Z}}$, substituted in (7.13), gives (7.12). The last
inequality in (7.12) follows from the fact that the $L_{\infty}$ norm of
$\mathfrak{b}$ majorates its $L_{p}$ norm, for any $p<\infty.$ ∎

###  7.3. Final estimates for $\gamma{\mathbf{Y}}^{*}$

We combine the results of the last two subsections in order to obtain the
final singular numbers estimate for ${\mathbf{Z}}=\gamma{\mathbf{Y}}^{*}.$

######  Proposition 7.3.

Let $\mathfrak{b}(y)$ be a matrix function satisfying the conditions of
Proposition 6.1 and $\tilde{\mathfrak{a}}(y)$ be a smooth matrix satisfying
$|\tilde{\mathfrak{a}}^{-1}(y)|,|\tilde{\mathfrak{a}}^{(}y)|\leq{\mathbf{C}}_{\tilde{\mathfrak{a}}}.$
Then for the singular numbers of the operator
$\gamma({\mathbf{Y}}[\mathfrak{b}])^{*}=\gamma(\mathfrak{b}\nabla(\tilde{{\mathbf{T}}}^{-1}))^{*}$
the singular numbers estimate holds

${\mathbf{n}}(2d,\gamma{\mathbf{Y}}[\mathfrak{b}]^{*})\leq
A({\mathbf{C}}_{\tilde{\mathfrak{a}}})\|\mathfrak{b}(.,1)\|_{L_{2d+2}(\Sigma)}^{{d}}.$
(7.14)

Note that, from the first glance, the statement of Propostion 7.3 coincides
with the one of Proposition 6.1. There is, however a critical improvement.
While the constant in (6.2) may depend, in an uncontrollable way, on the
matrix $\tilde{\mathfrak{a}}(y),$ the constant in (7.14) depends only on the
ellipticity constant $|\mathfrak{a}^{-1}|$ and the norm
$|\tilde{\mathfrak{a}}(y)|,$ but not on the matrix $\tilde{\mathfrak{a}}$
itself or its derivatives.

###### Proof.

For a given $\epsilon,$ we construct, following the procedure in Lemma 4.4,
the approximating _smooth_ matrix $\mathfrak{b}_{\epsilon}\in
C^{\infty}(\overline{\Omega})$ such that
$\|\mathfrak{b}_{\epsilon}\|_{L_{\infty}(\Sigma)}\leq
C_{0}\|\mathfrak{b}_{\epsilon}\|_{L_{\infty}(\Omega)}$ and

$\|\mathfrak{b}(.,1)-\mathfrak{b}_{\epsilon}(.,1)\|_{L_{2d+2}(\Sigma)}<\epsilon.$
(7.15)

By the Ky Fan inequality,
$n(\lambda_{1}+\lambda_{2},{\mathbf{K}}_{1}+{\mathbf{K}}_{2})\leq
n(\lambda_{1},{\mathbf{K}}_{1})+n(\lambda_{2},{\mathbf{K}}_{2}),$ it follows
in the usual way (see, e.g. [13], Sect.5), that

${\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{Y}}[\mathfrak{b}]^{*})^{(1+2d)^{-1}}\leq{\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{Y}}[\mathfrak{b}_{\epsilon}]^{*})^{(1+2d)^{-1}}+{\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{Y}}[\mathfrak{b}-\mathfrak{b}_{\epsilon}]^{*})^{(1+2d)^{-1}}.$
(7.16)

The last term in (7.16) tends to zero as $\epsilon\to 0$ by Proposition 6.1,
applied to $\mathfrak{b}-\mathfrak{b}_{\epsilon}$ in place of $\mathfrak{b}$.
Therefore, passing to the limit in (7.16), we obtain

${\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{Y}}[\mathfrak{b}]^{*})\leq\limsup_{\epsilon\to
0}{\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{Y}}[\mathfrak{b}_{\epsilon}]^{*}).$
(7.17)

Now we remember the estimate (7.12), which is valid, since
$\mathfrak{b}_{\epsilon}$ is smooth now. By the construction of the
approximating matrix $\mathfrak{b}_{\epsilon}$ in Lemma 4.4, the norm
$\|\mathfrak{b}_{\epsilon}(.,1)\|_{L_{2d+2}(\Sigma)}$ is controlled by the
same norm of $\mathfrak{b}$. Therefore we can pass to limit as $\epsilon\to 0$
in the inequality

${\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{Y}}[\mathfrak{b}]^{*})\leq
C\lim\sup_{\epsilon\to
0}\|\mathfrak{b}_{\epsilon}(.,1)\|^{d}_{L_{2d+2}(\Sigma)}\leq
C\|\mathfrak{b}(.,1)\|^{\frac{d}{2}}_{L_{2d+2}(\Sigma)},$ (7.18)

which gives the required estimate, since the norm
$\|\mathfrak{b}(.,1)\|_{L_{\infty}}\leq C\|\mathfrak{a}\|_{L_{\infty}}$ can be
absorbed in the coefficient $C.$ ∎

##  8\. Conclusion of the proof

The uniformity property of the estimate in (7.14) enables us to establish the
crucial approximation result for the operators
${\mathbf{H}}_{\tilde{\mathfrak{a}}}.$

######  Theorem 8.1.

Let $\mathfrak{a}$ be a matrix satisfying the conditions of Lemma 4.4, and for
a given $\epsilon>0$, $\tilde{\mathfrak{a}}=\tilde{\mathfrak{a}}_{\epsilon}$
be the approximating smooth matrix constructed in this Lemma. Then

${\mathbf{n}}^{\sup}(2d,\gamma({\mathbf{Y}}[\mathfrak{a}-\tilde{\mathfrak{a}}])^{*})\leq
C(\mathfrak{a})\epsilon^{d},$ (8.1)

and

$\lim_{\epsilon\to
0}{\mathbf{n}}^{\sup}(d,{\mathbf{H}}_{\mathfrak{a}}-{\mathbf{H}}_{\tilde{\mathfrak{a}}})\to
0.$ (8.2)

###### Proof.

We recall the representation
${\mathbf{H}}_{\mathfrak{a}}-{\mathbf{H}}_{\tilde{\mathfrak{a}}}=(\gamma({\mathbf{Y}}[\mathfrak{a}-\tilde{\mathfrak{a}}])^{*})^{*}(\gamma{\mathbf{X}}^{*}).$
Therefore, the second estimate will follow in the usual way, compare, e.g.,
[13], from the first one, by the Ky Fan inequality for the product of
operators,

$n(\lambda_{1}\lambda_{2},{\mathbf{K}}_{1}{\mathbf{K}}_{2})\leq
n(\lambda_{1},{\mathbf{K}}_{1})+n(\lambda_{2},{\mathbf{K}}_{2}),$ (8.3)

where
${\mathbf{K}}_{1}=\gamma{\mathbf{X}}^{*},\,{\mathbf{K}}_{2}=(\gamma({\mathbf{Y}}[\mathfrak{a}-\tilde{\mathfrak{a}}])^{*})^{*},$
since for ${\mathbf{K}}_{1}$ we already know the estimate
${\mathbf{n}}^{\sup}(2d,\gamma{\mathbf{X}}^{*})<\infty$. Namely, for a given
$\lambda>0,$ we set $\lambda_{1}=\epsilon^{-d/4}\lambda,$
$\lambda_{2}=\epsilon^{d/4}\lambda$ in (8.3),
$\lambda^{2}=\lambda_{1}\lambda_{2},$ then both terms on the right in (8.3)
get a small factor as $\epsilon\to 0.$

To prove (8.1), we apply Proposition 7.3 with
$\mathfrak{b}=\mathfrak{a}-\tilde{\mathfrak{a}}.$ ∎

Finally, using the basic asymptotic perturbation lemma, Lemma 2.2, we
establish our main result.

###### Proof.

Of Theorem 1.2. Let $\mathfrak{a}$ be the coefficient matrix in our P-S
problem in a smooth domain. Consider its smooth approximation
$\tilde{\mathfrak{a}}$ constructed according to Lemma 4.4, with $p=2d+2$. For
the operator ${\mathbf{H}}_{\tilde{\mathfrak{a}}}$ the asymptotic formula
(1.6) is known. For the difference,
${\mathbf{H}}_{\mathfrak{a}}-{\mathbf{H}}_{\tilde{\mathfrak{a}}}$ we have the
estimate (8.2). In these conditions, Lemma 2.2 grants that for the limit