Weyl asymptotics for Poincar\'{e}-Steklov eigenvalues in a domain with Lipschitz boundary

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

where ∂ ν is the derivative in the direction of the exterior normal ν(x) at x ∈ Σ.Further on, we assume that the boundary Σ is Lipschitz and we equip Σ with the surface measure induced by the Lebesgue measure in 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 • arding and Åke Plejel, with cooperation of Lars Hörmander, the author considered a moderately smooth domain Ω 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 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 ∞ cannot be approximated by smooth functions in L ∞ norm.We, however, construct an approximation of the coefficients in a certain L p -related norm, p < ∞, 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 ± (λ, K) we denote the distribution function of positive (negative) eigenvalues of the self-adjoint operator K; for an arbitrary compact operator K, n(λ, K) denotes the counting function of singular numbers of K.In the case when the operator is described by a quadratic form B[u] in the Hilbert space with norm defined by means of the quadratic form A[u], we replace K by the ratio of these quadratic forms in the notation n(λ, .),n ± (λ, .).For θ > 0, we denote by n sup ± (θ, K) the quantity lim sup λ→0 λ θ n ± (λ, K).Without the subscript ±, these notations concern the distribution of singular numbers of these operators.Without the superscript sup, these symbols denote the limits lim λ→0 λ θ n ± (λ, 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., n sup (θ, (X.Y )), etc.
We consider a bounded domain Ω ⊂ R d+1 with connected Lipschitz boundary Σ.(The connectedness condition can be easily removed, at the cost of certain notational complications.)This means that Σ = ∂Ω can be covered by a finite collection of co-ordinate neighborhoods U ι such that the portion of Σ in U ι can be, in a properly rotated Euclidean co-ordinate system, (x ′ , x d+1 ) ∈ U ι , represented by the equation x d+1 = ψ (ι) (x ′ ), x ′ ∈ V ⊂ R d with Lipschitz function ψ (ι) (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 L in divergence form, with the normal derivative ∂ ν replaced by the derivative along the conormal associated with 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 a (1. 2) The weighted N D problem for L is the eigenvalue problem with a real function ρ(x).The conormal differential operator ∂ a has symbol ı a(x)ξ, ν x .
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 Ω and 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 a(x) its sesquilinear form a a a x : a a a x (ξ) = j,k a jk (x)ξ j ξ k ; let a a a x (ξ, η) be the corresponding bilinear form.For x ∈ Σ, ξ ′ ∈ T * x Σ and the normal vector ν at x, we denote by this is a function positively homogeneous in ξ ′ of order 1.Then the density α ± (x) is defined as and the asymptotic formulas to be proved are where µ Σ is the surface measure on Σ generated by the embedding of Σ into R d+1 .For a continuous matrix-function 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 a is continuous only at the boundary Σ.In our case, the matrix 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 Σ, 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 Ω ⊂ R d+1 be a bounded domain; suppose that Σ = ∂Ω is a surface of class C 1 and ρ ∈ L ∞ (Σ).Suppose finally that the coefficients a jk are continuous in Ω.Then the asymptotics (1.6) holds.Moreover, the same conclusion is correct if Σ is Lipschitz but belongs to C 1 in a neighborhood of the support of ρ, outside a closed set of zero surface measure.The function ρ may belong to L d (Σ) for d > 1 and any L q (Σ) with q > 1, for d = 1.
The elementary formulation of the main result of our present paper is the following.
Theorem 1.2.Let Ω ⊂ R d+1 be a bounded domain with Lipschitz boundary Σ.Suppose that the operator L a is uniformly elliptic in Ω, the coefficients a jk are continuous in Ω. Suppose that ρ ∈ L d (Σ) for d > 1 and ρ belongs to the Orlicz class L log L(Σ) 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) ∈ L ∞ (Ω) is continuous at the boundary if there exists a continuous function a b ∈ C(Σ), Σ = ∂Ω such that for any where B(x ′ , δ) is the ball in R d+1 with radius δ, centered at x ′ .Since the functions in L ∞ (Ω) 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 Σ. Having this definition, the formulation of the more general result that we are going to prove here, is as follows.
Theorem 1.3.Let Ω ⊂ R d+1 be a bounded domain with Lipschitz boundary Σ.Suppose that the operator L a is uniformly elliptic in Ω, 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 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 L but with bounded coefficients, possibly, discontinuous at the boundary, in a domain with smooth boundary.This operator is approximated by operators L with smooth coefficients which converge to the coefficients of L in a certain L p -based norm, with some p < ∞.Then it is proved that the compact operators describing the P-S spectrum for L converge so strongly that it becomes possible to pass to the limit in the formulas for coefficients n(d, .) in the Steklov eigenvalue asymptotics for L.

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 v 0 (x) to the operator L, namely to consider the operator L+v 0 (x).Such perturbation is relatively compact with respect to self-adjoint realizations of 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 v 0 ≡ 1, but this does not simplify the matter.For economy of symbols, we use the notation L for the operator with v 0 already absorbed, further on.At a proper moment, we will use the freedom in choosing v 0 .
We denote by T ≡ T a the operator L with Neumann boundary conditions on Σ. Namely, it corresponds to the quadratic form with domain u ∈ H 1 (Ω); this form is equivalent to the square of the standard norm u H 1 (Ω) in the Sobolev space H 1 (Ω).Being considered in L 2 (Ω), this quadratic form defines the positive self-adjoint operator T satisfying 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 B[f ]  A[f ] with f ∈ X is a shorthand for the spectrum of the operator defined by the quadratic form B[f ] in the Hilbert space X equipped with the norm defined by the quadratic form A[f ].

2.2.
A freedom in the choice of the weight function ρ 0 .Estimates (2.3), (2.4) mean that the quadratic form (2.2) is bounded in H 1 (Ω), the equality in (2.2) can be therefore extended by continuity to all traces of functions in H 1 (Ω) and therefore the form ρ ρ ρ ρ ρ ρ ρ ρ ρ 0 defines a bounded operator S = S[ρ 0 ] in H 1 (Ω), with norm estimated by ρ 0 L d (Σ) , resp., ρ 0 L log L(Σ) .If the function ρ 0 is real-valued, this operator is self-adjoint.In fact, this operator is compact and its eigenvalues satisfy a power order estimate.(2.6) If the function ρ 0 is real-valued, estimates of the form (2.6) hold, separately, for positive and negative eigenvalues, (2.7) For d = 1 the above results hold, with L d (Σ) replaced by L log L(Σ).
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 Ω ′ containing Ω strictly inside and apply Theorem 3.3 in [38], for the particular choice of N = d + 1, l = 1, s = d, d > 1, µ being the surface measure on the surface Σ.This theorem states the following.Let Ω ′ ⊂ R N be a bounded open set, µ be a compactly supported measure in Ω ′ satisfying, for where B(X, r) is the ball with radius r centered at X.Then, for the particular 2), the estimates (2.6), (2.7) hold.We use this result in the following way.First, note that for a compact Lipschitz surface Σ of dimension d with µ being the d-dimensional Hausdorff measure on Σ, 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 Therefore, if on some subspace L ⊂ H 1 (Ω) of dimension n, we have By the variational principle, this implies the inequality for the counting functions of eigenvalues of operators, 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 (Ω), which we used, is valid under somewhat weaker restrictions than the Lipschitz property.Namely, in [25], the class of (ǫ − δ) -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 (ǫ − δ)-domain in 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 (ǫ − δ) domain, a bounded extension operator E for Sobolev spaces exists.Therefore, if the d-dimensional Hausdorff measure on ∂Ω satisfies (2.8) with s = d and Ω is an (ǫ − δ)-domain then for the operator S[ρ 0 , Ω] with a function ρ 0 ∈ L d (µ), which can be considered as a natural generalization of the N D operator, the estimate (2.6), (2.7) is valid.
The operator S ≡ S[ρ 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 K be a compact operator, θ > 0, and for sufficiently small ǫ one can split K into the sum, K = K ǫ + K ′ ǫ so that for K ǫ the singular numbers asymptotic formula holds, n(θ, K ǫ ) = A ǫ , (2.9) and the operators K ′ ǫ are asymptotically small in the sense as ǫ → 0. Then the limit lim ǫ→0 A ǫ = A exists and the asymptotics holds, An analogous statement is valid for the positive/negative eigenvalues of a selfadjoint operator K.One should replace in the above formulation n by n ± 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 S[ρ 0 ]; it follows from estimates (2.6), (2.7) Corollary 2.3.Suppose that for the weights ρ 0 in some set Y ⊂ L d (Σ), d > 1, the asymptotic formula (1.6) is established.Then this formula is valid for all ρ 0 ∈ Y, the closure of Y in the norm of L d (Σ).The same statement is valid for the dimension d = 1, with L d (Σ) replaced by the Orlicz space L log L(Σ).
In particular, this means that we can restrict ourselves to considering as Y, the set of continuous, and, later, smooth functions ρ 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 ρ 0 , we first approximate it in a proper integral norm by a function ρ which has 'separated' positive and negative parts, dist (supp ρ + , supp ρ − ) > 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 ρ 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 ∈ Σ, we can represent a a a x (ξ), a a a x (ξ, η) as matrix products (where ξ * denotes the row-vector, matrix-adjoint to the column-vector ξ.)In these notations, we have where Θ(x) = (ν * a(x)ν)a(x) − a(x)νν * a(x) (note that the last matrix here has rank 1, together with νν * .).Now let Π x be the embedding of T * x Σ into R d+1 , and Π * x be the adjoint operator, the projection of R d+1 to T * x Σ.Then (2.12) can be written as The meaning of the matrix Θ ′ x can be understood in the following way.We represent the (d + 1) × (d + 1) matrix Θ in an orthogonal frame where one of basis vectors is the normal vector ν x .Next, we strike out in this matrix the row and the column corresponding to ν x .What remains is just our matrix Θ ′ x .Using (2.14), we can establish certain estimates for the quantity n ± (d, 1.3), uniform in a class of elliptic operators.
It follows from the Cauchy-Schwartz inequality that

15)
For a positive matrix a, the exact equality in the Cauchy-Schwartz inequality (2.15) is possible only if the vectors ξ ′ and ν are parallel.The latter may never happen, therefore the inequality in (2.15) is strict; by homogeneity we obtain ) ) depends continuously on ξ ′ = 0; by the compactness of the unit sphere, we have where k(x) is the ellipticity constant of the matrix a(x) at the point x ∈ Σ.We suppose that the matrix a(x) is uniformly elliptic, a(x) (2.18) In the equivalent representation, the estimate for the matrix Θ ′ x also follows, det(Θ ′ x ) ≥ Ck −d 0 .The last estimates enable us to establish the following convergence property.
Lemma 2.4.Let a s (x), s ∈ [0, 1) be a family of Hermitian matrix functions on Σ, such that they satisfy ellipticity estimates a s (x) ≥ k 0 uniformly in x, s, and are also uniformly bounded, |a s (x)| ≤ k 1 .Suppose that the matrices a s converges as s → 0 to a 0 in L q (Σ), q < ∞.Denote by ̺ s (x) the function Proof.The statement follows easily from the fact that ̺ s (x) 2 is a rational function of the entries of the matrix 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.
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, (−∆ Σ ) 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 β(x, ξ ′ ) for the case of a smooth (or 'almost smooth') surface Σ.We discuss now the latter formula.Suppose first that the coefficients matrix a(x) does not depend on x ∈ Ω.Consider the fundamental solution R(x − y) for the operator L, so that L x R(x − y) = δ(x − y).With this fundamental solution the classical potential operators are associated, namely, the single layer potential operator S, and the 'direct value' of the conormal derivative of the single layer potential (often called the Neumann-Poincaré operator) the latter integral understood in the principal value sense.The adjoint to D ′ f is the integral operator D, called the direct value of the double layer potential, Then the N D operator is, up to lower order terms, the composition For a smooth boundary Σ, the operators S and D are order −1 pseudodifferential operators on Σ.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, in co-ordinates where x d+1 is directed along the normal vector ν x to the tangent plane at the point x ∈ Σ and the corresponding co-ordinates ξ ′ , ξ d+1 , with ξ ′ in the cotangent plane and ξ d+1 directed along ν x , see, e.g., [4], (4.28).Calculations by (2.20) give the expression for β(x, ξ ′ ) = s(x, ξ ′ ) −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), L x R(x, y) = δ(x − y) and in (2.20) one should change a(ξ ′ , ξ d+1 ) −1 to a(x ′ , 0, ξ ′ , ξ d+1 ) −1 , see [3], Sect.3.4.For a non-smooth, Lipschitz surface, the function β(x, ξ ′ ) 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.

