---
abstract: 'A plasmon of a bounded domain $\Omega\subset{\mathbb{R}}^n$ is a non-trivial bounded harmonic function on ${\mathbb{R}}^n\setminus\partial\Omega$ which is continuous at $\partial\Omega$ and whose exterior and interior normal derivatives at $\partial\Omega$ have a constant ratio.  We call this ratio a plasmonic eigenvalue of $\Omega$. Plasmons arise in the description of electromagnetic waves hitting a metallic particle $\Omega$. We investigate these eigenvalues and prove that they form a sequence of numbers converging to one.  Also, we prove regularity of plasmons, derive a variational characterization, and prove a second order perturbation formula.  The problem can be reformulated in terms of Dirichlet-Neumann operators, and as a side result we derive a formula for the shape derivative of these operators.' 
address: 'Institut für Mathematik, Carl-von-Ossietzky Universität Oldenburg, 26111 Oldenburg, Germany'
author:
- Daniel Grieser
title: The plasmonic eigenvalue problem
---

Problem and results
===================

Let $\Omega\subset{\mathbb{R}}^n$, $n\geq 3$ be a bounded domain with smooth boundary $\partial\Omega$ and connected complement ${\mathbb{R}}^n\setminus\Omega$. For functions $u_-$ on ${\overline{\Omega}}$ and $u_+$ on ${\mathbb{R}}^n\setminus\Omega$, assumed to be smooth up to the common boundary $\partial\Omega$, we consider the following problem: $$\begin{aligned}
\Delta u_\pm=0\quad &&\text{in}\quad \Omega \text{ resp. }{ \mathbb{R}}^n\setminus {\overline{\Omega}}\label{I1}\\
u_--u_+=0 \quad &&\text{on}\quad  \partial\Omega\label{I2}\\
\varepsilon\partial_n u_- +\partial_n u_+ =0\quad &&\text{on}\quad  \partial\Omega\label{I3},\\
u_+(x) = O(|x|^{2-n}) \quad &&\text{as}\quad |x|\to\infty \label{I4}\end{aligned}$$ where $\partial_n u_\pm=\boldsymbol{n}\cdot \nabla u_\pm$ is the derivative of $u_\pm$ in the direction of the unit normal $\boldsymbol{n}$ pointing out of $\Omega$. Here ${\varepsilon}\in{\mathbb{R}}$ is a constant.  Thus, the function $u$ on ${\mathbb{R}}^n$ defined by $u_-$ on $\Omega$ and by $u_+$ on ${\mathbb{R}}^n\setminus\Omega$ is continuous at ${{\partial\Omega}}$, but its normal derivative may have a jump, and the ratio of inner and outer normal derivative is required to be constant on ${{\partial\Omega}}$. It will be convenient to also allow ${\varepsilon}=\infty$, in which case is interpreted as $\partial_n u_- \equiv 0$, $\partial_n u_+\not\equiv 0$.

We call $\varepsilon\in{{\mathbb{R}}\cup\{\infty\}}$ a *plasmonic eigenvalue* of $\Omega$ if the system - has a solution $u\not\equiv 0$, and then $u$ is called a plasmonic eigenfunction or (surface) [*plasmon* ]{} for $\varepsilon$. This name is justified by the origin in the physics of plasmons, and by the fact, established below, that ${\varepsilon}^{-1}$ is indeed an eigenvalue of a certain operator on $\partial\Omega$. For a given plasmonic eigenvalue we call the linear space of solutions $u$ its plasmonic eigenspace.  What is usually denoted $\varepsilon$ in the physics literature is $-\varepsilon$ in our notation.  Our choice of sign is motivated by the fact that $\varepsilon>0$ then.  A corresponding problem may also be formulated for $n=2$, but this requires a few minor adjustments, see Section \[subsec:2d\]. 

Plasmons have been studied extensively in the physics literature, though principally for simple geometries, where explicit calculations are possible, and for more general domains by numerical methods.  See the end of the introduction for more on the physics background and references.  The purpose of this paper is to put the problem in its natural mathematical setting and to bring some more refined mathematical techniques to bear on the study of plasmonic eigenvalues.  This will show that also for general domains interesting conclusions can be drawn even without numerical calculation.  In particular, we will employ microlocal analysis to analyze the asymptotic behavior of the plasmonic eigenvalues, and differential geometry to calculate the perturbation of plasmonic eigenvalues under perturbation of the domain, to second order.  This is the order needed for investigating the behavior under random perturbations of the domain, a topic of physical interest.  Also, the perturbation results can be useful in obtaining good approximations for non-infinitesimal deformations of the sphere, for example, see [@BieEtAl09]. 

We now state our results.  First, we have basic results on existence and regularity of plasmons, and asymptotics of plasmonic eigenvalues.  Denote $${{\mathcal H}}= \{u \text{ satisfying } \eqref{I1}, \eqref{I2}, \eqref{I4}\},$$ and on ${{\mathcal H}}$ consider the scalar product (for simplicity we assume all functions to be real-valued) $$\label{eqn:scprod+} 
 (u,v)_+ := \int_{\Omega^c} \nabla u\cdot\nabla v$$ where $\Omega^c={\mathbb{R}}^n\setminus\Omega$. This is finite since harmonic functions on $\Omega^{c}$ satisfying have $\nabla u (x)= O(|x|^{1-n})$. Nondegeneracy follows from the unique solvability of the Dirichlet problem for $\Omega$. Thus ${{\mathcal H}}$ with $(\cdot,\cdot)_+$ is a pre-Hilbert space. 

\[thm:existence+asymp\] Let $\Omega\subset{\mathbb{R}}^n$, $n\geq3$ be a bounded domain with smooth boundary.  Then the finite plasmonic eigenvalues of $\Omega$ form a sequence ${\varepsilon}_1,{\varepsilon}_2,\dots$ of positive numbers.  For ${\varepsilon}\neq 1$ the corresponding eigenspaces are finite dimensional. 

Furthermore, ${\varepsilon}_k\to 1$ as $k\to \infty$. Also, if $E_{\varepsilon}$ is the space of plasmons $u$ with eigenvalue ${\varepsilon}$ then we have an orthogonal decomposition $${{\mathcal H}}= \bigoplus_{{\varepsilon}\in{{\mathbb{R}}\cup\{\infty\}}} E_{\varepsilon}$$ Also, ${\varepsilon}=\infty$ is a plasmonic eigenvalue, with one-dimensional eigenspace $$E_\infty = \{u\in{{\mathcal H}}:\, u_-\text{ is constant}\}$$

See Section \[subsec:2d\] for the case $n=2$.

One may also formulate the eigenvalue problem with weaker regularity assumptions.  We then have the following result. 

\[thm:regularity\] Let $u$ be a solution of $\eqref{I1}-\eqref{I4}$ with $u_-\in H^s(\Omega)$, $u_+\in H^{s}_{{{\rm loc}}}(\Omega^c)$ for some $s>3/2$.

If ${\varepsilon}\neq 1$ then $u$ is $C^\infty$ on ${\mathbb{R}}^n\setminus{{\partial\Omega}}$ up to the boundary from each side. 

Here $H^s$ is the Sobolev space of order $s$, and $H^s_{{\rm loc}}(\Omega^c)$ is the space of functions whose restriction to $K_R=\Omega^c\cap\{|x|<R\}$ is in $H^s(K_R)$ for all $R$. The range of $s$ is chosen so that the restrictions of $u$ and ${\partial_n}u$ to the boundary are well-defined by the Sobolev restriction theorem. 

Examples show that the statement of the theorem is not true in general for ${\varepsilon}=1$, see Section \[subsec:disk\].  Note that the only issue is boundary smoothness since harmonic functions are always smooth in the interior of the domain. 

The main ingredient in the proof of Theorems \[thm:existence+asymp\] and \[thm:regularity\] is a reformulation as a problem on ${{\partial\Omega}}$ in terms of Dirichlet-to-Neumann operators, see Section \[subsec:dir-neu\].  This is closely related to the fact that for a harmonic function $u_-$ on $\Omega$ we have $$\label{eq:ibp}
\int_{{{\partial\Omega}}} u_- \partial_n u_- = \int_\Omega |\nabla u_-|^2 =: \|u\|_-^2$$ which follows from Green’s formula.  Similarly, for a harmonic function $u_+$ on ${\mathbb{R}}^n\setminus\overline{\Omega}$ satisfying , we have $$\label{eq:ibp+}
-\int_{{{\partial\Omega}}} u_+ \partial_n u_+ =  \int_{\Omega^c} |\nabla u_+|^2
=: \|u\|_+^2$$ Theorem \[thm:existence+asymp\] was formulated in terms of the scalar product $(\cdot,\cdot)_+$ corresponding to the norm $\|\cdot\|_+$ on ${{\mathcal H}}$. The orthogonality there also holds with respect to the bilinear form $(u,v)_- := \int_\Omega \nabla u \cdot \nabla v$ corresponding to $\|\cdot\|_-$ which, however, is not positive definite – it is zero on $E_\infty$.

Equations and imply that for a plasmon $u$ the plasmonic eigenvalue equals $$\label{eq:el-functional}
{\varepsilon}= R_p (u),\qquad R_p (u) := \frac{\int_{\Omega^c} |\nabla u_+|^2}{\int_{\Omega} |\nabla u_-|^2}$$ Denote by ${\dot{H}}^1({\mathbb{R}}^n)=\{u\in L^2_{{\rm loc}}({\mathbb{R}}^n): \nabla u \in L^2({\mathbb{R}}^n)\}$ the first homogeneous Sobolev space.  Then we have:

\[thm:euler-lagrange\] - is the Euler-Lagrange equation of the functional $R_p$ on ${\dot{H}}^1({\mathbb{R}}^n)\setminus 0$. More precisely, the finite plasmonic eigenvalues are precisely the positive critical values of $R_p$ and the associated plasmons are the corresponding critical points of $R_p$.

Note that zero is also a critical value of $R_p$, with any smooth function with non-empty support contained in $\Omega$ as critical point.  Of course, these do not correspond to plasmons. 

In this theorem we use $H^1_{{\rm loc}}$ instead of $H^1$ since functions satisfying need not be in $L^2({\mathbb{R}}^n)$ for $n\leq 4$.

Next we consider the behavior of plasmonic eigenvalues under perturbation of the domain.  Thus, fix $\Omega$ and let $a:\partial\Omega\to{\mathbb{R}}$ be a smooth function, and for real numbers $h$ close to zero consider the bounded domain $\Omega(h)$ with smooth boundary $$\label{eq:def Omega_h}
 \partial\Omega(h) = \{x + h a(x)\,{{\boldsymbol{n}}}(x):\, x\in \partial\Omega\}.$$ Thus, $\Omega(0)=\Omega$ and $\Omega(h)$ is obtained by shifting the boundary in the normal direction by $ha$. $h$ is the order of magnitude of the shift and $a$ is a ’shape’ function. 

In the following theorem, we state the perturbation formulas only in the physically interesting case $n=3$. Generalization to any $n\geq2$ is straightforward. 

