\section{\bf Introduction}

Let $\Omega\subset\mathbb R^d$ be a bounded and connected domain with
Lipschitz boundary 
$\partial\Omega$. Consider on $\Omega$ the Neumann eigenvalue problem
\begin{gather}\label{Problem:Neumann}
\begin{cases}
-\Delta f= \mu f&\mbox{ in }\Omega,\\
\partial_\nu f= 0 &\mbox{ on }\partial\Omega,
\end{cases}
\end{gather}
as well as the Steklov eigenvalue problem
\begin{gather}\label{Problem:Steklov}
\begin{cases}
\Delta u=0&\mbox{ in }\Omega,\\
\partial_\nu u=\sigma u&\mbox{ on }\partial\Omega.
\end{cases}
\end{gather}
Here $\Delta$ is the Laplacian, and $\partial_\nu$ is the outward pointing normal derivative. Both problems consist in finding the eigenvalues $\mu$ and $\sigma$ such that there exist non-trivial smooth solutions to the boundary value problems \eqref{Problem:Neumann} and \eqref{Problem:Steklov}. For both problems, the spectra form discrete unbounded sequences
$$ 0 = \mu_0 < \mu_1 \le \mu_2 \le \dotso \nearrow \infty$$
and
$$0=\sigma_0<\sigma_1\leq\sigma_2\leq\cdots\nearrow\infty,$$
where each eigenvalue is repeated according to multiplicity. The corresponding
eigenfunctions $\set{f_k}$ and $\set{u_k}$ have
natural normalisations as orthonormal bases of $\RL^2(\Omega)$ and $\RL^2(\partial\Omega)$, respectively.

\subsection{From Steklov to Neumann : heuristics}
Let us start by painting with broad brushes the relationships between the Neumann
and Steklov eigenvalue problems. They exhibit many similar features, and it is not a surprise that they do so. Indeed, in both cases the eigenvalues are those of a differential or pseudo-differential operator, namely the Laplacian and the Dirichlet-to-Neumann map, whose kernels consist of constant functions. Moreover, in both cases, the natural isoperimetric type problem consists in maximizing $\mu_k$ and $\sigma_k$ (instead of minimizing it as is usual for the Dirichlet problem). The relation
between the two boundary value problem is not solely heuristic and incidental.
Indeed, it is known from the works of Arrieta--Jim\'enez-Casas--Rodriguez-Bernal \cite{Arrieta2008} and Lamberti--Provenzano \cite{LambertiProvenzano2015,LambertiProvenzano2017} that
one can recover the Steklov problem as a limit of weighted Neumann
problems
\begin{gather}\label{Problem:Neumanndensity}
\begin{cases}
-\Delta f= \mu \rho_\varepsilon f&\mbox{ in }\Omega,\\
\partial_\nu f= 0 &\mbox{ on }\partial\Omega,
\end{cases}
\end{gather}
where $\rho_\varepsilon$ is a density function whose support converges to the boundary
as $\varepsilon \to 0$. If we are to interpret the Neumann problem as finding the
frequencies and modes of vibrations of a free boundary membrane, this means that
the Steklov problem represents the frequencies and modes of a membrane whose
mass is concentrated at the boundary. The reader should also refer to the
work of Hassannezhad--Miclo~\cite[Section 4]{HassannezhadMiclo}, where a
similar construction is used to prove a Cheeger-type inequality for Steklov
eigenvalues of a compact Riemannian manifold with boundary.

Our primary goal in this paper is to establish a link in the reverse
direction, by realizing the Neumann problem as a limit of appropriate Steklov
problems. This is achieved in two steps. The first one is to accumulate uniformly
distributed boundary elements inside the domain 
$\Omega$. {This is done by perforating the interior of the domain with small holes that are uniformly distributed}. On these new boundary components, we consider the Steklov boundary
conditions; this step is known as
\emph{the homogenisation process}. We assume that the ratio of
the radii of the holes to the distance between them is at a Cioranescu--Murat type  critical regime. Then, the eigenvalues and eigenfunctions of the Steklov problem converge to those of a \emph{dynamical eigenvalue problem}
  \begin{gather}\label{problem:homo}
    \begin{cases}
      -\Delta U=  A_{d} \beta  \Sigma U&\mbox{ in }\Omega,\\
      \partial_\nu U=\Sigma U&\mbox{ on }\partial\Omega,
    \end{cases}
  \end{gather}
  where $\beta \ge 0$ is the critical regime parameter and $A_d$ is the area of the unit sphere
  in $\mathbb R^{d}$. Its eigenvalues form a discrete unbounded sequence: $$\Sigma_{0,\beta}<\Sigma_{1,\beta}\leq\Sigma_{2,\beta}\leq\cdots\nearrow\infty,$$
and once again the functions associated to the eigenvalue $\Sigma_{0,\beta} = 0$ are constant. 

\begin{remark}
In \cite{BelowFrancois}, Joachim von Below and Gilles François studied an eigenvalue
problem that is equivalent
  to Problem~\eqref{problem:homo}, which stems from a parabolic equation
  with dynamical boundary conditions. Indeed, they study the eigenvalue
  problem
  \begin{equation}
  \label{eq:vbf}
   \begin{cases}
   - \Delta u = \lambda u & \text{in }  \Omega,\\
    \partial_\nu u = \lambda \alpha u & \text{on } \partial \Omega.
   \end{cases}
  \end{equation}
If $\lambda$ is an eigenvalue of Problem~\eqref{eq:vbf}, then $\Sigma=\alpha^{-1} \lambda$ is an eigenvalue of Problem~\eqref{problem:homo} with parameter $\beta = \frac{1}{\alpha A_d}$.
\end{remark}

The parameter $\beta$ in \eqref{problem:homo} can be interpreted as a weight on the interior of the domain,
with the boundary $\partial \Omega$ having constant weight $1$. In order to recover
the Neumann problem, the second step will therefore be to send the parameter 
$\beta$ to $\infty$, putting all the weight inside the domain. Under an
appropriate normalisation, eigenvalues and eigenfunctions of Problem
\eqref{problem:homo} converge to those of Problem \eqref{Problem:Neumann},
completing the circle for the relation between the Steklov and the Neumann
problems.
\subsection{The homogenisation process}
Consider a family of problems obtained by
removing periodically placed balls from the domain $\Omega$.
More precisely, given $0 < \varepsilon < 1$, and $\bk\in\mathbb Z^d$, define the cube
$$Q_{\bk}^\varepsilon := \varepsilon \bk+\left[-\frac{\varepsilon}{2},\frac{\varepsilon}{2}\right]^d\subset\mathbb R^d,$$
and define the set of indices
\begin{gather*}
I^\varepsilon := \set{\bk\in\mathbb Z^d\,:\,Q_{\bk}^\varepsilon \subset\Omega}.
\end{gather*}
Let $r_\varepsilon$ be  an increasing positive function of $\varepsilon$ with $r_\varepsilon < \varepsilon/2$.
For $\bk\in \mathbb Z^d$, define 
$$T_{\bk}^\varepsilon:=B\left(\varepsilon\bk,r_\varepsilon\right)\subset Q_{\bk}^\varepsilon$$
and set
\begin{equation}
 T^\varepsilon := \bigcup_{\bk \in I^\varepsilon} T^\varepsilon_\bk \subset \Omega.
\end{equation}
Consider the family of perforated domains
$$\Omega^\varepsilon=\Omega\setminus \overline{T^\varepsilon}.$$
\begin{figure}\label{figure:Omegaeps}
  \includegraphics[width=9cm]{OmegaN.pdf}
  \caption{The domain $\Omega_\varepsilon$}
\end{figure}
The Steklov eigenvalues of $\Omega^\varepsilon$ are written as
$\sigma_k^\varepsilon:=\sigma_k(\Omega^\varepsilon),$ and we write $\set{u_k^{\varepsilon}}$ for a corresponding complete sequence of eigenfunctions, normalized by
$$\int_{\partial\Omega^\varepsilon}(u_k^{\varepsilon})^2\, \mathrm{d} A=1.$$
Our first main result is the following critical regime homogenisation theorem
for the Steklov problem.
\begin{theorem}\label{thm:homogenisation}
  Suppose that $r_\varepsilon^{d-1} \varepsilon^{-d} \to \beta$ for some $\beta\in [0,\infty)$, as $\varepsilon\searrow 0$.
    Then $\sigma_k^{\varepsilon}$ converges
    to the eigenvalue $\Sigma_{k,\beta}$ of \eqref{problem:homo}.
    The functions $U_k^{\varepsilon}\in \RH^1(\Omega)$ obtained by harmonic extension of a normalized Steklov
    eigenfunction $u_k^{\varepsilon}$ over the holes $T^\varepsilon\subset\Omega$ 
    form a sequence which weakly converges in $\RH^1(\Omega)$ to a solution $U_k$ associated with
    $\Sigma_{k,\beta}$ of \eqref{problem:homo}.
\end{theorem}

\begin{remark} \label{rem:multiple}
  If the eigenvalue $\Sigma_k=\Sigma_{k,\beta}$ is multiple of multiplicity $m$, \emph{i.e.}
  $$\Sigma_{k-1} < \Sigma_k = \dotso = \Sigma_{k+m-1} <
  \Sigma_{k+m},$$
  the convergence statement in the previous theorem is understood in the following
  sense. Given a basis $U_k,\dotsc,U_{k+m-1}$ for the eigenspace associated with $\Sigma_k$, there is a family of $m\times m$ orthogonal matrices $M(\varepsilon)$ such that
  \begin{equation}
   M(\varepsilon) \begin{pmatrix}
            U_k^\varepsilon \\
            \vdots \\
            U_{k+m-1}^\varepsilon
           \end{pmatrix}
\longrightarrow \begin{pmatrix}
            U_k \\
            \vdots \\
            U_{k+m-1}
           \end{pmatrix}
  \end{equation}
as $\varepsilon \to 0$. One could also be content with the weaker statement that if the eigenvalues are multiple, the convergence statement of theorem \ref{thm:homogenisation} is only true up to taking a subsequence.
\end{remark}


\begin{remark}
  Literature on homogenisation theory is often concerned with the situation where holes are proportional to their reference cell. That is, $r_\varepsilon=c\varepsilon$ for some constant $c\in (0,1/2)$. 
  In this case one has $r_\varepsilon^{d-1} \varepsilon^{-d}\to\infty$. It follows
  from~\cite{ceg2} that $\sigma_k^\varepsilon\to 0$.
  Indeed it is proved there that any bounded domain $\Omega \subset \mathbb R^d$ satisfies 
\begin{equation}
  \sigma_k(\Omega)\abs{\partial\Omega}^{\frac 1{d-1}} \le C_{d,k},
\end{equation}
where the number $C_{d,k}>0$ depend only on the dimension $d$ and index $k$.
The hypothesis that $r_\varepsilon^{d-1} \varepsilon^{-d}\to\infty$ implies that
$\abs{\partial\Omega^\varepsilon} \to \infty$, which forces $\sigma_k^\varepsilon \to 0$, as claimed. Note that this also corresponds to the homogenisation regime which was studied by Vanninathan in~\cite{vanninathan} for a slightly different problem, for which the Dirichlet boundary condition was imposed on $\partial\Omega$ and the Steklov condition on $\partial T^\varepsilon$.

The regime that we consider in Theorem \ref{thm:homogenisation} is the critical
regime for the Steklov problem, where we observe a change of behaviour in the limiting problem. This is akin to the situation studied by Rauch--Taylor~\cite{RauchTaylor} and Cioranescu--Murat \cite{CioranescuMurat}.
\end{remark}





\subsection{Convergence to the Neumann problem and spectral comparison theorems}
The $\beta$ parameter in Problem \eqref{problem:homo} can be interpreted as a relative weight between the interior of $\Omega$ and its boundary $\partial \Omega$ for the behaviour of that problem, see Section \ref{sec:dynamical} for details on this interpretation. Our second main result is the following theorem, describing
the specific dependence on $\beta$ in \eqref{problem:homo}.

\begin{theorem}\label{thm:limitbeta}
  For each $k\in\mathbb N$, the eigenvalue $\Sigma_{k,\beta}$ depends continuously on $\beta\in[0,\infty)$ and satisfies
	\begin{gather*}
	\lim_{\beta\to\infty}{A_d}\beta \Sigma_{k,\beta} = \mu_k.
	\end{gather*}
	The eigenfunctions $\set{U_{k,\beta}}$ satisfy
\begin{equation}
 \beta^{1/2} U_{k,\beta} \to f_k
\end{equation}
weakly in $\RH^1(\Omega)$ as $\beta \to \infty$, where $f_k$ is the $k$th non--trivial Neumann eigenfunction.
\end{theorem}
\begin{remark}
 We make the second observation that this convergence \emph{cannot} be uniform in $k$, as that would contradict \cite[Theorem 4.4]{BelowFrancois}. 
\end{remark}
The relationships between isoperimetric type problems for the Neumann and Steklov eigenvalue problems have been investigated for the first few eigenvalues in~\cite{FreiLauW2018,FreiLauS2018} from the point of view of the Robin problem. Our methods also allow us to investigate the relationship between these isoperimetric problems for every eigenvalue rank $k$. 




The combination of Theorem~\ref{thm:limitbeta} and Theorem~\ref{thm:homogenisation} allows the transfer of known bounds for Steklov eigenvalues to bounds for Neumann eigenvalues. For instance, we can combine these two theorems with~\cite[Theorem 1.3]{ceg2}, asserting that for bounded Euclidean domains with smooth boundary
\begin{equation} \label{eq:boundceg2}
\sigma_k(\Omega)|\partial\Omega|^{\frac{2}{d-1}} \le C(d) k^{2/d}.
\end{equation}
This leads to the following.
\begin{cor} \label{cor:boundd}
  The Neumann eigenvalues of a bounded domain $\Omega\subset\mathbb R^d$ satisfy
  \begin{gather} \label{eq:dimbound}
    \mu_k(\Omega)|\Omega|^{2/d}\leq C(d)k^{2/d},
  \end{gather}
  where the constant $C(d)$ is exactly that of\, \cite[Theorem 1.3]{ceg2}.
\end{cor}
\begin{remark}
  The existence of a constant depending only on the dimension in inequality \eqref{eq:dimbound} is already known. In fact, Kröger obtained a better constant in \cite{Kroger1992}. However, it follows from Corollary \ref{cor:boundd} that any improvement to the bound  
  \eqref{eq:boundceg2} will transfer to bounds on Neumann eigenvalues.
  \end{remark}
  One of the original motivation for this project was the study of the following quantity:
$$\hat \sigma_k^*:=\sup\left\{\sigma_k(\Omega)|\partial\Omega|\,:\,\Omega\subset\mathbb R^2 \text{ bounded with smooth boundary}\right\}.$$
  In dimension $d=2$, we are able to get a stronger version of Corollary
  \ref{cor:boundd} in the sense that we obtain a direct link between $\hat \sigma_k^*$ and
  $$\hat \mu_k^*:=\sup\left\{\mu_k(\Omega)|\Omega|\,:\,\Omega\subset\mathbb R^2 \text{ bounded with smooth boundary}\right\}.$$
  In that case, we obtain.
  \begin{theorem} \label{thm:neumannd=2}
   For $d = 2$ and every $k \in \mathbb N$,
   \begin{equation}
   \hat \mu_k^* \le \hat \sigma_k^*.
   \end{equation}
  \end{theorem}
  \begin{remark}
  From~\cite{fraschoen3} we have that
  $$\sigma_1(\Omega)|\partial\Omega|<4\pi.$$
  It follows from Theorem~\ref{thm:homogenisation} and Theorem~\ref{thm:limitbeta} that
  $$\mu_1(\Omega)|\Omega|\leq 4\pi.$$
  Of course this bound is already known. Indeed the optimal upper bound is given by the famous Szeg\H{o}--Weinberger theorem:
  $$\mu_1(\mathbb{D})\pi\approx 3.39\pi.$$
  One can remark, however, that $4\pi$ corresponds to the P\'olya bound for the first Neumann eigenvalue. We can now also see that obtaining bounds of a similar type for $\sigma_k$ would also transfer to $\mu_k$.
\end{remark}

The previous discussion also yields the following corollary.
\begin{cor}
  $$\hat \sigma_1^*\geq \mu_1(\mathbb{D})\times\pi\approx 3.39\pi.$$
\end{cor}

Indeed, by the Szego-Weinberger inequality, we have that
\begin{equation}
 \hat \sigma^*_1 \ge \pi \mu_1(\mathbb{D}) \approx 3.39 \pi.
\end{equation}
Furthermore, this number will be approached as close as desired in homogenisation sequence of pierced unit disks with large enough parameter $\beta$.  Of course, it is known from \cite{fraschoen2} that if one is allowed to optimise amongst surfaces rather than Euclidean domains, then there is a sequence of surfaces such that $\sigma_1(\Omega_n)\abs{\partial \Omega_n} \to 4\pi$.
This leads to the natural conjecture
\begin{conjecture}
  $$\hat \sigma_1^* = \mu_1(\mathbb{D})\times\pi\approx 3.39\pi.$$
\end{conjecture}
Note that the previous best known lower bound for $\hat \sigma_1^*$ was 
attained on some concentric annulus, whose first normalised Steklov eigenvalue is approximately $2.17 \pi$, see \cite{gpsurvey}. We also observe that a similar analysis yields that for any $\Omega \subset \mathbb R^2$ bounded with smooth boundary,
$$ \hat \sigma_k^* \ge \mu_k(\Omega) \abs{\Omega}.$$ In particular, it follows from Weyl's law for Neumann eigenvalues that there exists a sequence $a_k \sim 4\pi k$ such that
$$\hat \sigma_k^* \ge a_k.$$