Some geometry considerations
In this Section we introduce a global co-ordinate system in an interior collar neighborhood of the boundary Σ = ∂Ω.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 Σ ⊂ R d+1 has a tangent plane almost everywhere.This fact enables one to describe explicitly the surface measure on Σ generated by the Lebesgue measure on R d+1 .Namely, on the local Lipschitz graph where πE is the projection of the set E ⊂ Σ to R d .It is important to note that for a Lipschitz function ψ, the integrand in (3.1) is a bounded function, an algebraic expression of the partial derivatives of ψ.
A detailed study of geometric properties of Lipschitz domains is presented in [23].In particular, for a domain Ω ⊂ R d+1 with Lipschitz boundary, the distributional gradient grad (χ Ω ) of the characteristic function of Ω is a vector measure which is absolutely continuous with respect to the measure µ Σ in (3.1), 2) with the vector-valued density ν, |ν(x)| = 1 coinciding with the normal at x almost everywhere on Σ. (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 Ω to the ones in Lipschitz cylinders.In our study, it is more convenient to define the surface Σ 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 ∞ -vector field Γ Γ Γ(x), x ∈ R d+1 such that at the points x ∈ Σ = ∂Ω, those points where the tangent plane exists, |Γ Γ Γ(x)| = 1 and Γ Γ Γ(x) forms an acute angle with the normal ν(x) to Σ at x, not exceeding some constant which is strictly smaller than π/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 Ω ⊂ R d+1 be a bounded domain with Lipschitz boundary Σ.Let Γ Γ Γ be the transversal smooth vector field, as above.In a neighborhood of Σ denote by ℓ x the integral curve of this vector field passing through x.Then, for a certain interior collar neighborhood U of Σ there exists a smooth surface Σ Σ Σ Σ Σ Σ Σ Σ Σ ⊂ U, which has exactly one intersection point with every integral curve ℓ of the field Γ Γ Γ in U, and this intersection is transversal.
Proof.In some two-sided collar neighborhood U of Σ, the integral curves ℓ x ′ , x ′ ∈ Σ, of the vector field Γ Γ Γ form a one-dimensional foliation.We define the function f(x), x ∈ U, as the distance of the point x to the intersection point x ′ ∈ Σ along the integral curve ℓ x ′ passing through x (with proper sign, positive inside Σ.) Then x = (x ′ , t), t = f(x) form a new co-ordinate system in U.In these co-ordinates the function f is continuous in U, differentiable almost everywhere and has derivatives of any order in t variable, in other words, along ℓ x ′ .Moreover, f is strictly decaying in t variable, i.e., along the curves ℓ x ′ .
We take a smooth non-negative cut-off function χ ǫ (x), diam supp χ ǫ < ǫ, χ ǫ (x)dx = 1, with sufficiently small ǫ, and consider the convolution, the smooth function f ǫ = f * χ ǫ .Sufficiently close to Σ, this function is smooth, growing in t variable, and it has the same sign as f outside the ǫ-neighborhood of Σ.We consider the level surface Σ Σ Σ Σ Σ Σ Σ Σ Σ for f ǫ , f ǫ (x) = δ for some δ > 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 ℓ x ′ and is uniformly nontangent with these trajectories.Due to the monotonicity of f, the level δ can be chosen in such way that the surface Σ Σ Σ Σ Σ Σ Σ Σ Σ lies strictly inside Ω; we may also suppose that dist (Σ,Σ Σ Σ Σ Σ Σ Σ Σ Σ) > r 0 > 0.

Change of variables.
We introduce now new co-ordinates in Ω near Σ.Namely, for the smooth surface Σ Σ Σ Σ Σ Σ Σ Σ Σ, just constructed, for a point x ∈ U, we take now x = (x ′ , x d+1 ), where According to the properties of the surface Σ Σ Σ Σ Σ Σ Σ Σ Σ, this is a smooth co-ordinates change from the initial Euclidean co-ordinates.
In these co-ordinates, the initial surface Σ is described globally by the equation We denote by C 0 the truncated cylinder with base Σ Σ Σ Σ Σ Σ Σ Σ Σ and top Σ: which is diffeomorphic to a collar neighborhood of Σ in Ω, and denote by Ω 0 the complement of this cylinder in Ω, In this setting the domain Ω can be considered as a smooth manifold with Lipschitz boundary Σ, Our next aim now is to perform a change of variables Φ Φ Φ of Ω to the standard This change of variables is the identity mapping, Φ Φ Φ(x) = x on Ω 0 and an adjoining part of the cylinder C 0 [ψ ψ ψ], 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 Γ Γ Γ namely, with the mapping so, at the top, Σ, we have y d+1 = 1.Thus, in new co-ordinates, the Lipschitz boundary Σ of the domain Ω is transformed to the smooth boundary of the domain Ω Ω Ω: y ′ ∈ Σ Σ Σ Σ Σ Σ Σ Σ Σ, y d+1 = 1.We denote by x = Ψ Ψ Ψ(y) the mapping inverse to Φ Φ Φ (it, obviously, exists and it is the inverse stretching along the y d+1 variable near Σ Σ Σ Σ Σ Σ Σ Σ Σ × {1}).The mapping Φ Φ Φ is is a bilipschitz mapping of Ω to Ω Ω Ω.Its restriction to Σ is a bilipschitz mapping to Σ Σ Σ Σ Σ Σ Σ Σ Σ.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.
(1) The components of the mappings Φ Φ Φ, Ψ Ψ Ψ are continuous, their first order derivatives in x ′ , resp., y ′ 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 The Jacobian matrices of the mappings Ψ Ψ Ψ and Φ Φ Φ are bounded; their entries are continuous in y d+1 , resp., x d+1 , variables as functions with values in All these properties follow automatically from the definition of the mapping Ψ Ψ Ψ and the chain rule for derivatives, taking into account that the first order derivatives of a Lipschitz function are bounded.

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 Ω with cylindrical exit, as presented in the previous section.Having the quadratic forms ρ ρ ρ ρ ρ ρ ρ ρ ρ 0 , a 0 , we find out what happens with them under the transformation Ψ Ψ Ψ.
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 = Φ Φ Φ(x), First, we consider the quadratic form where J(y ′ , 1) is the Jacobian matrix of the mapping Ψ Ψ Ψ and . This quadratic form can be written as The quadratic form a 0 [u], see (2.1), is transformed to v ∈ H 1 (Ω Ω Ω).The transformed coefficients ǎj,k (y), 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 Σ in the sense of (1.7).This property persists after the continuous change of variables x = Ψ Ψ Ψ(y), so the functions a j,k (Ψ Ψ Ψ(y)) are continuous at the new, smooth, boundary Σ Σ Σ Σ Σ Σ Σ Σ Σ in the sense of (1.7).Next, to obtain ǎj,k (y), we need to multiply a j,k (Ψ Ψ Ψ(y)) by some algebraic combinations of derivatives of the mapping Ψ Ψ Ψ, namely, the derivatives of Ψ Ψ Ψ, appearing in the chain rule when passing from x-derivatives to y-derivatives and also by the Jacobian of the mapping Ψ Ψ Ψ 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 (Ψ Ψ Ψ(y)) for y ∈ Σ Σ Σ Σ Σ Σ Σ Σ Σ.However, these derivatives are still continuous in y d+1 variable as functions with values in L ∞ (Σ Σ Σ Σ Σ Σ Σ Σ Σ).We introduce the following definition.
Definition 4.1.Let q(y) = q(y ′ , y d+1 ) be a bounded measurable function in the variables Denote by q ⋄ (y) ≡ q ⋄ (y ′ , y d+1 ) the continuation of q(y ′ ) as a function not depending on y d+1 , q ⋄ (y ′ , y d+1 ) := q(y ′ ), we say that q(y) has If this is the case, we redefine the initial function q(y) on Σ Σ Σ Σ Σ Σ Σ Σ Σ, by setting q(y ′ , 1) = q(y ′ ), thus changing q on a set of d + 1-dimensional measure 0. After this, we say that One can easily see that the product of a function This 'very week' continuity property proves to be sufficient for our approximation construction.
We recall that the ratio of quadratic forms ρ ρ ρ ρ ρ ρ ρ ρ ρ 0 /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 (Ω Ω Ω).Therefore, our transformations can be summarized in the following way.
in a domain Ω Ω Ω with smooth boundary Σ Σ Σ Σ Σ Σ Σ Σ Σ, and where the quadratic form a[v] is elliptic and has bounded coefficients, which are 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 Ω and its, now smooth, boundary by Σ (= Σ × {1} in our y coordinates).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 ∞ (Σ) continuity.In the notations of the quadratic forms ρ ρ ρ ρ ρ ρ ρ ρ ρ 0 , 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 ρ.
The quadratic form a[v] with domain H 1 (Ω) defines the self-adjoint positive 'Neumann' operator T = T a in L 2 (Ω).This elliptic operator acts, formally, as 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 T is a positive self-adjoint elliptic operator with bounded coefficients, the domain of the square root T 1 2 can be described explicitly: it coincides with H 1 (Ω) and It was already stated that we are free in the choice of the function v. Now we fix this choice, supposing that v ∈ C ∞ 0 (Ω) (after the above co-ordinates change).Since a and therefore the ratio (4.3) transforms to Here and further on, we denote by γ : f → γf the trace operator acting as the restriction to Σ of a function f defined on the domain Ω with cylindric exit C having the boundary lid Σ identified with Σ×{1}.We recall that the boundary Σ is smooth.Therefore, this restriction is known to act as a bounded operator from H s (Ω) to H s−1/2 (Σ) for all positive s except the half-integer ones (s ∈ N + 1 2 ).In particular, this implies that the operator γT − 1 2 is bounded as acting from L 2 (Ω) to L 2 (Σ) (and even to H 1 2 (Σ)).Following Sect.2, we restrict ourselves to ρ ≥ 0 and ρ ∈ L ∞ , therefore, the operator ρ As a result, the numerator in (4.4) can be written as .
In this way, we have reduced our spectral problem to the study of the spectrum of the self-adjoint operator acting in L 2 (Ω).Since the nonzero discrete spectrum of the product does not change under the cyclic permutation of the factors, the nonzero eigenvalues of G a are the same as the ones of the self-adjoint operator acting in L 2 (Σ).
It is more convenient to have a somewhat different representation of the operator H a .
Lemma 4.3.The following equality holds Proof.We start with the equality [ρ where we have the product of two bounded operators on the left.Next we take adjoint of both parts in (4.7), This is a bounded operator acting from L 2 (Σ) to H 1 (Ω), again a product of two bounded operators.Therefore, the restriction operator γ is well defined on the range of the operator in (4.8) and we can apply γ on the left, and then multiply by the bounded function ρ 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 a.The condition (3) means that the L ∞ norm for the approximating matrix can be arranged to be not depending on ǫ, being controlled, by the L ∞ norm of the initial matrix 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 Σ.
In the case when the matrix a is not supposed to be uniformly positive definite, the simplified reasoning succeeds: namely we just skip all estimates concerning a −1 , therefore, we may simply set a 1 (y) = 0.
We will use Lemma 4.4 twice; first time when approximating the elliptic coefficients matrix a(y), where the control of the ellipticity constant is important, and the second time for the approximation of the matrix b = a − ã when proving the required operator convergence of H ã to H a .4.4.The approximating operator.Consider, for a given ǫ > 0, the approximating coefficient matrix ã constructed in Section 4.3 for the elliptic coefficient matrix a and the corresponding elliptic differential operator L = L ã, with the same zero order term v ∈ C ∞ (Ω).With the coefficient matrix ã we associate the Neumann operator T = T ã, by means of the quadratic form ã ã ã[v], v ∈ H 1 (Ω).
Since the boundary of Ω is smooth, as well as the coefficients matrix ã, some more can be said about the properties of the operator (4.5) for this matrix.The corresponding compact operator, which we denote by H ã, has, by Lemma 4.3, the form 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 ã replacing a, due to the existing results, see Theorem 1.1.Next, since the coefficients of T ã are smooth, the domain and action of the operator T ã can be described explicitly.Namely, by the standard elliptic regularity results, the domain D( T) of T is the subspace in the Sobolev space H 2 (Ω), consisting of functions v satisfying the classical Neumann boundary conditions, k ãj,k ν k ∂ j v| ∂Ω = 0, where ν k are the components of the unit normal vector to Σ = ∂Ω.
Lemma 4.5.For a smooth matrix ã, for any f ∈ L 2 (Σ), the function v = (ρ Proof.Let g ∈ C ∞ 0 (Ω) be a smooth function with compact support in Ω.Consider the sesquilinear form I[f, g] = ( L(ρ . The function g belongs to the domain of the self-adjoint operator T, therefore, (4.12) Since g = 0 on Σ, the last expression in (4.12) is zero.Therefore, the function L(ρ , and this implies that L(ρ 4.5.H a as a pseudodifferential operator.Here we establish an important property of the approximating operator H ã, needed further for evaluating the rate of approximation of spectra.Proposition 4.6.For a smooth approximating matrix ã and a smooth weight ρ, the operator H ã is an order −1 pseudodifferential operator on the boundary Σ. Proof.When finding, in the smooth case, an expression for the operator H ã, we study its component (ρ 1 2 γ T−1 ) * first.Recall that T−1 is the inverse of the realization of the operator L with Neumann boundary conditions, γ 1 u = 0, where γ 1 is the conormal derivative.Continue the coefficients ã outside Ω in a smooth way.We denote by R the fundamental solution for L; so that LR − 1, R L − 1 are infinitely smoothing operators; R is the order −2 pseudodifferential operator with symbol r(y, η) having the principal term r −2 (y, η) = ã ã ãy (η) −1 , where ã ã ãy (η) = j,k ãj,k (y)η j η k is the principal symbol of the operator L. Let d > 1 first.Then the operator R acts on functions on Ω 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 r −2 (y, η), 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 (Ω).In order to find a representation for u = T−1 f, we set u 0 = Ω R(y, y − z)f (z)dz.This function satisfies the equation Lu 0 = f in Ω but not the Neumann boundary condition γ 1 u = 0, where γ 1 is the conormal derivative at the boundary, corresponding to the elliptic operator L. We construct the correction u 1 , such that Lu 1 = 0, γ 1 u 1 = −γ 1 u 0 .To do this, we consider the The operator , it gives the identity operator.Therefore the solution of the boundary problem Lu = f in Ω, γ 1 u = 0 on Σ can be expressed as As a result, the operator T−1 can be represented as 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 R is treated as a truncated pseudodifferential operator, where e Ω is the operator of extension by zero of a function in Ω to the whole R d+1 , R is a pseudodifferential operator in R d+1 with symbol r(y, η) and χ Ω is the restriction of functions in R d+1 to Ω.We also need an expression for the operator T−1 0 , the resolvent of the Dirichlet problem for the equation Lu = f in Ω.In a similar way, where G G G 0 is the Poisson operator (the Green function) solving the nonhomogeneous Dirichlet boundary problem for L in Ω: Lu = 0, γu = h; this means that γG G G 0 = 1 1 1.
We pass now to γ T−1 .Since γ T0 −1 = 0, we can write We recall now that γG G G 0 = 1 1 1.Further on, the function u 1 = G G G 1 γ 1 Rf satisfies the equation Lu 1 = 0 in Ω, therefore γu 1 and γ 1 u 1 are connected by the Neumann-to Dirichlet operator N D, γu 1 = N Dγ 1 u 1 .We set all this into (4.17) and obtain Consequently, Finally, we have 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, R is a truncated pseudodifferential operator of order −2, γ is a trace operator of zero order, further, (γR) * is a Poisson operator, and as a result, γ(γR) * is a pseudodifferential operator of order −1.Similarly, γ 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 R and N D, as we will see in Sect.7.Finally, the multiplication by the smooth function ρ 1 2 produces, again, a pseudodifferential operator.

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 a(y) is continuous, it can be approximated in C(Ω), both from above and from below, by smooth elliptic matrices ã± (y), and the required closeness of spectra of operators T and 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 ψ 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.

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 a, ã be the matrices of coefficients of the operators L, L, described in Sect.4,so that a, ã−1 belong to L ∞ (Ω), ã, ã−1 ∈ C ∞ (Ω) and ã| Σ − a| Σ is small in the sense of Lemma 4.4.
Consider the operators T, T, the Neumann operators for L, L.
Lemma 5.1.Under the above conditions, the following factorization is valid: where X, Y are bounded operators acting from ) 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 (5.3) where the angle brackets denote the scalar product in C d+1 .
In our conditions, the domains of the operators T 1 2 , T1 2 coincide, both are equal to the Sobolev space H 1 (Ω).It follows, in particular, that We set here u = T−1 f, v = T −1 g, where f, g are arbitrary elements in L 2 (Ω).These functions u, v belong to H 1 (Ω), therefore (5.4) is satisfied.Such substitution leads to (5.3): and, similarly, for (T Using (5.2), we arrive at the representation of the difference of the operators H a and H ã, see (4.5), (4.10): ( We need a somewhat different representation for the operator in (5.5).
Lemma 5.2.The following equality is valid , and obtain (5.6).We will use expression (5.6) to evaluate the singular numbers of the difference H a − H ã. Since the function ρ is bounded, the multiplication by ρ 1 2 preserves spectral estimates, therefore, it suffices to drop this weight in further estimates.5.2.Spectral estimates for the operator γX * .We represent the operator X = ∇T −1 in the following way as a product of two bounded operators, 2 is controlled by the ellipticity constant of L and does not depend on the approximation ã.Therefore, s-numbers of the operator γX * are majorated by s-numbers of γT − 1 2 , n sup (2d, γX * ) ≤ Cn sup (2d, γT − 1 2 ).Next we have n(λ, γT 7) The expression on the right in (5.7) characterizes the counting function for the eigenvalues of the ratio We set f = T ) and obtain the ratio since a[u] ≥ C u 2 H 1 (Ω) , with constant determined by ellipticity constant of 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 with constant C depending only on the L ∞ norm of the matrix a −1 , therefore, returning to the original problem, only on the ellipticity constant of the initial operator and the Lipshitz norm of the function ψ.
6. Spectral estimates for the operator γY * .Rough estimates In the study of the second term in the factorization (5.5), (5.6), namely, γY[b] * , where we set Y[b] = b∇ T−1 , we will need to obtain an asymptotic singular numbers estimate for n sup (2d, γY * ) = lim sup λ→0 λ 2d n(λ, γY * ).This will be done in two steps.First, in this section, for an arbitrary matrix function b, L ∞ (Σ)continuous at Σ × {1}, we find an estimate of n sup (2d, γY * ) in terms of certain integral norm of b(., 1).The constant in this estimate will, unfortunately, depend, in an uncontrollable manner, on the approximating matrix ã.However, this estimate enables us to restrict ourselves to smooth matrices b, using Lemma 2.2.Then, in the next section, for a, now smooth, matrix b, we can use the pseudodifferential calculus to establish the asymptotic estimate for singular numbers of γY * , containing the integral norm of b but depending now only on the ellipticity constant and some algebraic norm of the approximating matrix ã and its inverse, or, what is the same, on these bounds for the initial matrix a. Finally, we collect our estimates and establish the asymptotic singular numbers bound for b = 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 X, therefore our reasoning uses essentially the smoothness of the coefficients ã(y).
As a preparation, we note that we can suppose that b = 0 outside an (arbitrarily small) neighborhood of the boundary Σ of the domain Ω.This observation has been used many times in papers on our topic, including [39], [2], [41].A simple explanation is that if dist (supp b, ∂Ω) > 0 then the integral kernel of T−1 is a smooth function on the set supp b × Σ (which is separated from the diagonal) and the eigenvalues of the corresponding operator decay faster than any power, Therefore, we may suppose that b is supported in a conveniently small neighborhood of Σ where it possesses the properties discussed in Section 3.
Proposition 6.1.Let ã(y), y ∈ Ω, be a smooth elliptic matrix function and b(y In fact, starting with the identity .
And now we take adjoints, since b is a bounded function.We consider now the function w(y) = (γ T−1 ) * g ∈ L 2 (Ω) for g ∈ L 2 (Σ).By Lemma 4.5, this function satisfies the second order elliptic equation Lw = 0 in Ω.Its restriction to Σ is γw = γ(γ T−1 ) * g, and, since by Proposition 4.6 the operator γ(γ T−1 ) * is an order −1 pseudodifferential operator, the inequality holds By the elliptic regularity property, for the solution w of the second order elliptic equation Lw = 0 in Ω with the Dirichlet boundary condition γw, the Sobolev spaces estimate holds: now for all s ∈ R 1 .Note here, that the constants concealed in the ′ ≍ ′ symbol in (6.5) depend on the particular value of s and, what is important, on the operator L, and may deteriorate when L is changing, while the derivatives of the coefficients of L grow in the process of approximation, even with the ellipticity constant preserved.Now we pass to estimating the singular numbers of the operator (γY * ) * , whose squares, by (6.3), are described by the ratio where, recall, w = (γ T−1 ) * g, g ∈ L 2 (Σ).
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 , Lw = 0 in Ω. ( 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 ∈ H s (Ω), s > 0, be a solution of an elliptic equation Lw = 0 with smooth coefficients in a bounded domain Ω and let r(y) be the distance from the point y ∈ Ω to the boundary of Ω.Then, for any number κ ≥ 0 such that s + κ is an integer, for a constant C = C(s, κ, Ω) not depending on w.
The constant in (6.8) depends, of course, on the operator 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, for s ≥ 1.
We consider the case d > 1 first, thus excluding temporarily the case of a two-dimensional domain Ω.We apply Lemma 6.2 for the particular values s = 3  2 , κ = 1 2 , s + κ = 2, and thus reduce our task to estimating the eigenvalues of the ratio Ω |b(y)∇w(y)| 2 dy Ω r(y)(|∇ 2 w| 2 + |∇w| 2 dy , Lw = 0; (6.9) where, in our case, r(y) = 1 − y d+1 If we drop the condition Lw = 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 b(y) by its matrix norm b(y) = |b(y)| = ( |b j,k (y)| 2 ) 1 2 .After this, we replace the gradient of w, ∇w = (∂ 1 w, . . ., ∂ d+1 w) by an arbitrary vector function with d + 1 components w ι .This widens the set of functions where the variational ratio is considered, therefore, we arrive at the ratio .10)In this way, the spectral problem (6.9) splits into the direct sum of d + 1 identical scalar spectral problems, 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 Σ corresponds to y d+1 = 1, r(y) = 1 − y d+1 , and the corresponding change is made in the formulation.
Set , w ∈ H 1 (Ω), (6.14) the estimate holds The crucial importance of this estimate is that it involves the integral norm of 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 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.
In fact, a stronger statement is formulated in Lemma 5.4 in [46], namely, that the asymptotic bounds for n ± for the functions b and 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 ∈ C(Ω) was imposed there.In fact, only the condition (6.16) was actually used in [41].A short but, hopefully, sufficient explanation is the following.
We apply Lemma 6.4 for (6.11) and this gives us the desired estimate for the singular numbers of γY * The above reasoning breaks down in the two-dimensional case, m = d + 1 = 2, since here k = 1 and we have the equality k = 2 m instead of the required inequality k > 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 κ, 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 κ. (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 κ in (6.8) to be not 1 2 but 3 2 , so, s + κ = 3, and the weighted inequality (6.8) takes the form 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].
L 4 (Σ) b(., 1) Formally, the proof of Lemma 4.1 in [41] requires b to be continuous, however the passage to the discontinuous b which is L ∞ (Σ)-continuous at Σ 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 Z[b], therefore, the spectral estimate of order d for (6.6) produces the singular numbers estimate of order 2d for γY * 7. Sharp estimates for Z = γY * We recall that the constant C(ã) in (6.2) depends in a non-controllable way on the approximating matrix ã (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 ã.) Now we, for the case of a smooth matrix b(y), establish an estimate for n sup (2d, γY * ) in terms of the norm of b in an integral metric and the bounds for the principal symbol of the operator L.
In this section we establish asymptotic bounds for singular numbers of the operator Z with smooth matrix b ∈ C ∞ (Ω), which depend only on the ellipticity constant of the matrix ã (in other words, on the ellipticity constant of the operator L), its norm and certain integral norm of the matrix b(., 1).7.1.The structure of Z * Z.The squares of singular numbers of Z are eigenvalues of the self-adjoint operator Z * Z acting in L 2 (Σ).We discuss its structure in more detail now.We have Since all coefficient functions entering in 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 γ T−1 and its adjoint were described in detail in Sect.4.4 and 4.5.Since R is an integral operator with Hermitian kernel R(y, z) the operator in (4.18) is an integral operator acting from Ω to Σ with kernel where N D y , γ y , γ 1,y denote the N-to-D operator N D, the trace γ and the conormal derivative γ 1 acting upon the y variable at y ∈ Σ, γ 1,y = γ y ∂ ν ã(y) .Recall that N D is an order −1 pseudodifferential operator at the boundary with symbol β(y ′ , η ′ ), see (1.4).
The adjoint operator (γ T−1 ) * is, therefore, an integral operator acting from Σ to Ω with the adjoint kernel, is the integral operator with kernel with kernel Z given by (7.2), acting on Σ.
As elements in the Boutet-de-Montvel algebra, (γ −N Dγ 1 )R is a trace operator acting from Ω to Σ, its adjoint ((γ − N Dγ 1 )R) * is a Poisson operator acting from Σ to Ω, and the whole composition, W, has integral kernel of the form By the Boutet de Monvel calculus, W is a pseudodifferential operator of order −1 on Σ.We need to find the dependence of the principal symbol of this operator on the matrices b(y) and a(y), y ∈ Σ.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 (b(σ)∇ σ ) * (b(σ)∇ σ ).This is a pseudodifferential operator of order 2 in Ω with principal symbol satisfying Further on, the truncated self-adjoint pseudodifferential operator M b = -(b∇R) * (b∇R) of order −2 has principal symbol ) where, recall, r −2 (σ, ς) = ã(σ)ς, ς −1 is the principal symbol of the fundamental solution R, so M b is an operator of order −2.Next, we need to make the restriction of M b to the boundary, y ∈ Σ.
The trace operators γ − N Dγ 1 acting from both sides upon M ≡ M b produce an order −1 pseudodifferential operator on the boundary with principal symbol, again, containing the matrix b * b and a homogeneous symbol depending algebraically on the matrix ã and its inverse ã−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, M (γ) b := γ(γM b ) * , which is the restriction of M b to the boundary.Its symbol can be evaluated as the expression (7.9) is calculated in local co-ordinates at the point y ∈ Σ, where y d+1 axis is directed along the normal to Σ.In a similar way, other terms in the symbol of Z * 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, ı ã(y)η, ν and the symbol β(y, η ′ ) of the N D operator, see (1.4).All 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 Z. Proposition 7.2.Let L be an elliptic operator with leading coefficients matrix ã(y) ∈ C ∞ (Ω).Suppose that for all y ∈ Ω, Let b(y) be a smooth symmetric matrix in C ∞ (Ω).Then for the singular numbers of the operator Z = b∇(γ T −1 ) * the estimate holds with constant A(C a ) depending only on the L ∞ norms of ã−1 and ã.
Proof.As we just found, the operator Z * Z is an order −1 pseudodifferential operator on the boundary Σ, with leading symbol W(y ′ , η ′ ), (y By the formula ( 23) in [14], for the eigenvalue counting function n(λ, Z * Z) the asymptotics holds (in our notations) and ω is the standard volume form on the unit sphere S d−1 .Thus, our estimate for the principal symbol W of the operator Z * Z, substituted in (7.13), gives (7.12).The last inequality in (7.12) follows from the fact that the L ∞ norm of b majorates its L p norm, for any p < ∞.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 ã(y), the constant in (7.14) depends only on the ellipticity constant |a −1 | and the norm |ã(y)|, but not on the matrix ã itself or its derivatives.
Proof.For a given ǫ, we construct, following the procedure in Lemma 4. L 2d+2 (Σ) , (7.18) which gives the required estimate, since the norm b(., 1) L∞ ≤ C a L∞ can be absorbed in the coefficient C.
Finally, using the basic asymptotic perturbation lemma, Lemma 2.2, we establish our main result.
Proof.Of Theorem 1.2.Let a be the coefficient matrix in our P-S problem in a smooth domain.Consider its smooth approximation ã constructed according to Lemma 4.4, with p = 2d + 2. For the operator H ã the asymptotic formula (1.6) is known.For the difference, H a −H ã we have the estimate (8.2).In these conditions, Lemma 2.2 grants that for the limit operator H a the eigenvalue asymptotics is valid and the coefficient in the asymptotics is given by the limit in the formula (1.6), using Lemma 2.5.

Appendix A. Potential theory approach
We consider here the Poincaré eigenvalue asymptotics.However, in [6], M.S.Agranovich and B.A.Amosov succeeded in proving order sharp two-sided eigenvalue estimates for potential type operators on Lipschitz surfaces, which, in their turn, produce two-sided eigenvalue estimates for N D.Moreover, a localization technic was developed in [6].It turned out that if the non-smooth singularities of the surface Σ are localized in the sense that the surface is smooth outside a closed set of zero measure, then the contribution of these singularities of Σ both to the behavior of eigenvalues of S and to the non-compactness of D is negligible, and the Weyl asymptotic formula for the P-S problem was proved in this setting.
In that period, after M.S.Agranovich started popularizing the problem on Steklov eigenvalue asymptotics for Lipschitz surfaces without additional restrictions, Prof. Grigory Tashchiyan and the author of this paper managed in 2006 to dispose of one of complications.Namely, in [37], we used the approximation of Lipschitz surfaces by smooth ones, in the sense of a certain integral norm, to prove that the corresponding potential type integral operators, properly relocated, so that they act on one and the same surface, can be compared, their difference becomes small in the spectral sense (see Lemma 6.3), and the passage to limit in formulas for spectral asymptotics becomes possible.Although this result did not solve the Lipschitz Steklov eigenvalue asymptotics problem, it solved a closely related one.Namely, in the interface spectral problem, where the spectral parameter is placed, for solutions of the (Laplace, in the simplest case) equation at the jump of the normal derivative across the Lipschitz surface S ⊂ Ω : the solution involves only the single layer potential operator S and the Weyl asymptotics for the Lipschitz surface is justified.Rather optimistically, we hoped that the Steklov problem will succumb as well.
One of the approaches we tried to explore was the following.Although we know that the double layer operator D is not expected to be compact, probably, the composition DS can be shown to be weaker than S itself.The notion 'weaker' should imply that the singular numbers of DS decay faster than the ones of S, In the eigenvalue studies, there are a lot of occasions when the multiplication by a non-compact operator nevertheless improves the rate of decay of the eigenvalues of the product.Hypothesis.Consider the composition DS of single and double layer potentials.Recall that S is the integral operator with kernel C|x − y| 1−d over a Lipschitz surface of dimension d > 1.Let T (x, y) be the integral kernel of the domains contains all C 1 domains, but does not contain all Lipschitz domains.By Theorem 4.35 in [24], singular integral operators of the type of double layer potential D, are compact in L 2 (Σ).By a localization of this property, as it was done, e.g., in (see, e.g., [6], [5], [2]) we arrive at the following result.where ρ e is the characteristic function of a small neighborhood of the set E, having surface measure less than ǫ.This leads to the corresponding splitting of the operator H in (4.5) into H = H ρǫ + H 1−ρǫ .For the first of these operators, the smallness for n sup (d, H ρǫ ) can be proved on the base of Theorem 2.1.For the remaining operator, H 1−ρǫ , it follows from the V MO 1 compactness property for the double layer potential that the n sup (d, H ρǫ ) = 0.One may hope that the approximation of a Lipschitz surface by V MO 1 ones may lead to a more straightforward proof of the Weyl asymptotics for the P-S eigenvalues.

1 . Introduction 1 . 1 .
The Poincaré-Steklov problem.Let Ω ⊂ R d+1 be a connected bounded open set.We suppose that the boundary Σ = ∂Ω is connected.The classical Poincaré-Steklov (P-S) eigenvalue problem consists in the study of the spectral properties for the Laplacian in Ω with the spectral parameter entering in the boundary condition,

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

G 1
Green function G 1 (y, z) for the Neumann problem, in other words, the integral kernel of the solution operator G G G 1 of the Neumann problem for L .This means that for a smooth function h on Σ,LG G G 1 h := Ly Σ G 1 (y, z)h(z)dz = 0, y ∈ Ω, (y, z)h(z)dz = h(y 0 ).(4.14)

7 . 1 .
these terms are majorated by |b(y)| 2 , with bounded dependence on |ã(y) −1 | and |ã(y)|, y ∈ Σ. 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 The principal symbol W(y, η ′ ), (y, η ′ ) ∈ T * Σ of the order −1 pseudodifferential operator W = Z * Z admits the estimate

Proposition B. 1 .
Let the Ω be a Lipschitz domain and let the boundary Σ belong to VMO 1 outside a closed set E of the surface measure zero.Then for the N D operator the Weyl asymptotic formula holds.The proof, leading, again, to (A.7), consists in splitting the Steklov problem ∂ ν u = λ −1 u on Σ into two weighted ones, ∂ ν u = λρ ǫ u, and ∂ ν u = λ(1 − ρ ǫ )u, (B.1)