[**Notation from differential geometry:**]{} $H$ and $K$ denote the mean and Gau[ß]{} curvature of ${{\partial\Omega}}$, $W_0=W-HI$ the trace free part of the Weingarten map $W$, and ${\operatorname{div}}_\partial$ and $\nabla_\partial$ the divergence and gradient operator in the surface ${{\partial\Omega}}$. The definitions of these quantities are recalled in Section \[secperturb\].  Also, $$\label{eq:def scalar product dO}
 \langle g,g' \rangle := \int_{{{\partial\Omega}}} g g' \,dS,\qquad g,g'\in {C^\infty}({{\partial\Omega}})$$ where $dS$ is surface measure (all functions are assumed real-valued), and for an operator $A$ on ${C^\infty}({{\partial\Omega}})$ we let $${\langle g \, |\, A \, | \, g' \, \rangle } := \langle g, Ag' \rangle$$ Similarly, for vector fields $V,V'$ on ${{\partial\Omega}}$ and operator $A$ acting on vector fields, ${\langle V \, |\, A \, | \, V' \, \rangle } := \langle V, AV' \rangle :=\int_{{{\partial\Omega}}} V\cdot AV'\,dS$. All operators appearing below are self-adjoint. 

\[thm:perturbation\] Let ${\varepsilon}\neq 1$ be a finite plasmonic eigenvalue of $\Omega\subset{\mathbb{R}}^3$ with eigenspace $E$. Then there are $h_0>0$ and real analytic functions $h\mapsto {\varepsilon}^{(i)}(h)$, $h\mapsto u^{(i)}(h)$ defined for $|h|<h_0$, $i=1,\dots,\dim E$, such that for each $h$ the numbers ${\varepsilon}^{(i)}(h)$, $i=1,\dots,\dim E$ are plasmonic eigenvalues of $\Omega(h)$ with eigenfunctions $u^{(i)}(h)$, and $u^{(1)}(0),\dots,u^{(\dim E)}(0)$ are a basis of $E$.

For a fixed analytic branch ${\varepsilon}(h)$, $u(h)$ of eigenvalue and eigenfunction with ${\varepsilon}(0)={\varepsilon}$ and $\|u(h)\|_-=1$ for each $h$ we have for $\dot{\varepsilon}:=\frac{d}{dh}{\varepsilon}(0), \ \dot{u}(x):=\frac{d}{dh}{u}(x,0),\
\ddot{\varepsilon}:=\left(\frac{d}{dh}\right)^2{\varepsilon}(0)
$ the formulas:

$$\begin{aligned}
{\dot{{\varepsilon}}}& = q_1(u) := ({\varepsilon}+1) \big[
- {\langle {\nabla_\partial}u_- \, |\, a \, | \, {\nabla_\partial}u_- \, \rangle } +
{\varepsilon}{\langle {\partial_n}u_- \, |\, a \, | \, {\partial_n}u_- \, \rangle }\big] \label{eq:epsdot}\\
&= ({\varepsilon}+ 1) \int_{{{\partial\Omega}}} a \left[-|{\nabla_\partial}u_-|^2 + {\varepsilon}({\partial_n}u_-)^2 \right]\, dS \notag \\
\label{eq:epsddot}
\tfrac12{\ddot{{\varepsilon}}}& =
\begin{array}[t]{c@{\,\langle\,}l@{\,|\,}c@{\,|\,}l@{\,\rangle}}
& {\nabla_\partial}u_-  & {-{\dot{{\varepsilon}}}a - ({\varepsilon}+1)a^2 W_0}  &
{{\nabla_\partial}u_-} \\ + &  {\nabla_\partial}u_- & ({\varepsilon}^2-1) a & {{\nabla_\partial}(a{\partial_n}u_-)} \\
+ & {{\nabla_\partial}u_-} & {-({\varepsilon}+1)a} & {{\nabla_\partial}{\dot{u}}_-}\\
+ & u_- & -{\dot{{\varepsilon}}}& {\partial_n}{\dot{u}}_-\\
+  & {{\partial_n}u_-} & {a{\varepsilon}({\dot{{\varepsilon}}}+ ({\varepsilon}+1)aH)} & {{\partial_n}u_-} \\
+ &  {{\partial_n}u_-} & {{\varepsilon}({\varepsilon}+1)a} & {{\partial_n}{\dot{u}}_-}
\end{array}\end{aligned}$$

Here, ${\dot{u}}$ is a solution of the inhomogeneous system $$\label{eq:udot}
\begin{aligned}
\Delta \dot{u} &=0\quad &&\text{in}\quad {\mathbb{R}}^3\setminus
\partial\Omega,\\
\dot{u}_--\dot{u}_+ &= -(\varepsilon+1) a\partial_n u_-\quad &&\text{on}\quad \partial\Omega,\\
\varepsilon \partial_n \dot{u}_- + \partial_n
\dot{u}_+&=-\dot{\varepsilon}\partial_n u_-+(\varepsilon+1)
{{\operatorname{div}}_\partial}(a{\nabla_\partial}u_-)\quad&&\text{on}\quad
\partial\Omega
\end{aligned}$$

Also, $u(0)$ is a critical point of the functional $R_p'(u) = \frac{q_1(u)}{\|u\|_-^2}$ on $E$, with $q_1$ defined in , with critical value ${\dot{{\varepsilon}}}$.

The meaning of analyticity of $h\mapsto u^{(i)}(h)$ will be explained in Section \[subsec:analytic\].  Note that ${\dot{u}}$ in is not unique, but the value of ${\ddot{{\varepsilon}}}$ calculated from is independent of the choice of ${\dot{u}}$. There is a similar inhomogeneous system determining ${\ddot{u}}$, see -.  In the proof of the theorem, the formulas for ${\dot{{\varepsilon}}}$, ${\ddot{{\varepsilon}}}$ are obtained from the requirement that these systems have a solution.  The last statement of the theorem helps to identify which $u$ to use in in case of degeneracy, i.e. $\dim E>1$.

As a side result we also obtain a perturbation formula for the Dirichlet-Neumann operator which may be of independent interest, see Theorem \[thm:DN-perturb\].  This is sometimes called the shape derivative of the Dirichlet-Neumann operator.  Recently there has been some interest in such results in the context of water waves, see [@lannes] for example. 

Physics background and previous work {#physics-background-and-previous-work .unnumbered}
------------------------------------

The problem - with $n=3$ arises in the physics of surface plasmon polaritons.  These are waves which couple electromagnetic fields to the electron gas of a metal, and thus propagate along the interface of a metal and a dielectric material or vacuum.  As outlined in several recent reviews, the particular interest in these excitations stems from the fact that by altering the structure of the metal’s surface, the properties of these plasmons, including their interaction with light, can be specifically tailored [@Atwater07; @BarnesEtAl03; @Maier07; @PitarkeEtAl07; @ZiaEtAl06].  This offers a high potential for developing new types of photonic devices, combining the speed of photonics with the nanoscale dimensions of electronics.  Possible applications include, among others, plasmonic chips, data storage, nanolithography, and highly sensitive molecular detectors. 

Our setting models a metallic particle whose shape is given by $\Omega$, surrounded by vacuum.  The function $u$ is the potential of the electric field, and the equations arise from Maxwell’s equations in the quasistatic approximation, which assumes constancy of the magnetic field, and is justified if the wavelength is much larger than the dimensions of the particle, which is the case for particles on the nano scale.  See [@landau-lifschitz Chapter 7] for details.  In the context of surface plasmons the quasistatic approximation has often been used in the literature, see [@KlimovGuzatov07b; @KlimovGuzatov07a] for example, and its validity has been verified experimentally. 

The boundary conditions , are obtained from continuity of the tangential component of the electric field and of the normal component of the electric displacement, see [@griffiths Chapter 4.4].  The number $-{\varepsilon}$ is the relative permittivity of the material of the particle, which is coupled to the frequency by a law specific to the material.  Thus, plasmonic eigenvalues correspond to frequencies of electromagnetic waves at which there is a strong excitation of the electron plasma in the metallic particle $\Omega$. The name plasmonic eigenvalue was coined in [@GrieserRuting09], where also the asymptotics ${\varepsilon}_k\to1$ was proved.  The paper [@BieEtAl09] contains the first order perturbation formula . 

The spectrum of the Laplacian on a bounded domain, and its relation to the geometry of the domain, has been studied intensively for many years, so it may be worth while to compare it to our problem.  There are similarities, for example the existence of discrete spectrum and the presence of the low energy and the high energy regime, which require very different mathematical techniques.  There are also clear differences, for example the asymptotic behavior and the fact that plasmonic eigenvalues are scale invariant while the Laplacian eigenvalues are not. 

Acknowledgement {#acknowledgement .unnumbered}
---------------

I am grateful to Gohar Harutyunyan for checking the perturbation formula independently. 

Basic observations {#sec:basics}
==================

Relation to Dirichlet-Neumann operators; proofs of Theorems \[thm:existence+asymp\] and \[thm:regularity\] {#subsec:dir-neu}
----------------------------------------------------------------------------------------------------------

Theorems \[thm:existence+asymp\] and \[thm:regularity\] are proved by reformulating the plasmonic eigenvalue problem as an ’honest’ eigenvalue problem for an operator on $C^\infty(\partial\Omega)$ which is built from Dirichlet-to-Neumann operators.  For basics on harmonic functions in $\Omega$ and ${\mathbb{R}}^n\setminus{\overline{\Omega}}$, and for proofs of some of the facts stated below, see for example [@folland].  The behavior at infinity of harmonic functions may, for example, be analyzed by using a Fourier decomposition on large circles. 

The (interior) Dirichlet-to-Neumann operator of $\Omega$ maps Dirichlet to Neumann data of a harmonic function.  That is, $${{\mathcal N}}_- : {C^\infty(\partial\Omega)}\to {C^\infty(\partial\Omega)},\quad g \mapsto \partial_n u_-$$ where $u_-$ is the unique solution of the Dirichlet problem $\Delta u_-=0$ in $\Omega$, $u_{-|\partial\Omega} = g$. Similarly, the exterior Dirichlet-to-Neumann operator is $${{\mathcal N}}_+ : {C^\infty(\partial\Omega)}\to {C^\infty(\partial\Omega)},\quad g \mapsto \partial_n u_+$$ where $u_+$ is the unique solution of the Dirichlet problem $\Delta u_+=0$ in ${\mathbb{R}}^n\setminus{\overline{\Omega}}$, $u_{+|\partial\Omega} = g$, $u_+(x)= O(|x|^{2-n})$ as $|x|\to\infty$. It is standard that this decay requirement already follows from boundedness of $u_+$.

${{\mathcal N}}_{\pm}$ are self-adjoint operators on $L^2(\partial\Omega)$ with the standard scalar product .  The symmetry of ${{\mathcal N}}_-$ follows immediately from Green’s formula, and for ${{\mathcal N}}_+$ can be reduced to that of ${{\mathcal N}}_-$ by use of the Kelvin transform. 

Note that while ${{\mathcal N}}_-$ is a positive operator by , ${{\mathcal N}}_+$ is negative because the normal ${{\boldsymbol{n}}}$ points into ${{\mathbb{R}}^n\setminus\overline{\Omega}}$ rather than out. 

\[prop:dirichlet-neumann\] The plasmonic eigenvalue problem - is equivalent to the generalized eigenvalue problem $$\label{eq:gen ev pr}
({\varepsilon}{{\mathcal N}}_- + {{\mathcal N}}_+) g = 0,$$ in the sense that smooth solutions $({\varepsilon},u)$ of - are in 1-1 correspondence with smooth solutions $({\varepsilon},g)$ of , via $g = u_{|\partial\Omega}$.

By the previous discussion we have a 1-1 correspondance $$\label{eq:1-1 corr}
\begin{aligned}
{{\mathcal H}}& \longleftrightarrow
C^\infty({{\partial\Omega}}) \\
u &\mapsto u_{|{{\partial\Omega}}} \notag
\end{aligned}$$ Under this correspondence is equivalent to . 

For $g,g' \in C^\infty({{\partial\Omega}})$ we denote $$(g,g')_+ := - \langle g, {{\mathcal N}}_+ g' \rangle$$ If $u,u'\in{{\mathcal H}}$ are the functions on ${\mathbb{R}}^n$ corresponding to $g,g'$ under the correspondence then Green’s formula, cf. ,  gives $$\label{eqn:isometry}
 (g,g')_+  = (u,u')_+$$ where the latter is defined in .  This justifies the notation and shows that $(\cdot,\cdot)_+$ is a scalar product on $C^\infty({{\partial\Omega}})$. The latter fact also follows from the well-known fact that ${{\mathcal N}}_+$ is invertible. 

The problem can easily be transformed into a standard eigenvalue problem: The plasmonic eigenvalues are the reciprocals of the eigenvalues of the operator $$A=-{{\mathcal N}}_+^{-1}{{\mathcal N}}_-.$$ This operator is symmetric with respect to $(\cdot,\cdot)_+$. Introduce the Hilbert space $${\overline{{{\mathcal H}}}}:= \text{ completion of $C^\infty(\partial\Omega)$ with respect to $(\cdot,\cdot)_+$.}$$

\[Proof of Theorem \[thm:existence+asymp\]\] Positivity of plasmonic eigenvalues follows from : $R_p(u)$ must be non-zero since otherwise $u_+\equiv 0$, so $u_{-|{{\partial\Omega}}}=0$, hence $u\equiv 0$.

We use some concepts and facts from microlocal analysis.  See [@folland] for background.  It is well-known that ${{\mathcal N}}_-,{{\mathcal N}}_+$ are first order classical pseudodifferential operators on $C^\infty(\partial\Omega)$ with principal symbols $$\sigma({{\mathcal N}}_-) = |\xi|,\quad  \sigma({{\mathcal N}}_+) = -|\xi|.$$ The fact that ${{\mathcal N}}_+$ is an invertible elliptic first order operator implies that ${\overline{{{\mathcal H}}}}=H^{1/2}(\partial\Omega)$, the Sobolev space of order $1/2$. By the standard theory of pseudodifferential operators, $A=-{{\mathcal N}}_+^{-1}{{\mathcal N}}_-$ is a zeroth order pseudodifferential operator with principal symbol $-(-|\xi|^{-1}) \cdot |\xi|=1$, hence $A = I + R$ where $R$ is a pseudodifferential operator of order $-1$. Therefore, $R$ is a compact operator on ${\overline{{{\mathcal H}}}}$, and since $A$ is selfadjoint, so is $R$, so $R$ has real nonzero eigenvalues $r_1,r_2,\dots \to 0$ of finite multiplicity, and may have a zero eigenvalue $r_0=0$ of any multiplicity.  The eigenspaces span ${\overline{{{\mathcal H}}}}$ and are orthogonal with respect to $(\cdot,\cdot)_+$. Then $a_k=1+r_k$ are the eigenvalues of $A$, so $a_k^{-1}$ are the plasmonic eigenvalues and $a_k^{-1}\to 1$ as $k\to\infty$. By the plasmonic eigenspaces are orthogonal with respect to the scalar product $(\cdot,\cdot)_+$ on ${{\mathcal H}}$. There is an eigenvalue $a_k=0$ which corresponds to the kernel of ${{\mathcal N}}_-$, that is, the constants in ${\overline{{{\mathcal H}}}}$, and these correspond to the space $E_\infty$.

Clearly, the correspondence in Proposition \[prop:dirichlet-neumann\] also holds for $H^s$-plasmons.  The operator ${\varepsilon}{{\mathcal N}}_-+{{\mathcal N}}_+$ has principal symbol $({\varepsilon}-1)|\xi|$. If ${\varepsilon}\neq 1$ then this is invertible for each $\xi\neq0$, hence the operator ${\varepsilon}{{\mathcal N}}_-+{{\mathcal N}}_+$ is invertible then.  So by elliptic regularity its kernel consists of smooth functions.  These functions are the boundary values of the corresponding plasmons, so these are smooth up to the boundary by elliptic boundary regularity. 

Variational principle {#subsec:var princ}
---------------------

We calculate the first variation of $R_p$. Let $u\in {\dot{H}}^1({\mathbb{R}}^n)$ be such that $\int_\Omega |\nabla u|^2 >0$ and let ${\varepsilon}=R_p(u)$. For $v\in C_0^\infty({\mathbb{R}}^n)$, the space of smooth, compactly supported functions, we have $$\begin{gathered}
\frac{\delta R_p}{\delta u}(v)
= \frac 2 {\int_\Omega |\nabla u_-|^2} \left( \int_{\Omega^c}\nabla u_+\,\nabla v_+ - {\varepsilon}\int_\Omega \nabla u_-\, \nabla v_- \right) 
\label{eqn:ibp varpr}\\
\qquad\qquad
= \frac 2 {\int_\Omega |\nabla u_-|^2} \left( \int_{{{\partial\Omega}}} (\partial_n u_+ + {\varepsilon}\partial_n u_-) v - \int_{\Omega^c} (\Delta u_+) v_+ + {\varepsilon}\int_{\Omega} (\Delta u_-) v_- \right) 
\label{eqn:ibp varpr2}\end{gathered}$$ using integration by parts.  If $u$ satisfies - then it follows that $\frac{\delta R_p}{\delta u}(v) =0$ for all $v\in C_0^\infty({\mathbb{R}}^n)$. Now $C_0^\infty({\mathbb{R}}^n)$ is dense in ${\dot{H}}^1({\mathbb{R}}^n)$, and for $u$ satisfying in addition the right hand side of depends continuously on $v$, so $u$ is a critical point of $R_p$, and $R_p(u)={\varepsilon}>0$ by definition and Theorem \[thm:existence+asymp\]. 

Conversely, assume $u$ is a critical point of $R_p$ with ${\varepsilon}=R_p(u)>0$. Using smooth test functions $v$ localized near a point in $\Omega$ or the interior of $\Omega^c$ one sees from , that $\Delta u=0$ on ${\mathbb{R}}^n\setminus \partial \Omega$, and then using $v$ localized near points on ${{\partial\Omega}}$ one gets ${\varepsilon}\partial_n u_- + \partial_n u_+ =0$ at ${{\partial\Omega}}$. By the Sobolev restriction theorem we have $u_-=u_+$ at ${{\partial\Omega}}$. Finally, a harmonic function $u$ on $\Omega^c$ with $\nabla u\in L^2(\Omega^c)$ satisfies , so $u$ is a plasmon with eigenvalue ${\varepsilon}$.

Scale invariance {#subsec:scale inv}
----------------

As opposed to the eigenvalue problem for the Laplacian, the plasmonic eigenvalue problem is scale invariant: If $t>0$ then the domains $\Omega$ and $t\Omega = \{tx:\, x\in\Omega\}$ have the same plasmonic eigenvalues, since an eigenfunction $u$ for $\Omega$ turns into an eigenfunction $u_t(x) := u(x/t)$ for $t\Omega$ with the same ${\varepsilon}$.

The two-dimensional case {#subsec:2d}
------------------------

The condition on the behavior of $u$ at infinity is natural for all dimensions $n\geq 2$ in that it gives unique solvability of the Dirichlet problem on $\Omega^c$ for any boundary data on ${{\partial\Omega}}$. However, this condition is not suitable for formulating an analogue of the plasmonic eigenvalue problem for $n=2$ because any constant function $u$ will satisfy - for [*any*]{} value of ${\varepsilon}$. On the other hand, adding a constant to a plasmon $u$ with arbitrary plasmonic eigenvalue will yield a plasmon with the same eigenvalue.  Therefore, it is natural to consider plasmons as elements of the quotient space ${{\mathcal H}}/\{\text{constants}\}$. Since any bounded harmonic function on $\Omega^c\subset{\mathbb{R}}^2$ is of the form $\,\text{constant}+ O(|x|^{-1})$, this is equivalent to replace the decay condition by the condition $$\label{I4 n=2}
u_+(x) = O(|x|^{-1}) \quad \text{as}\quad |x|\to\infty \tag{$\ref{I4}_{n=2}$}$$ Denote the space of plasmons satisfying , and by ${{\mathcal H}}'$. On this space $(\cdot,\cdot)_+$ is a scalar product, and Theorem \[thm:existence+asymp\] holds, except that now we have the decomposition $${{\mathcal H}}' = \bigoplus_{{\varepsilon}\in{\mathbb{R}}} E_{\varepsilon}$$ and there is no space $E_\infty$. Theorems \[thm:regularity\] and \[thm:euler-lagrange\] hold as in $n\geq3$, with the same proofs. 

In Section \[subsec:dir-neu\] the definition and discussion of the Dirichlet-Neumann operators ${{\mathcal N}}_\pm$ before Proposition \[prop:dirichlet-neumann\], and also Proposition \[prop:dirichlet-neumann\], hold for $n=2$ (where condition is used to define ${{\mathcal N}}_+$).  However, the operator ${{\mathcal N}}_+$ is not invertible if $n=2$. More precisely, we have $$\ker{{\mathcal N}}_+=\ker{{\mathcal N}}_-=\{\text{constants}\}$$ and by self-adjointness both ranges are equal to the orthogonal complement of the constants (with respect to the standard scalar product ), which we denote by ${{\mathcal R}}$. Therefore, both ${{\mathcal N}}_+$ and ${{\mathcal N}}_-$ descend to invertible operators on $C^\infty({{\partial\Omega}})/\{\text{constants}\}$, which may be identified with ${{\mathcal R}}$. Then the plasmonic eigenvalues are again the reciprocals of the eigenvalues of $A=-{{\mathcal N}}_+^{-1}{{\mathcal N}}_-$ considered as operator ${{\mathcal R}}\to{{\mathcal R}}$ since the isomorphism ${{\mathcal H}}\longleftrightarrow C^\infty({{\partial\Omega}})$ induces an isomorphism ${{\mathcal H}}/\{\text{constants}\} \longleftrightarrow C^\infty({{\partial\Omega}})/\{\text{constants}\}$. The operator $A$ is invertible, hence $\infty$ is not a plasmonic eigenvalue.  The proof of Theorem \[thm:existence+asymp\] carries over if one extends $A$ to ${\overline{{{\mathcal H}}}}$ by continuity and by setting it zero on constants. 

There is a subtlety, which is the reason why quotient spaces are more natural in this discussion than their replacements ${{\mathcal H}}'$ and ${{\mathcal R}}$: The boundary values of harmonic functions on $\Omega^c$ satisfying condition are [*not*]{} the functions in ${{\mathcal R}}$. Rather, they are characterized as follows.  There is a unique harmonic function $v$ on $\Omega^c$ having boundary values $v_{|{{\partial\Omega}}}=0$ and satisfying $v(x) \sim \frac1{2\pi} \log |x|$ as $|x|\to\infty$. This function may be constructed by choosing any point $y\in\Omega$ and letting $v=G-G'$ where $G(x)=\frac1{2\pi}\log|x-y|$ is the Newton potential centered at $y$ and $G'$ is the unique bounded harmonic function on $\Omega^c$ having the same boundary values on ${{\partial\Omega}}$ as $G$. Let $g_0=\partial_\nu v\in C^\infty({{\partial\Omega}})$. Then we have:

> $g\in C^\infty({{\partial\Omega}})$ is the boundary value of a harmonic function on $\Omega^c$ satisfying if and only if $\langle g,g_0\rangle=0$.

Necessity follows from Green’s formula, and then sufficiency from a dimension argument.  The function $g_0$ is not constant, except when $\Omega$ is a disk.  This means that the harmonic extension $u$ of $g\in{{\mathcal R}}$ will not satisfy the bound in general.  In order to obtain a plasmon satisfying one has to add a constant to $u$.

Perturbation of plasmonic eigenvalues {#secperturb}
=====================================

In this section we prove Theorem \[thm:perturbation\]. 

Analyticity {#subsec:analytic}
-----------

The analytic dependence of the plasmonic eigenvalues on $h$ is proved by standard methods: First, translate into a problem on a fixed domain, then show that the operators involved depend analytically on $h$, then apply standard perturbation theory.  We provide the details: In a first step, we translate the problem for $\Omega(h)$ to a ($h$-dependent) problem for the $h$-independent domain $\Omega$. To do this, we extend the vector field $a{{\boldsymbol{n}}}$ on ${{\partial\Omega}}$ to a smooth vector field $V$ on ${\mathbb{R}}^n$, which vanishes outside some small neighborhood $U$ of ${{\partial\Omega}}$. This can be done by first extending $a{{\boldsymbol{n}}}$ to a neighborhood of ${{\partial\Omega}}$ and then multiplying by a smooth function on ${\mathbb{R}}^n$ supported in this neighborhood and equal to one on ${{\partial\Omega}}$. Then define maps $T_h:{\mathbb{R}}^n\to{\mathbb{R}}^n, x\mapsto x+ h V(x)$. For any $h$, $T_h$ is the identity outside $U$, and by compactness of ${{\partial\Omega}}$ and the inverse function theorem it is a diffeomorphism on ${\mathbb{R}}^n$ for $h$ sufficiently close to zero.  By definition, $T_h$ is a diffeomorphism $\Omega\to\Omega(h)$.

Now the map $T_h$ is used to ’transport’ data related to $\Omega(h)$ back to $\Omega$: Let ${{{\bf g}}_{\text{eucl}}}$ be the Euclidean metric on ${\mathbb{R}}^n$ and $\bf{{\bf g}}_h = T_h^* {{{\bf g}}_{\text{eucl}}}$ its pull-back.  Denote by $\Delta_h$, $\nabla_h$ the Laplace-Beltrami and gradient operators for the Riemannian metric $\bf{{\bf g}}_h$ on ${\mathbb{R}}^n$ and by ${{\boldsymbol{n}}}_h$ the normal vector field with respect to $\bf{{\bf g}}_h$ on ${{\partial\Omega}}$. Then the plasmonic eigenvalue problem for $\Omega(h)$ is equivalent to the problem for $\Omega$ with $\Delta, \partial_n$ replaced by $\Delta_h, \partial_{n_h}$. More precisely, $u$ solves the former if and only if $T_h^*u:=u\circ T_h$ solves the latter.  Also, in the discussion of the variational principle one needs to replace $\nabla$ by $\nabla_h$, and in Proposition \[prop:dirichlet-neumann\] the Dirichlet-Neumann operators need to be replaced by those for the metric $\bf{{\bf g}}_h$, denoted ${{\mathcal N}}_{\pm, h}$.

We can now state the precise meaning of analyticity in Theorem \[thm:perturbation\]. 

\[prop:analytic\] Let ${\varepsilon}\neq 1$ be a finite plasmonic eigenvalue of $\Omega$ with eigenspace $E$. Then there are $h_0>0$ and real analytic functions $${\varepsilon}^{(i)}: (-h_0,h_0)\to {\mathbb{R}},\quad u^{(i)}:(-h_0,h_0) \to \{u\in C({\mathbb{R}}^n): \, u_{|{\overline{\Omega}}}\text{ and }u_{|{\mathbb{R}}^n\setminus \Omega}\text{ are smooth}\}$$ for $i=1,\dots,\dim E$ such that $u^{(1)}(0),\dots,u^{(\dim E)}(0)$ form a basis of $E$ and the pairs $((T_h^{-1})^*u^{(i)}(h), {\varepsilon}^{(i)}(h))$ satisfy the plasmonic eigenvalue problem for the domain $\Omega(h)$, for each $h\in(-h_0,h_0)$.

The analyticity of $u^{(i)}$ in the given space is to be understood as analyticity of the restrictions in the Hilbert spaces $H^s(\Omega)$ and $H^s(\Omega^c\cap\{|x|<R\})$ for each $R>0$ and $s\in{\mathbb{R}}$.

We first consider the reduction of the problem given in and after Proposition \[prop:dirichlet-neumann\].  Thus, we are dealing with an eigenvalue problem for the operator $A_h=-{{\mathcal N}}_{+,h}^{-1}{{\mathcal N}}_{-,h}$ on $L^2({{\partial\Omega}})$. The analytic dependence of plasmonic eigenvalues and eigenfunctions, as elements of $L^2({{\partial\Omega}})$, follows from standard perturbation theory, see [@Kat:PTLO], if we can show that $h\mapsto {{\mathcal N}}_{\pm,h}$ are analytic families of operators.  We show this for ${{\mathcal N}}_h:={{\mathcal N}}_{-,h}$, the argument for ${{\mathcal N}}_{+,h}$ is analogous.  By [@Kat:PTLO Theorem VII.4.2] we only need to show that the associated family of quadratic forms, ${{\boldsymbol{q}}}_{h}(g) = \langle g, {{\mathcal N}}_{h} g \rangle$, is analytic in $h$, in the sense that the form domain is independent of $h$ and that $h\to {{\boldsymbol{q}}}_{h}(g)$ is analytic for each $g$ in this domain.  Now since ${{\mathcal N}}_{h}$ is an elliptic pseudodifferential operator of order one, the domain of ${{\boldsymbol{q}}}_{h}$ is the Sobolev space $H^{1/2}({{\partial\Omega}})$ for each $h$. Also, for each $g\in H^{1/2}({{\partial\Omega}})$ we have (compare ) $ \langle g,{{\mathcal N}}_{h}g \rangle =  \int_{\Omega} |\nabla_h v_{h} |^2 \, {\text{dvol}}_h
$ where $v_{h}\in H^1(\Omega)$ is the solution of the Dirichlet problem on $\Omega$ with boundary value $g$ (here we need to impose also the decay condition in the case of ${{\mathcal N}}_{+,h}$).  So we need to prove analyticity of $h\mapsto  \langle g,{{\mathcal N}}_{h}g \rangle =  \int_{\Omega} |\nabla_h v_{h} |^2 \, {\text{dvol}}_h
$.  We first prove analyticity of $h\mapsto v_h$: Since the metric $\bf{{\bf g}}_h$ depends analytically on $h$ (in fact, ${{\bf g}}_h(x)$ is given by the matrix $(I+h\, DV(x))^t (I+h\,DV(x))$ in standard coordinates), so does $\Delta_h$. Let $G_{h}$ be the inverse of $\Delta_h$ on $\Omega$ with Dirichlet boundary conditions at ${{\partial\Omega}}$. Then $G_{h}$ is analytic in $h$ by [@Kat:PTLO Theorem VII.1.3].  We have $v_{h} = v - G_h \Delta_h v$ where $v=v_{0}$ is the unperturbed solution.  This shows that $v_h$ is analytic in $h$. Finally, $\nabla_h v_h$ and ${\text{dvol}}_h = \sqrt{\det \bf{{\bf g}}_h}\,dx$ are analytic in $h$, so we have proved the analytic dependence of ${{\mathcal N}}_h$ on $h$.

Now the cited theorems give the existence of ${\varepsilon}^{(i)}(h)$ and of branches $g^{(i)}(h)\in L^2({{\partial\Omega}})$, depending analytically on $h$. Then since ${\varepsilon}^{(i)}(0)={\varepsilon}\neq 1$ we can choose $h_0$ small enough so that ${\varepsilon}^{(i)}(h)\neq1$ for all $h\in(-h_0,h_0)$. Then the standard parametrix construction gives, for any $s\in{\mathbb{N}}$, a left parametrix $Q_h$ of $P_h={\varepsilon}^{(i)}(h){{\mathcal N}}_{-,h} + {{\mathcal N}}_{+,h}$ with error term $R_h=Q_hP_h-I\in\Psi^{-s}({{\partial\Omega}})$, with $Q_h$ depending on $h$ analytically.  Then $0 = Q_hP_hg^{(i)}(h) = g^{(i)}(h) + R_hg^{(i)}(h)$, and this shows that $g^{(i)}(h)$ is analytic with values in $H^s({{\partial\Omega}})$. Finally, this implies that $u^{(i)}(h)$, defined as the solution of the Dirichlet problem with respect to $\Delta_h$ with boundary values $g^{(i)}(h)$, depends analytically on $h$ in any Sobolev space. 

The perturbation calculation {#subsec:perturb calc}
----------------------------

### Procedure

We first explain our procedure.  Consider - with all quantities depending on $h$. Differentiate these equations in $h$ and set $h=0$. This yields a system of equations $$\begin{aligned}
\Delta {\dot{u}}& = 0 \quad\text{in }{\mathbb{R}}^3\setminus\partial\Omega
\label{eqn1a}\\
{\dot{u}}_- - {\dot{u}}_+ & = F_1 \quad\text{on }\partial\Omega
\label{eqn2a}\\
{\varepsilon}\partial_n {\dot{u}}_- + \partial_n {\dot{u}}_+ & = -{\dot{{\varepsilon}}}\partial_n u_- + G_1 \quad\text{on } \partial\Omega
\label{eqn3a}\end{aligned}$$ satisfied by ${\dot{u}}$, where $F_1,G_1$ are explicit expressions in $u$ and ${\varepsilon}$, see ,.  This is an inhomogeneous version of the system -.  Since the homogeneous system has a nontrivial solution (the given $u$), solvability of the inhomogeneous system implies certain compatibility conditions for the right hand sides $F_1$, $-{\dot{{\varepsilon}}}\partial_n + G_1 $, see Lemma \[lem:compat\].  This yields a formula for ${\dot{{\varepsilon}}}$. Similarly, differentiating twice in $h$ before setting $h=0$ gives $$\begin{aligned}
\Delta {\ddot{u}}& = 0 \quad\text{in }{\mathbb{R}}^3\setminus\partial\Omega
\label{eqn1b}\\
{\ddot{u}}_-- {\ddot{u}}_+ & = F_2 \quad\text{on }\partial\Omega
\label{eqn2b}\\
{\varepsilon}\partial_n {\ddot{u}}_- + \partial_n {\ddot{u}}_+ & = -{\ddot{{\varepsilon}}}\partial_n u_- + G_2 \quad\text{on } \partial\Omega
\label{eqn3b}\end{aligned}$$ with $F_2,G_2$ given in , , and in the same way as before this yields the formula for ${\ddot{{\varepsilon}}}$.

### Calculation of derivatives at ${{\partial\Omega}}$

In order to determine $F_1,G_1,F_2,G_2$ in -, we need to calculate certain derivatives at the boundary.  It is important to express these in invariant differential geometric terms.  For a point $p\in{{\partial\Omega}}$ we decompose orthogonally ${\mathbb{R}}^3=T_p{\mathbb{R}}^3 = T_p {{\partial\Omega}}\oplus {\operatorname{span}}{{{\boldsymbol{n}}}(p)}$ and express all quantities with respect to this decomposition.  For example, for a function $u$ we have $$\nabla u = {\nabla_\partial}u + (\partial_n u ) {{\boldsymbol{n}}}.$$

Denote by ${{\boldsymbol{n}}}(q,h)$ the unit normal to ${{\partial\Omega}}_h$ at $q\in{{\partial\Omega}}_h$. Recall that $K$, $H$, $W$ denote the Gauss curvature, mean curvature and Weingarten map of ${{\partial\Omega}}$, respectively.  Let ${\Delta_\partial}= {{\operatorname{div}}_\partial}{\nabla_\partial}$ be the Laplace-Beltrami operator on ${{\partial\Omega}}$. The definitions of these and the other quantities is recalled in the proof of Lemma \[lem:geom dO in Koords u,v\]. 

\[prop:deriv at dO\] Fix $p\in{{\partial\Omega}}$. Denote by ${{\boldsymbol{\dot{n}}}}(p)$, ${{\boldsymbol{\ddot{n}}}}(p)$ the first and second derivative in $h$, evaluated at $h=0$, of the unit normal ${{\boldsymbol{n}}}(q_h,h)$ at $q_h= p + h a(p){{\boldsymbol{n}}}(p) \in{{\partial\Omega}}_h$. Then we have at $p$ $$\begin{aligned}
\dot{\boldsymbol{n}}  & =-{\nabla_\partial}a  \\
\ddot{\boldsymbol{n}} & =-W({\nabla_\partial}a^2) -|{\nabla_\partial}a |^2 \boldsymbol{n} \end{aligned}$$ Furthermore, if $x\mapsto u(x,h)$ is harmonic on $\Omega_h$ for each $h$, smooth in $(x,h)$ up to the boundary, then $h\mapsto \nabla u(q_h,h)$ has first and second derivatives at $h=0$ $$\begin{aligned}
\big(\nabla u\big)^{\cdot} &= a(T_1 u + T_1' \partial_n u) && + a(N_1 u + N_1' {\partial_n}u)\,{{\boldsymbol{n}}}&& + \nabla {\dot{u}}\\
\big(\nabla u\big)^{\cdot\cdot}&= a^2(T_2 u + T_2' {\partial_n}u) && + a^2(N_2 u + N_2' {\partial_n}u)\,{{\boldsymbol{n}}}&& + 2 (\nabla {\dot{u}})^{\cdot} - \nabla {\ddot{u}}\end{aligned}$$ where $$\begin{aligned}
T_1 &=  W {\nabla_\partial}, & T_1' &=  {\nabla_\partial}, \\
N_1 &= -{\Delta_\partial}, &   N_1' &= 2H, \\
N_2 &= 2 ( -H{\Delta_\partial}-{{\operatorname{div}}_\partial}W {\nabla_\partial}+ {\nabla_\partial}H\cdot {\nabla_\partial}),
& N_2' &= -{\Delta_\partial}+8H^2-2K.\end{aligned}$$

Here the operators $T_i,T_i', N_i, N_i'$ act on functions defined on ${{\partial\Omega}}$.

The tangential component of $\big(\nabla u\big)^{\cdot\cdot}$ is not needed, so we don’t calculate it.  The term $(\nabla{\dot{u}})^{\cdot}$ can be calculated using the same formula as for $(\nabla´u)^{\cdot}$ since ${\dot{u}}(h)$ is also harmonic for each $h$. The same conclusion holds if $u(\cdot,h)$ is harmonic on ${\mathbb{R}}^n\setminus\overline{\Omega_h}$ for each $h$.

The main point in the formulas for the variation of $\nabla u$ is that they involve at most first derivatives in the normal direction.  The second and third normal derivatives, which occur in calculating the left hand side, are expressed in terms of tangential derivatives and curvature terms, using the harmonicity of $u$. See Lemma \[lemma2\]. 

Fix $p\in \partial\Omega.$ We introduce coordinates in ${\mathbb{R}}^3$ so that $p=0$ and that the tangent plane to $\partial\Omega$ in $p$, $T_p\partial\Omega$, is horizontal and $\Omega$ lies below it.  Then we can parametrize the surface ${{\partial\Omega}}$ near $p$ as a graph: $$\label{1.1}
x(v,w)=(v,w,f(v,w)),$$ so that, near $p$, $\partial\Omega=\{x(v,w): v,w\quad \text{near}\quad 0\}$, and $f=0$, $\nabla f=0$ at $(0,0)$. We may rotate coordinates so that $f$ has the form $$\label{1.2}
f(v,w)=\frac{\lambda}{2}v^2+\frac{\mu}{2}w^2+{{{\mathcal O}}}^3$$ for some $\lambda,\mu\in{\mathbb{R}},$ where here und in the sequel ${{{\mathcal O}}}^k$ denotes any quantity vanishing at least to order $k$ at $v=w=0.$ Geometrically, $\lambda,\mu$ are the principal curvatures of $\partial\Omega$ at $p$.

Let $\boldsymbol{n}=\boldsymbol{n}(v,w)$ be the upward unit normal: $$\label{eq:normal}
\boldsymbol{n}(v,w)=\frac{(-f_v(v,w),
-f_w(v,w),1)}{\sqrt{1+f_v^2(v,w)+f_w^2(v,w)}}.$$ Expanding this to second order we get $$\begin{aligned}
\boldsymbol{n}(0,0)&=&(0,0,1),\nonumber\\
\boldsymbol{n}_v(0,0)&=&(-\lambda,0,0),\label{1.3}\\
\boldsymbol{n}_{vv}(0,0)&=&(-f_{vvv},-f_{vvw},-\lambda^2)\nonumber.\end{aligned}$$ and analogously with $v,w$ interchanged and $\lambda$ replaced by $\mu$. Here and in the sequel derivatives of $f$ are evaluated at $(0,0)$ if no arguments are given. 

We now express geometric quantities of ${{\partial\Omega}}$ in the coordinates $v,w$.

\[lem:geom dO in Koords u,v\] For the metric tensor $g$, second fundamental form $II$ and Weingarten map $W$ we have in the coordinates $v,w$ $$\begin{aligned}
g &= {{\rm Id}}+ {{\mathcal O}}^2
\label{1.5} \\
\label{1.6}
II &=W+{{\mathcal O}}^2 =
\left( \begin{array}{cc} \displaystyle{\lambda}& 0\\
 0 & \displaystyle{\mu}
\end{array}\right)
+
\left( \begin{array}{cc} \displaystyle{v f_{vvv}+w
f_{vvw}}& \displaystyle{v f_{vvw}+w f_{vww}}\\
\displaystyle{v f_{vvw}+w f_{vww}}& \displaystyle{v
f_{wwv}+w f_{www}}
\end{array}\right)+{{{\mathcal O}}}^2\end{aligned}$$ Here ${{\rm Id}}$ is the identity matrix.  The derivatives of $f$ are evaluated at the origin. 

The intrinsic differential operators in the surface are $$\label{eq:div nabla Laplace}
{\nabla_\partial}=\nabla+ {{\mathcal O}}^2,\ {{\operatorname{div}}_\partial}= {\operatorname{div}}+ {{\mathcal O}}^1,\ {\Delta_\partial}= \Delta + {{\mathcal O}}^1$$ where $\nabla,{\operatorname{div}},\Delta$ are the standard (Euclidean) operators in $v,w$.

We note the following consequences of the lemma, needed below: The mean curvature is $$H=\frac{1}{2}\text{trace}\,W=\frac{1}{2}(\lambda+\mu+v\Delta f_v+w\Delta f_w)+{{{\mathcal O}}}^2.$$ and then simple calculations give $$\label{1.10}
{{\operatorname{div}}_\partial}W {\nabla_\partial}= \lambda \partial_v^2+\mu \partial_w^2 +\nabla (\Delta f)\cdot \nabla +{{\mathcal O}}^1$$ and $$\label{1.11}
{\nabla_\partial}H= \frac{1}{2}\nabla\Delta f + {{\mathcal O}}^1.$$

For the metric tensor $g$ we have $$\begin{aligned}
g(v,w)&=& \left(
\begin{array}{cc}
\displaystyle{g_{11}}&\displaystyle{g_{12}}\\
\displaystyle{g_{21}}&\displaystyle{g_{22}}
\end{array}
\right)
 =\left(
\begin{array}{cc}
\displaystyle{\partial_v x\cdot\partial_v x}&\displaystyle{\partial_v x\cdot\partial_w x}\\
\displaystyle{ \partial_w x\cdot\partial_v
x}&\displaystyle{ \partial_w x\cdot\partial_w x}
\end{array}
\right)\nonumber\\
&=&\left(
\begin{array}{cc}
\displaystyle{1+f_v^2(v,w)}&\displaystyle{f_v(v,w) f_w(v,w)}\\
\displaystyle{f_v(v,w) f_w(v,w)}&\displaystyle{1+f_w^2(v,w)}
\end{array}\right) = {{\rm Id}}+ {{\mathcal O}}^2\nonumber\end{aligned}$$

The second fundamental form, expressing curvature, is defined by $$II=\left(
\begin{array}{cc}
\displaystyle{II_{11}}&\displaystyle{II_{12}}\\
\displaystyle{II_{21}}&\displaystyle{II_{22}}
\end{array}
\right)=\left(\begin{array}{cc}
\displaystyle{\boldsymbol{n}(v,w)\cdot \partial_v^2 x(v,w)}&\displaystyle{\boldsymbol{n}(v,w)\cdot \partial_{vw}^2 x(v,w)}\\
\displaystyle{\boldsymbol{n}(v,w)\cdot \partial_{vw}^2
x(v,w)}&\displaystyle{\boldsymbol{n}(v,w)\cdot \partial_w^2 x(v,w)}
\end{array}\right),$$ where $$\begin{aligned}
\boldsymbol{n}(v,w)\cdot \partial_v^2
x(v,w)&=&\frac{(-f_v(v,w),-f_w(v,w),1)}{\sqrt{1+f_v^2(v,w)+f^2_w(v,w)}}\cdot
(0,0,f_{vv}(v,w))=\frac{f_{vv}(v,w)}{\sqrt{1+f_v^2+f_w^2}}\\
&=&f_{vv}(v,w)+{{{\mathcal O}}}^2,\\
\boldsymbol{n}(v,w)\cdot \partial_{vw}^2 x(v,w)&=&f_{vw}(v,w)+{{{\mathcal O}}}^2,\\
\boldsymbol{n}(v,w)\cdot \partial_w^2 x(v,w)&=&f_{ww}(v,w)+{{{\mathcal O}}}^2.\end{aligned}$$ Expanding to first order in Taylor series gives the formula for $II$ in . 

The Weingarten map (or shape operator), again expressing curvature, has components $$W_i^j=\sum\limits_{k=1}^2 II_{ik}g^{kj},\quad i,j=1,2,$$ where $\Big(g^{kj}\Big)_{k,j=1,2}$ is the inverse matrix of $g$. From (\[1.5\]) it follows that $g^{-1}={{\rm Id}}+{{{\mathcal O}}}^2$ and hence $W=II\cdot g^{-1}=II+{{{\mathcal O}}}^2$, so is proved. 

The gradient of $u$ in the surface is $${\nabla_\partial}u=\left(
\begin{array}{cc}
\displaystyle{g^{11}}&\displaystyle{g^{12}}\\
\displaystyle{g^{21}}&\displaystyle{g^{22}}
\end{array}
\right)\left(
\begin{array}{c}
u_v\\
u_w
\end{array}
\right)=\Big( {{\rm Id}}+{{{\mathcal O}}}^2\Big)\left( \begin{array}{c}
u_v\\
u_w
\end{array}\right)=\left( \begin{array}{c}
u_v\\
u_w
\end{array}\right)+{{{\mathcal O}}}^2.$$

Next, $g={{\rm Id}}+{{\mathcal O}}^2$ implies $\sqrt{\det g}=1+{{\mathcal O}}^2$, so for a vector field $Z=(Z_v, Z_w)$ the divergence in the surface is $$\begin{aligned}
{{\operatorname{div}}_\partial}Z&=&\frac{1}{\sqrt{\text{det}\, g}}\Big(\partial_v(\sqrt{\text{det}\,g}Z_v)+\partial_w(\sqrt{\text{det}\,g}Z_w)\Big)\\
&=&\partial_v Z_v+\partial_w Z_w+{{{\mathcal O}}}^1,\end{aligned}$$ and this gives also the result for the Laplace-Beltrami operator, ${\Delta_\partial}={{\operatorname{div}}_\partial}{\nabla_\partial}$.

We now prove the formulas for ${{\boldsymbol{\dot{n}}}},{{\boldsymbol{\ddot{n}}}}$ in Proposition \[prop:deriv at dO\].  We parametrize $\partial\Omega_h$ by $$Y(v,w,h)=x(v,w)+ha(v,w)\boldsymbol{n}(v,w)$$ with $x(v,w)$ as in .  For $p=0$ we have $q_h=Y(0,0,h)$, and $N(h)=Y_v(0,0,h)\times Y_w(0,0,h)$ is a non-normalized normal to ${{\partial\Omega}}_h$ at this point.  From we have $Y_v(0,0,h)=\big(1-a\lambda h,0,h a_v\big) $ and analogously $Y_w(0,0,h)=\big(0,1-a\mu h,h a_w\big)$, so (with $a$ and its derivatives taken at $(0,0)$) $$\begin{aligned}
N(h) &= (-a_v h+aa_v\mu h^2,\, -a_w h+aa_w\lambda  h^2,
1-a(\mu+\lambda)h+a^2\lambda\mu h^2) \\
&=A+Bh+Ch^2\end{aligned}$$ where $A=(0,0,1),\ B=(-a_v,-a_w,-2H a),\ C=(a a_v\mu,aa_w\lambda,a^2 K)$. Expanding ${{\boldsymbol{n}}}(q_h,h)=\frac{N(h)}{||N(h)||}$ in Taylor series results in ${{\boldsymbol{n}}}(q_h,h) =$ $$A+\Big(B-A(A\cdot B)\Big)h+\left( C-B(A\cdot B)-A\Big(\frac{||B||^2}{2}+A\cdot C-\frac{3}{2}(A\cdot B)^2\Big)\right)h^2+O(h^3).$$ The coefficient of $h$ is $(-a_v,-a_w,0)=-{\nabla_\partial}a $ by , and the coefficient of $h^2/2$ is $(-2aa_v\lambda,-2a a_w\mu,-{|\nabla a|^2}) =  -W({\nabla_\partial}a^2)-|{\nabla_\partial}a|^2\boldsymbol{n}$ by and , which proves the first part of Proposition \[prop:deriv at dO\]. 

To prove the second part of Proposition \[prop:deriv at dO\] it is convenient to use coordinates normal with respect to ${{\partial\Omega}}$, near $p=0$. That is, we parametrize a point $X=(X_1,X_2,X_3)$ near $p$ by coordinates $(v,w,s)$ where $s$ is the signed distance from $X$ to ${{\partial\Omega}}$ and $x(v,w)$ is the point on ${{\partial\Omega}}$ closest to $X$. That is, $$X(v,w,s)=x(v,w)+s\boldsymbol{n}(v,w).$$

\[lemma1\] The derivatives in the coordinates $X_1, X_2, X_3$ and $v,w,s$ at the point $p=0$ are related as follows: $$\begin{aligned}
&&\partial_{X_1}=\partial_v,\quad \partial_{X_2}=\partial_w,\quad
\partial_{X_3}=\partial_s;\nonumber\\
&&\partial^2_{X_1}=\partial_v^2-\lambda\partial_s,\quad
\partial_{X_2}^2=\partial_w^2-\mu\partial_s,\quad
\partial^2_{X_3}=\partial_s^2;\nonumber\\
&&\partial_{X_1}\partial_{X_3}=\partial_v\partial_s+\lambda\partial_v,\quad
\partial_{X_2}\partial_{X_3}=\partial_w\partial_s+\mu\partial_w;\label{1.4}\\
&&\partial^2_{X_1}\partial_{X_3}=\partial^2_v\partial_s-\lambda
\partial_s^2+2\lambda\partial^2_v-\lambda^2\partial_s+f_{vvv}\cdot
\partial_v+f_{vvw}\cdot \partial_w;\nonumber\\
&&\partial^2_{X_2}\partial_{X_3}=\partial^2_w\partial_s-\mu
\partial_s^2+2\mu\partial^2_w-\mu^2\partial_s+f_{www}\cdot
\partial_w+f_{wwv}\cdot \partial_v\nonumber.\end{aligned}$$

By differentiating and using the chain rule repeatedly, one gets $$\begin{aligned}
&&\partial_s=\frac{\partial X}{\partial x}\partial_X = \boldsymbol{n}\cdot \partial_X, \quad \partial_s^2=\boldsymbol{n}\cdot\boldsymbol{n}\cdot
\partial_X^2;\\
&&\partial_v=X_v\cdot \partial_X,\quad
\partial_v^2=\partial_v(X_v\cdot \partial_X)=X_v\cdot X_v\cdot
\partial_X^2+X_{vv}\cdot \partial_X;\\
&&\partial_v \partial_s=\boldsymbol{n}\cdot X_v\cdot
\partial_X^2+\boldsymbol{n}_v\cdot \partial_X,\\
&&\partial_v^2\partial_s=\boldsymbol{n}\cdot
X_v\cdot X_v\cdot
\partial_X^3+\boldsymbol{n}\cdot X_{vv}\cdot
\partial^2_X+2\boldsymbol{n}_v\cdot X_v\cdot
\partial_X^2+\boldsymbol{n}_{vv}\cdot\partial_X,\end{aligned}$$ and analogous formulas for the $w$-derivatives[^1]. 

We now evaluate all quantities at $p=0$, i.e. $v=w=s=0$. Using and $
X_{v}(0,0,0)=(1,0,0),\ X_{vv}(0,0,0)=(0,0,\lambda)
$ we get at $p$ $$\begin{aligned}
&&\partial_s=\partial_{X_3},\quad
\partial_s^2=\partial_{X_3}^2;\\
&&\partial_v=\partial_{X_1},\quad
\partial_v^2=\partial^2_{X_1}+\lambda\partial_{X_3};\\
&&\partial_v
\partial_s=\partial_{X_1}\partial_{X_3}-\lambda\partial_{X_1},\\
&&\partial_v^2 \partial_s=\partial^2_{X_1}\partial_{X_3}+\lambda
\partial^2_{X_3}-2\lambda\partial^2_{X_1}-f_{vvv}\partial_{X_1}-f_{vvw}\partial_{X_2}-\lambda^2
\partial_{X_3}.\end{aligned}$$

Analogous formulas hold with $v$ replaced by $w$, $\lambda$ by $\mu$ and with $\partial_{X_1},\partial_{X_2}$ interchanged.  Solving these equations for the derivatives $\partial_{X_i}$ one gets the desired relations. 

\[lemma2\] Let $u$ be harmonic on $\Omega$, smooth up to ${{\partial\Omega}}$. Then $$\begin{aligned}
({{\boldsymbol{n}}}\cdot\nabla)^2 u&=&N_1 u + N_1' {\partial_n}u\label{1.12},\\
({{\boldsymbol{n}}}\cdot\nabla)^3 u&=&N_2 u + N_2' {\partial_n}u \label{1.13}\end{aligned}$$ at $p=0$, with $N_i, N_i'$ as given in Proposition \[prop:deriv at dO\]. 

Here $\left[({{\boldsymbol{n}}}\cdot\nabla)^2 u\right](p)$ is understood as $\left({{\boldsymbol{n}}}(p)\nabla\right)^2 u$, evaluated at $p$; that is, ${{\boldsymbol{n}}}$ is not differentiated[^2].  In the coordinates above, we have at $p=0$ $$({{\boldsymbol{n}}}\cdot\nabla)^2 u=\partial^2_{X_3}u,\quad ({{\boldsymbol{n}}}\cdot\nabla)^3 u =
\partial^3_{X_3}u$$

A short calculation, using (\[1.4\]), and $\partial_s=\partial_{X_3}=\partial_{n}$, gives, at $0$, $$\begin{aligned}
\partial^2_{X_3}u &=&-(\partial^2_{X_1}+\partial^2_{X_2})u\\
&=&-(\partial_v^2-\lambda \partial_s)u-(\partial_w^2-\mu\partial_s)u\\
&=&-(\partial_v^2+\partial_w^2)u+(\lambda+\mu)\partial_s u\\
&=&-{\Delta_\partial}u+2H\partial_{n} u.\end{aligned}$$ Next we evaluate $\partial^3_{X_3}u$ at $0$. By the last formulas in (\[1.4\]) we get $$\begin{aligned}
\partial^3_{X_3}u &=& -(\partial^2_{X_1}+\partial^2_{X_2})\partial_{X_3}u\\
&=&-(\partial_v^2 \partial_s+\partial_w^2 \partial_s)u+(\lambda+\mu)\partial^2_s u\\
&-&2(\lambda\partial_v^2+\mu\partial_w^2)u+(\lambda^2+\mu^2)\partial_s u\\
&-&(\underbrace{f_{vvv}\partial_v+f_{www}\partial_w+f_{vvw}\partial_w+f_{wwv}\partial_v}_{\nabla(\Delta f)\cdot\nabla})u.\end{aligned}$$

Using , (\[1.10\]), (\[1.11\]) and then $\partial_s=\partial_{n}$ (everywhere), $\partial^2_s u =\partial^2_{X_3}u=(-{\Delta_\partial}+2H\partial_{n})u$ (at $p$ only) we obtain $$\begin{aligned}
\partial^3_{X_3}u &=-{\Delta_\partial}\partial_s u+2H\partial_s^2 u+(4H^2-2K)\partial_s u
-2{{\operatorname{div}}_\partial}(W({\nabla_\partial}u))+2{\nabla_\partial}H\cdot {\nabla_\partial}u\\
&= \Big(-{\Delta_\partial}+8H^2-2K\Big)\partial_{n}u+
2\Big(-H{\Delta_\partial}-{{\operatorname{div}}_\partial}(W({\nabla_\partial}))+{\nabla_\partial}H\cdot{\nabla_\partial}\Big)u\end{aligned}$$ as claimed. 

Now we prove the second part of Proposition \[prop:deriv at dO\].  We have $\frac{d}{dh} \left[(\nabla u)(q_h,h)\right] = (\dot{q}_h\cdot\nabla)(\nabla u)(q_h,h) + (\nabla {\dot{u}})(q_h,h)$. With $q_h=p+ha(p){{\boldsymbol{n}}}(p)$ we get $\dot{q}_h=a(p){{\boldsymbol{n}}}(p)$. At $p=0$, ${{\boldsymbol{n}}}(p)\cdot\nabla = \partial_{X_3}$. So we get, using (\[1.4\]) and (\[1.12\]), with everything evaluated at $p=0$, $h=0$, $$\begin{aligned}
\big(\nabla u\big)^{\cdot}&=&a\big(\boldsymbol{n}\cdot \nabla\big)\nabla u+\nabla \dot{u}=
a\left(\begin{array}{c}\partial_{X_3}\partial_{X_1}u\\
\partial_{X_3}\partial_{X_2}u\\
\partial^2_{X_3}u\end{array}\right)+\nabla \dot{u}\\
&=&a\left\{\left(\begin{array}{c}
\partial_v \partial_n u\\
\partial_w \partial_n u\\
0
\end{array} \right)+\left(\begin{array}{c}
\lambda \partial_v u\\
\mu \partial_w u\\
0
\end{array}\right)+\left(\begin{array}{c}
0\\
0\\
N_1 u + N_1' {\partial_n}u
\end{array}\right)\right\}+\nabla\dot{u}\\\end{aligned}$$ The first term is $a{\nabla_\partial}\partial_n u$ by , and the second is $aW({\nabla_\partial}u)$ by , . 

Finally, $$\big(\nabla u\big)^{\cdot\cdot}=a^2 \big(\boldsymbol{n}\cdot \nabla\big)^2 \nabla u +2a\big(\boldsymbol{n}\cdot \nabla\big)\nabla \dot{u}+\nabla\ddot{u}$$ and the normal component of the first term is $a^2(\partial_{X_3}^3 u) {{\boldsymbol{n}}}= a^2(N_2 u + N_2' {\partial_n}u){{\boldsymbol{n}}}$ by (\[1.13\]) while the last two terms can be rewritten as $2(\nabla{\dot{u}})^{\cdot} - \nabla {\ddot{u}}$. This finishes the proof of Proposition \[prop:deriv at dO\]. 

\[lem:partial n u dot\] Let $u$ be as in Proposition \[prop:deriv at dO\].  The derivatives of $h\mapsto u(q_h,h)$ at $h=0$ are $$\label{eq:udot uddot}
(u)^{\cdot} = a\,{\partial_n}u + {\dot{u}},\qquad (u)^{\cdot\cdot} = -a^2{\Delta_\partial}u + 2a^2 H {\partial_n}u + 2a\, {\partial_n}{\dot{u}}+ {\ddot{u}}$$ The derivatives of $h\mapsto(\partial_n u)(q_h,h)$ at $h=0$ are $$\begin{aligned}
\label{2.3}
\big(\partial_n u\big)^{\cdot} &= P_1(u) + Q_1(\partial_n u) + \partial_n\dot{u}\\
\label{2.7}
\big(\partial_n u\big)^{\cdot\cdot} &=P_2(u)+Q_2\big(\partial_n u\big)+2P_1(\dot{u})+2Q_1 (\partial_n \dot{u}) +\partial_n \ddot{u}\end{aligned}$$ with $$\begin{aligned}
 P_1&=&-{{\operatorname{div}}_\partial}a{\nabla_\partial},\\
 Q_1&=&2aH,\\
 P_2&=&-2\Big(a^2 H {\Delta_\partial}+
{{\operatorname{div}}_\partial}(a^2 W) {\nabla_\partial}- a^2({\nabla_\partial}H)\cdot{\nabla_\partial}\Big),\\
 Q_2&=&-{{\operatorname{div}}_\partial}a^2{\nabla_\partial}+ a^2(8H^2-2K)-|{\nabla_\partial}a|^2,\\
 \end{aligned}$$

The first equation of is obvious.  The second one follows from $\left(u(q_h,h)\right)^{\cdot\cdot}= a^2 ({{\boldsymbol{n}}}\cdot\nabla)^2 u  + 2a({{\boldsymbol{n}}}\cdot\nabla) {\dot{u}}+ {\ddot{u}}$ and . 

Equation follows from $ \big(\partial_n u\big)^{\cdot}=\big(\boldsymbol{n}(q_h,h)\cdot\nabla u(q_h,h)\big)^{\cdot}=\dot{\boldsymbol{n}}\cdot\nabla u+
 \boldsymbol{n}\cdot \big(\nabla u\big)^{\cdot},
 $ and Proposition \[prop:deriv at dO\]: The terms acting on $u$ are $$-({\nabla_\partial}a)\cdot{\nabla_\partial}+ a N_1 = -({\nabla_\partial}a)\cdot{\nabla_\partial}- a{\Delta_\partial}= -{{\operatorname{div}}_\partial}a {\nabla_\partial}=P_1$$ and there is one term acting on ${\partial_n}u$, which is $a N_1'=Q_1$.

For we write $
 \big(\partial_n u\big)^{\cdot\cdot}=\big( \boldsymbol{n}(q_h,h)\cdot \nabla u(q_h,h)\big)^{\cdot\cdot} = \ddot{\boldsymbol{n}}\cdot \nabla u+2\dot{\boldsymbol{n}}\cdot \big(\nabla u\big)^{\cdot}+\boldsymbol{n}\cdot \big(\nabla u\big)^{\cdot\cdot}
$, use Proposition \[prop:deriv at dO\] and sort by terms involving $u$, ${\partial_n}u$, ${\dot{u}}$, ${\partial_n}{\dot{u}}$. The terms acting on $u$ are $$\begin{aligned}
 & -W({\nabla_\partial}a^2)\cdot{\nabla_\partial}-2({\nabla_\partial}a) a T_1 + a^2 N_2  \\
=& -({\nabla_\partial}a^2) \cdot W {\nabla_\partial}- ({\nabla_\partial}a^2)\cdot W{\nabla_\partial}+a^2 N_2\qquad\text{ (since $W$ is selfadjoint)} \\
=& -2 \left[ ({\nabla_\partial}a^2)\cdot W{\nabla_\partial}+ a^2{{\operatorname{div}}_\partial}W {\nabla_\partial}+ a^2H{\Delta_\partial}- a^2({\nabla_\partial}H){\nabla_\partial}\right]= P_2\end{aligned}$$ The terms acting on ${\partial_n}u$ are $$\begin{aligned}
& -|{\nabla_\partial}a|^2 - 2({\nabla_\partial}a)aT_1'+ a^2 N_2' \\
=&
-|{\nabla_\partial}a|^2 - ({\nabla_\partial}a^2)\cdot{\nabla_\partial}- a^2{\Delta_\partial}+ a^2(8H^2-2K) = Q_2\end{aligned}$$ The terms involving ${\dot{u}}$ and ${\partial_n}{\dot{u}}$ arise from one of the $h$-derivatives in $\big( \boldsymbol{n}(q_h,h)\cdot \nabla u(q_h,h)\big)^{\cdot\cdot}$ falling on the second argument of $u$, hence are twice the corresponding terms in $({\partial_n}{\dot{u}})^{\cdot}$.

### End of proof of the perturbation formulas

Now let $u(\cdot,h)$, ${\varepsilon}(h)$ be a solution of - for $\Omega_h$, for each $h$, smoothly depending on all variables (up to the boundary of $\Omega_h$ from either side). 

We first prove - and calculate $F_1,G_1,F_2,G_2$. If $x\not\in{{\partial\Omega}}$ then $x\not\in \partial\Omega_h$ for sufficiently small $h$, and , follow.  The boundary conditions are $$\begin{aligned}
\label{eq:u(h) eqn}
u_-(q_h,h) - u_+(q_h,h) &=0 \\
\label{eq:dn u(h) eqn}
{\varepsilon}(h)\partial_n u_- (q_h,h) + \partial_n u_+ (q_h,h) &=0\end{aligned}$$ where $q_h= p +ha(p) {{\boldsymbol{n}}}(p)$, for each $p\in{{\partial\Omega}}$ and $h$. We differentiate and use the first equation in for $u_-$ and $u_+$, at $h=0$. This gives ${\dot{u}}_--{\dot{u}}_+=-a({\partial_n}u_- - {\partial_n}u_+)$. Using at $h=0$ we obtain with $$\label{eq:F}
 F_1 = -  ({\varepsilon}+1) a{\partial_n}u_-$$ Next, we differentiate at $h=0$ and use for $u_-$ and $u_+$. We obtain $$\begin{aligned}
0 &= {\dot{{\varepsilon}}}{\partial_n}u_- + {\varepsilon}({\partial_n}u_-)^{\cdot} + ({\partial_n}u_+)^{\cdot}  \\
&=
{\dot{{\varepsilon}}}{\partial_n}u_- + P_1 ({\varepsilon}u_- + u_+) + Q_1 ({\varepsilon}{\partial_n}u_- + {\partial_n}u_+)
+ {\varepsilon}{\partial_n}{\dot{u}}_- + {\partial_n}{\dot{u}}_+ \\
&= {\dot{{\varepsilon}}}{\partial_n}u_-  + ({\varepsilon}+1) P_1 u_- + {\varepsilon}{\partial_n}{\dot{u}}_- + {\partial_n}{\dot{u}}_+\end{aligned}$$ using , at $h=0$, hence with $$\label{eq:G}
 G_1 = -({\varepsilon}+ 1) P_1 u_-$$

Differentiating twice and using for $u_-$ and $u_+$ we obtain $$\begin{aligned}
0
&= \left(u_-\right)^{\cdot\cdot} - \left(u_+\right)^{\cdot\cdot} \\
&= - a^2 {\Delta_\partial}(u_- - u_+) + a Q_1 ({\partial_n}u_- - {\partial_n}u_+)  + 2a ({\partial_n}{\dot{u}}_- - {\partial_n}{\dot{u}}_+) + {\ddot{u}}_- - {\ddot{u}}_+\end{aligned}$$ By , at $h=0$ the first term vanishes and the second equals $({\varepsilon}+1) a Q_1 {\partial_n}u_-$. By the third term is $({\varepsilon}+1)2a {\partial_n}{\dot{u}}_- + {\dot{{\varepsilon}}}2a{\partial_n}u_- - 2a G_1$. Rearranging and using we obtain with $$\label{eq:F_2}
 F_2 = -({\varepsilon}+ 1)D - 2a{\dot{{\varepsilon}}}{\partial_n}u_-$$ $$D = a\big( 2 P_1 u_- +  Q_1  {\partial_n}u_- + 2{\partial_n}{\dot{u}}_-  \big)$$

Finally, we differentiate twice and get, using and , $$\begin{aligned}
0
&=  \ddot{\varepsilon}\partial_{n} u_{-}+2\dot{\varepsilon}(\partial_{n} u_{-})^\cdot+\varepsilon (\partial_{n} u_{-})^
 {\cdot\cdot}+(\partial_{n} u_{+})^{\cdot\cdot}  \\
&= {\ddot{{\varepsilon}}}{\partial_n}u_- + 2{\dot{{\varepsilon}}}(P_1 u_- + Q_1 {\partial_n}u_- + {\partial_n}{\dot{u}}_-) \\
&+ P_2 ({\varepsilon}u_- + u_+) + Q_2 ({\varepsilon}{\partial_n}u_- + {\partial_n}u_+)
+ 2 P_1 ({\varepsilon}{\dot{u}}_- + {\dot{u}}_+) + 2 Q_1 ({\varepsilon}{\partial_n}{\dot{u}}_- + {\partial_n}{\dot{u}}_+) \\
&+ {\varepsilon}{\partial_n}{\ddot{u}}_- + {\partial_n}{\ddot{u}}_+  \\
\intertext{The terms in the third line of this equation are, using \eqref{eq:u(h) eqn}, \eqref{eq:dn u(h) eqn} at $h=0$ and \eqref{eqn2a}, \eqref{eqn3a},}
& ({\varepsilon}+1) P_2 u_- + 0 + 2P_1 (({\varepsilon}+1) {\dot{u}}_- - F_1) +
2 Q_1 (-{\dot{{\varepsilon}}}{\partial_n}u_- + G_1) \\
\intertext{
Using $F_1 = -  ({\varepsilon}+1) a{\partial_n}u_-$, $G_1 = -({\varepsilon}+ 1) P_1 u_-$ we obtain }
0 &= {\ddot{{\varepsilon}}}{\partial_n}u_- + 2{\dot{{\varepsilon}}}A + ({\varepsilon}+1)B +  {\varepsilon}{\partial_n}{\ddot{u}}_- + {\partial_n}{\ddot{u}}_+\end{aligned}$$ with $$\begin{aligned}
A &= P_1 u_- + {\partial_n}{\dot{u}}_- \\
B &= (P_2  - 2 Q_1 P_1) u_- + 2 P_1 a {\partial_n}u_- + 2 P_1 {\dot{u}}_-\end{aligned}$$ Rearranging, we get with $$\label{eq:G_2}
G_2 = -2{\dot{{\varepsilon}}}A - ({\varepsilon}+1) B$$ Before proceeding we simplify $B$. Note that $$\begin{aligned}
P_2 - 2Q_1 P_1
&=&-2a^2 H{\Delta_\partial}+2a^2 \big({\nabla_\partial}H\big)\cdot {\nabla_\partial}-2{{\operatorname{div}}_\partial}(a^2 W) {\nabla_\partial}+4aH\big( a {\Delta_\partial}+({\nabla_\partial}a)\cdot {\nabla_\partial}\big)\\
&=&-2{{\operatorname{div}}_\partial}\big(a^2 W\big){\nabla_\partial}+2{{\operatorname{div}}_\partial}\big(a^2 H\big){\nabla_\partial}\\
&=&-2{{\operatorname{div}}_\partial}\big(a^2 W_0){\nabla_\partial}\end{aligned}$$ where $W_0=W-HI$. Thus $$B=-2{{\operatorname{div}}_\partial}\big(a^2 W_0\big){\nabla_\partial}u_- + 2 P_1\big( a {\partial_n}u_- + \dot{u}_-\big)$$

We have proved .  In order to calculate ${\dot{{\varepsilon}}}$, ${\ddot{{\varepsilon}}}$, we use the following lemma. 

\[lem:compat\] Let ${\varepsilon}$ be a plasmonic eigenvalue with eigenspace $E$. Given functions $F,{\tilde{G}}$ on ${{\partial\Omega}}$, the inhomogeneous system $$\begin{aligned}
\Delta v & = 0 \quad\text{in }{\mathbb{R}}^3\setminus\partial\Omega
\\
v_- - v_+ & = F \quad\text{on }\partial\Omega
\\
{\varepsilon}\partial_n v_- + \partial_n v_+ & = {\tilde{G}}\quad\text{on } \partial\Omega\end{aligned}$$ has a solution $v$ satisfying the decay condition iff $$\langle {\tilde{G}},w_{-} \rangle = {\varepsilon}\langle F,\partial_n w_{-} \rangle$$ for all $w\in E$. Here $\langle\cdot,\cdot\rangle$ is the scalar product on $L^2({{\partial\Omega}})$.

We reformulate this in terms of Dirichlet-Neumann operators, compare Proposition \[prop:dirichlet-neumann\].  Let $g=v_{+}$ at ${{\partial\Omega}}$. Then the second equation says $v_- = F+g$, and the first and third that ${\varepsilon}{{\mathcal N}}_- (F+g) + {{\mathcal N}}_+ g = {\tilde{G}}$. Therefore, a solution $v$ exists if and only if the equation $$({\varepsilon}{{\mathcal N}}_- + {{\mathcal N}}_+) g = {\tilde{G}}-{\varepsilon}{{\mathcal N}}_- F$$ has a solution $g$. Since the operator ${\varepsilon}{{\mathcal N}}_- + {{\mathcal N}}_+$ is selfadjoint, this is equivalent to $\langle {\tilde{G}}-{\varepsilon}{{\mathcal N}}_- F,w_{-}\rangle = 0$ for all $w\in E$. This is equivalent to the claim. 

We apply the lemma with $v={\dot{u}}$ and $w=u$ and $F=F_1$, see , ${\tilde{G}}=-{\dot{{\varepsilon}}}{\partial_n}u_- + G_1$, see .  Using the assumption $\langle u_- ,{\partial_n}u_- \rangle = \|u_-\|_-^2 = 1$ we get $${\dot{{\varepsilon}}}= \langle G_1,u_- \rangle - {\varepsilon}\langle F_1, {\partial_n}u_- \rangle$$ Integrating by parts $\langle G_1,u_- \rangle = ({\varepsilon}+1) \langle {{\operatorname{div}}_\partial}a {\nabla_\partial}u_-, u_- \rangle = -({\varepsilon}+1) \langle a{\nabla_\partial}u_-,{\nabla_\partial}u_-\rangle$ we obtain . 

Now we apply the lemma with $v={\ddot{u}}$ and $w=u$ and $F=F_2$, see , ${\tilde{G}}=-{\ddot{{\varepsilon}}}{\partial_n}u_- + G_2$, see .  We obtain $${\ddot{{\varepsilon}}}= \langle G_2,u_- \rangle - {\varepsilon}\langle F_2, {\partial_n}u_- \rangle$$ Now we sort by like terms.  For example, the terms involving $u_-,u_-$ are $$-2{\dot{{\varepsilon}}}\langle P_1 u_-,u_- \rangle - ({\varepsilon}+1)(-2) \langle {{\operatorname{div}}_\partial}a^2 W_0 {\nabla_\partial}u_-,u_- \rangle$$ where $P_1 = -{{\operatorname{div}}_\partial}a{\nabla_\partial}$. Integrating by parts yields the first line in , and the other terms are obtained in a similar way. 

Shape derivative of Dirichlet-Neumann operator
----------------------------------------------

As a side result of our calculations above we can calculate the derivative of the Dirichlet-Neumann operator ${{\mathcal N}}:={{\mathcal N}}_-$ when changing the domain as in .  Since the domain ${{\partial\Omega}}_h$ varies, we have to state precisely what we mean by this.  Recall from Section \[subsec:analytic\] the diffeomorphism $T_h:\Omega\to\Omega_h,\ x\mapsto x+h V(x)$ where $V=a{{\boldsymbol{n}}}$ on ${{\partial\Omega}}$. This induces an isomorphism $T_h^*:{C^\infty}({{\partial\Omega}}_h) \to {C^\infty}({{\partial\Omega}}),\ g\mapsto T_h^*g = g\circ T_h$. Then we consider $${{\mathcal N}}_h := T_h^* \circ {{\mathcal N}}_{\Omega_h} \circ (T_h^*)^{-1} : {C^\infty}({{\partial\Omega}}) \to {C^\infty}({{\partial\Omega}})$$ where ${{\mathcal N}}_{\Omega_h}:{C^\infty}({{\partial\Omega}}_h) \to {C^\infty}({{\partial\Omega}}_h)$ is the Dirichlet-Neumann operator for $\Omega_h$. The operator ${{\mathcal N}}_h$ is the same as the one introduced in Section \[subsec:analytic\].  Explicitly, it is given as follows: Let $g\in{C^\infty}({{\partial\Omega}})$. Let $u(\cdot,h)$ be the solution of $$\begin{aligned}
\Delta u(\cdot,h) &= 0 \quad\text{ in }\Omega_h \label{eq:DN1}\\
u(T_h(x),h) &= g(x) \quad\text{ for all }x\in {{\partial\Omega}}\label{eq:DN2}\\
\text{and then set } ({{\mathcal N}}_h g)(x) &= ({\partial_n}u) (T_h(x),h).  \label{eq:DN3}\end{aligned}$$

\[thm:DN-perturb\] The derivative at $h=0$ of the familiy of Dirichlet-Neumann operators defined above is given by $$\label{eq:DN perturb}
{\dot{{{\mathcal N}}}}= -{{\operatorname{div}}_\partial}a \nabla + 2aH {{\mathcal N}}- {{\mathcal N}}a {{\mathcal N}}$$

Differentiate , , at $h=0$ and use (where $q_h=T_h(p)$) for the right side of .  This gives $$\begin{aligned}
\Delta {\dot{u}}&= 0 \quad\text{ in }\Omega\\
a {\partial_n}u + {\dot{u}}&= 0 \quad\text{ at } {{\partial\Omega}}\\
{\dot{{{\mathcal N}}}}g &= -{{\operatorname{div}}_\partial}a {\nabla_\partial}u + 2aH {\partial_n}u + {\partial_n}{\dot{u}}\end{aligned}$$ Rewrite the second equation as ${\dot{u}}= -a{\partial_n}u$ at ${{\partial\Omega}}$, then together with the first equation this gives ${\partial_n}{\dot{u}}= {{\mathcal N}}(-a{\partial_n}u) = -{{\mathcal N}}(a {{\mathcal N}}g)$, and then the third equation gives the claim. 

For the exterior Dirichlet-Neumann operator ${{\mathcal N}}_+$ one has the same formula, replacing ${{\mathcal N}}$ by ${{\mathcal N}}_+$ everywhere.  Same proof. 

Using the second derivative formulas in Section \[subsec:perturb calc\] one could also derive a formula for the second derivative of ${{\mathcal N}}_h$ in $h$.

Note that the principal symbol of the first and third term in is $a|\xi|^2$, hence ${\dot{{{\mathcal N}}}}$ is a pseudodifferential operator of order at most one, which is to be expected since each ${{\mathcal N}}_h$ is a pseudodifferential operator of order one. 

Splitting at first order
------------------------

We now prove the last statement of Theorem \[thm:perturbation\], which gives also another (though similar) proof of the formula for ${\dot{{\varepsilon}}}$. We use the characterization in terms of Dirichlet-Neumann operators, Proposition \[prop:dirichlet-neumann\], where we use the operators ${{\mathcal N}}_h$ transplanted from $\Omega_h$ to $\Omega$ as in the previous section.  If $g(x,h) = u(T_h(x),h)$ for $x\in{{\partial\Omega}}$ then $({\varepsilon}(h){{\mathcal N}}_{-,h} + {{\mathcal N}}_{+,h}) g(h) = 0$ for each $h$. Differentiating in $h$ at $h=0$ gives $$({\dot{{\varepsilon}}}{{\mathcal N}}_- + {\varepsilon}{\dot{{{\mathcal N}}}}_- + {\dot{{{\mathcal N}}}}_+)g + ({\varepsilon}{{\mathcal N}}_- + {{\mathcal N}}_+){\dot{g}}=0$$ We take the scalar product of this with any $g'\in E_\partial :=\{u_{|{{\partial\Omega}}}:u\in E\}$. Using selfadjointness of ${{\mathcal N}}_\pm$ and $({\varepsilon}{{\mathcal N}}_- + {{\mathcal N}}_+)g'=0$ we get $$\langle ({\dot{{\varepsilon}}}{{\mathcal N}}_- + {\varepsilon}{\dot{{{\mathcal N}}}}_- + {\dot{{{\mathcal N}}}}_+)g, g' \rangle = 0
\quad \forall g'\in E_\partial$$ or equivalently $$- \langle ({\varepsilon}{\dot{{{\mathcal N}}}}_- + {\dot{{{\mathcal N}}}}_+)g, g' \rangle = {\dot{{\varepsilon}}}(g,g')_-$$ A simple calculation using for ${{\mathcal N}}={{\mathcal N}}_\pm$ reveals that the left hand side equals $q_1(u,u')$, where $u'$ is the harmonic function on $\Omega$ with boundary values $g'$, and the claim follows. 

Examples {#sec:examples}
========

Planar disk {#subsec:disk}
-----------

For the disk $\Omega=\{x\in{\mathbb{R}}^2:\, |x| < 1\}$ the only plasmonic eigenvalue is ${\varepsilon}=1$, and any harmonic function in $\Omega$ whose value and normal derivative at ${{\partial\Omega}}$ are defined (e.g. any function in $H^s(\Omega)$ with $s>3/2$) is a plasmon.  This follows from the fact that if $u_-$ is harmonic in $\Omega$ then its Kelvin transform $u_+ (x) = u_-(\frac x{|x|^2})$ is harmonic in $|x|>1$. Clearly, $u_-$ and $u_+$ have the same boundary values at $|x|=1$ and their normal derivatives add to zero, so ${\varepsilon}=1$. The orthogonality in Theorem \[thm:existence+asymp\] then shows that there can be no other solutions.  Choosing $u$ with $u_{|{{\partial\Omega}}}$ non-smooth we see that the regularity statement in Theorem \[thm:regularity\] does not hold. 

Perturbed ball {#subsec:ex ball}
--------------

For the unit ball, $\Omega=\{x\in{\mathbb{R}}^3:\, |x| < 1\}$, the solutions of the plasmonic eigenvalue problem are given by $$u_{-,k} (r,\omega) = r^k Y_k (\omega),\ u_{+,k} (r,\omega) = r^{-k-1} Y_k (\omega),\ {\varepsilon}_k = \frac{k+1}k$$ for $k=0,1,\dots$. Here $r>0$, $\omega$ with $|\omega|=1$ are polar coordinates (i.e. $x=r\omega$, $r=|x|$, $\omega=x/|x|$).  The function $Y_k$ varies over the space of spherical harmonics of degree $k$, i.e. $r^k Y_k(\omega)$ is a $k$-homogeneous harmonic polynomial on ${\mathbb{R}}^3$. Therefore, ${\varepsilon}_k$ has multiplicity $(2k+1)$.

It is clear that these functions are solutions, and since the $Y_k$ span $L^2(S^2)$, there can be no further solutions. 

Since the problem is scale invariant, one gets the same ${\varepsilon}$ for spheres of arbitrary radius.  This implies that for $a\equiv 1$ we should get ${\dot{{\varepsilon}}}_k = 0$. As a check on our calculations we can see this directly:

From $-{\Delta_\partial}Y_k = k(k+1) Y_k$ we have for $u_-=u_{-,k}$ $${\int_{{\partial\Omega}}}|{\nabla_\partial}u_-|^2 = -{\int_{{\partial\Omega}}}u_-{\Delta_\partial}u_- = k(k+1) {\int_{{\partial\Omega}}}|Y_k|^2.$$ Also, $\partial_n u_- = \partial_r u_- = k Y_k$ at $r=1$, so the right hand side of equals for $a\equiv 1$ $$\left[ k(k+1) + {\varepsilon}_k k^2 \right] {\int_{{\partial\Omega}}}|Y_k|^2 = 0.$$

For arbitrary $a$ we calculate ${\dot{{\varepsilon}}}$ in case $k=1$ and $u_-(x,y,z) = z$, where ${\varepsilon}=2$. This choice of $u_-$ (out of the three-dimensional space, spanned by the functions $x,y,z$, of solutions for $k=1$) corresponds to an exterior electric field pointing in the $z$-direction, in the physical problem leading to -.  We parametrize the sphere as graph $z=\pm \sqrt{1-x^2-y^2}$ .  Since $u_-$ is 1-homogeneous, we have $\partial_r u_- = u_-=z$ at $r=1$. Then $|{\nabla_\partial}u_-|^2 = |\nabla u_-|^2 - |\partial_r u_-|^2 = 1-z^2$. So the bracket in equals $ (1-z^2) -2z^2 = 1-3z^2 = 3(x^2+y^2) -2$. The volume element of the sphere is $dx\,dy/\sqrt{1-x^2-y^2}$. To normalize $u$ we divide by $\|u\|_-^2=\int_\Omega |\nabla u_-|^2 = \int_\Omega 1 = {{\operatorname{vol}}}(\Omega) = \frac43\pi$. Finally, ${\varepsilon}+1 = 3$. Therefore, we get $$\begin{gathered}
{\dot{{\varepsilon}}}_1 = \frac9{4\pi} \int_{x^2+y^2<1} A(x,y) \frac{3(x^2+y^2)-2}{\sqrt{1-x^2-y^2}}\, dx dy ,\\
 A(x,y) := a(x,y,\sqrt{1-x^2-y^2}) + a(x,y,-\sqrt{1-x^2-y^2}) .\end{gathered}$$

Example: Half space {#subsec:ex plane}
-------------------

We always assumed $\Omega$ to be compact.  However, a similar problem can also be formulated for non-compact $\Omega$. For example, consider the half space $$\Omega = \{(x',x_n):\ x'\in {\mathbb{R}}^{n-1},\ x_n<0\}$$ We then consider the problem -, but we need to replace the decay condition .  In view of the regularity discussion it is natural to consider $H^{1/2}$ boundary data, hence we demand a uniform $H^{1/2}$ condition: there is a constant $C$ so that $$\label{eqn:L2 bounded}
 \int_{{\mathbb{R}}^{n-1}}|u(x',x_n)|^2\,dx' + 
 \int_{{\mathbb{R}}^{n-1}}|u(x',x_n)|\, |\nabla_{x'} u(x',x_n)|\,dx'
 \leq C \ \text{ for all } x_n$$ In particular, we can take the Fourier transform in the $x'$-variables, ${\hat{u}}(\xi',x_n) = \int_{{\mathbb{R}}^{n-1}} e^{-ix'\cdot \xi'} u(x',x_n)\,dx'$, and $\Delta u=0$ yields $(\partial_{x_n}^2 - |\xi'|^2){\hat{u}}= 0$ for $x_n\neq0$, with solutions ${\hat{u}}_\pm = a_\pm e^{-|\xi'| x_n} + b_\pm e^{|\xi'| x_n}$ (if $\xi'\neq0$).  Then implies $b_+=a_-=0$, and gives $b_-=a_+$, so $$u(x',x_n) = \int_{{\mathbb{R}}^{n-1}} e^{ix'\cdot\xi'} h(\xi') e^{-|\xi'|\, |x_n|} \,d\xi'$$ for a function $h$ on ${\mathbb{R}}^{n-1}$ satisfying $\int_{{\mathbb{R}}^{n-1}} |h(\xi')|^2 (1+|\xi'|)\, d\xi' < \infty$. We have $h=(2\pi)^{-n} {\hat{g}}$ for the boundary value $g(x')=u(x',0)$. Also $\partial_n u_\pm (x',0) = \mp\int e^{ix'\cdot\xi'} |\xi'| h(\xi')\, d\xi'$ (distributional Fourier transform), which shows that ${\varepsilon}=1$ for any $g$.

Therefore, ${\varepsilon}=1$ is the only plasmonic eigenvalue for the half space.  Incidentally, the explicit formulas show that the Dirichlet-Neumann operator is precisely the Fourier multiplier by $|\xi'|$, i.e. ${{\mathcal N}}_-g = (|\xi'|{\hat{g}})\check{ }$. 

[10]{} H. A. Atwater, [*The Promise of Plasmonics*]{}, Scientific American, April 2007, 56. 

W. L. Barnes, A. Dereux und T. W. Ebbesen, [*Surface plasmon subwavelength optics*]{}, Nature [**424**]{}, 824 (2003). 

S.-A. Biehs, D. Grieser, M. Holthaus, O. Huth, F. Rüting, and H. Uecker, [*Perturbation theory for plasmonic eigenvalues*]{}, Phys.  Rev. B [**80**]{}, 245405 (2009). 

Gerald B. Folland, [*Introduction to partial differential equations.  2nd ed.*]{},  Princeton University Press, 1995. 

D. Grieser and F. Rüting, [*Surface plasmon resonances of an arbitrarily shaped nanoparticle: high-frequency asymptotics via pseudo-differential operators*]{}, J. Phys.  A: Math.  Theor. [** 42**]{}, 135204 (2009). 

D. J. Griffiths, [*Introduction to Electrodynamics*]{}, 3rd ed.,  Prentice Hall, 1999. 

T. Kato, [*Perturbation theory for linear operators.*]{},  Grundlehren der Mathematischen Wissenschaften, 132.  Berlin etc.:  Springer-Verlag, 1966. 

V. V. Klimov and D. V. Guzatov, [*Strongly localized plasmon oscillations in a cluster of two metallic nanospheres and their influence on spontaeneous emission of an atom*]{}, Phys.  Rev. B [**75**]{}, 024303 (2007). 

V. V. Klimov and D. V. Guzatov, [*Plasmonic atoms and plasmonic molecules*]{}, Appl.  Phys.  A [**89**]{}, 305 (2007). 

L. D. Landau and E. M. Lifschitz, [ *Lehrbuch der theoretischen Physik.  Vol.  VIII: Elektrodynamik der Kontinua.*]{}  Akademie Verlag, 1974. 

D. Lannes, [*Well-posedness of the water-waves equations*]{}, J. AMS [**18**]{}, 605-654 (2005). 

S. A. Maier, [*Plasmonics: Fundamentals and Applications*]{}, Springer Verlag, Berlin, 2007. 

J. M. Pitarke, V.M. Silkin, E. V. Chulkov, and P. M. Echenique, [*Theory of surface plasmons and surface-plasmon polaritons*]{}, Rep. Prog.  Phys. [** 70**]{}, 1 (2007). 

R. Zia, J. A. Schuller, A. Chandran und M. L. Brongersma, [*Plasmonics: the next chip-scale technology*]{}, materials today [**9**]{}, 20 (2006). 

[^1]: For vectors $V,W$ we write $V\cdot W\cdot\partial_X^2 =\sum\limits_{i,j=1}^3 V_i W_j \partial_{X_i}\partial_{X_j}$, for short, and similar for third derivatives in $X$.

[^2]: Of course it would not make sense to differentiate it since it is only defined at ${{\partial\Omega}}$. But this is in contrast to the consideration of $({\partial_n}u)^{\cdot}$ below, where the normal depends on $h$, hence must be differentiated.