\subsection{Discussion}
Homogenisation theory is a young branch of mathematics which star\-ted around
the 1960's. Its general goal is to describe macroscopic properties of materials
through their microscopic structure. To the best of our knowledge, the first
papers to study periodically perforated domains from a rigourous mathematical
point of view are those of Marchenko and Khruslov from the early 1960's
(e.g.~\cite{MarchenkoKhruslov1964}) leading to their influential
book~\cite{MarchenkoKhruslov} in 1974. The topic became widely known in the West
with the work of  Rauch and Taylor on the \emph{crushed ice
problem}~\cite{RauchTaylor} in 1975 and then with the publication in 1982
of~\cite{CioranescuMurat} by Cioranescu and Murat. Many of these early results
were concerned with the Poisson problem $\Delta u_\varepsilon=f$ under Dirichlet
boundary conditions $u_\epsilon=0$ on $\partial\Omega_\epsilon$. The limiting
behaviour of the solution $u_\varepsilon$ depends on the rate at which $r_\varepsilon\searrow
0$. Three regimes are considered. If the size of the holes $r_\epsilon$ tends to
zero very fast, then in the limit the solutions tend to those of the Poisson
problem on the original domain $\Omega$, while if the size of the holes are big
enough the solutions tend to zero. The main interest comes from the critical
regime, in which case the solutions tend to solutions of a new elliptic
problem.

In this paper we are concerned with the much less studied Steklov spectral
problem \eqref{Problem:Steklov}. From the point of view of homogenisation
theory, this problem is atypical in the sense that the function spaces which
occur for different values of the parameter $\varepsilon>0$ are not naturally related.
This means in particular that this problem does not yield to any of the usual
general frameworks used in homogenisation theory (such as \cite[Chapter
11]{JikovKozlovOlenik}). Nevertheless, several authors have considered
homogenisation for this problem, using ad hoc methods depending on the specifc
situation considered. The behaviour of Steklov eigenvalues under singular
perturbations such as the perforation of a single hole has also been studied
in~\cite{GryCristo2014, Nazarov2014}.


Our main inspiration for this work is the paper \cite{vanninathan}. Several
papers have also considered homogenisation in the situation where the Dirichlet
condition is imposed on the outer boundary while the Steklov condition is
considered only on the boundary of the holes
\cite{ChiadoPiatNazarovPiatnitski,Douanla}. In these papers the holes are
proportional to the size of the reference cell. The novelty of our
homogenisation result in the case of the Steklov problem is that consider holes
that are shrinking much faster than that, in a critical regime where the
limiting problem
is fundamentally different. They are in fact shrinking at the precise rate which
makes their total surface area (or perimeter in dimension $2$) is asymptotically
comparable with the volume of the domain. This is similar to the work of
Rauch--Taylor \cite{RauchTaylor} for the Neumann problem, Cioranescu--Murat
\cite{CioranescuMurat} for the Dirichlet problem and Kaizu for the Robin problem
\cite{Kaizu}.

The \emph{energy method} of Tartar (see
\cite[Section 1.3]{allaire} for an exposition) has been used extensively in the
study of homogenisation problems at critical regimes, see
\cite{CioranescuMurat,Kaizu} and more recently \cite{CheredDondlRosler}, where
they obtain norm-resolvant convergence for the Dirichlet, Neumann and Robin
problem. That method uses an auxiliary function, satisfying some
energy-minimising PDE in the fundamental cells, in order to derive convergence
of the problem in the weak formulation. The method, in its Robin or Neumann
form, is boundary-condition agnostic and as such is ill-suited for the Steklov
problem, where the normalisation is with respect to $\RL^2(\partial \Omega^\varepsilon)$.
Indeed, while the technique could be used to obtain some form of convergence, it
will not be able to transfer boundary estimates to $\RL^2(\Omega)$, and ensure
that the limit solution doesn't degenerate to the trivial one.

Nevertheless, one can interpret our technique as a variation on the energy
method, adapted for problems defined on the boundary. We are also using an auxiliary PDE in order
to derive convergence, but it does not stem from compensated compactness. The
main difference, however, is that we can deduce interior estimates from those on
the boundary of the periodic holes from our auxiliary problem, see Lemma
\ref{lem:normconv}.

 We note that in this paper we have chosen to consider only spherical holes.
In fact, it should be possible to consider more general convex holes
obtained as scaled copies $r_\varepsilon\omega$ of a fixed convex set $\omega$, as
it is done in classical homogenization literature. Here, convexity of the holes
would be required for $\RL^\infty$ estimates of the Steklov eigenfunctions, see
Lemma \ref{lem:tracebv}.

Since we are motivated by spectral questions, namely to explore a new link
between the Steklov problem and the Neumann problem, we decided to avoid the 
technical difficulties which would occur by considering
more general inclusions (which would probably provide a more complicated limit problem
with some kind of capacity of the holes and a polarization matrix).
This allows us to get the simpler dynamical eigenvalue problem (3) and emphasise
clearly the link with the the classical
Neumann eigenvalue problem by letting the parameter $\beta \to +\infty$.


\subsection{Structure of the proof and plan of the paper}

In Section~\ref{section:NotationsFunctionSpaces} we formally describe properties
of the various eigenvalue problems that we study as well as the functions spaces over which they are defined. While they are well-known, the notation used for all of them often collides. In that section, we fix notation once and for all for the remainder of the paper for definedness and ease of references. 

In Section~\ref{section:ComparisonTheorems} we study one of the main technical tools in this paper, properties of harmonic extensions of functions on annuli to the interior disk. Results are separated in two categories: those that rely on the fact that the functions satisfy a Robin-type boundary condition, and more 
general results that do not rely on such a thing. Most of our results will be obtained by considering the Fourier expansion of a function
$$
u(\rho,\boldsymbol{\theta}) = \sum_{\ell,m} a_{\ell}^m(\rho) Y_\ell^m(\boldsymbol{\theta})
$$
in spherical harmonics, and obtaining our inequalities term by term for every $a_\ell^m$. 

Section \ref{sec:homogenisation} is the \emph{pièce de résistance} of this paper.
It is where we show Theorem \ref{thm:homogenisation} and the proof proceeds in many steps. We first prove that the family of harmonic extensions $U^\varepsilon:=U_k^\varepsilon$ \/ is bounded in the Sobolev space $H^1(\Omega)$ hence there
   exists a subsequence $\varepsilon_n\searrow 0$ such that $U^{\varepsilon_n}$ weakly converges to a function $U\in H^1(\Omega)$. This allows us to consider properties of the weak limit $U$, and of an associated limit $\Sigma$ to the eigenvalue sequence $\sigma_k^{\varepsilon_n}$.
   
   It is then not so hard to show that, using the weak formulations of Problems \eqref{Problem:Steklov} and \eqref{problem:homo} that the limit of the homogenised Steklov problem contains terms corresponding to the limit dynamical eigenvalue problem, plus some spurious terms that must be shown to
   converge to zero. This is done by studying two representations of the weak formulation using Green's identity either towards the inside of the holes or the inside of the domains $\Omega^\varepsilon$. The first one is used to show that
   the functionals that arise in the study of the homogenisation problem are
   uniformly bounded, hence we can use smooth test functions. In the second representation, we are therefore allowed to use smooth test functions, which allows us to recover convergence to zero of the spurious terms. A key step in this argument is to  understand the limit behaviour of an auxiliary homogenisation problem for the transmission eigenvalue problem (see Proposition~\ref{prop:limitproblem}).
   
   Once we have established convergence to a solution of Problem \ref{problem:homo}, we end up showing that the limit eigenpair $(\Sigma,U)$ does not degenerate to the trivial function. Using variational characterisation of eigenvalues and eigenfunctions, we can also show that we get complete spectral convergence, and that subsequences are not needed.
   
  Finally, in Section \ref{Section:LargeBeta} we show convergence to the Neumann
  problem as $\beta \to \infty$. The method is similar to the one used at the end
  of the previous section, but many of the inequalities are more subtle. We also
  show the comparison theorems between Steklov and Neumann eigenvalues in this section.
  
  In this paper, we use $c$ and $C$ to mean constants whose precise value is not important to our argument, and whose exact value may change from line to line. We use the notations $f = \bigo g$ and $f \ll g$ interchangeably to mean that there exists a constant $C$ such that $\abs{f(x)} \le Cg(x)$.
  
\section{\bf Notation and function spaces}
\label{section:NotationsFunctionSpaces}
Four different eigenvalue problems will be used. The goal of this section is to introduce them and fix the relevant notation.
Throughout the paper, we use real valued functions.

\subsection{The Steklov problem on $\Omega$ and $\Omega^\varepsilon$}

Given a bounded domain $\Omega$ whose boundary $\partial \Omega$ is smooth, the Dirichlet-to-Neumann 
operator (DtN map) $\Lambda$ acts on $C^\infty(\partial\Omega)$ as
\begin{equation}
  \Lambda f = \partial_\nu \hat f,
\end{equation}
where $\hat f$ is the harmonic extension of $f$ to the interior of $\Omega$. The DtN map is an
elliptic, positive, self--adjoint pseudodifferential operator of order $1$. Because $\partial\Omega$ is compact, it follows from
standard theory of such operators, see \emph{e.g} \cite{shubin}, that the eigenvalues form a non--negative unbounded sequence
$\set{\sigma_k : k \in \mathbb N_0}$ and that there exists an orthonormal basis of $(f_k)$ of $\RL^2(\partial\Omega)$ such that
$\Lambda f_k=\sigma_kf_k$.
The  harmonic extensions $u_k=\hat f_k$ satisfy the Steklov problem
\begin{equation}
  \begin{cases}
   - \Delta u_k = 0 & \text{in } \Omega,\\
    \partial_\nu u_k = \sigma_k u_k &\text{on } \partial\Omega.
  \end{cases}
\end{equation}
In general, we use the same symbol $u_k$ for the function on $\Omega$ and for its trace on the boundary $\partial\Omega$.
The eigenvalue sequence for $\Omega^\varepsilon$ is
denoted $\set{\sigma_k^\varepsilon : k \in \mathbb N_0)}$, with corresponding eigenfunctions $u_k^\varepsilon$.
The
eigenfunctions $u_k$, respectively $u_k^\varepsilon$, form an orthonormal basis with
respect to the inner products
\begin{equation}
  (f,g)_\partial := \int_{\partial \Omega} f g \, \mathrm{d} A,
\end{equation}
respectively
\begin{equation}
  (f,g)_{\partial^\varepsilon} := \int_{\partial \Omega^\varepsilon} f g \, \mathrm{d} A.
\end{equation}
The $k$-th nonzero eigenvalue $\sigma_k$ is characterised by
\begin{equation}
  \label{eq:varcharsigm}
  \sigma_k = \inf\set{\frac{\int_\Omega \abs{\nabla u}^2 \, \mathrm{d} \bx}{\int_{\partial \Omega}u^2
  \, \mathrm{d} A}:u \in \RH^1(\Omega) \text{ and }  (u,u_j)_{\partial}=0 \mbox{ for } 0 \le j < k}.
\end{equation}
The eigenvalues
$\sigma_k^\varepsilon$ have the same characterisation, integrating over $\Omega^\varepsilon$
and $\partial \Omega^\varepsilon$ respectively instead, and with the orthogonality being
with respect to $(\cdot,\cdot)_{\partial^\varepsilon}$.

\subsection{Dynamical eigenvalue problem} \label{sec:dynamical}
For $\beta \in (0,\infty)$, consider
the eigenvalue problem:
\begin{equation}\label{eq:dyneigenprob}
\begin{cases}-\Delta U = A_d \beta \Sigma U & \text{in } \Omega,\\
    \partial_\nu U = \Sigma U & \text{on } \partial \Omega.
\end{cases}
\end{equation}
where $A_d$ is the area of the unit sphere in $\mathbb R^d$.
Problem \eqref{eq:dyneigenprob} was introduced with a slightly different normalization in \cite{BelowFrancois}, where it is called a \emph{dynamical eigenvalue problem}.
The eigenvalues and
eigenfunctions are those of the operator
\begin{equation}
  P := \begin{pmatrix} -(A_d\beta)^{-1} \Delta & 0 \\ 0 & \partial_\nu \end{pmatrix}.
\end{equation}
This unbounded operator is 
defined on an appropriate domain in the space $\RL^2_{A_d \beta}(\Omega) \times \RL^2 (\partial \Omega)$ which consists simply of
$\RL^2(\Omega) \times \RL^2 (\partial \Omega)$ equipped
with the inner product defined by
\begin{equation}
  (f,g)_\beta := A_d \beta\int_\Omega f g \, \mathrm{d} \bx + \int_{\partial \Omega} f
  g \, \mathrm{d} A,
\end{equation}
The dynamical eigenvalue problem~\eqref{eq:dyneigenprob}
has a discrete sequence of eigenvalues
\begin{equation}
  0 = \Sigma_{0,\beta} < \Sigma_{1,\beta} \le \Sigma_{2,\beta} \le \dotso \nearrow
  \infty.
\end{equation}
Let $X\subset \RL^2_{A_d \beta}(\Omega) \times \RL^2 (\partial \Omega)$ be the subspace
defined by
\begin{equation}
  X:= \set{U = (u,\tau u) : u \in \RH^1(\Omega)},
\end{equation}
where $\tau:\RH^1(\Omega) \to \RL^{2}(\partial\Omega)$ is the trace operator.
The eigenfunctions $U_{k,\beta}$ associated to $\Sigma_{k,\beta}$ form a basis of $X$,
but not of $\RL^2_{A_d\beta}(\Omega) \times \RL^2(\partial\Omega)$. The eigenvalues
$\Sigma_{k,\beta}$ are characterised by

\begin{equation}
  \Sigma_{k,\beta} = \inf\set{\frac{\int_\Omega \abs{\nabla U}^2 \, \mathrm{d}
    \bx}{A_d \beta\int_\Omega U^2 \, \mathrm{d} \bx + \int_{\partial \Omega} U^2
      \, \mathrm{d} A}:U \in \RH^1(\Omega) \text{ and }  (U,U_{j,\beta})_{\partial}=0 \mbox{ for }0 \le j < k
    }.
\end{equation}

\subsection{The Neumann eigenvalue problem}
We will also make use of the classical Neumann eigenvalue problem
\begin{equation}
  \begin{cases}-\Delta f = \mu f & \text{in } \Omega,\\
    \partial_\nu f = 0 &
    \text{on } \partial \Omega.
  \end{cases}
\end{equation}
The Neumann eigenvalues form an increasing sequence
\begin{equation}
  0 = \mu_0 < \mu_1 \le \mu_2 \le \dotso \infty,
\end{equation}
with associated eigenfunctions $f_k$, orthonormal with respect to the
$\RL^2(\Omega)$ inner product
\begin{equation}
  (f,g) := \int_\Omega f g \, \mathrm{d}\bx.
\end{equation}
The eigenvalues are characterised by
\begin{equation}
  \mu_k = \inf\set{\frac{\int_\Omega \abs{\nabla f}^2\, \mathrm{d} \bx }{\int_{\Omega} f^2
  \, \mathrm{d} \bx}: f \in \RH^1(\Omega) \text{ and } (f,f_j)=0 \mbox{ for }j=0,1,\cdots,k-1}. 
\end{equation}


\section{\bf Comparison theorems}\label{section:ComparisonTheorems}
In this section, we derive comparison inequalities that will be used repeatedly.
For $0 < r < R$, set $$A_{r,R} = B(0,R) \setminus\overline{B(0,r)}.$$
For any set $X$ where it is defined, the Dirichlet energy of a function $f:X\rightarrow\mathbb R$ is
$$\mathcal{D}(f) :=\int_{X}|\nabla f|^2\,\, \mathrm{d}\bx.$$
 The Dirichlet energy of $f$ on a subset $Y\subset X$ is written
$\mathcal{D}(f;Y)$. Note that since every problem under consideration is self-adjoint no
generality is lost by studying real-valued functions.

Lemmas \ref{lemma:localEstimateFourier} and \ref{lemma:Sobolevconst} are proved by
  using a Fourier series decomposition for functions in an annulus.
  We note that even though the proofs rely on estimates on the radial part of
  these functions, we do not claim that the inequalities proved therein are
  realised by radial optimisers. The standard Schwarz symmetrisation arguments do not
apply here because the domain of the functions are annuli rather than balls.

\subsection{Comparison theorems for functions satisfying a Steklov boundary
condition}

\begin{lemma}\label{lemma:localEstimateFourier}
  Fix a postive real number $\sigma>0$. For any $0 < r < R \le 1$,
  let $u\in C^\infty(\overline{A_{r,R}})$ be such that
  \begin{equation}
  \begin{cases}
  \Delta u=0&\mbox{ in }A_{r,R},\\
  \partial_\nu u=\sigma u&\mbox{ on } \partial B(0,r).
  \end{cases}
  \end{equation}
  Consider the function $h:B(0,r)\rightarrow\mathbb R$ defined by
  $$\begin{cases}
  h=u&\mbox{ on }\partial B(0,r),\\
  \Delta h=0&\mbox{ in }B(0,r).
  \end{cases}
  $$
  Then as the ratio $r/R$ goes to $0$,
  \begin{equation}
    \label{eq:dircomp}
    \mathcal{D}(h) \le 5 \mathcal{D}(u) \left(\frac r R \right)^d \left( 1 + \bigo{\left(\frac r
    R \right)^d} \right).
  \end{equation}
  \end{lemma}

\begin{proof}
  For every $\ell \ge 0$, we denote by $N_\ell$ the dimension of the space
  $H_\ell$ of spherical harmonics of order $\ell$ and denote by $Y_\ell^m(\boldsymbol{\theta})$, 
  $1 \le m \le N_\ell$ the standard orthonormal basis of spherical harmonics on the unit sphere.
  On $A_{r,R}$, the function $u$ admits a Fourier decomposition in spherical harmonics
\begin{equation}
  u(\rho,\boldsymbol{\theta}) = \sum_{\substack{\ell \ge 0 \\1 \le m \le N_\ell}}
  a_\ell^m(\rho)
  Y_\ell^m(\boldsymbol{\theta}).
  \label{eq:fourierspherical}
\end{equation}
We start by studying the form of the coefficients $a_\ell^m$.
\newline
\noindent
\textbf{Case $d>2$}.
The harmonicity condition on $u$ implies that the radial parts
$a_\ell^m(\rho)$ are given by
\begin{equation}
  a_\ell^m(\rho) = c_\ell^m \rho^{\ell} + c_{-\ell}^m \rho^{-\ell + 2 - d}.
  \label{eq:radial}
\end{equation}
By convention the coefficients $c_0^1$ and $c_{-0}^1$ are assumed to be different, the minus sign referring as for the other coefficients to the solution blowing up at the origin. The Steklov condition for $u$ on $\partial B(0,r)$ along with the orthogonality of
the spherical harmonics $Y_\ell^m$ imply
\begin{equation}
  - \partial_\rho a_\ell^m(r) = \sigma a_\ell^m(r),
  \label{eq:steklovint}
\end{equation}
which yields the relations
\begin{equation}
  c_{-\ell}^m = \overbrace{\left(\frac{\ell + r \sigma}{\ell - 2 + d - r
  \sigma}\right)}^{M:=}r^{2\ell+d-2}c_\ell^m.
  \label{eq:relcoeff}
\end{equation}
In turns this yields the following explicit expression for the radial functions
\begin{equation}
  \label{eq:radexplicit}
  a_\ell^m(\rho) = c_\ell^m \rho^{\ell}\left(1+ M
  \left(\frac{r}{\rho}\right)^{2\ell+d-2}\right).
\end{equation}
For $r \le \frac{d - 2}{2\sigma}$, it follows that
\begin{equation}
\label{eq:Mdg2}
  \frac{\ell}{\ell - 2 + d}\le  M \le 1.
\end{equation}


  \noindent
\textbf{Case $d=2$}.
Equations~\eqref{eq:radial} holds except at $\ell = 0$,
in which case, 
$$a_0^1(\rho) = c_0^1 + c_{-0}^1 \log \rho.$$
In that case, \eqref{eq:radexplicit} becomes
\begin{align*}
 a_0^1 &= c_0^1 \left(1 + \frac{r\sigma}{1 + r \log r} \log(1/\rho) \right) \\
 &= c_0^1 \left(1 + M \log(1/\rho)\right).
 \end{align*}
Observe for the sequel that when $d = 2$ and $\ell = 0$, then $M = \bigo{r\sigma}$, and if $\ell > 0$, then
\begin{equation} \label{eq:Md2}
  M = 1 + \bigo{r \sigma}.
\end{equation}



 Because inequality~\eqref{eq:dircomp} is invariant
under scaling, it is sufficient to prove the case $R= 1$ and to let $r \to 0$.
The harmonic extension of $u$ to $B(0,r)$ is given by
\begin{equation}
  h(\rho,\boldsymbol{\theta}) = \sum_{\substack{\ell \ge 0 \\ 1 \le m \le N_\ell}}
  a_\ell^m(r) \frac{\rho^\ell}{r^\ell} Y_\ell^m(\boldsymbol{\theta}).
  \label{eq:extension}
\end{equation}
The Dirichlet energy of $h$ is
\begin{equation}
  \mathcal{D}(h) = \sum_{\substack{\ell \ge 1 \\ 1 \le m \le N_\ell}} \ell a_\ell^m(r)^2
  r^{d-2}.
  \label{eq:dirichh}
\end{equation}
On the other hand, the Dirichlet energy of $u$ is given, from Green's identity,
and the Steklov condition on $\partial B(0,r)$
\begin{equation}
  \mathcal{D}(u) = \sum_{\ell,m} \sigma a_\ell^m(r)^2r^{d-1} + \sum_{\ell,m}
  a_\ell^m(1)\partial_\rho a_\ell^m(1).
  \label{eq:dirichu}
\end{equation}
Our goal is now to find a bound on \eqref{eq:dirichh} in terms of
\eqref{eq:dirichu}.
It is sufficient to show that for each $\ell \ge 1$,
\begin{equation}
  \ell a_\ell^m(r)^2 \le 4 r^2 a_\ell^m(1) \partial_\rho a_\ell^m(1)\left( 1 +
  \bigo{r^d}    \right).
  \label{eq:sufficientbound}
\end{equation}
Suppose without loss of generality that $c_\ell^m = 1$. The substitution of \eqref{eq:relcoeff} into \eqref{eq:radial} imply that
\begin{equation}
   a_\ell^m(r)^2 =  r^{2 \ell} \left( 1 + M \right)^2 
  \le 5r^2,
  \label{eq:coeffh}
\end{equation}
and that
\begin{equation}
  a_\ell^m(1)\partial_\rho a_\ell^m(1) = \ell \left( 1 + \frac{(2-d)
    M}{\ell}
r^{2\ell + d - 2} + \frac{( 2 - d - \ell) M^2}{\ell} r^{2(2 \ell + d-2)}\right)
= \ell + \bigo{r^d}.
  \label{eq:coeffu}
\end{equation}
Hence, dividing \eqref{eq:coeffh} by \eqref{eq:coeffu}, and using the bounds
on $M$ from \eqref{eq:Mdg2} and \eqref{eq:Md2} we have that for $\ell\ge 1$,
\begin{equation}
\ell a_\ell^m(r)^2 \le 5 a_\ell^m(1) \partial_\rho
a_\ell^m(1)r^2 \left( 1 + \bigo{r^d} \right).
  \label{eq:boundcoeffs}
\end{equation}
This proves  our claim.

\end{proof}

\subsection{General $\RH^1$ comparison theorems on annuli and balls}

The next two lemmas do not depend on any  specific boundary condition.
The first one gives bounds for Sobolev constants of annuli.
\begin{lemma}
  \label{lemma:Sobolevconst}
  For \/ $0 < r < R < 1$, define
  \begin{equation}
    \label{eq:rayleighgamma}
    \gamma(r,R) := \inf\set{
      \frac{\int_{A_{r,R}} \abs{\nabla u}^2 + u^2\,\, \mathrm{d} \bx}{\int_{\partial
    B(0,r)} u^2\,dA}:u \in \RH^1(A_{r,R}), u\big|_{\partial B(0,r)} \not \equiv 0}.
  \end{equation}
 Suppose that $R \ge c r^{\frac{d-1}{d}} \ge 2r$ for some $c > 0$. Then, there
  is a constant $C$ depending only on the dimension and on $c$ such that
  \begin{equation}
    \gamma(r,R) \ge C \min\set{R^d r^{1-d},r^{\frac 1 d - 1} }.
  \end{equation}
\end{lemma}

\begin{proof}
 
  As earlier, write a function $u \in \RH^1(A_{r,R})$ as
\begin{equation}
  u(\rho,\boldsymbol{\theta}) = \sum_{\ell,m} a_\ell^m(\rho) Y_\ell^m(\boldsymbol{\theta}).
\end{equation}
Using the notation $u_{\boldsymbol{\theta}}$ for the tangential gradient, the Dirichlet energy of $u$ is expressed as
\begin{equation}
  \begin{aligned}
    \mathcal{D}(u) &= \int_r^R \int_{\mathbb S^{d-1}} \left( u_\rho^2 + \rho^{-2} u_{\boldsymbol{\theta}}^2 \right)
    \rho^{d-1} \, \mathrm{d} \boldsymbol{\theta} \, \mathrm{d} \rho \\
    &\geq
    \int_r^R \int_{\mathbb S^{d-1}} \left( u_\rho^2\right)
    \rho^{d-1} \, \mathrm{d} \boldsymbol{\theta} \, \mathrm{d} \rho \\
    &= \sum_{\ell,m} \int_r^R \left(a_\ell^m(\rho)'\right)^2 \rho^{d-1} \, \mathrm{d}
    \rho.
  \end{aligned}
\end{equation}
On the other hand, the denominator in \eqref{eq:rayleighgamma} is given by
    \begin{equation}
      \int_{\partial B(0,r)} u^2 \, \mathrm{d} A = r^{d-1} \sum_{\ell,m} a_\ell^m(r)^2.
    \end{equation}
    Combining these last two expressions in \eqref{eq:rayleighgamma} and defining the
    density $w(\rho) = \left(\frac \rho r \right)^{d-1}$, we see it is enough to prove that
\begin{equation} 
  \label{eq:ineqsum}
  \sum_{\ell,m} \int_r^R \left(\left( \partial_\rho a_\ell^m(\rho) \right)^2 +
  a_\ell^m(\rho)^2\right)
  w(\rho) \, \mathrm{d} \rho \ge  C \min\set{R^d r^{1-d},r^{\frac 1 d - 1} } \sum_{\ell,m}
  a_\ell^m(r)^2.
\end{equation}
Indeed, working term by term, we prove that any smooth function
$f : [r,R] \to \mathbb R$ satisfies
\begin{equation}
  \int_r^R \left(f'(\rho)^2 + f(\rho)^2\right)w(\rho) \, \mathrm{d} \rho \ge
  C \min\set{R^d r^{1-d},r^{\frac 1 d - 1} } f(r)^2.
\end{equation}
To this end, assume without loss of generality that $f(r) = 1$.
Following a strategy that was used in \cite{CEG3} and in \cite{CianciGirouard}, consider the two following situations.

Let $t\in (r,R)$, to be fixed later.

\noindent
\textbf{Case a.} Suppose first that for all $\rho \in (t,R)$,
$\abs{f(\rho)} \ge 1/2$. It follows from monotonicity and explicit integration that
\begin{equation}
  \int_r^R \abs{f(\rho)}^2 w(\rho) \, \mathrm{d} \rho \ge 
  \frac{r^{1-d}}{4d}\left(R^d - t^d \right).
\end{equation}

\noindent
\textbf{Case b.} Suppose there exists $\rho_0 \in (t,R)$ such that
$f(\rho_0) < 1/2$. Splitting the integral, using the fact that $w(\rho) \ge 1$ and is increasing
for all $\rho$ together with the Cauchy--Schwarz inequality leads to
\begin{equation}
  \label{eq:sumsquares}
  \begin{aligned}
    \int_r^R f'(\rho)^2 w(\rho) \, \mathrm{d} \rho &\ge \int_r^t f'(\rho)^2 +
    \int_t^{\rho_0} f'(\rho)^2w(\rho) \, \mathrm{d} \rho \\
    &\ge \frac{1}{t - r}\left( \int_r^t f'(\rho) \, \mathrm{d} \rho \right)^2 +
    \left(\frac t r \right)^{d-1}\frac{1}{\rho_0 - t}\left(\int_t^{\rho_0} f'(\rho) \, \mathrm{d}
    \rho\right)^2.
  \end{aligned}
\end{equation}
  By hypothesis $R<1$ so that $\frac{1}{\rho_0-t}>1$. This leads to
\begin{equation}
  \begin{aligned}
    \int_r^R f'(\rho)^2 w(\rho) \, \mathrm{d} \rho 
    \ge \frac 1 2 \min\set{\frac{1}{t - r},\left( \frac t r
    \right)^{d-1}}\left(\int_r^{\rho_0}f'(\rho) \, \mathrm{d} \rho\right)^2 ,
  \end{aligned}
\end{equation}
Choosing $t = \min\set{r^{\frac{d-1}{d}},R/2}$ guarantees that
  $\min\set{\frac{1}{t - r},\left( \frac t r
    \right)^{d-1}}=\left( \frac t r
  \right)^{d-1}$ so that
\begin{equation}
  \begin{aligned}
    \int_r^R \abs{f'(\rho)}^2 w(\rho) \, \mathrm{d} \rho 
    \ge \frac{1}{2}\left( \frac t r
    \right)^{d-1}\left(\int_r^{\rho_0}f'(\rho) \, \mathrm{d} \rho\right)^2.
  \end{aligned}
\end{equation}
It follows from the definition of $t$ that
$$\left(\frac{t}{r}\right)^{d-1}\geq\overbrace{\min\left\{1,\left(\frac{c}{2}\right)^{d-1}\right\}}^{A:=}r^{\frac{1-d}{d}}.$$
We can bound asymptotically the last integral in \eqref{eq:sumsquares} as
\begin{equation}
  \begin{aligned}
    \int_r^R f'(\rho)^2 w(\rho) \, \mathrm{d} \rho
    & \geq Cr^{\frac{1-d}{d}}
  \left(\int_r^{\rho_0}f'(\rho)\, \mathrm{d} \rho\right)^2\\
  &\geq \frac{A}{4}r^{\frac{1-d}{d}}.
\end{aligned}
\end{equation}
This also ensures that $R^d - t^d \ge (1-2^{-d})R^d$.
Since both situations are exclusive, inequality \eqref{eq:ineqsum} holds, finishing the proof.
\end{proof}
The next lemma compares $\RL^2$ norms on $B(0,r)$ with $\RH^1$ norms on
$B(0,R)$. Here, for any $\Omega' \subset \Omega$, the norm on $\RH^1(\Omega')$ is given by
\begin{equation}
  \norm{u}_{\RH^1(\Omega')}^2 = \norm{u}_{\RL^2(\Omega')}^2 + \norm{\nabla
  u}^2_{\RL^2(\Omega')^2}.
\end{equation}
\begin{lemma} \label{lemma:smalltolarge}
 For \/ $0 < r < R\leq 1$, if $R \ge c r^{\frac{d-1}{d}}$ for some $c>0$, there is a constant $C$ depending only on $c$ and on the dimension such that for all $u \in \RH^1(B(0,R))$, 
\begin{equation} \label{eq:boundrR}
 \norm{u}_{\RL^2(B(0,r))} \le C r^{1/2} \norm{u}_{\RH^1(B(0,R))}.
\end{equation}
 \end{lemma}
\begin{proof}
   Let $u\in H^1(B(0,R))$. Given $r\in (0,R)$,
   $$\int_{B(0,r)}u^2\,\, \mathrm{d} \bx=\int_{0}^r\rho^{d-1}\int_{S^{d-1}}u^2\,\, \mathrm{d}\theta\,\, \mathrm{d}\rho
   =\int_{0}^r\|u\|_{L^2(\partial B(0,\rho))}^2\,\, \mathrm{d}\rho
   $$
   It follows from the definition of $\gamma$ in Lemma~\ref{lemma:Sobolevconst} that
   $$
   \|u\|_{\partial B(0,\rho)}^2\leq\frac{1}{\gamma(\rho,R)}\|u\|_{H^1(A_{\rho,R})}^2
   \leq\frac{1}{\gamma(\rho,R)}\|u\|_{H^1(B(0,R))}^2.
   $$
   Substitution in the above leads to
   $$\int_{B(0,r)}u^2\,\, \mathrm{d} \bx\leq
   \|u\|_{H^1(B(0,R))}^2\int_{0}^r\frac{1}{\gamma(\rho,R)} \, \mathrm{d} \rho.
   $$
   It follows from Lemma~\ref{lemma:Sobolevconst} that
     \begin{align*}
       \int_{0}^r\frac{1}{\gamma(\rho,R)}\,\, \mathrm{d}\rho&\leq\frac{1}{C}\int_{0}^r\left(\frac{\rho^{d-1}}{R^d}+\rho^{1-\frac{1}{d}}\right)\, \mathrm{d}\rho\\
       &=\frac{1}{C}\left(\frac{r^d}{dR^d}+\frac{r^{2-\frac{1}{d}}}{2-\frac{1}{d}}\right)\\
       &\leq
       \frac{1}{C}\left(\frac{r}{cd}+\frac{r^{2-\frac{1}{d}}}{2-\frac{1}{d}}\right)\qquad\longleftarrow\text{ since }R^d\geq cr^{d-1}\\
       &\leq \tilde{C}r.
     \end{align*}
\end{proof} 
Finally, we require the following lemma about the behaviour of the boundary trace operator as a domain gets shrunk.

\begin{lemma}
 \label{lem:traceconstant}
 Let $\Omega \subset \mathbb R^d$ be a bounded open set and denote
 \begin{equation}
  \gamma(\Omega):= \inf\set{\frac{\int_\Omega \abs{\nabla u}^2 + u^2 \, \mathrm{d} \bx}{\int_{\partial \Omega} u^2 \, \mathrm{d} A} : u \in \RH^1(\Omega), u\big|_{\partial \Omega} \not \equiv 0}.
 \end{equation}
 Then, the following inequality holds as $\varepsilon \to 0$:
 \begin{equation}
  \gamma(\varepsilon\Omega) \ge \varepsilon \gamma(\Omega).
 \end{equation}
\end{lemma}
\begin{proof}
 Consider $u \in \RH^1(\varepsilon\Omega)$, $u\big|_{\partial \varepsilon\Omega} \not \equiv 0$. Then, with the change of variable $\by = \varepsilon^{-1} \bx$,
 \begin{align}
  \int_{\varepsilon \Omega} \abs{\nabla u(\bx)}^2 + u(\bx)^2 \, \mathrm{d} \bx &=
  \varepsilon^{d} \int_{\Omega} \varepsilon^{-2} \abs{\nabla u(\by)}^2 + u(\by)^2 \, \mathrm{d} \by \\
  &\ge \varepsilon^{d} \int_{\Omega} \gamma(\Omega) \int_{\partial \Omega} u(\by)^2 \, \mathrm{d} A 
  \\
&= \varepsilon \gamma(\Omega) \int_{\partial \varepsilon \Omega} u(\bx)^2 \, \mathrm{d} A.
 \end{align}

\end{proof}

\subsection{Uniform bounds on Steklov eigenfunctions}

In order to obtain convergence of the eigenfunctions we will need that their
$\RL^\infty$ norm stays bounded. In this subsection, since we will need to
understand specific interplay between boundary surface area and volume, we will diverge
from our convention and denote the area of a boundary $\partial \Omega$ by
$\mathcal{H}^{d-1}(\partial \Omega)$. In \cite{bgt}[Theorem 3.1], the authors prove that for any
Steklov eigenfunction $u$ of a domain $\Omega$ with eigenvalue $\sigma$,
\begin{equation}
  \norm{u}_{\RL^\infty(\Omega)} \le C \norm{u}_{\RL^2(\partial \Omega)}
\end{equation}
where $C$ depends continuously on the dimension $d$, $\sigma$, $\abs{\Omega}$
and the norm of the trace application
\begin{equation}
  T : \RB\RV(\Omega) \to \RL^1(\partial \Omega).
\end{equation}
It is clear that in our situation, the dimension and $\abs{\Omega}$ stay
bounded. The eigenvalues will be shown later to also stay bounded, but we use
to control the norm of $T$. In \cite{bgt}[Proposition 5.1], they give the
following
condition under which said norm stays bounded for a family of domains. If
$\set{\Omega^\varepsilon}$ is a family of open bounded domains with
$\mathcal{H}^{d-1}(\partial\Omega^\varepsilon) < \infty$ and such that $\Omega^\varepsilon \subset K$ for some
bounded set $K$, and if there exists constants $Q$ and $\delta$ such that for every
$x \in \partial \Omega^\varepsilon$
\begin{equation}
  \label{eq:tracebound}
  \sup \set{\frac{\mathcal{H}^{d-1}(\partial^* E \cap \partial^* \Omega^\varepsilon)}{\mathcal{H}^{d-1}(\partial^*
    E \cap \Omega^\varepsilon)} : E \subset \Omega^\varepsilon \cap B_\delta(x),
  \operatorname{Per}(E,\Omega_n) < \infty} \le Q,
\end{equation}
then the norm of $T : \RB\RV(\Omega^\varepsilon) \to \RL^1(\partial \Omega^\varepsilon)$ is uniformly bounded in $\varepsilon$. Here, $\partial^* E$
denotes the reduced bondary of $E$, which in general may be much smaller than
the topological boundary.

The next lemma is inspired by \cite{bgt}[Example 2]. Note that in their example,
$r_\varepsilon = \smallo{\varepsilon^{\frac{2d-1}{d-1}}}$, which means that the radius of the
holes if one order of magnitude in $\varepsilon$ smaller than
the critical level at which our holes are going to $0$. Parts of the proof are in a similar spirit but we need more
precise estimates separately around every boundary component.

\begin{lemma}
  \label{lem:tracebv}
  For all $\varepsilon>0$ sufficiently small, the norm of $T : \RB\RV(\Omega^\varepsilon) \to
  \RL^1(\partial\Omega^\varepsilon)$ is uniformly bounded in $\varepsilon$. 
\end{lemma}

\begin{proof}
  Following \eqref{eq:tracebound}, for any $x_\varepsilon \in \partial \Omega^\varepsilon$, we
  want to give a uniform upper bound for the ratio
  \begin{equation}
    \label{eq:ratiobounded}
    \frac{\mathcal{H}^{d-1}(\partial^* E \cap \partial \Omega^\varepsilon)}{\mathcal{H}^{d-1}(\partial^* E \cap
    \Omega^\varepsilon)}
  \end{equation}
  for sets $E \subset B_\delta(x_\varepsilon) \cap \Omega^\varepsilon$. Let us make the
  observation that for $\varepsilon>0$ small enough,
  \begin{equation}
    \#\set{\bn \in I^\varepsilon : Q_\bn^\varepsilon \cap B_\delta(x_\varepsilon) \ne \emptyset} \le
    2\omega_d\left(\frac \delta \varepsilon\right)^d
  \end{equation}
  and for some $M$, $\mathcal{H}^{d-1}(\partial \Omega \cap \overline{B_\delta(x_\varepsilon)}) \le
  M\delta^{d-1}$. On one hand, setting $\tilde \beta = \max(\beta,1)$, for any $E \subset \Omega^\varepsilon \cap B_\delta(x_\varepsilon)$
    of finite perimeter,
  \begin{equation}
    \label{eq:choicedelta}
    \mathcal{H}^{d-1}(\partial^* E \cap \partial \Omega^\varepsilon) \le \mathcal{H}^{d-1}(\partial \Omega^\varepsilon
    \cap \overline{B_\delta(x_\varepsilon)}) \le \left(M + 2d\tilde \beta
      \delta\right) \delta^{d-1}.
  \end{equation}
  On the other hand, we now need to find a lower bound for the denominator in
  \eqref{eq:ratiobounded} terms of the numerator. 
  By definition of $I^\varepsilon$, for all $\bn \in I^\varepsilon$ we have that $\partial
  \Omega^\varepsilon \cap Q_\bn^\varepsilon = \partial T_\bn^\varepsilon$. From this, we can decompose
  \begin{equation}
    \label{eq:decompositionboundary}
    \mathcal{H}^{d-1}(\partial^*E \cap \partial \Omega^\varepsilon) = \mathcal{H}^{d-1}(\partial^* E \cap \partial
    \Omega) + \sum_{\bn \in I^\varepsilon} \mathcal{H}^{d-1}(\partial^*E \cap \partial T_\bn^\varepsilon).
  \end{equation}
  We first observe that if $\delta$ is chosen such that $\omega_d \delta^d \le
  \abs{\Omega}/2$ and $\mathcal{H}^d(E)\leq 1$, the trace inequality for $BV(\Omega)\rightarrow L^1(\partial\Omega)$ and the relative isoperimetric inequality relative to $\Omega$
  \cite{bgt}[Inequality (2.2)]  applied to the characteristic function
  $\chi_E\in BV(\Omega)$  implies that there is a constant $C$, whose precise
  value may change from line to line but which depens only on $\Omega$, such that
  \begin{equation}
    \label{eq:extbdry}
    \begin{aligned}
      \mathcal{H}^{d-1}(\partial^*E \cap \partial \Omega)
      &\le C \mathcal{H}^{d-1}(\partial^*E \cap \Omega)\\& = C
      (\mathcal{H}^{d-1}(\partial^*E \cap (\partial \Omega^\varepsilon \setminus \partial \Omega)) +
      \mathcal{H}^{d-1}(\partial^*E \cap \Omega^\varepsilon))\\
      &=
      C(\sum_{\bn \in I^\varepsilon} \mathcal{H}^{d-1}(\partial^*E \cap \partial T_\bn^\varepsilon) +
      \mathcal{H}^{d-1}(\partial^*E \cap \Omega^\varepsilon))
    \end{aligned}
  \end{equation}
  Substitution in~\eqref{eq:decompositionboundary} implies
  \begin{equation}\label{eq:extbdry2}
    \mathcal{H}^{d-1}(\partial^*E \cap \partial \Omega^\varepsilon) \leq C\left(\mathcal{H}^{d-1}(\partial^* E \cap \Omega^\varepsilon) + \sum_{\bn \in I^\varepsilon} \mathcal{H}^{d-1}(\partial^*E \cap \partial T_\bn^\varepsilon)\right).
  \end{equation}  
  For $\bn \in I^\varepsilon$ and $t\in (0,\varepsilon/4)$, define
  \begin{equation}
    F_{\bn,t} := \{x\in E\,:\,\dist_{\partial T_\bn^\varepsilon}(x) \leq t\}\subset Q_\bn^\varepsilon.
  \end{equation}  
  Assume that there is $t \in (0,\varepsilon/4)$ such that
  \begin{equation}
    \label{eq:localisedboundary}
    \mathcal{H}^{d-1}(\partial^* F_{\bn,t} \cap Q_{\bn}^\varepsilon) \le 2\mathcal{H}^{d-1}(\partial^* E \cap
    Q_\bn^\varepsilon).
  \end{equation}
  Since projections on convex sets are nonexpansive and $F_{\bn,t} \subset
  Q_{\bn}^\varepsilon$, we have that
  \begin{equation}
    \label{eq:locbdry}
    \begin{aligned}
      \mathcal{H}^{d-1}(\partial^* E \cap \partial T_\bn^\varepsilon) &= \mathcal{H}^{d-1}(\partial^* F_{\bn,t} \cap \partial
      T_\bn^\varepsilon) \\
      &\le \mathcal{H}^{d-1}(\partial^* F_{\bn,t} \cap Q_{\bn}^\varepsilon) \\
      &\le 2\mathcal{H}^{d-1}(\partial^*E \cap  Q_{\bn}^\varepsilon).
    \end{aligned}
  \end{equation}
  If \eqref{eq:localisedboundary} holds for each $\bn\in I^\varepsilon$ then it follows from \eqref{eq:extbdry2} that
  \begin{equation}
    \begin{aligned}
      \mathcal{H}^{d-1}(\partial^*E \cap \partial \Omega^\varepsilon) &\leq C\left(\mathcal{H}^{d-1}(\partial^* E \cap \Omega^\varepsilon) + \sum_{\bn \in I^\varepsilon}\mathcal{H}^{d-1}(\partial^*E \cap  Q_{\bn}^\varepsilon) \right)\\
      &\leq C\mathcal{H}^{d-1}(\partial^* E \cap \Omega^\varepsilon),
    \end{aligned}
  \end{equation}
  which completes the proof in this case.
  
  Otherwise, let $J^\varepsilon := \set{\bn \in I^\varepsilon : \text{equation
      \eqref{eq:localisedboundary} does not hold}}$.
  For $\bn\in J^\varepsilon$, set
  \begin{equation}
    h_\bn(t) := \mathcal{H}^{d-1}\left(\{x\in \partial^*F_{\bn,t}\,:\, \dist_{\partial T_\bn^\varepsilon}(x) =
    t\}\right).
  \end{equation}
  Since \eqref{eq:localisedboundary} does not hold,
  $$2h_\bn(t)\geq \mathcal{H}^{d-1}(\partial^*F_{\bn,t} \cap Q_\bn^\varepsilon).$$
  It follows from the
  relative isoperimetric inequality with respect to $Q_\bn^\varepsilon$  that
  \begin{equation}
      \begin{aligned}
        c\mathcal{H}^d(F_{\bn,t})^{\frac{d-1}{d}}
        &\le \mathcal{H}^{d-1}(\partial^*F_{\bn,t} \cap Q_\bn^\varepsilon) \le 2 h_\bn(t).
      \end{aligned}
    \end{equation}
  The coarea formula gives $\partial_t \mathcal{H}^{d}(F_{\bn,t}) = h_\bn(t)$, and it follows by integration and the relative isoperimetric inequality with respect to $Q_\bn^\varepsilon$ that
  \begin{equation}
    \label{eq:lowerbound}
    \begin{aligned}
      d\mathcal{H}^d(F_{\bn,\varepsilon/4})^{1/d} &= \int_0^{\varepsilon/4} \frac{h_\bn(t)}{\mathcal{H}^{d}(F_{\bn,t})^{\frac{d-1}{d}}} \, \mathrm{d} t 
      \ge C \varepsilon.
    \end{aligned}
  \end{equation}
  That is for $\bn\in J^\varepsilon$,
  \begin{equation}
    \label{eq:firstpart}
  \begin{aligned}
    \mathcal{H}^d(E\cap Q_\bn^\varepsilon)&\geq\mathcal{H}^d(F_{\bn,\varepsilon/4})\\
    &\geq C\varepsilon^d\\
    &=\mathcal{H}^d(Q_\bn^\varepsilon).
  \end{aligned}
  \end{equation}
  This implies that
  for $\varepsilon>0$ small enough one has
  \begin{equation}
    \begin{aligned}
      \mathcal{H}^d(E\cap Q_\bn^\varepsilon)
      &\geq\frac{C}{\tilde \beta}\mathcal{H}^{d-1}(\partial T_\bn^\varepsilon).
    \end{aligned}
  \end{equation}  
  The isoperimetric inequality for $E\subset\mathbb R^d$ gives
  \begin{equation}
    \begin{aligned}
      \sum_{\bn \in J^\varepsilon}\mathcal{H}^{d-1}(\partial T_\bn^\varepsilon)&\leq\tilde \beta\sum_{\bn \in J^\varepsilon}\mathcal{H}^d(E \cap Q_\bn^\varepsilon)\\
      &\leq \tilde \beta\mathcal{H}^d(E)\\
      &\le C\tilde \beta
      \mathcal{H}^{d-1}(\partial^*E)^{\frac{d}{d-1}}\\
      &=C\tilde \beta\left(\mathcal{H}^{d-1}(\partial^*E \cap \partial \Omega) + \mathcal{H}^{d-1}(\partial^*E \cap
    \partial T^\varepsilon) + \mathcal{H}^{d-1}(\partial^*E \cap \Omega^\varepsilon)\right)^{\frac{d}{d-1}}.
    \end{aligned}
  \end{equation}
  Together with \eqref{eq:extbdry} and \eqref{eq:locbdry} this leads to the existence of a constant $C$ which can depend on $\beta$, such that
    \begin{equation}
      \label{eq:2ndpart}
      \begin{aligned}
        \sum_{\bn \in J^\varepsilon}\mathcal{H}^{d-1}(\partial T_\bn^\varepsilon)\le C \left(\mathcal{H}^{d-1}(\partial^*E \cap \Omega^\varepsilon) + \sum_{\bn \in J^\varepsilon}\mathcal{H}^{d-1}(\partial T_\bn^\varepsilon)\right)^{\frac{d}{d-1}}.
  \end{aligned}
    \end{equation}
    Because $J^\varepsilon\neq\emptyset$, dividing and factoring, this leads to
  \begin{align}
    1\le C\left(\sum_{\bn \in J^\varepsilon}\mathcal{H}^{d-1}(\partial T_\bn^\varepsilon)\right)^{\frac{1}{d-1}}\left(\frac{\mathcal{H}^{d-1}(\partial^*E \cap \Omega^\varepsilon)}{\sum_{\bn \in J^\varepsilon}\mathcal{H}^{d-1}(\partial T_\bn^\varepsilon)} + 1\right)^{\frac{d}{d-1}}.
  \end{align}
  For each $\bn\in J^\varepsilon$, $Q_\bn^\varepsilon\subset B_{2\delta}(x_\varepsilon)$ and it follows
  from \eqref{eq:choicedelta} that we can choose $\delta$ small
  enough depending on $\Omega$, the dimension, and $\beta$ but not on $\varepsilon$
  such that
  \begin{equation}
     C\left(\sum_{\bn \in J^\varepsilon}\mathcal{H}^{d-1}(\partial
     T_\bn^\varepsilon)\right)^{\frac{1}{d-1}}<\frac{1}{4}.
  \end{equation}
  This implies that there is a constant $c>0$ such that for all $\varepsilon > 0$,
  \begin{equation}
    \label{eq:globbdry}
    1 \le \frac{\mathcal{H}^{d-1}(\partial^*E \cap \Omega^\varepsilon)}{\sum_{\bn \in J^\varepsilon} \mathcal{H}^{d-1}(\partial T_\bn^\varepsilon)}.
  \end{equation}
  Combining \eqref{eq:extbdry2}, \eqref{eq:locbdry} and \eqref{eq:globbdry}, provides a constant $C$ such that
  \begin{equation}
    \mathcal{H}^{d-1}(\partial^*E \cap \partial \Omega^\varepsilon) \le C \mathcal{H}^{d-1}(\partial^*E \cap
    \Omega^\varepsilon).
  \end{equation}
\end{proof}


\begin{remark}
  When $r_\varepsilon = \smallo{\varepsilon^{\frac{d}{d-1}}}$, (i.e. in the subcritical
  regime), the previous result along with \cite{bgt}[Theorem 4.1] implies
  convergence of the eigenvalues of the Steklov problem on $\Omega^\varepsilon$ to the
  eigenvalues of the Steklov problem on $\Omega$. 
\end{remark}




\section{\bf Homogenisation of the Steklov problem} \label{sec:homogenisation}
 Let us first establish some basic facts related to the geometry of the
 homogenisation problem, under the assumption that $r_\varepsilon^{d-1} \varepsilon^{-d}
 \to \beta \in [0,\infty)$ as $\varepsilon \to 0$. The number of holes $N(\varepsilon)=\#I^\varepsilon$ satisfies
$$N(\varepsilon) \sim |\Omega|\varepsilon^{-d}  \qquad\mbox{ as }\varepsilon \to 0.$$
This implies,
\begin{gather}\label{eq:sizeholes}
  \abs{T^\varepsilon} = \sum_{\bk\in I^\varepsilon}|T_{\bk}^\varepsilon| = \bigo{r_\varepsilon}\qquad\mbox{and}\qquad
 \abs{\partial T^\varepsilon}=  \sum_{\bk\in I^\varepsilon}|\partial T_{\bk}^\varepsilon|\sim A_d \beta |\Omega|.
\end{gather}

The remainder of the section is split into three parts. In the first one we
extend the functions $u_k^\varepsilon$ to the whole of $\Omega$, in order to obtain
weak $\RH^1$ convergence, up to taking subsequences. In the second part, we
prove that those converging subsequences converge to solutions of Problem
\eqref{problem:homo}. Finally, we prove in the third part that the only
functions they can converge to are the corresponding eigenfunction in
\eqref{problem:homo}, implying convergence as $\varepsilon \to 0$, with this understood
in the sense of Remark \ref{rem:multiple} if the limit problem has eigenvalues
that are not simple.
\subsection{Extension of eigenfunctions} 
For $k \ge 1$, recall that $u_k^{\varepsilon}:\Omega^\varepsilon\rightarrow\mathbb R$ is the $k$'th Steklov
eigenfunction on $\Omega^\varepsilon$:
\begin{gather*}
  \begin{cases}
    \Delta u_k^{\varepsilon}=0&\mbox{ in }\Omega^\varepsilon,\\
    \partial_{\nu}u_k^{\varepsilon}=\sigma_k^{\varepsilon}u_k^{\varepsilon}&\mbox{ on
    }\partial\Omega^\varepsilon.
  \end{cases}
\end{gather*}
Recall also that the eigenfunctions $u_k^{\varepsilon}$ are normalized by requiring that
$$\int_{\partial\Omega^\varepsilon}\left(u_k^{\varepsilon}\right)^2 \, \mathrm{d} A=1.$$
Define the function $U_k^{\varepsilon}\in \RH^1(\Omega)$ to be the harmonic extension of
$u_k^{\varepsilon}$ to the interior of the holes:
\begin{gather*}
  \begin{cases}
    U_k^{\varepsilon}=u_k^{\varepsilon}&\mbox{ in }\overline{\Omega^\varepsilon},\\
    \Delta U_k^{\varepsilon}=0&\mbox{ in } T^\varepsilon.
  \end{cases}
\end{gather*}



  
\begin{lemma}\label{lemma:boundedUn}
  There is a sequence $\varepsilon_n \to 0$ such that $U_k^{\varepsilon_n}$ has a weak
  limit in $\RH^1(\Omega)$.
\end{lemma}
\begin{proof}

  It suffices to show that $\set{U_k^\varepsilon : 0 <\varepsilon \le 1}$ is bounded in
  $\RH^1(\Omega)$.   Recall that 
  \begin{equation}
  \norm{U_k^{\varepsilon}}^2_{H^1} = \norm{U_k^{\varepsilon}}^2_{\RL^2} +
  \mathcal{D}\left(U_k^{\varepsilon}\right).
\end{equation}
We first bound the $\RL^2$ norm of $U_{k}^{\varepsilon}$.
Let $\lambda$ be the first eigenvalue of the following Robin problem on $\Omega$:
\begin{gather*}
\begin{cases}
-\Delta u= \lambda u&\mbox{ in }\Omega,\\
\partial_\nu u= - u &\mbox{ on }\partial\Omega.
\end{cases}
\end{gather*}
It is well known (see \emph{e.g.} \cite{bucurfreitaskennedy}) that $\lambda>0$ and that it admits the following characterization:
$$\lambda=\inf_{v \in \RH^1(\Omega)} \frac{\int_\Omega|\nabla v|^2 \, \mathrm{d}
\bx+\int_{\partial\Omega}v^2\, \mathrm{d} S}{\int_{\Omega}v^2\, \mathrm{d} \bx}.$$
Applying this to $v=U_k^{\varepsilon}$ leads to
\begin{equation}
  \label{eq:L2bound}
\begin{aligned}\int_\Omega \left(U_k^{\varepsilon}\right)^2 \, \mathrm{d} \bx
  &\leq\frac{1}{\lambda}\left(\int_\Omega\left|\nabla U_k^{\varepsilon}\right|^2 \, \mathrm{d}
  \bx+\int_{\partial\Omega}\left(U_k^{\varepsilon}\right)^2 \, \mathrm{d} A\right)\\
  &\leq \frac{1}{\lambda}\left(\mathcal{D}\left(U_k^{\varepsilon}\right)+1\right).
  \end{aligned}
\end{equation}
  It is therefore sufficient to bound the Dirichlet energy. We first see that
\begin{equation}
  \label{eq:dirusep}
  \mathcal{D}\left(U_k^{\varepsilon};\Omega\right) = \mathcal{D}\left(u_k^{\varepsilon};\Omega^\varepsilon\right) +
  \mathcal{D}\left(U_k^{\varepsilon};T^\varepsilon\right)
  =\sigma_k^\varepsilon+ \mathcal{D}\left(U_k^{\varepsilon};T^\varepsilon\right).
\end{equation}
It follows from Lemma \ref{lemma:localEstimateFourier} and monotonicity of the Dirichlet energy,
that the contribution from the holes is
\begin{equation}
  \label{eq:dirutrous}
\begin{aligned}
\mathcal{D}\left(U_k^{\varepsilon};T^\varepsilon\right)&=\sum_{\bk\in
I^\varepsilon}\mathcal{D}\left(U_{k}^{\varepsilon};T_{\bk}^\varepsilon\right)\\
&\leq
\sum_{\bk\in I^\varepsilon} 5\left(\frac{r_\varepsilon}{\varepsilon}\right)^d
\mathcal{D}\left(u_k^{\varepsilon};Q_{\bk}^\varepsilon\right)
\left( 1 + \bigo{\left(\frac{ r_\varepsilon}{\varepsilon} \right)^d} \right)
\\
&\leq 
5 \left(\frac{r_\varepsilon}{\varepsilon}\right)^d \mathcal{D}\left(u_k^{\varepsilon};\Omega^\varepsilon\right) 
\left( 1 + \bigo{\left( \frac{r_\varepsilon}{\varepsilon} \right)^d} \right) \\
     &\leq C \sigma_k(\Omega^\varepsilon),
\end{aligned}
\end{equation}
for some constant $C$.  Combining \eqref{eq:dirutrous} and
  \eqref{eq:L2bound}, we see that, to bound $\norm{U_k^\varepsilon}_{\RH^1(\Omega)}$,
  it is sufficient to find a bound for 
  $\sigma_k^\varepsilon$ independent of $\varepsilon$. The variational characterisation for
  Steklov eigenvalues can be rewritten as
  \begin{equation}
    \sigma_k^\varepsilon = \min_{\substack{E \subset \RL^2(\partial \Omega^\varepsilon) \\ \dim(E)
    = k+1}} \max_{u \in E} \frac{\mathcal{D}(u)}{\norm{u}^2_{\RL^2(\partial\Omega^\varepsilon)}}.
  \end{equation}
  We use eigenfunctions of the dynamical eigenvalue problem as a test subspace
  for $\sigma_k^\varepsilon$. Namely, setting $E = \operatorname{span}(U_0,\dotsc,U_{k})$
  we see that for $\varepsilon$ small enough it spans a $k+1$ dimensional subspace of
  $\RL^2(\partial\Omega^\varepsilon)$. Indeed, they are an orthonormal set with respect to
  $(\cdot,\cdot)_\beta$, and for every $0 \le j,\ell \le k$,
  \begin{equation}
    (U_j,U_\ell)_{\partial^\varepsilon} \xrightarrow{\varepsilon \to 0} (U_j,U_\ell)_{\beta}.
  \end{equation}
  Therefore, from the characterisation of the eigenvalues $\Sigma_{k,\beta}$,
  \begin{equation}
    \begin{aligned}
    \sigma_k^\varepsilon &\le \max_{U \in E} \frac{\mathcal{D}(U)}{(U,U)_{\partial^\varepsilon}} \\
    &\le \Sigma_{k,\beta} + \smallo{1}.
  \end{aligned}
  \end{equation}
  This completes the proof that $\set{U_k^\varepsilon}$ is bounded so that there is
  a converging subsequence as $\varepsilon \to 0$.

\end{proof}

From now on we will abuse notation and relabel that sequence $\varepsilon_n \to 0$
along which $U_k^\varepsilon$ has a weak limit as $\varepsilon \to 0$ again.


\subsection{Establishing the limit problem}
Our aim by the end of this subsection is to prove the following weaker version
of Theorem \ref{thm:homogenisation}.
\begin{prop}\label{prop:weak}
  Let $k\in\mathbb N$.
  As $\varepsilon\to 0$, the pairs $\left(\sigma_k^\varepsilon,U_k^\varepsilon\right)$ converge to a solution
  $(\Sigma,U)$ of \eqref{problem:homo}, the convergence of the functions $U_k^\varepsilon$
  being weak in $\RH^1(\Omega)$. 
\end{prop}

Up to choosing a subsequence, we assume that $\sigma_k^{\varepsilon}$ converges to some
number $\Sigma$
and also that
$\set{U_k^{\varepsilon}}\subset \RH^1(\Omega)$ is weakly converging in $\RH^1(\Omega)$
to some $U \in \RH^1(\Omega)$, from which we
also get strong convergence to $U$ in $\RL^2(\Omega)$.
Considering the real-valued test function $V\in\RH^1(\Omega)$, we see that
\begin{align*}
  \int_\Omega \nabla U_k^{\varepsilon}\cdot\nabla V \, \mathrm{d} \bx
&=
\int_{\Omega^\varepsilon} \nabla u_k^{\varepsilon}\cdot\nabla V \, \mathrm{d} \bx+\int_{T^\varepsilon} \nabla
U_k^{\varepsilon}\cdot\nabla V \, \mathrm{d} \bx\\
&=\sigma_k^{\varepsilon}\int_{\partial\Omega^\varepsilon}u_k^{\varepsilon}V \, \mathrm{d} \bx+\int_{T^\varepsilon} \nabla
U_k^{\varepsilon}\cdot\nabla V \, \mathrm{d} \bx\\
&=\sigma_k^{\varepsilon}\int_{\partial\Omega}u_k^{\varepsilon}V \, \mathrm{d}
A+\sigma_k^{\varepsilon}\int_{\partial T^\varepsilon}u_k^{\varepsilon}V \, \mathrm{d} A +\int_{T^\varepsilon} \nabla
U_k^{\varepsilon}\cdot\nabla V \, \mathrm{d} \bx.
\end{align*}
Letting $\varepsilon\to0$ leads, if the limits exist, to
\begin{gather*}
\int_\Omega \nabla U\cdot\nabla V \, \mathrm{d} \bx-\Sigma\int_{\partial\Omega}UV
\, \mathrm{d} A=
\Sigma\lim_{\varepsilon \to 0}
\int_{\partial T^\varepsilon}u_k^{\varepsilon} V \, \mathrm{d} A+\lim_{\varepsilon \to
0}\int_{T^\varepsilon} \nabla
  U_k^{\varepsilon}\cdot\nabla V \, \mathrm{d}\bx.
\end{gather*}
It follows from the Cauchy--Schwarz inequality that
$$\int_{T^\varepsilon} \nabla U_k^{\varepsilon}\cdot\nabla V \, \mathrm{d} \bx\leq\left(\int_{T^\varepsilon} \left|\nabla
U_k^{\varepsilon}\right|^2 \, \mathrm{d} \bx\int_{T^\varepsilon}|\nabla V|^2 \, \mathrm{d} \bx\right)^{1/2},$$
which tends to 0 according to Lemma \ref{lemma:localEstimateFourier}.
It follows that
\begin{gather*}
\int_\Omega \nabla U\cdot\nabla V \, \mathrm{d} \bx-\Sigma\int_{\partial\Omega}UV
\, \mathrm{d} A=
\Sigma\lim_{\varepsilon \to 0}\int_{\partial T^\varepsilon}u_k^{\varepsilon}V \, \mathrm{d} A,
\end{gather*}
and all that is left to do is to analyse the last term.

\begin{prop}\label{prop:limitproblem}
  Suppose that $\varepsilon^{-d}r_\varepsilon^{d-1}\to\beta\geq 0$. Then,
	for each $V\in \RH^1(\Omega)$ the following holds,
    \begin{equation}
      \label{eq:toprove}
      \lim_{\varepsilon \to 0}\int_{\partial T^\varepsilon}u_k^{\varepsilon}V \, \mathrm{d} A=
      A_d\beta\int_{\Omega}UV \, \mathrm{d} \bx.
    \end{equation}
\end{prop}


\begin{remark}
    The functional $V\mapsto\int_{\partial T^\varepsilon}u_k^{\varepsilon}V$ is bounded on $H^1(\Omega)$. By the 
    Riesz--Fr\'echet representation theorem, there exists
    a function $\xi^\varepsilon\in H^1(\Omega)$ such that
    $$\int_{\partial T^\varepsilon}u_k^{\varepsilon}V =\int_{\Omega}\nabla\xi^\varepsilon\cdot
    \nabla V+\xi^\varepsilon V\,dx\qquad\forall V\in H^1(\Omega).$$
    Using appropriate test functions shows that $\xi^\varepsilon$ is the weak solution of the following transmission problem:
    \begin{gather}\label{Problem:transmission}
      \begin{cases}
        \Delta\xi^\varepsilon=0&\mbox{ in }\Omega^\varepsilon\cup T^\varepsilon,\\
        \partial_\nu \xi^\varepsilon_++\partial_\nu \xi^\varepsilon_-=u_k^\varepsilon&\mbox{ on }\partial T^\varepsilon,\\
        \partial_\nu \xi^\varepsilon=0&\mbox{ on }\partial\Omega.
      \end{cases}
    \end{gather}
    Proposition~\ref{prop:limitproblem} is an homogenisation result for this problem. It means that in the limit as $\varepsilon\to 0$, the solution converges to that of the following problem:
    \begin{gather}
      \begin{cases}
        -\Delta \Xi+(1-A_d\beta)\Xi=0&\mbox{ in }\Omega,\\
        \partial_\nu \Xi=0&\mbox{ on }\partial\Omega.
      \end{cases}
    \end{gather}
    Transmission problems have recently been the subject of investigation through means of homogenisation, see for example \cite{dmfz}.
  \end{remark}

   The proof of Proposition \ref{prop:limitproblem} is divided in three main
  steps. In the first step, we justify that we can use smooth test functions in
the limit \eqref{eq:toprove}. In order to do so, we will use an inner
representation of the lefthandside in \eqref{eq:toprove}, defined in terms of
extensions to the holes, to show that it is bounded, uniformly in
$\varepsilon$. In this representation, however, it is hard to explicitly compute the
limit problem. 

In our second step, we introduce an outer representation of the
lefthandside in \eqref{eq:toprove}, through integration
on the outside of the holes. This representation is given in terms of an
auxilliary function $\Psi$, for which we derive some regularity properties.

In the final step, we use this latter representation to show that the limit
\eqref{eq:toprove} indeed holds. Here, we reap rewards from the previous steps
and use explicitly the properties of the
auxilliary function $\Psi$, as well as better estimates awarded from the fact
that we can test against smooth functions.

\begin{proof}[Proof of Proposition \ref{prop:limitproblem}]

  Define the family of bounded functionals $L_\varepsilon:\RH^1(\Omega)\rightarrow\mathbb R$ by
  $$L_\varepsilon(V):=\int_{\partial T_\varepsilon} V \, \mathrm{d} A.$$
  
	\noindent\textbf{Step 1}: Inner representation of $L_\varepsilon$.
        
        Define $\phi_\varepsilon:\mathbb R^d\rightarrow\mathbb R$ by
        \begin{equation}\label{eq:explicitphi}
          \phi_\varepsilon(\bx) =
          \begin{cases}
            \frac{ \abs \bx^2}{2r_\varepsilon}&\text{for }x\in B(0,r_\varepsilon),\\
            0&\text{elsewhere}.
          \end{cases}
        \end{equation}
        By periodizing along $\varepsilon\mathbb Z^d$ we obtain the function
        $\Phi_\varepsilon : \mathbb R^d \to \mathbb R$ given by
	\begin{equation}
	  \Phi_\varepsilon(\bx) := \sum_{\bk \in I^\varepsilon} \phi_\varepsilon(\bx - \varepsilon \bk).
	\end{equation}
        \begin{lemma}\label{lemma:innerrep}
          The functional $L_\varepsilon:\RH^1(\Omega)\rightarrow\mathbb R$
          admits the following representation:
          \begin{equation}\label{eq:innerrep}
            L_\varepsilon(V) = \frac{d}{r_\varepsilon} \int_{T^\varepsilon} V \, \mathrm{d} \bx + \int_{T^\varepsilon} \nabla \Phi_\varepsilon \cdot \nabla V \, \mathrm{d} \bx.
          \end{equation}
        \end{lemma}
        \begin{proof}
          It is straightforward to check that
          \begin{equation}
	    \begin{cases}
	      \Delta \Phi_\varepsilon = \frac{d}{r_\varepsilon} & \text{in } T^\varepsilon,\\
	      \partial_\nu \Phi_\varepsilon = 1 & \text{on } \partial T^\varepsilon,\\
	      \Phi = 0 & \text{in } \Omega^\varepsilon.
	    \end{cases}
	  \end{equation}
          The function $\Phi_\varepsilon$ therefore satisfies the weak identity
          \begin{equation}
            \int_{T^\varepsilon} \nabla \Phi_\varepsilon \cdot V = - \frac{d}{r_\varepsilon} \int_{T^\varepsilon} V + \int_{\partial T^\varepsilon} V,\qquad\forall V\in H^1(\Omega).
          \end{equation}
        \end{proof}
    
For each $k\in\mathbb N$ and $\varepsilon>0$, the functional $\tilde{L_\varepsilon}:H^1(\Omega)\rightarrow\mathbb R$ is defined by
$$\tilde{L_\varepsilon}(V):=L_\varepsilon(u_k^{\varepsilon}V).$$
\begin{lemma} \label{lem:banstein}
There is an $\varepsilon_0>0$ such that the family
$\left\{\tilde{L_\varepsilon}\right\}_{\varepsilon>0}\subset \left(H^1(\Omega)\right)^*$
is uniformly bounded for $0 < \varepsilon \le \varepsilon_0$.
\end{lemma}
\begin{proof}
Given $V \in \RH^1(\Omega)$,
  it follows from the Lemma~\ref{lemma:innerrep} that
  \begin{equation}
	\label{eq:innerrepbound}
	\begin{aligned}
		L_\varepsilon(U_k^\varepsilon V) = \frac{d}{r_\varepsilon} \int_{T^\varepsilon} U_k^\varepsilon V \, \mathrm{d} \bx + \int_{T^\varepsilon} \nabla \Phi_\varepsilon \cdot \nabla(U_k^\varepsilon V) \, \mathrm{d} \bx.
	\end{aligned}
\end{equation}
To bound the first term, start by using the Cauchy--Schwarz inequality:
\begin{equation}\label{eq:innerrepbound2}
	\abs{\frac{d}{r_\varepsilon}\int_{T^\varepsilon} U_k^\varepsilon V \, \mathrm{d} \bx} \le \frac{d}{r_\varepsilon}\norm{U_k^\varepsilon}_{\RL^2(T^\varepsilon)} \norm{V}_{\RL^2(T^\varepsilon)}.
\end{equation}
It follows from Lemma~\ref{lemma:smalltolarge} that there is $C > 0$ depending
only on $\beta$ such that
  \begin{equation}\label{ineq:demiborne1}
    \begin{aligned}
      \norm{V}_{\RL^2(T^\varepsilon)} &\le C r_\varepsilon^{1/2} \norm{V}_{\RH^1(\Omega)},\\
      \norm{U_{k}^\varepsilon}_{\RL^2(T^\varepsilon)} &\le C r_\varepsilon^{1/2}
      \norm{U_k^\varepsilon}_{\RH^1(\Omega)},
  \end{aligned}
  \end{equation}
  so that
  \begin{equation}
    \label{ineq:demiborne2}
    \sup_{\varepsilon \in (0,1)} \abs{\frac{d}{r_\varepsilon} \int_{T^\varepsilon}U_k^\varepsilon V \, \mathrm{d} \bx
    } < \infty.
  \end{equation}

To bound the second term in~\eqref{eq:innerrepbound}, the generalised Hölder inequality leads to
\begin{equation*}
  \begin{aligned}
    \abs{\int_{T^\varepsilon} \nabla \Phi_\varepsilon \cdot \nabla(U_k^\varepsilon V) \, \mathrm{d} \bx}
    &\le \abs{\int_{T^\varepsilon} U_k^\varepsilon \nabla \Phi_\varepsilon \cdot \nabla V + V \nabla \Phi_\varepsilon \cdot \nabla U_k^\varepsilon \, \mathrm{d} \bx} \\
    &\le \norm{U_k^\varepsilon}_{\RL^2(T^\varepsilon)} \norm{\nabla V}_{\RL^2(T^\varepsilon)}\norm{\nabla \Phi_\varepsilon}_{\RL^\infty(T^\varepsilon)}\\
    &\qquad+ \norm{\nabla U_k^\varepsilon}_{\RL^2(T^\varepsilon)} \norm{V}_{\RL^2(T^\varepsilon)} \norm{ \nabla \Phi_\varepsilon}_{\RL^\infty(T^\varepsilon)}\\
    &\leq\left(\norm{U_k^\varepsilon}_{\RL^2(T^\varepsilon)}+\norm{\nabla U_k^\varepsilon}_{\RL^2(T^\varepsilon)}\right)\norm{\nabla V}_{H^1(\Omega)}\norm{\nabla \Phi_\varepsilon}_{\RL^\infty(T^\varepsilon)}\\
    &\leq\norm{U_k^\varepsilon}_{H^1(T^\varepsilon)}\norm{\nabla V}_{H^1(\Omega)}.
  \end{aligned}
\end{equation*}
In the last inequality we have used $\norm{\nabla \Phi_\varepsilon}_{\RL^\infty(T^\varepsilon)}=1$, which follows from~\eqref{eq:explicitphi}.
This quantity is uniformly bounded as $\epsilon\searrow 0$ since we have shown in the proof of Lemma~\ref{lemma:boundedUn} that  $U_k^\varepsilon$ is bounded in $H^1(\Omega)$.
Together with~\eqref{ineq:demiborne2} this proves for each $V\in H^1(\Omega)$ the existence of a constant $C$ such that
$|L_\varepsilon(U_k^\varepsilon V)|\leq C\|V\|$ for each $\varepsilon$,
and the conclusion follows from the Banach--Steinhaus theorem.
\end{proof}


\medskip
\noindent\textbf{Step 2}: Outer representation of $L_\varepsilon$.
Consider the torus $\mathcal{C}=\mathbb{T}^d = \mathbb R^d / \mathbb Z^d$ and introduce the fundamental cell $\mathcal{C}^\varepsilon$ as the perforated torus
$$\mathcal{C}^\varepsilon:=\mathcal{C}\setminus \overline{B\left(0,\rho_\varepsilon\right)},$$
where $\rho_\varepsilon:= \varepsilon^{-1} r_\varepsilon$ is the renormalised radius.
Following \cite{vanninathan}, we define the function 
$\psi_\varepsilon\in \RH^1(\mathcal{C}^\varepsilon)$ through the weak
variational problem:
\begin{equation}
  \label{eq:psivar}
  \int_{\mathcal{C}^\varepsilon}\nabla\psi_\varepsilon\cdot\nabla
  V=-c_\varepsilon\int_{\mathcal{C}^\varepsilon}V+\int_{\partial B(0,\rho_\varepsilon)}V.
\end{equation}
By taking $V\equiv 1$, one sees that the necessary and sufficient condition for existence
of a solution (see e.g. \cite[Theorem 5.7.7]{TaylorI}) is 
$$c_\varepsilon=\frac{A_d \rho_\varepsilon^{d-1}}{|\mathcal{C}^\varepsilon|}\sim A_d (\rho_\varepsilon)^{d-1}.$$
Uniqueness of the solution is guaranteed by requiring that $\psi_\varepsilon$
be orthogonal to constants on $\mathcal{C}^\varepsilon$. Therefore, $\psi_\varepsilon$ is the
unique function such that
\begin{gather*}
  \begin{cases}
    \Delta\psi_\varepsilon=c_\varepsilon&\mbox{ in }\mathcal{C}^\varepsilon\\
    \partial_\nu\psi_\varepsilon=1&\mbox{ on }\partial B(0,\rho_\varepsilon)
  \end{cases}
  \qquad\mbox{ and }\qquad
  \int_{\mathcal{C}^\varepsilon}\psi_\varepsilon=0.
\end{gather*}

Consider the union of all cells strictly contained in $\Omega$,
$$\tilde\Omega^\varepsilon:=\bigcup_{\bk\in I^\varepsilon}Q_{\bk}^\varepsilon\subset\Omega^\varepsilon.$$
Define the function $\Psi_\varepsilon: \mathbb R^d \setminus \bigcup_{\bk \in \mathbb Z^d}
T_{\bk}^\varepsilon \rightarrow\mathbb R$
as the scaled lift of $\psi_\varepsilon$. That is, if $q : \mathbb R^d \to \mathbb R^d/\mathbb Z^d$ is the
covering map, then 
$$\Psi_\varepsilon(\varepsilon \bx):= \psi_\varepsilon(q(\bx)).$$
This function satisfies
\begin{gather*}
\begin{cases}
\Delta\Psi_\varepsilon = \varepsilon^{-2} c_\varepsilon &\mbox{ in }\tilde\Omega^\varepsilon,\\
\partial_\nu\Psi_\varepsilon= \varepsilon^{-1} &\mbox{ on }\partial T^\varepsilon.
\end{cases}
\end{gather*}
\begin{lemma}The functional  $L_\varepsilon:\RH^1(\Omega)\rightarrow\mathbb R$ 
admits the following representation:
$$L_\varepsilon(V)= \varepsilon\int_{\tilde\Omega^\varepsilon}\nabla\Psi_\varepsilon\cdot\nabla V \, \mathrm{d} \bx
+ \varepsilon \int_{\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon} V \partial_\nu \Psi_\varepsilon
\, \mathrm{d} A
+\varepsilon^{-1}c_\varepsilon\int_{\tilde\Omega^\varepsilon}V \, \mathrm{d} A.$$
\end{lemma}

The proof is immediate since for $V\in \RH^1(\Omega)$ the following holds:
\begin{equation}
\int_{\tilde\Omega^\varepsilon}\nabla\Psi_\varepsilon\cdot\nabla V \, \mathrm{d} \bx=
-\varepsilon^{-2}c_\varepsilon\int_{\tilde\Omega^\varepsilon}V \, \mathrm{d} \bx+\varepsilon^{-1}\int_{\partial
T^\varepsilon}V \, \mathrm{d} A
+ \int_{\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon} V \partial_\nu \Psi_\varepsilon \, \mathrm{d}
A.
\end{equation}







We establish the following claim concerning $\psi_\varepsilon$. 
\begin{lemma}
  \label{lemma:sobolevpsieps} 
  There is a constant $C$, depending only on the dimension and on $\beta$, such that
  \begin{equation}
    \norm{\psi_\varepsilon}_{\RH^1(\mathcal{C}^\varepsilon)} \le C \varepsilon^{\frac 1 2 + \frac 1 d}.
  \end{equation}
  Furthermore, for any $s > 1$, any compact set $K \subset \mathcal{C}$, containing the origin in its
  interior, there is a constant $C'$ depending only on the dimension,
  $\beta$, $s$,
  and on $K$ such that
  \begin{equation}
    \label{eq:hspsi}
      \norm{\psi_\varepsilon}_{\RH^s(\mathcal{C}^\varepsilon \setminus K)} \le  C' \varepsilon^{\frac 1
        2 + \frac 1 d}.
  \end{equation}
  In particular, this implies $\norm{D^\alpha \psi_\varepsilon}_{\RL^\infty(\mathcal{C}^\varepsilon \setminus
    K)}$ decays
  as $\varepsilon^{\frac 1 2 + \frac 1 d}$ for any multi-index $\alpha$. 
\end{lemma}
\begin{proof}
  Observe that since $\psi_\varepsilon$ has mean $0$ on $\mathcal{C}^\varepsilon$, 
  the Poincaré--Wirtinger inequality implies
  \begin{equation}
    \label{eq:pw}
    \mu_1(\mathcal{C}^\varepsilon) \norm{\psi_\varepsilon}^2_{\RL^2(\mathcal{C}^\varepsilon)} \le \norm{\nabla
      \psi_\varepsilon}^2_{\RL^2(\mathcal{C}^\varepsilon)^d},
  \end{equation}
  where $\mu_1(\mathcal{C}^\varepsilon)$ is the first non-zero Neumann eigenvalue of $\mathcal{C}^\varepsilon$. Observe that
  $$\mu_1(\mathcal{C}^\varepsilon) \to 4\pi^2\qquad\text{as }\varepsilon \to 0.$$
  Indeed, $\mu_1(\mathcal{C}^\varepsilon)$ 
  is the first non-zero Neumann eigenvalue of a punctured $d$-dimensional
  torus, which is known to converge to the first nonzero eigenvalue of the torus
  itself as $\varepsilon\searrow 0$. See~\cite[Chapter IX]{chaveleig} for instance.
  
  Take $V = \psi_\varepsilon$ in the variational characterisation
  \eqref{eq:psivar} of $\psi_\varepsilon$, and consider $\varepsilon$ to be small enough
  that $\mu_1(\mathcal{C}^\varepsilon) \ge 1$. Using the Cauchy--Schwarz inequality
  yields
  \begin{equation}
    \label{eq:gradpsi}
    \begin{aligned}
      \norm{\psi_\varepsilon}_{\RH^1(\mathcal{C}^\varepsilon)}^2 &\le 2\int_{\mathcal{C}^\varepsilon}|\nabla\psi_\varepsilon|^2\, \mathrm{d} \bx  \\
      &\le 2\int_{\partial
        B(0,\rho_\varepsilon)}\psi_\varepsilon \, \mathrm{d} A\\
      &\leq 2\sqrt{A_d \rho_\varepsilon^{d-1}}\|\psi_\varepsilon\|_{\RL^2(\partial
        B(0,\rho_\varepsilon))}\\
      &\leq 2\sqrt{A_d\rho_\varepsilon^{d-1}}\|\tau_\varepsilon\|\|\psi_\varepsilon\|_{\RH^1(\mathcal{C}^\varepsilon)},
    \end{aligned}
  \end{equation}
  where $\tau_\varepsilon$ is the trace operator $\RH^1(\mathcal{C}^\varepsilon) \to \RL^2(\partial
  B(0,\rho_\varepsilon))$. From the definition of $\tau_\varepsilon$ and monotonicity of the
  involved integrals, we have that
\begin{equation}
\norm{\tau_\varepsilon} \le \gamma\left(\varepsilon^{\frac{1}{d-1}},1\right)^{-1},
\end{equation}
where $\gamma$ is defined in Lemma \ref{lemma:Sobolevconst}. We therefore deduce
that $\norm{\tau_\varepsilon} \ll \varepsilon^{1/d}$, the implicit constant depending only on
the dimension and on $\beta$. Finally, dividing both sides in \eqref{eq:gradpsi}
by $\norm{\psi_\varepsilon}_{\RH^1}$, observing that $\rho_\varepsilon^{d-1} \sim \beta\varepsilon$ and
inserting the bound for $\norm{\tau_\varepsilon}$ in \eqref{eq:gradpsi} finishes the
proof of the first inequality.

For the second one, we have from \cite[Theorem 8.10]{gilbargtrudinger} that there is
a constant $\tilde C$ depending only on the dimension, on $K$, and on $s$ such that
\begin{equation}
  \norm{\psi_\varepsilon}_{\RH^s(\mathcal{C}^\varepsilon \setminus K)} \le \tilde C \left(
  \norm{\psi_\varepsilon}_{\RH^1(\mathcal{C}^\varepsilon)} + \norm{c_\varepsilon}_{\RH^s(\mathcal{C}^\varepsilon)}
\right). 
\end{equation}
The second term is of order $\rho_\varepsilon^{d-1} \sim \varepsilon$, and bounds for the
$\RH^1$ norm of $\psi_\varepsilon$ that were obtained in \eqref{eq:gradpsi} conclude
the proof of the second inequality.

As for the remark considering the $\RL^\infty$ bounds on derivatives of
$\psi_\varepsilon$, we have from inequality \eqref{eq:hspsi} that
\begin{equation}
  \norm{D^\alpha \psi_\varepsilon}_{\RH^{s-\abs \alpha}(\mathcal{C}^\varepsilon \setminus K)} 
  \le \norm{\psi_\varepsilon}_{\RH^s(\mathcal{C}^\varepsilon \setminus K)} \le
  C'\varepsilon^{\frac 1 2 + \frac 1 d}.
\end{equation}
Choosing $s > \abs \alpha + \frac{d+1}{2}$, and $\varepsilon$ small enough that
$\mathcal{C}^\varepsilon \setminus K = \mathcal{C} \setminus K$, we have that
$$\RH^{s- \abs \alpha}(\mathcal{C}^\varepsilon\setminus K) \hookrightarrow \RL^\infty(\mathcal{C}^\varepsilon
\setminus K)$$
with norm  independent of $\varepsilon$, concluding the proof.
    \end{proof}


\noindent
\textbf{Step 3}: Computing the limit problem.
It follows from Lemma \ref{lem:banstein} that it is sufficient to verify convergence for a test function $V \in C^\infty(\overline\Omega)$. From the outer representation for $L_\varepsilon$ we have that
\begin{equation}
  \label{eq:outerreptilde}
  \begin{aligned}
    \tilde L_\varepsilon(V) &= \int_{\partial T^\varepsilon}u_k^{\varepsilon}V \\
    &= \varepsilon \int_{\tilde\Omega^\varepsilon}\nabla\Psi_\varepsilon\cdot\nabla (u_k^{\varepsilon}V) +
    \int_{\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon} u_k^\varepsilon V \partial_\nu
    \Psi_\varepsilon
    +\underbrace{\varepsilon^{-1}c_\varepsilon\int_{\tilde\Omega^\varepsilon}u_k^{\varepsilon}V}_{\to
    A_d\beta\int_\Omega UV},
\end{aligned}
\end{equation}
where convergence of the last term stems from strong $\RL^2$ convergence.
We now show that the two other terms converge to $0$. For the first one, 
the Cauchy--Schwarz inequality gives
\begin{equation}
  \label{eq:afterCS}
  \int_{\tilde\Omega^\varepsilon} \varepsilon \nabla\Psi_\varepsilon\cdot\nabla (u_k^{\varepsilon}V)\leq
  \left(\int_{\tilde\Omega^\varepsilon} \varepsilon^{2}|\nabla\Psi_\varepsilon|^2
  \int_{\tilde\Omega^\varepsilon}|\nabla (u_k^{\varepsilon}V)|^2\right)^{1/2}.
\end{equation}
Let us first observe that since $V$ is in $ C^\infty(\overline \Omega)$, we have that 
\begin{equation}
  \begin{aligned}
    \int_{\tilde \Omega^\varepsilon} \abs{\nabla(u_k^{\varepsilon} V)}^2 
    &\le 2 \sup_{\bx \in \tilde \Omega^\varepsilon} \left( \abs{V(\bx)}^2 + \abs{\nabla V(\bx)}^2
    \right) \|U_k^{\varepsilon}\|_{\RH_1(\Omega)} \\
    &= \bigo 1.
  \end{aligned}
  \label{eq:boundgraduv}
\end{equation}
On to the other term, note that
\begin{equation}
  \label{eq:Psitopsi}
\begin{aligned}
\int_{\tilde\Omega^\varepsilon} \varepsilon^2 |\nabla\Psi_\varepsilon|^2&=
\varepsilon^2 \sum_{\bk\in I^\varepsilon}\int_{Q^\varepsilon_{\bk}} |\nabla\Psi_\varepsilon|^2 \\
&=\sum_{\bk\in I^\varepsilon}\varepsilon^{d}\int_{\mathcal{C}^\varepsilon}|\nabla\psi_\varepsilon|^2 \\
&\sim |\Omega|\int_{\mathcal{C}^\varepsilon}|\nabla\psi_\varepsilon|^2.
\end{aligned}
\end{equation}
It follows from Lemma \ref{lemma:sobolevpsieps} that 
\begin{equation}
  \int_{\mathcal{C}^\varepsilon} \abs{\nabla \psi_\varepsilon}^2 \, \mathrm{d} \bx = \bigo{\varepsilon^{1 + \frac 2 d}},
\end{equation}
hence that term indeed goes to $0$. For the second term in
\eqref{eq:outerreptilde}, we have from the generalised Hölder inequality that
\begin{equation}
  \varepsilon \int_{\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon} u_k^\varepsilon V \partial_\nu
  \Psi_\varepsilon \, \mathrm{d} A \le \varepsilon \norm{V}_{\RL^2(\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon)}
  \norm{u_k^\varepsilon}_{\RL^2(\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon)}
\norm{\partial_\nu \Psi_\varepsilon}_{\RL^\infty(\partial \tilde \Omega^\varepsilon \setminus \partial
T^\varepsilon)},
\end{equation}
we now analyse each of those norms. 

First, by scaling of derivatives we have
\begin{equation}
\begin{aligned}
 \norm{\varepsilon \partial_\nu \Psi_\varepsilon}_{\RL^\infty(\partial \Omega^\varepsilon \setminus \partial T^\varepsilon)} &\le 
 \sup_{\abs \alpha = 1}\norm{D^\alpha \psi_\varepsilon}_{\RL^\infty(\mathcal{C}^\varepsilon \setminus B(0,1/4))} \\
 &\ll \varepsilon^{\frac 1 2 + \frac 1 d},
 \end{aligned}
\end{equation}
where the last bound holds by Lemma \ref{lemma:sobolevpsieps}. We move on to
$\norm{u_k^\varepsilon}_{\RL^2(\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon)}$. Denote
  by $\tilde I^\varepsilon$ the set of indices $\bn \in \mathbb Z^d $ such that $\partial
  Q_\bn^\varepsilon \cap \partial \tilde \Omega^\varepsilon \ne \varnothing$. One can see that
\begin{equation}
  \label{eq:lastestimate}
  \norm{u_k^\varepsilon}_{\RL^2(\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon)} \le \sum_{\bn \in \tilde I^\varepsilon}\norm{U_k^\varepsilon}_{\RL^2(\partial Q_\bn^\varepsilon)}.
\end{equation}
On the other hand, it follows from Lemma \ref{lem:traceconstant} that there is a constant $C$ (in fact, the trace constant of the unit cube) such that
\begin{equation}
\begin{aligned}
 \sum_{\bn \in \tilde I^\varepsilon}\norm{U_k^\varepsilon}_{\RL^2(\partial Q_\bn^\varepsilon)} &\le C\varepsilon^{-1/2} \sum_{\bn \in \tilde I^\varepsilon}\norm{U_k^\varepsilon}_{\RH^1(Q_{\bn}^\varepsilon)} \\
 &\le C\varepsilon^{-1/2} \norm{U_k^\varepsilon}_{\RH^1(\Omega)}.
 \end{aligned}
\end{equation}
Finally, since $V$ is fixed and smooth we have that
\begin{equation}
\begin{aligned}
 \norm{V}_{\RL^2(\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon)} &\le \abs{\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon}^{1/2} \sup_{\bx} V(\bx) \\
 &\le C \abs{\partial \Omega}^{1/2} \sup_\bx V(\bx).
 \end{aligned}
\end{equation}
All in all, this implies that the second term in \eqref{eq:outerreptilde} is bounded by a constant times $\varepsilon^{1/d}$, hence goes to $0$
as $\varepsilon \to 0$, finishing the proof of Lemma \ref{prop:limitproblem}.
\end{proof}



We are now ready to complete the proof of the main result of this subsection.
\begin{proof}[Proof of Proposition~\ref{prop:weak}]
  We know from \cite{BelowFrancois} that we only need to show that the solutions
converge to a solution $(\Sigma,U)$ of the weak formulation of Problem \ref{problem:homo},1
  which is that for any test function $V$,
  \begin{equation}
    \int_\Omega \nabla U \cdot \nabla V \, \mathrm{d} \bx = \Sigma\left(A_d \beta
    \int_\Omega U V \, \mathrm{d} \bx + \int_{\partial \Omega} U V \, \mathrm{d} A\right).
  \end{equation}
  The weak formulation of the Steklov problem on $\Omega^\varepsilon$ is that for all test functions $V$,
  \begin{equation}
    \int_{\Omega^\varepsilon} \nabla u_k^\varepsilon \cdot \nabla V \, \mathrm{d} \bx = \sigma_k^\varepsilon
    \left(\int_{\partial \Omega} u_k^\varepsilon V \, \mathrm{d} A + \int_{\partial T^\varepsilon} u_k^\varepsilon V
    \, \mathrm{d} A\right).
  \end{equation}
  The convergence of the gradient terms follows from weak convergence in $\RH^1$
  of $U_k^\varepsilon$ to $U$ and Lemma \ref{lemma:localEstimateFourier}. We already
  have that $\sigma_k^\varepsilon \to \Sigma$. The integrals on $\partial \Omega$
  converge by weak convergence in $\RH^1$ and compactness of the trace operator
  on $\partial \Omega$. Finally, convergence of the interior term comes from
  Proposition
  \ref{prop:limitproblem}.
\end{proof}


\subsection{Spectral convergence of the problem}

We need the following technical lemma.
\begin{lemma} \label{lem:normconv}
 As $\varepsilon \to 0$, we have that
  \begin{equation}
    \tilde L_\varepsilon(u_k^\varepsilon) \to A_d \beta \int_{\Omega} U^2 \, \mathrm{d} \bx.
  \end{equation}
\end{lemma}  
\begin{remark}
  Observe that this situation is specific to the sequence $u_k^\varepsilon$.
  Indeed, there are sequences $\set{v_\varepsilon}$ converging
  weakly to some $v \in \RH^1(\Omega)$ such that
  \begin{equation}
  \lim_{\varepsilon \to 0} \tilde L_\varepsilon (v_\varepsilon) \ne A_d \beta \int_\Omega U v \, \mathrm{d}
  \bx.
  \end{equation}
\end{remark}
\begin{proof}

  From the outer representation \eqref{eq:outerreptilde}, we have that
\begin{equation}
  \begin{aligned}
    \tilde L_\varepsilon(u_k^\varepsilon) 
    &= \varepsilon \int_{\tilde\Omega^\varepsilon}\nabla\Psi_\varepsilon\cdot\nabla (u_k^{\varepsilon})^2
    \, \mathrm{d} \bx +
\varepsilon \int_{\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon} (u_k^\varepsilon)^2 \partial_\nu
\Psi_\varepsilon \, \mathrm{d} \bx +
    \varepsilon^{-1}c_\varepsilon\int_{\tilde\Omega^\varepsilon}(u_k^{\varepsilon})^2 \, \mathrm{d}
  \bx, \\
    &= \varepsilon \int_{\tilde\Omega^\varepsilon}2\nabla\Psi_\varepsilon\cdot u_k^\varepsilon\nabla
    u_k^\varepsilon \, \mathrm{d} \bx
+\varepsilon \int_{\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon} (u_k^\varepsilon)^2 \partial_\nu \Psi_\varepsilon \, \mathrm{d} \bx
    +\varepsilon^{-1}c_\varepsilon\int_{\tilde\Omega^\varepsilon}(u_k^{\varepsilon})^2 \, \mathrm{d}
  \bx. \\
  \end{aligned}
\end{equation}
The last term converges towards the desired $A_d \beta \int_\Omega U^2 \, \mathrm{d} \bx$
once again by strong $L^2$ convergence of the sequence $U_k^\varepsilon$. 
To study the first term, let us now introduce the sets
\begin{equation}
 \omega^\varepsilon := \bigcup_{\bk \in I^\varepsilon} B\left(\varepsilon \bk,\frac \varepsilon 4\right) \setminus B(\varepsilon\bk,r_\varepsilon) \subset \tilde \Omega^\varepsilon,
\end{equation}
and decompose
\begin{equation}
  \varepsilon \int_{\tilde \Omega^\varepsilon} \nabla \Psi_\varepsilon \cdot \nabla(u_k^\varepsilon)^2
  \, \mathrm{d} \bx =
  2 \left(\int_{\omega^\varepsilon} + \int_{\tilde \Omega^\varepsilon \setminus \omega^\varepsilon} \right) \varepsilon \nabla \Psi_\varepsilon \cdot u_k^\varepsilon \nabla u_k^\varepsilon \, \mathrm{d}
  \bx.
\end{equation}
 Let us first consider the integral over $\omega^\varepsilon$. It follows from \cite{bgt} that the
$\RL^\infty$ norms of the Steklov eigenfunctions $u_k^\varepsilon$ is bounded, uniformly in terms of
$\sigma_k^\varepsilon$, and the norm of the trace $T : \RB\RV(\Omega^\varepsilon) \to
\RL^1(\partial \Omega^\varepsilon)$, which we have shown in Lemma \ref{lem:tracebv} to be
bounded. This, along with the
scaled $\RH^1$ norm estimate for $\varepsilon \nabla \Psi_\varepsilon$ from Lemma
\ref{lemma:sobolevpsieps}, the $\RL^2$ boundedness of $\nabla u_k^\varepsilon$ from
Lemma \ref{lemma:boundedUn} and the generalised Hölder inequality yields

\begin{equation}
   \label{eq:gpsigu2}
   \begin{aligned}
    \int_{\omega^\varepsilon}2 \varepsilon \nabla \Psi_\varepsilon \cdot u_k^\varepsilon \nabla
    u_k^\varepsilon \, \mathrm{d} \bx &\le 2 \norm{\varepsilon \nabla \Psi_\varepsilon}_{\RL^2(\Omega^\varepsilon)^d}
    \norm{\nabla u_k^\varepsilon}_{\RL^2(\Omega^\varepsilon)^d}
    \norm{u_k^\varepsilon}_{\RL^\infty(\omega^\varepsilon)} \\
    &\ll \varepsilon^{\frac 1 2 + \frac 1 d}.
    \end{aligned}
\end{equation}

For the integral over $\tilde \Omega^\varepsilon \setminus \omega^\varepsilon$, observe that scaling Lemma \ref{lemma:sobolevpsieps} yields
\begin{equation}
\begin{aligned}
 \norm{\varepsilon\nabla \Psi_\varepsilon}_{\RL^\infty(\tilde \Omega^\varepsilon \setminus \omega^\varepsilon)^d} &= \sup_{\abs \alpha = 1} \norm{D^\alpha \psi_\varepsilon}_{\RL^\infty(\mathcal{C}^\varepsilon \setminus B(0,1/4))} \\
 &\ll \varepsilon^{\frac 1 2 + \frac 1 d}.
 \end{aligned}
\end{equation}

Inserting that bound into the generalised Hölder inequality, along with $\RH^1$ boundedness of $u_k^\varepsilon$ yields
\begin{equation}
   \begin{aligned}
    \int_{\omega^\varepsilon}2 \varepsilon \nabla \Psi_\varepsilon \cdot u_k^\varepsilon \nabla
    u_k^\varepsilon \, \mathrm{d} \bx &\le 2 \norm{\varepsilon \nabla \Psi_\varepsilon}_{\RL^\infty(\tilde \Omega^\varepsilon \setminus \omega^\varepsilon)^d}
    \norm{\nabla u_k^\varepsilon}_{\RL^2(\Omega^\varepsilon)^d}
    \norm{u_k^\varepsilon}_{\RL^2(\Omega^\varepsilon)} \\
    &\ll \varepsilon^{\frac 1 2 + \frac 1 d}.
    \end{aligned}
\end{equation}

Finally, for the integral on the boundary $\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon$, we have
from Hölder that
\begin{equation}
  \label{eq:delnupsiuk2}
\varepsilon \int_{\partial\tilde \Omega^\varepsilon \setminus \partial T^\varepsilon} (u_k^\varepsilon)^2 \partial_\nu
\Psi_\varepsilon \, \mathrm{d} A \le \varepsilon\abs{\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon}\norm{\partial_\nu \Psi_\varepsilon}_{\RL^\infty(\partial \tilde \Omega^\varepsilon \setminus \partial
T^\varepsilon)}\norm{(u_k^\varepsilon)^2}_{\RL^\infty(\partial \tilde \Omega^\varepsilon \setminus \partial T^\varepsilon)}
\end{equation}
We have as in the proof of Lemma \ref{prop:limitproblem} that the $\RL^2$ norm of $u_k^\varepsilon$ is uniformly bounded, from uniform boundedness of the trace operator. From Lemma \ref{lemma:sobolevpsieps} we have that $$\varepsilon \norm{\partial_\nu
\Psi_\varepsilon}_{\RL^\infty(\partial\tilde \Omega^\varepsilon \setminus \partial T^\varepsilon)} =
\bigo{\varepsilon^{\frac 1 2 + \frac 1 d}},$$ and as earlier, we have
$$ \norm{(u_k^\varepsilon)^2}_{\RL^\infty(\partial \tilde \Omega^\varepsilon \setminus \partial
T^\varepsilon)} \le C,$$
from Lemma \ref{lem:tracebv}. Finally, it follows from standard lattice packing
theory that $\# \tilde I^\varepsilon \le C \varepsilon^{1-d}$
for some $C$ depending only on $\Omega$, so that $\abs{\partial \tilde \Omega^\varepsilon
\setminus \partial T^\varepsilon} \le C$. Combining these three estimates, we have indeed
that the product in
\eqref{eq:delnupsiuk2} is going to $0$ as $\varepsilon \to 0$, concluding the proof.
\end{proof} 


Until now, we have shown that the harmonic extensions to the holes in
$\Omega^\varepsilon$
of Steklov eigenpairs $(\sigma_k^{\varepsilon},u_k^{\varepsilon})$ converge weakly in $\RH^1$
and strongly in $\RL^2$ to a solution $(\Sigma,U)$ to the
problem
\begin{equation}
  \begin{cases}
   - \Delta U = A_d \beta \Sigma U & \text{in } \Omega, \\
    \partial_\nu U = \Sigma U & \text{on } \partial \Omega.
  \end{cases}
\end{equation}
It remains to be shown that the convergence is to the ``right'' eigenpair
$(\Sigma_{k,\beta},U_k)$. Spectral convergence of this type is a staple of
  homogenisation theory, see for example \cite{ACG,NPP} in a standard setting or
\cite{FL} using the theory of $E$-convergence. In both cases,
the general theory cannot be directly applied here since the Hilbert space
$\RL^2_{A_d\beta}(\Omega) \times \RL^2(\partial \Omega)$ on which the
limit problem is self-adjoint has no natural embedding to the Hilbert spaces
$\RL^2(\partial \Omega^\varepsilon)$, compare e.g. with \cite[Definitions 1--3]{FL}. Our methods use instead the quadratic forms
associated with the eigenproblems directly, similar methods were used in
spectral prescription, see e.g. \cite{CdV87}.

For the remainder of this section, as $\beta$ is
fixed, we will write simply $\Sigma_k$.
We first start by the following lemma, showing that the limit function $U$ does
not degenerate to the $0$ function.
\begin{lemma}
  Let $U$ be such that $U_k^{\varepsilon} \to U$ weakly in $\RH^1(\Omega)$. Then, 
  \begin{equation}
    (U,U)_\beta = 1.
  \end{equation}
\end{lemma}
\begin{proof}
 By compactness of the
  trace operator on $\partial \Omega$, we have that
  \begin{equation}
    \int_{\partial \Omega} (U_k^{\varepsilon})^2 \, \mathrm{d} \bx \to \int_{\partial \Omega} U^2 \, \mathrm{d}
    \bx.
  \end{equation}
  as $\varepsilon \to 0$. Moreover, by Lemma \ref{lem:normconv}, we have that
  \begin{equation}
    \int_{\partial T^\varepsilon} \left( U_k^\varepsilon \right)^2 \, \mathrm{d} \bx \to A_d\beta\int_\Omega U^2\, \mathrm{d}\bx 
  \end{equation}
  as $\varepsilon \to 0$. Hence $(U_k^\varepsilon,U_k^\varepsilon)_{\partial^\varepsilon} \to (U,U)_\beta$.
  Since $U_k^\varepsilon$ has been normalised to $\RL^2(\partial \Omega^\varepsilon)$ norm $1$, this concludes the
  proof.
\end{proof}
We are now ready to complete the proof of our first main result.
\begin{proof}[Proof of Theorem \ref{thm:homogenisation}]
We first show that all the eigenvalues converge. We proceed by induction on the eigenvalue rank $k$. The case $k = 0$ is
trivial. Indeed, we then have that the eigenvalues $\sigma_0^{\varepsilon} \equiv 0$
obviously converge to $\Sigma_0 = 0$ and the normalised eigenfunctions $U_0^{\varepsilon}(\bx) = \abs{\partial
\Omega^\varepsilon}^{-1/2}$, which converges to the constant fonction 
$$U_0(\bx) =
\left(\abs{\partial \Omega} + A_d\beta\abs{\Omega}\right)^{-1/2}.$$
Suppose now that for all $0 \le j \le k-1$, we have that $U_j^{\varepsilon}$ converges to $U_j$ weakly in $\RH^1(\Omega)$.
We first show that 
\begin{equation}
  \label{eq:asympboundsigstar}
  \Sigma_k \ge \sigma_k^{\varepsilon} + \smallo 1.
\end{equation}
In order to do so, we will show that the eigenfunctions $U_k$ are good
approximations to appropriate test functions for the variational
characterisation \eqref{eq:varcharsigm} of $\sigma_k^\varepsilon$.

Observe that by compactness of the trace operator on $\partial \Omega$ and by Lemma
\ref{prop:limitproblem}
\begin{equation}
  \lim_{\varepsilon \to 0} \int_{\partial \Omega^\varepsilon} u_j^{\varepsilon} U_k \, \mathrm{d} A = \int_{\partial
  \Omega} U_k U_j\, \mathrm{d} A + A_d\beta \int_\Omega U_k U_j \, \mathrm{d} \bx= 0
\end{equation}
for all $0 \le j \le k-1$. Hence, we can write
\begin{equation}
  \label{eq:ukalmostortho}
  U_k = V^{\varepsilon} + \sum_{j=0}^{k-1} \eta_j^{\varepsilon} u_j^{\varepsilon} \, \mathrm{d} A
\end{equation}
where for all $0 \le j \le k-1$ and all $\varepsilon > 0$, 
\begin{equation}
 \int_{\partial \Omega^\varepsilon} V^\varepsilon u_j^\varepsilon = 0
\end{equation}
and $\eta_j^{\varepsilon} \to 0$ as $\varepsilon \to 0$.
 Now, we have that
\begin{equation}
  \begin{aligned}
    \Sigma_k &= \int_\Omega \abs{\nabla U_k}^2 \, \mathrm{d} \bx \\
    &\ge \int_{\Omega^\varepsilon} \abs{\nabla V^{\varepsilon}}^2  + \sum_{j=0}^{k-1}(\eta_j^{\varepsilon}) \nabla
    u_j^{\varepsilon} \cdot \nabla V^\varepsilon + \sum_{j,l = 0}^{k-1} \eta_j^\varepsilon
    \eta_k^\varepsilon \nabla u_j^\varepsilon \cdot \nabla u_k^\varepsilon \, \mathrm{d} \bx. \\
\end{aligned}
\end{equation}
Since the $u_j^\varepsilon$ and $V^\varepsilon$ are bounded in $\RH^1(\Omega)$, the integral
of the two sums in the previous equation go to $0$ as $\varepsilon \to 0$. On the other
hand, we have that for all $\varepsilon > 0$, $V^\varepsilon$ is an appropriate test function
for $\sigma_k^\varepsilon$. Hence,
\begin{equation}
  \int_{\Omega^\varepsilon} \abs{\nabla V^\varepsilon}^2 \, \mathrm{d} \bx \ge \sigma_k^\varepsilon \int_{\partial
  \Omega^\varepsilon} (V^\varepsilon)^2 \, \mathrm{d} A.
\end{equation}
It follows from the decomposition \eqref{eq:ukalmostortho} and
the fact that by Lemma \ref{prop:limitproblem},
\begin{equation}
  \int_{\partial \Omega^\varepsilon} U_k^2 \, \mathrm{d} A \to \int_{\partial \Omega} U_k^2 \, \mathrm{d} A + A_d\beta
  \int_\Omega U_k^2 \, \mathrm{d} \bx = 1
\end{equation}
that
\begin{equation}
  \lim_{\varepsilon\to 0} \int_{\partial \Omega^\varepsilon} (V^\varepsilon)^2 \, \mathrm{d} \bx  = 1,
\end{equation}
implying that indeed $\Sigma_k \ge \sigma_k^{\varepsilon} + \smallo 1$. We now show that 
$$
\Sigma_k \le \sigma_k^{\varepsilon}.
$$
Let $(\Sigma,U)$ be the limit eigenpair for
$(\sigma_k^{\varepsilon},u_k^{\varepsilon})$ and suppose that 
$$(\Sigma,U) =
(\Sigma_j,U_j)$$
for some $0 \le j \le k-1$. We have
that  
\begin{equation}
  \begin{aligned}
  0 &=   \lim_{\varepsilon \to 0} \int_{\partial \Omega^\varepsilon} u_k^{\varepsilon} \, \mathrm{d} A \\&= 
  \lim_{\varepsilon \to 0}\int_{\partial \Omega^\varepsilon}
  u_k^\varepsilon U_j \, \mathrm{d} A + \int_{\partial \Omega^\varepsilon} u_k^\varepsilon( u_j^\varepsilon - U_j) \, \mathrm{d} A
\end{aligned}
\end{equation}
The first term converges to $1$ by the assumption that $(U,U_j)_\beta = 1$ and Lemma
\ref{prop:limitproblem}. As for the second term, we have by the Cauchy--Schwarz
inequality and the normalisation $\norm{u_k^\varepsilon}_{\RL^2(\partial \Omega^\varepsilon)}=1$ that
\begin{equation}
  \begin{aligned}
  \int_{\partial \Omega^\varepsilon} u_k^\varepsilon( u_j^\varepsilon - U_j) \, \mathrm{d} A &\le \norm{u_j^\varepsilon -
  U_j}_{\RL^2(\partial \Omega^\varepsilon)} \\
  &\to 0
\end{aligned}
\end{equation}
as $\varepsilon \to 0$ by Lemma \ref{lem:normconv}. This results in a contradiction. We therefore deduce that $\sigma_k^{\varepsilon}$
converges to some eigenvalue $\Sigma$ of problem \ref{problem:homo} that
is larger than $\Sigma_j$ for all $j \le k-1$. Combining this with the
bound \eqref{eq:asympboundsigstar} implies that $\sigma_k^{\varepsilon}$ converges to
$\Sigma_k$, and convergence of the eigenfunction therefore follows, when the
eigenvalues are simple. 

Otherwise, if $\Sigma_k$ has multiplicity $m$, \emph{i.e.} 
\begin{equation}
  \Sigma_{k-1} < \Sigma_k = \dotso = \Sigma_{k+m-1} < \Sigma_{k+m},
\end{equation}
observe that the above argument still yields convergence of $\sigma_j^\varepsilon$ to
$\Sigma_j$ for all $k \le j < k+m$. For the eigenfunctions, start by fixing a
basis $U_k,\dotsc,U_{k+m - 1}$ of the eigenspace associated with $\Sigma_k$. Observe that along
any subsequence there is a further subsequence such that all the eigenfunctions
$U_j^\varepsilon$, converge simultaneously to solutions of Problem \ref{problem:homo}.
Since for all $\varepsilon$, the functions $U_j^\varepsilon$ were $\RL^2(\partial
\Omega^\varepsilon)$-orthogonal, in the limit they are still orthogonal, this time with respect
to $(\cdot,\cdot)_\beta$. This implies that in the limit they span the
eigenspace associated with $\Sigma_k$. As such, the projection on the span of
$\set{U_j^\varepsilon : k \le j < k+m}$ converges to the projections on the span of
$\set{U_j : k \le j < k+m}$. Since this was
true along any subsequence, it is also true for the whole sequence, proving
convergence of the projections in the sense alluded to in Remark \ref{rem:multiple}.
\end{proof}
\section{\bf Dynamical boundary conditions with large parameter}\label{Section:LargeBeta}
The goal of this final section is to understand the limit as $\beta$ becomes large of the eigenvalues $\Sigma_{k,\beta}$ and of the corresponding eigenfunctions $U_{k,\beta}$, normalized by
\begin{equation*}
 1 = (U_{k,\beta},U_{k,\beta})_\beta = \int_{\partial \Omega} U_{k,\beta}^2 \, \mathrm{d} A + A_d \beta \int_{\Omega} U_{k,\beta}^2 \, \mathrm{d} \bx.
\end{equation*}
Recall that the Neumann eigenvalues of $\Omega$ are
$$0=\mu_0\leq\mu_1\leq\mu_2\leq\cdots\nearrow\infty.$$
We are now ready to prove our second main result.

\begin{proof}[Proof of Theorem \ref{thm:limitbeta}]
For $k$ fixed, we start by showing that $\beta \Sigma_{k,\beta}$ is bounded. Consider the min--max characterisation of $\Sigma_{k,\beta}$ to obtain
\begin{equation}
\label{eq:minmaxbSig}
 \tilde \Sigma_{k,\beta}:=\beta \Sigma_{k,\beta} = \min_{\substack{E \subset \RH^1(\Omega) \\ \operatorname{dim}(E) = k+1}} \max_{f \in E \setminus \set 0} \frac{\int_{\Omega} \abs{\nabla f}^2 \, \mathrm{d} \bx}{\frac 1 \beta \int_{\partial \Omega} f^2 \, \mathrm{d} A + A_d \int_\Omega f^2 \, \mathrm{d} \bx}.
\end{equation}
The quotient on the righthand side of \eqref{eq:minmaxbSig} is clearly bounded uniformly in $\beta$ for any $k+1$ dimensional subspace of smooth functions on $\Omega$. 
We can therefore suppose that a subsequence in $\beta$ of $\tilde \Sigma_{k,\beta}$ converges, say to $\tilde \Sigma_{k,\infty}$. Let us now prove that $\tilde U_{k,\beta} := \beta^{1/2} U_{k,\beta}$ is a bounded family in $\RH^1(\Omega)$. The normalisation on $U_{k,\beta}$ implies that
\begin{equation*}
 1 \ge A_d \beta\int_\Omega U_{k,\beta}^2 \, \mathrm{d} \bx 
  = A_d \int_\Omega \tilde U_{k,\beta}^2 \, \mathrm{d} \bx.
\end{equation*}
For the Dirichlet energy, we have that
\begin{equation}
 \int_{\Omega} \abs{\nabla\tilde U_{k,\beta}}^2 \, \mathrm{d} \bx = \beta \Sigma_{k,\beta},
\end{equation}
which was already shown to be bounded. Therefore there is also a weakly convergent subsequence in in $\RH^1(\Omega)$ as $\beta \to \infty$, converging to say $\tilde U_{k,\infty}$. We take the subsequence to coincide with the one for $\tilde \Sigma_{k,\beta}$. Observe as well that the normalisation condition on $U_{k,\beta}$ prevents the limit $\tilde U_{k,\infty}$ from vanishing identically.

The functions $\tilde U_{k,\beta}$ satisfy the following weak variational characterisation for any element $V \in \RH^1(\Omega)$
\begin{equation}
\int_{\Omega} \nabla \tilde U_{k,\beta} \cdot \nabla V\,\, \mathrm{d}\bx = \tilde\Sigma_{k,\beta}\left(\beta^{-1}\int_{\partial \Omega} \tilde U_{k,\beta} V \, \mathrm{d} A + A_d \int_\Omega \tilde U_{k,\beta} V \, \mathrm{d} \bx  \right).
\end{equation}
Letting $\beta \to \infty$, weak convergence of $\tilde U_{k,\beta}$ in $\RH^1(\Omega)$ implies that the limit satisfies the weak identity
\begin{equation}
  \int_\Omega \nabla \tilde U_{k,\infty} \cdot\nabla V\,\, \mathrm{d}\bx = \tilde \Sigma_{k,\infty} A_d \int_\Omega \tilde U_{k,\infty} V \, \mathrm{d} \bx. 
\end{equation}
In other words, $\tilde U_{k,\infty}$ is a solution to the Neumann eigenvalue problem with eigenvalue $\mu = \tilde \Sigma_{k,\infty} A_d$.

We now proceed by recursion on $k$ to show convergence to the right eigenpair. Once again, the statement is trivial for $k = 0$ and the constant eigenfunction. Assume that we have convergence for the first $k-1$ eigenpairs.
We now proceed in a similar fashion as in the proof of the spectral convergence to Problem \eqref{problem:homo}. We repeat the argument because the inequalities are more subtle.
 We first show that
 \begin{equation}
 \label{eq:neumannupper}
  \tilde \mu_k(\Omega):= \frac{\mu_k(\Omega)}{A_d} \ge \tilde \Sigma_{k,\beta}(1 + \smallo 1).
 \end{equation}

 Write
 \begin{equation}
  f_k = F_\beta + \sum_{j=0}^{k-1} (f_k, U_{j,\beta})_\beta U_{j,\beta},
 \end{equation}

with $F_\beta \perp_\beta U_{j,\beta}$ for $0 \le j < k$. We have that
\begin{equation}
\label{eq:scalardecomp}
 (f_k,U_{j,\beta})_\beta = (f_k,\beta^{-1/2}f_j)_\beta + (f_k,U_{j,\beta}-\beta^{-1/2}f_j)_\beta .
\end{equation}
The first inner product develops as
\begin{equation}
 (f_k,\beta^{-1/2}f_j)_\beta = \beta^{-1/2} \int_{\partial\Omega} f_k f_j \, \mathrm{d} A +A_d\beta^{1/2}\int_{\Omega} f_k f_j \, \mathrm{d} \bx. 
\end{equation}
The first term clearly goes to $0$ as $\beta \to \infty$, and the second one vanishes by orthogonality of the Neumann eigenfunctions in $\RL^2(\Omega)$. We now turn our attention to the second inner product in \eqref{eq:scalardecomp}. We have that
\begin{equation}
 \lim_{\beta \to \infty} \beta^{-1/2} \int_{\partial \Omega}f_k (\tilde U_{j,\beta} - f_j) \, \mathrm{d} A = 0
\end{equation}
by weak $\RH^1$ convergence of $\tilde U_{j,\beta}$ and compactness of the trace operator. On the other hand, strong $\RL^2$ convergence implies that
\begin{equation}
\label{eq:tightcontrol}
 A_d \beta^{1/2} \int_{\Omega} f_k (\tilde U_{j,\beta} - f_j) \, \mathrm{d} \bx = \smallo{\beta^{1/2}}.
\end{equation}
All in all, this implies that
\begin{equation}
 (f_k,U_{j,\beta})_\beta = \smallo{\beta^{1/2}}
\end{equation}
for all $0 \le j < k$. 
We now write
\begin{equation}
\label{eq:firstapproxneumann}
 \begin{aligned}
  \tilde \mu_k &= \frac{1}{A_d} \int_\Omega \abs{\nabla f_k}^2 \, \mathrm{d} \bx \\
  &\ge \frac{1}{A_d} \int_{\Omega} \abs{\nabla F_\beta}^2 \, \mathrm{d} \bx - \frac{1}{A_d}\sum_{j=0}^{k-1} (f_k,U_{j,\beta})_\beta^2 \Sigma_{j,\beta}.
 \end{aligned}
\end{equation}
Since $\beta \Sigma_{j,\beta}$ is bounded and equation \eqref{eq:tightcontrol} implies that $(f_k,U_{j,\beta})_\beta^2 = \smallo{\beta}$, we deduce that the last term in \eqref{eq:firstapproxneumann} goes to $0$. We now observe that by the variational characterisation of $\Sigma_{k,\beta}$,
\begin{equation}
\begin{aligned}
  \frac{1}{A_d}\int_{\Omega} \abs{\nabla F_\beta}^2 \, \mathrm{d} \bx &\ge \frac{\Sigma_{k,\beta}}{A_d}\left(\int_{\partial \Omega} F_\beta^2 + A_d \beta \int_\Omega F_\beta^2\right), \\
 &= \tilde \Sigma_{k,\beta} \left(\frac{1}{A_d\beta}\int_{\partial\Omega} F_\beta^2 \, \mathrm{d} A  + \int_\Omega F_\beta^2 \, \mathrm{d} \bx\right).
 \end{aligned}
\end{equation}
In the same way as we obtained the bound on the last term in \eqref{eq:firstapproxneumann}, the first integral is $\smallo \beta$. As for the second one, we write
\begin{equation*}
 \begin{aligned}
  \int_\Omega F_\beta^2  \, \mathrm{d} \bx &= \int_{\Omega} \left(f_k - \sum_{j=0}^{k-1} (f_k, U_{j,\beta})_\beta U_{j,\beta}^2\right)^2 \, \mathrm{d} \bx \\
  &= 1 - \sum_{j=0}^{k-1} \int_\Omega (f_k,U_{j,\beta})_\beta f_j U_{j,\beta}\, \mathrm{d} \bx + \sum_{j,\ell = 0}^{k-1} \int_\Omega (f_k,U_{j,\beta})_\beta(f_k,U_{\ell,\beta})_\beta U_{j,\beta}U_{\ell,\beta} \, \mathrm{d} \bx.
 \end{aligned}
\end{equation*}
Strong $\RL^2(\Omega)$ convergence of $\beta^{1/2} U_{j,\beta}$ and the fact that $(f_k,U_{j,\beta}) = \smallo{\beta^{1/2}}$ imply that the last two integrals converge to $0$ as $\beta \to \infty$. 

We have therefore obtained that
\begin{equation}
 \tilde \mu_k(\Omega) \ge \tilde \Sigma_{k,\beta}\left(1 + \smallo 1\right),
\end{equation}
\emph{i.e.} we have indeed proven assertion \eqref{eq:neumannupper}.

Suppose now that $(\tilde \Sigma_{k,\beta}, \tilde U_{k,\beta})$ converge to a Neumann eigenpair $(\tilde \mu_j, f_j)$ for some $j < k$. Then, we have that 
\begin{equation}
 \begin{aligned}
  1 &= \lim_{\beta \to \infty} \int_\Omega f_j \tilde U_{k,\beta} \, \mathrm{d} \bx \\
  &= \lim_{\beta \to \infty} \int_{\Omega} (f_j - \tilde U_{j,\beta})\tilde U_{k,\beta} \, \mathrm{d} x + \int_\Omega \tilde U_{k,\beta} U_{j,\beta} \, \mathrm{d} \bx \\
  &\le \lim_{\beta \to \infty}\norm{f_j - \tilde U_{j,\beta}}_{\RL^2(\Omega}\norm{\tilde U_{k,\beta}}_{\RL^2(\Omega)} + \abs{\int_\Omega U_{k,\beta} U_{j,\beta} \, \mathrm{d} \bx} \\
  &= 0,
 \end{aligned}
\end{equation}
where the limit comes from strong convergence in $\RL^2(\partial \Omega)$ to $0$ of $U_{j,\beta}$ and our recursion hypothesis. This is a contradiction, hence the convergence is to the correct eigenpair. This implies convergence of the whole sequence if the Neumann spectrum of $\Omega$ is simple. If there are eigenvalues with multiplicity, the same procedure as for the homogenisation problem yields once again convergence.

As for continuity in $\beta$, the same proof goes through in exactly the same way, except for the fact that we do not need to show the boundedness results in $\beta$.
\end{proof}

We can now prove the comparison results between Steklov and Neumann eigenvalues.


\begin{proof}[Proof of Corollary \ref{cor:boundd}]

It is proved in~\cite[Theorem 1.4]{ceg2} that any bounded domain $\Omega\subset\mathbb R^d$ satisfies
\begin{gather}
  \sigma_k(\Omega)|\partial\Omega|\leq C(d)|\Omega|^{\frac{d-2}{d}}k^{2/d},
\end{gather}
where $C(d)$ is a constant which depends only on the dimension.
When applied to $\Omega^\varepsilon$ this leads, after taking $\varepsilon\to 0$, to
\begin{gather}
  \Sigma_{k,\beta}(|\partial\Omega|+A_d\beta|\Omega|)\leq C(d)|\Omega|^{\frac{d-2}{d}}k^{2/d}.
\end{gather}
Taking the limit $\beta\to\infty$ leads to
\begin{gather}
  \mu_k|\Omega|\leq C(d)|\Omega|^{\frac{d-2}{d}}k^{2/d}.
\end{gather}
This is equivalent to
\begin{gather}
  \mu_k|\Omega|^{2/d}\leq C(d)k^{2/d}.
\end{gather}
\end{proof}

Finally, all is left to do is to prove 
the previous theorem has the following corollary in dimension $d = 2$. It allows one to transform universal bounds for Steklov eigenvalues into universal bounds for Neumann eigenvalues.

We write $\Omega^\varepsilon_\beta$ for a domain $\Omega^\varepsilon$ as constructed earlier whose holes are exactly of radius $r_\varepsilon^{d-1} = \beta \varepsilon^{d}$.


\begin{proof}[Proof of Theorem \ref{thm:neumannd=2}]
 We have from Theorem \ref{thm:homogenisation} that
 \begin{equation}
 \begin{aligned}
  \lim_{\varepsilon \to 0} \sigma_k\left(\Omega_\beta^\varepsilon\right)\abs{\partial \Omega_\beta^\varepsilon} &= \Sigma_{k,\beta}(\Omega)\left(\abs{\partial \Omega} + A_d \beta \abs \Omega\right) \\
  &= \frac{\abs{\partial \Omega}} \beta \tilde \Sigma_{k,\beta} + A_d \abs \Omega \tilde \Sigma_{k,\beta}. 
  \end{aligned}
 \end{equation}
Now, the first term clearly goes to $0$ as $\beta \to \infty$ while, by Theorem \ref{thm:limitbeta}, we have that
\begin{equation}
 \lim_{\beta \to \infty} A_d \abs \Omega \tilde \Sigma_{k,\beta} = \mu_k(\Omega) \abs \Omega.
\end{equation}
\end{proof}









\subsection*{Acknowledgements}
The authors would like to thank Iosif Polterovich for several very useful
conversations. We would also like to thank Valeri Smyshlyaev and Ilya Kamotski
for pointing us towards many classical results in homogenisation theory. We
would also like to thank Almut Burchard for possible further related questions.
The authors also would like to thank Gérard Philippin, who was helpful in the
early stage of this project. 
We thank both referees for a careful reading and
many comments which greatly helped in improving the paper. In particular, we
thank one of them for comments that led to Lemma \ref{lem:tracebv}, and the
other for pointing out reference \cite{HassannezhadMiclo}.
A.G. acknowledges support from NSERC. A.H. was partially supported by the ANR
project ANR-18-CE40-0013 -- SHAPO on Shape Optimisation. He wants
also to thank the Centre de Recherche Math\'ematiques de Montr\'eal and the
Département de Mathématiques of Université Laval for grants and fruitful
atmosphere during his stay there in 2018. 
J.L. was supported by the NSERC postdoctoral fellowship, and EPSRC grant number EP/P024793/1.

\bibliographystyle{plain}
