Finite element approximation for the fractional eigenvalue problem
The purpose of this work is to study a finite element method for finding solutions to the eigenvalue problem for the fractional Laplacian. We prove that the discrete eigenvalue problem converges to the continuous one and we show the order of such convergence. Finally, we perform some numerical experiments and compare our results with previous work by other authors.
\section{Introduction and Main Results}

	Anomalous diffusion phenomena are ubiquitous in nature \cite{Klafter, MetzlerKlafter}, and
	the study of nonlocal operators has been an active area of research in different branches of 
	mathematics. Such operators arise in applications as image processing~\cite{BuadesColl, GattoHesthaven,  GilboaOsher, 
	YifeiZhang}, finance~\cite{CarrHelyette, RamaTankov}, electromagnetic fluids~\cite{McCayNarasimhan}, peridynamics 
	\cite{Silling}, porous media flow~\cite{BensonWheatcraft, CushmanGinn}, among others. 
	
	\medskip
	
	One striking example of a nonlocal operator is the fractional Laplacian $(-\Delta)^s$, defined by
	\[
		(-\Delta)^su(x) \coloneqq  2 C(n,s)
		\int_{\mathbb{R}^n}\dfrac{u(x)-u(y)}{|x-y|^{n+2s}}\,dy,
		\quad x\in\mathbb{R}^n.
	\]
	Here the integral is understood in the principal value sense and the normalization constant $C(n,s)$ is given by
	\[
		C(n,s) \coloneqq \frac{2^{2s-1} s 
		\Gamma(s+\frac{n}{2})}{\pi^{n/2} 
		\Gamma(1-s)} .
	\]
	In the theory of stochastic processes, this operator appears as the infinitesimal generator of a stable L\'evy process, see 
	for instance \cite{Bertoin, Valdinoci}. Moreover, the fractional Laplacian is also one of the simplest examples of a 
	pseudo-differential operator, because its symbol is just $P(\xi)=|\xi|^{2s}.$

	\medskip
	
	An interesting problem concerning the fractional Laplacian is to find its eigenvalues on bounded domains. 
	Namely, to find a positive number $\lambda$ (eigenvalue) and a function $u\not\equiv0$ (eigenfunction) such that
	\begin{equation}\label{eq:int.autovalores}
		\begin{cases}
			(-\Delta)^s u = \lambda  u &\mbox{ in }\Omega, \\
				u = 0 &\mbox{ in }\Omega^c=\mathbb{R}^n\setminus\Omega,
		\end{cases}
	\end{equation}
	where $\Omega$ is a bounded smooth domain in $\mathbb{R}^n$ and $s\in(0,1)$.
	Observe that, due to the fact that pointwise values of $(-\Delta)^s u$ depend on the value of $u$ over the 
	whole space, the ``boundary'' conditions need to be imposed on the complement of $\Omega$. 
	
	Even if $\Omega$ is an interval, it is very challenging to obtain closed analytical expressions for the eigenvalues and eigenfunctions of the fractional Laplacian. This motivates the utilization of discrete approximations of this problem (see, for example, \cite{DuoZhang, Kwasnicki, ZoiaRossoKardar}); in this work we consider a finite element method. In first place, we prove that the discrete eigenvalue problem converges to the continuous one. Then, we show the order of convergence for eigenvalues and eigenfunctions, both in the energy and the $L^2$-norm. Finally, we perform some numerical experiments and compare our results with previous work by other authors. These results are in good agreement with our theory. 
	
The finite element method is flexible enough to deal with non-convex domains, and enables us to provide estimates and sharp upper bounds for eigenvalues even in this context. Moreover, as a consequence of our numerical experiments in the $L$-shaped domain $\Omega = [-1,1]^2\setminus [0,1]^2$, we conjecture that the first eigenfunction for this domain is as regular as the first one in any smooth domain.

	
	\subsection*{Main Results} In order to state our results we need some notations and definitions. 
	The natural functional space for the eigenvalue problem \eqref{eq:int.autovalores} is
	\[
		\widetilde{H}^s (\Omega) \coloneqq \left\{ v \in H^s({\mathbb{R}^n}) \colon \text{ supp } v \subset \bar{\Omega} \right\},
	\]
	where $H^s({\mathbb {R}}^n)$ is the space of all functions $u\in L^2({\mathbb {R}}^n)$ such that
	\[
		|v|_{H^s({\mathbb {R}}^n)}^2 \coloneqq \iint_{{\mathbb {R}}^{2n}} 
				\frac{|v(x)-v(y)|^2}{|x-y|^{n+2s}} \, dx \, dy<\infty. 
	\]
	Moreover, $\left({\mathbb {V}},\|\cdot\|_{\mathbb {V}}\right)\coloneqq \left(\widetilde{H}^s (\Omega), \sqrt{C(n,s)} \, |\cdot|_{H^s({\mathbb {R}}^n)} \right)$
	is a Hilbert space with the inner product
	\[
			\left\langle u, v\right\rangle \coloneqq 
			C(n,s) \iint_{\mathbb{R}^{2n}} 
			\frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} \, dx \, dy.
	\]
	Obviously, the constant $\sqrt{C(n,s)}$ has no effect on the definition of the space; it is included in order to make the notation simpler in the rest of the paper.
	The fractional space $H^s({\mathbb {R}}^n)$ can also be defined for any $s>1.$
	If $s = m + \sigma$, where $m \in\mathbb{N}$ and $\sigma \in (0,1)$, 
	$H^s({\mathbb {R}}^n)$ is the space of all functions $v \in H^m({\mathbb {R}}^m)$ such that its weak derivatives of order $m$ belong to $H^{\sigma}({\mathbb {R}}^n)$. For more details, see Section \ref{sec.preliminares}.
	
	\medskip
	
	In this context, the eigenvalue problem \eqref{eq:int.autovalores} has the following variational formulation:
	find $\lambda\in (0,+\infty)$ and 
	$u\in {\mathbb {V}} $ such that $u\not\equiv 0$ and  
	\begin{equation} \label{eq:int.autovalores.de.intro}
		\left\langle u, v \right\rangle = 
		\lambda (u, v)
		\quad \mbox{ for all } v \in {\mathbb {V}}, 
	\end{equation}
	where  $( \cdot, \cdot )\colon  L^2(\Omega) \times L^2(\Omega) \to \mathbb{R}$ is the bilinear form 
	\[
		( u, v) \coloneqq \int_{\Omega} u(x)v(x) \, dx.  
	\]
		
	In \cite{MR3002745},
	the authors show that there is an infinite 
	sequence of eigenvalues $\{\lambda^{(k)}\}_{k\in\mathbb{N}}$, 
	\[
		0<\lambda^{(1)}< \lambda^{(2)}\leq \cdots \leq 
		\lambda^{(k)}\leq\cdots,\quad  \lambda^{(k)}\to\infty 
			 \mbox{ as }  k\to \infty,
	\] 
	where the same eigenvalue can be repeated several
	times according to its multiplicity. The corresponding eigenfunctions $u^{(k)}$ (normalized by $\|u^{(k)}\|_{L^2(\Omega)}=1$) form a complete orthonormal set in $L^2(\Omega).$ See also \cite{MR3089742,ServadeiValdinoci, MR3271254} and the
	references therein for further details on the fractional eigenvalue problem.
		
	\medskip
		
	Let us now introduce the discrete space. Let  ${\mathcal {T}}_h$ be a family of triangulations of $\Omega$ satisfying
		\begin{align}
		& \exists \sigma > 0 \mbox{ s.t. } h_T \leq \sigma \rho_T,  \tag{Regularity} 
			\label{eq:regularity.intro}
	\end{align}
for any element $T \in \mathcal{T}_h$, where $h_T$ is the diameter of $T$ and $\rho_T$ is the radius of the largest ball contained in $T$.
	This is the only requirement we need to impose to our family of triangulations.
	
	We consider continuous piecewise linear functions on $\mathcal{T}_h$, namely
	$$
		{\mathbb {V}}_h \coloneqq \left\{ v \in {\mathbb {V}} \colon v \big|_T \in \mathcal{P}_1(T) \ \forall T \in {\mathcal {T}}_h \right\} .
	$$
		
	Our Galerkin approximation consists in looking for discrete 
	eigenvalues $\lambda_{h}\in \mathbb{R}$ and $u_h\in{\mathbb {V}}_h$ such that
	$u_h\not\equiv0$ and
	\begin{equation} \label{eq:int.autovalores.di.intro}
		\left\langle u_h, v\right\rangle=\lambda_h(u_h,v)\quad\forall
		v\in{\mathbb {V}}_h.
	\end{equation}
	We can order  the discrete eigenvalues of \eqref{eq:int.autovalores.di.intro} as follows
	\[
		0<\lambda_h^{(1)}\leq \lambda_h^{(2)}\leq \cdots \leq 
		\lambda_h^{(k)}\leq \cdots \le\lambda_h^{(\text{dim} {\mathbb {V}}_h )},
	\]
	where the same eigenvalue is repeated according to its multiplicity. The corresponding eigenfunctions 
	$u^{(k)}_h$  (normalized by $\|u^{(k)}_h\|_{L^2(\Omega)}=1$) form an 
	orthonormal set in $L^2(\Omega).$ 
		
	\medskip
				
	The purpose of this work is to prove the convergence of this 
	discrete problem to the continuous one.  
	To this end, we provide a Sobolev regularity result. 
	Namely, using the regularity results of Grubb \cite{Grubb} and a bootstrap argument, we prove 
	
	\begin{proposition}\label{prop:regintro}
		Let $u$ be an eigenfunction of $(-\Delta)^s$ in $\Omega$ with homogeneous Dirichlet boundary conditions. Then, 
		$u \in \widetilde{H}^{s+\nicefrac12-\varepsilon}(\Omega)$ for any $\varepsilon > 0.$
	\end{proposition}
	
	Another important tool is given by considering approximation properties of
	$\Pi_h \colon {\mathbb {V}} \to {\mathbb {V}}_h$, the projection with respect to 
	the $\| \cdot \|_{{\mathbb {V}}}$ norm. To be specific,
	
	\begin{proposition}\label{estimacionlinf0}
		Let $u$ be an eigenfunction of \eqref{eq:int.autovalores.de.intro}. Then, 
		for any $\varepsilon>0$ there exists a positive constant $C$ independent of $h$ 
		such that
		\begin{equation}
			\label{estimacionlinf1}
				\|u-\Pi_h u\|_{{\mathbb {V}}} \le C h^{\nicefrac12-\varepsilon}|u|_{H^{s+\nicefrac12-\varepsilon}(\Omega)},
		\end{equation}
		and
		\begin{equation}\label{eq:aproxL2}
			\|u-\Pi_h u\|_{L^2(\Omega)} \le	C h^{\nicefrac12+\alpha-\varepsilon} |u|_{H^{s+\nicefrac12-\varepsilon}(\Omega)}.
		\end{equation}
		Here, $\alpha = s$ if $s < \nicefrac12$ and $\alpha = \nicefrac12 - \varepsilon$ if $s\ge 
		\nicefrac12.$
	\end{proposition}	
	

	Our first result concerning the convergence is to determine it in a ``gap" distance,
	that is, we prove that the discrete eigenvalue problem \eqref{eq:int.autovalores.di.intro}  converges to the continuous one 
	\eqref{eq:int.autovalores.de.intro}, following the definition of convergence given in 
	\cite{Kato,Boffia} (see also \cite{Boffi}). More precisely, we define 
	the gap between Hilbert spaces $E, F \subset H$ by
	\[
		\delta(E,F)=\sup_{u\in E,\|u\|_{H}=1} \inf_{v\in F} \|u-v\|_{H},
		\quad 
		\hat{\delta}(E,F)=\max(\delta(E,F),\delta(F,E)).
	\]
	Then, taking $H={\mathbb {V}},$ we say that the discrete eigenvalue problem 
	\eqref{eq:int.autovalores.di.intro} converges to the continuous one 
	\eqref{eq:int.autovalores.de.intro} if, for any $\varepsilon>0$ and $k>0$, there is 
	$h_0>0$ such that 
	\[
		\max_{1\leq i \leq m(k)} |\lambda^{(i)}-\lambda_h^{(i)}|
		\leq \varepsilon,\qquad
		\hat{\delta}\left( \bigoplus_{i=1}^{m(k)}  E^{(i)},
		 \bigoplus_{i=1}^{m(k)}  E_h^{(i)}\right)\leq \varepsilon,
	\]
	for all $h<h_0$. 
	Here, $m(k)$ is the dimension of the space spanned by the first distinct $k$ eigenspaces
	and $E^{(i)}$ and $E^{(i)}_h$ are the eigenspace and the discrete eigenspace associated to $\lambda^{(i)}$
	and $\lambda^{(i)}_h,$ respectively.
	
	\begin{theorem}\label{teo:gap}
		The discrete eigenvalue problem 
		\eqref{eq:int.autovalores.di.intro} converges to the continuous one 
		\eqref{eq:int.autovalores.de.intro}.
	\end{theorem}	
	
	Having established the convergence of the discrete problem to the continuous one, we aim to study the order of such convergence. Our procedure is in the same spirit as in \cite{Boffi}, and we prove the following results
	
	\begin{theorem}\label{teo:ordensimple}
		Let $(u^{(k)}, \lambda^{(k)})$ be an eigenpair of multiplicity one.
		For any $\varepsilon>0,$ there is a positive constant $C$ independent of $h$ such that
		\[
			\|u^{(k)}-u^{(k)}_h\|_{L^2(\Omega)}\leq C h^{\alpha+\nicefrac12-\varepsilon},
			\quad		
			0\leq \lambda_h^{(k)}-\lambda^{(k)}\leq C h^{1-\varepsilon},
		\]
		and
		\[
			\|u^{(k)}-u^{(k)}_h\|_{{\mathbb {V}}}\leq C h^{\nicefrac12-\varepsilon}.
		\]
		Here  $u_h^{(k)}$ is chosen in such a way that $\left(\Pi_h u^{(k)}, u_h^{(k)}\right) \ge 0$
		and $\alpha$ is given by Proposition \ref{estimacionlinf0}.
	\end{theorem}
	
	\begin{theorem}\label{teo:ordenmultiple}
		Let $\lambda^{(k)}$ be an eigenvalue of multiplicity $m$ (that is, $\lambda^{(k)}=\lambda^{(k+1)}=\cdots=\lambda^{(k+m-1)}$ 
		and $\lambda^{(i)} \neq \lambda ^{(k)}$ for $i\neq k,\dots,k+m-1 $).
		For any $\varepsilon>0,$ there exists a positive constant $C$ independent of $h$ 
		such that
		\[
			0\leq \lambda_h^{(j)}-\lambda^{(k)}\leq C h^{1-\varepsilon} \quad \forall k\leq j\leq k+m-1.
		\]
		Moreover, if $u^{(k)}$ is an eigenfunction associated to $\lambda^{(k)}$, there is
		\[
		\{ w_h^{(k)} \} \subset E_h^{(k)} \oplus \ldots \oplus E_h^{(k+m-1)}
		\]
		such that
		\[
			\|u^{(k)}-w^{(k)}_h\|_{L^2(\Omega)}\leq C h^{\alpha+\nicefrac12-\varepsilon}
		\]
	and	
		\[
			\|u^{(k)}-w^{(k)}_h\|_{{\mathbb {V}}}\leq C h^{\nicefrac12-\varepsilon},
		\]
		where $\alpha$ is given by Proposition \ref{estimacionlinf0}.
	\end{theorem}
	
	Propositions \ref{prop:regintro} and \ref{estimacionlinf0} and Theorem \ref{teo:gap} 	
	are the key ingredients in the proof of the above two theorems. 
	
	Let us mention at this point that in \cite{Boffi}, where the author studies the eigenvalue problem for the laplacian, the order of convergence of the eigenvalues is proved in first place, by using the min-max formulation and estimate \eqref{estimacionlinf1}. If we followed the sketch of this proof for the fractional laplacian case we would obtain order $\nicefrac12+s$ when $s<\nicefrac12$, which is less than $1$. 
	So, in order to obtain the optimal order we need to prove first the convergence of the eigenvalues and afterwards the order of convergence of the eigenfunctions. In this manner, using both results and  estimate \eqref{estimacionlinf1}, we obtain optimal order for the eigenvalues for any $0<s<1$.
	
	 \subsection*{The paper is organized as follows}  Section \ref{sec.preliminares} collects the notation we employ, and reviews some previous works on the problem under consideration. In particular, regularity of eigenfunctions is proved. The section concludes with a discussion of certain aspects of finite element approximations of the fractional Laplacian. Afterwards, Section \ref{sec.conv} deals with the convergence of the discrete eigenvalue problem to the continuous one in gap distance. In Sections \ref{sec.orden} and \ref{sec.orden_multiple} the study of the order of convergence is developed, first for simple eigenvalues and later for eigenspaces of dimension greater than one. Finally, numerical experiments are discussed in Section \ref{sec:numerico}.
	
\subsection*{Acknowledgments} The authors would like to thank G. Grubb and R. Rodr\'iguez for the  valuable help they provided through clarifying discussions on the topic of this paper, and to F. Bersetche for improving the efficiency of the Matlab code employed. 
		
		
\section{Preliminaries and Definitions}
\label{sec.preliminares}
	In this section we review the basic aspects of the problem under consideration. In first place, we set notation regarding 
	Sobolev spaces. Afterwards, we analyze theoretical properties of the eigenvalue problem \eqref{eq:int.autovalores}. Regularity 
	results for weak solutions of fractional Laplace equations are recalled next. The section concludes with the introduction of 
	the finite element spaces we work with.


\subsection{Sobolev spaces}	
	Given an open set $\Omega \subset {\mathbb{R}^n}$ and $s \in(0,1)$, define the fractional Sobolev space $H^s(\Omega)$ as
	\[
		H^s(\Omega) \coloneqq \left\{ v \in L^2(\Omega) \colon |v|_{H^s(\Omega)} < \infty \right\},
	\]
	where $|\cdot|_{H^s(\Omega)}$ is the seminorm
	\[
		|v|_{H^s(\Omega)}^2 \coloneqq \iint_{\Omega^2} 
		\frac{|v(x)-v(y)|^2}{|x-y|^{n+2s}} \, dx \, dy. 
	\]
	Obviously, $H^s(\Omega)$ is a Hilbert space endowed with the norm $\|\cdot\|_{H^s(\Omega)} = \|\cdot\|_{L^2(\Omega)} + |\cdot|_{H^s(\Omega)} .$
	If $s>1$ and it is not an integer, the decomposition $s = m + \sigma$, where $m \in \mathbb{N}$ and $\sigma \in (0,1)$, 
	allows to define $H^s(\Omega)$ by setting
	\[
		H^s(\Omega) \coloneqq \left\{ v \in H^m(\Omega) \colon |D^\alpha v|_{H^{\sigma}(\Omega)} < \infty \text{ for all } \alpha \text{ s.t } 
		|\alpha| = m \right\}.
	\]
 
	Let us also define the space of functions supported in $\Omega$, 
	\[
		\widetilde{H}^s (\Omega) \coloneqq \left\{ v \in H^s({\mathbb{R}^n}) \colon \text{ supp } v \subset \bar{\Omega} \right\}.
	\]	
	For $0\le s \le 1$, this space may be defined through interpolation,
	\[ 
		\widetilde{H}^s (\Omega) = \left[L^2(\Omega), H^1_0(\Omega) \right]_{s}.
	\]
	Moreover, depending on the value of $s$, different characterizations of this space are available (see, for example 
	\cite[Chapter 11]{LionsMagenes}). If $s<\nicefrac12$ the space $\widetilde{H}^s (\Omega)$ coincides with $H^s (\Omega)$, and if $s>\nicefrac12$ it may 
	be characterized as the closure of $C^\infty_0(\Omega)$ with respect to the $|\cdot|_{H^s (\Omega)}$ norm. In the latter case, it is also customary to denote it by $H^s_0(\Omega)$. The particular case of $s=\nicefrac12$ gives raise to the Lions-Magenes space $H^{\nicefrac12}_{00}(\Omega)$, 
	which can be characterized by
	\[
		H^{\nicefrac12}_{00} (\Omega) \coloneqq \left\{ v \in H^{\nicefrac12}(\Omega) \colon \int_\Omega \frac{v(x)^2}{\text{dist}(x,\partial \Omega)} \, dx < 
		\infty \right\}.
	\]
	Note that the inclusions $H^{\nicefrac12}_{00}(\Omega) \subset H^{\nicefrac12}_{0}(\Omega) \subset H^{\nicefrac12}(\Omega)$ are strict.


	It is apparent that $\left\langle \cdot, \cdot\right\rangle$
	defines an inner product on $\widetilde{H}^s(\Omega).$ In addition,
	the norm induced by it, which is just 
	the $H^s({\mathbb{R}^n})$ seminorm, is equivalent to the 
	full $H^s({\mathbb{R}^n})$ norm on this space, because of the following well known result.

	\begin{proposition}[Poincar\'e inequality]
		There is a constant $c=c(\Omega,n,s)$ such that 
		\begin{equation} \label{eq:poincare}
			\| v \|_{L^2(\Omega)} \leq c |v|_{H^s({\mathbb{R}^n})}  \quad \forall v \in \widetilde{H}^s(\Omega) .
		\end{equation}
	\end{proposition}


\subsection{Eigenvalue Problem}
	In the sequel, we will work within the Hilbert space
	$$
		({\mathbb {V}}, \| \cdot \|_{\mathbb {V}}) \coloneqq (\widetilde{H}^s(\Omega), \sqrt{C(n,s)} \, |\cdot|_{H^s({\mathbb{R}^n})}).
	$$ 
	
	In \cite{MR3002745}, the authors prove that for any $k\in\mathbb{N}$
	the eigenvalues of \eqref{eq:int.autovalores.de.intro} can be characterized as follows:
	\[
		\lambda^{(k)}=
		\min\left\{\dfrac{\|u\|_{{\mathbb {V}}}^2}{\|u\|_{L^2(\Omega)}^2}
		\colon u\in {\mathbb {V}}^{(k)}\setminus
		\{0\}
		\right\},
	\]
	where ${\mathbb {V}}^{(1)}={\mathbb {V}}$ and
	\[
		{\mathbb {V}}^{(k)}\coloneqq\left\{u\in{\mathbb {V}}\colon
		\left\langle u,u^{(j)}\right\rangle=0\quad\forall j=1,\dots,k-1\right\}
	\]
	for all $k\ge2.$ Therefore, by the min-max theorem,
	\[
		\lambda^{(k)} = \min_{E \in S^{(k)}} \max_{u \in E} 
		\frac{\| u \|^2_{\mathbb {V}}}{\| u \|^2_{L^2(\Omega)}}
	\]
	where $S^{(k)}$ denotes the set of all $k-$dimensional subspaces of 
	${\mathbb {V}}.$
	
	In \cite{MR3002745}, the authors show that 
	$\lambda^{(1)}$ is simple, see also
	\cite{2014arXiv1412.4722D,MR3148135}. 
	Whereas in \cite{MR3307955} it is proved that 
	$u^{(k)}\in L^{\infty}(\Omega)$ for all $k\in\mathbb{N}.$ 
	Therefore, by \cite[Proposition 1.1]{RosOtonSerra},
	$u^{(k)}\in C^s({\mathbb {R}}^n)$ for all $k\in\mathbb{N}.$


\subsection{Regularity results}
	Given a function $f \in H^r(\Omega)$ ($r \ge -s$), let us consider the homogeneous Dirichlet problem for the fractional
	Laplacian,		
		\begin{equation}\label{eq:fuente}
			\begin{cases}
				(-\Delta)^s u = f  &\mbox{ in }\Omega, \\
				u = 0 &\mbox{ in }\Omega^c .
			\end{cases}	
		\end{equation}
	Existence and uniqueness of a weak solution $u \in \widetilde{H}^s(\Omega)$ of the above equation is an immediate consequence of the Lax-Milgram lemma. 
	Following Grubb \cite{Grubb}, it is possible to obtain Sobolev regularity results for the solution. In that paper, the author deals with
	H\"ormander $\mu-$spaces $H^{\mu(s)}_p$; see that work for a definition and
	further details. The following result is a particular case of Theorem 7.1 therein:

	\begin{theorem}\label{teo:Grubb}
		Let $\ell>s-\nicefrac12$ and assume 
		$u \in \widetilde{H}^{\sigma}(\Omega)$ for some $\sigma>s-\nicefrac12$ and consider a right hand side function
		$f\in H^{\ell-2s} (\Omega)$. Then, it holds that $u \in H^{s(\ell)}(\Omega)$.
	\end{theorem}
 
	In particular, considering $\ell = r+2s$ in the previous theorem and taking into account that 
	\[ 
		H^{s(r+2s)}(\overline \Omega)
		\begin{cases}
				= \widetilde{H}^{2s+r}(\Omega)  &\mbox{ if } 0 < s+r < \nicefrac12, \\
				\subset  \widetilde{H}^{s+\nicefrac12-\varepsilon}(\Omega) \ \forall \varepsilon > 0, &\mbox{ if } 
				\nicefrac12 \le s+r < 1,
			\end{cases}	
	\]
	(see \cite[Theorem 5.4]{Grubb}), we obtain the following:

	\begin{proposition}\label{prop:regHr}
		Let $f\in H^r(\Omega)$ for $r\geq -s$ and $u\in \widetilde{H}^s(\Omega)$ be the solution of the Dirichlet problem
		\eqref{eq:fuente}. Then, the following regularity estimate holds 
		$$
			|u|_{H^{s+\alpha}({\mathbb{R}^n})} \leq C(n, \alpha) \|f\|_{H^r(\Omega)}.
		$$
		Here, $\alpha = s+r$ if $s+r < \nicefrac12$ and $\alpha = 
		\nicefrac12 - \varepsilon$ if $s+r \ge \nicefrac12$, with $\varepsilon > 0$ arbitrarily small.
	\end{proposition}


	\begin{remark}
		Assuming further Sobolev regularity in the right hand side function does not imply that the solution will 
		be any smoother than what is given by the previous proposition. Indeed, if $f \in H^r(\Omega)$, then 
		Theorem \ref{teo:Grubb} gives $u \in H^{s(r+2s)}(\Omega)$, which can not be embedded 
		in any space sharper than $H^{s+\nicefrac12-\varepsilon}(\Omega)$ if $r+s\ge \nicefrac12$.

		Moreover, regularity estimates for eigenfunctions of the fractional Laplacian are derived in 
		\cite{Grubb_autovalores, MR3294242}. 
		These estimates are formulated in terms of H\"older norms. Letting $d$ be a smooth function that behaves like
		$\text{dist}(x,\partial \Omega)$ near the boundary of $\Omega$, it is shown that any eigenfunction $u$ of 
		\eqref{eq:int.autovalores} lies in the space $d^sC^{2s(-\varepsilon)}(\Omega)$, where the $\varepsilon$ is active only if 
		$s=\nicefrac12$ and that $\nicefrac{u}{d^s}$ does not vanish near $\partial\Omega$. 
		This shows that no further regularity than $H^{s+\nicefrac12-\varepsilon}(\Omega)$ should be expected for eigenfunctions. 
	\end{remark}

	Sobolev regularity of eigenfunctions is a consequence of Proposition \ref{prop:regHr}.


	\begin{proof}[Proof of Proposition \ref{prop:regintro}]
		Since $f=\lambda u\in H^s(\Omega)$ the case $s\ge \nicefrac14$ is a direct consequence of Proposition \ref{prop:regHr}. When $s<\nicefrac14$ 
		the result follows using  a bootstrap argument.
	\end{proof}

	\begin{remark}
		As in Grubb's work the domain $\Omega$ is required to be smooth, we assumed that hypothesis in the two previous 
		propositions. However, numerical evidence indicates that this assumption is not necessary (see Subsection 
		\ref{sub:2d}). Observe that the theory of Section \ref{sec.conv} does 
		not necessitate the smoothness of $\Omega$, and that in sections \ref{sec.orden} and \ref{sec.orden_multiple} it is only required in order to apply Theorem \ref{teo:Grubb}.
	\end{remark}

\subsection{Finite element approximations}

	Let ${\mathcal {T}}_h$ be a family of shape-regular triangulations  of $\Omega$ (see introduction).
	Observe that for all $h>0$ the discrete space ${\mathbb {V}}_h$ is a subset of the continuous one.

	We have the analogous  min-max characterization for eigenvalues of the discrete problem, 
	\[
		\lambda_h^{(k)} = \min_{E \in S_h^{(k)}} \max_{u \in E} 
		\frac{\| u \|^2_{\mathbb {V}}}{\| u \|^2_{L^2(\Omega)}},
	\]
	where $S_h^{(k)}$ denotes the set of all
	$k$ dimensional subspaces of ${\mathbb {V}}_h.$

	\medskip

	\begin{remark}\label{re:compara}
		It follows from ${\mathbb {V}}_h \subset {\mathbb {V}}$ that
		\[
			\lambda^{(k)}\leq \lambda ^{(k)}_h .
		\] 
	\end{remark}

	In the remainder of this section, we analyze approximation properties of finite element spaces. Consider a 
	quasi-interpolation operator
	$I_h\colon H^l(\Omega) \to {\mathbb {V}}_h $ satisfying the following estimate: there exists $C>0$ such that 
	for any $w\in H^{l}(\Omega)$,
	\begin{equation}
		\label{estiminterpol}
		\|w - I_h w\|_{{\mathbb {V}}} \le C h^{l-s}  |w|_{H^{l}(\Omega)} .
	\end{equation}
	Quasi-interpolation operators were introduced in \cite{clement} (see also \cite{ScottZhang}), and an estimate like 
	\eqref{estiminterpol} is derived, for example, in \cite{AcostaBorthagaray}.

 
	Given $u \in {\mathbb {V}}$, recall that by $\Pi_h u$ we denote the projection with respect to the $\|\cdot\|_{\mathbb {V}}$ norm. Clearly, this is the only function in ${\mathbb {V}}_h$ such that the Galerkin orthogonality
	\[
			\left\langle u - \Pi_h u , \, v_h \right\rangle = 0 \quad\forall v_h \in {\mathbb {V}}_h
	\]
	holds, or equivalently, 
	\begin{equation}\label{eq:cea}
		\| u - \Pi_h u \|_{{\mathbb {V}}} = \inf_{v_h \in {\mathbb {V}}_h} \|u - v_h \|_{{\mathbb {V}}}.
	\end{equation}
	Observe that if $u$ is the solution of \eqref{eq:fuente}, 
	then $\Pi_h u$ corresponds to the solution of the corresponding
	discrete problem on ${\mathbb {V}}_h$.

	Now, we can prove our estimates for the error. 

	\begin{proof}[Proof of Proposition \ref{estimacionlinf0}]
		The first inequality follows by combining estimate \eqref{eq:cea}, interpolation inequality \eqref{estiminterpol} and 
		Proposition \ref{prop:regintro}.

		The second one follows from the following  Aubin--Nitsche argument. 	
		Let $w\in{\mathbb {V}}$ be the weak solution of the boundary value problem
			\[
			      \begin{cases}
				    (-\Delta)^sw=u-\Pi_h u &\mbox{ in }\Omega,\\
				    w=0&\mbox{ in }\Omega^c.\\
			      \end{cases}
			\]
		Then, resorting to Galerkin orthogonality again we obtain
			\[  		 
			      \|u-\Pi_h u\|_{L^2(\Omega)}^2 = 
			      \left\langle  w, u-\Pi_h u \right\rangle 
			      \leq \|w - I_h w\|_{{\mathbb {V}}} \|u-\Pi_h u\|_{{\mathbb {V}}},
			\]
		where $I_h w \in {\mathbb {V}}_h$ is the interpolator of
			$w.$
			 
		Taking into account the regularity given by Proposition \ref{prop:regHr} with $r=0$, interpolation estimate \eqref{estiminterpol} gives
			\[
			      \|w - I_h w\|_{{\mathbb {V}}} \le C h^{\alpha}  |w|_{H^{s+\alpha}(\Omega)} 
			      \le C  h^{\alpha} \|u-\Pi_h u\|_{L^2(\Omega)}.
			\]
		Finally, using the error estimate \eqref{estimacionlinf1} we obtain 
		\begin{equation}\label{estimacionlinf} 		 
			 \|u-\Pi_h u\|_{L^2(\Omega)}^2 \le C h^{\nicefrac12+\alpha-\varepsilon}|u|_{H^{s+\nicefrac12-\varepsilon}(\Omega)}\|u-\Pi_h u\|_{L^2(\Omega)},
			 \end{equation}
		and then estimate \eqref{eq:aproxL2} follows.
	\end{proof}

\section{Convergence of eigenvalues and eigenfunctions in gap distance}
\label{sec.conv}

	In this section we prove Theorem \ref{teo:gap}, stating that the discrete eigenvalue problem 
	\eqref{eq:int.autovalores.di.intro}  converges to the continuous one 
	\eqref{eq:int.autovalores.de.intro} in the gap distance (recall the definition of convergence
	given in the introduction).
	Let us start by defining the solution operators of the continuous 
	and discrete problems, $T:L^2(\Omega)\to {\mathbb {V}}$ and 
	$T_h:L^2(\Omega)\to {\mathbb {V}}_h$. Given  $f\in L^2(\Omega)$, 
	we define $Tf\in {\mathbb {V}}$ as the  unique solution of 
	\begin{equation}\label{ecuT} 
		\left\langle Tf,v \right\rangle =( f, v) \quad \forall v\in {\mathbb {V}},
	\end{equation}
 	and    $T_h f\in {\mathbb {V}}_h$ as the  unique solution of 
 	\[
 		\left\langle T_h f,v_h \right\rangle =(f,v_h) \quad \forall v_h\in {\mathbb {V}}_h. 
	\]
	Observe that if $(u,\lambda)$ is an eigenpair, then $T(\lambda u) = u$ and $T_h (\lambda u) = \Pi_h u$. 
	
	To prove Theorem \ref{teo:gap}, by \cite[Proposition 7.4 and Remark 7.5]{Boffi}, we only need to show that
	the operators $T$ and $T_h$ are compact and
	\[
		\| T- T_h\|_{\mathcal{L}(L^2(\Omega),{\mathbb {V}})} \to 0 \mbox{ as } h \to 0.
	\]
	
	\medskip 
	
	\begin{lemma}\label{compacto}
		The operators $T$ and $T_h$ are compact.
	\end{lemma}
	\begin{proof}
		Let $\{f_k\}_{k\in\mathbb{N}}$ 
		be a bounded sequence in $L^2(\Omega).$  
		Then, there exists a subsequence of  $\{f_k\}_{k\in\mathbb{N}}$ 
		(still denoted by $\{f_k\}$) and $f\in L^2(\Omega)$ such that  
		$f_k \rightharpoonup f$ weakly in $L^2(\Omega).$ 
		Taking $v=T f_k$ in \eqref{ecuT},  we get
		\[ 
			\|T f_k\|^2_{{\mathbb {V}}}=(f_k, T f_k)\leq C \|T f_k\|_{L^2(\Omega)}
			\quad\forall k\in\mathbb{N}.
		\]
		Therefore, by Poincar\'e inequality \eqref{eq:poincare},  we have that 
		$\{T f_k\}_{k\in\mathbb{N}}$ is bounded in ${\mathbb {V}}$. 
		Thus, there exists a subsequence of  $\{f_k\}_{k\in\mathbb{N}}$
		(still denoted by $\{f_k\}_{k\in\mathbb{N}}$) and $u\in{\mathbb {V}}$
		such that  $T f_k \rightharpoonup u$ weakly in ${\mathbb {V}}.$ Hence
		\[
			\left\langle u,v \right\rangle =\lim_{k\to\infty} \left\langle Tf_k,v \right\rangle =\lim_{k\to\infty}( f_k, v)
			=( f, v)\quad\forall v\in{\mathbb {V}},
		\]
		that is, $u=Tf.$
		
		On the other hand, since the inclusion ${\mathbb {V}}\hookrightarrow L^2({\mathbb {R}}^n)$ is compact,
		passing if necessary to a subsequence we may assume
		\begin{align*}
			T f_k \rightharpoonup Tf &\mbox{ weakly in } {\mathbb {V}},\\
			T f_k \rightarrow Tf & \mbox{ strongly in } L^2({\mathbb {R}}^n).
		\end{align*}
		Then,
		\[ 
			\|T f_k\|_{{\mathbb {V}}}^2=(f_k, T f_k) \to (f, T f)=\|T f\|^2_{{\mathbb {V}}}
		\]
		as $n\to\infty.$ Since the space ${\mathbb {V}}$ is uniformly convex, it follows that
		\[
			T f_k \rightarrow Tf \mbox{ strongly in } {\mathbb {V}}.
		\]
		
		The proof for the operator $T_h$ is similar.
	\end{proof}
 
	 \begin{lemma}\label{convop} 
		The following norm convergence holds true:
		\[
			\| T- T_h\|_{\mathcal{L}(L^2(\Omega),{\mathbb {V}})} \to 0 \mbox{ as } h \to 0.
		\]
 	\end{lemma}

	\begin{proof}
 		For each $h$, take $f_h\in L^2(\Omega) $ 
		such that $\|f_h\|_{L^2(\Omega)}=1$ and
		\[
			\sup_{\|f\|_{L^2(\Omega)}=1} \| T f - T_h f\|_{{\mathbb {V}}}
			= \| T f_h - T_h f_h\|_{{\mathbb {V}}}.
		\]
		
		Then, to prove the result, it is enough to show that for any sequence $h_k\to0$
		there is a subsequence $\{h_{k_j}\}_{j\in{\mathbb {N}}}$ such that 
		\[
			\| T f_{h_{k_j}} - T_{h_{k_j}} f_{h_{k_j}}\|_{{\mathbb {V}}}\to 0 \quad\text{ as } j\to\infty.
		\]
		
		Let $\{h_k\}_{k\in{\mathbb {N}}}$ be a sequence sucht that $h_{k}\to0.$
 		It follows from $\|f_{h_k}\|_{L^2(\Omega)}=1$ for all $k\in{\mathbb {N}}$ that  
 		there exist a subsequence $\{f_{h_{k_j}}\}_{j\in\mathbb{N}}$ 
		of $\{f_{h_k}\}_{k\in{\mathbb {N}}}$ and $f\in L^2(\Omega)$ such that  
		$ f_ {h_{k_j}}  \rightharpoonup  f$ weakly in $L^2(\Omega)$. 
		Proceeding as in the proof of Lemma \ref{compacto} 
		(and passing if necessary to a subsequence), we may assume
 		\begin{align*}
			T_{h_{k_j}} f_{h_{k_j}} \rightharpoonup v &\mbox{ weakly in } {\mathbb {V}}, \\
			T_{h_{k_j}} f_{h_{k_j}} \rightarrow v &\mbox{ strongly in } L^2({\mathbb {R}}^n).
		\end{align*}
		
		On the other hand, it follows from \eqref{estiminterpol} that 
		\begin{align*}
			I_h \varphi  \rightarrow \varphi &\mbox{ strongly in } {\mathbb {V}}, \\
			I_h \varphi  \rightarrow \varphi &\mbox{ strongly in }  L^2({\mathbb {R}}^n),
		\end{align*}
		for any $\varphi\in C^{\infty}_0(\Omega).$ Therefore,  
		\[
			\left\langle v,\varphi \right\rangle =\lim_{j\to \infty} \left\langle T_{h_{k_j}} f_{h_{k_j}}, 
				I_{h_{k_j}} \varphi \right\rangle =\lim_{j \to \infty }
				( f_{h_{k_j}} ,I_{h_{k_j}} \varphi)
				=( f , \varphi)\quad\forall\varphi\in C^{\infty}_0(\Omega),
		\] 
		which means that $v=T f$. 
		Then,
		\[
			\| T f_{h_{k_j}} - T_{h_{k_j}} f_{h_{k_j}} \|^2_{{\mathbb {V}}}=( f_{h_{k_j}}, T f_{h_{k_j}} - T_{h_{k_j}} f_{h_{k_j}}) \to 0.
		\]
	 \end{proof}
	Now we conclude the convergence of the discrete eigenvalue problem to the continuous one in the gap distance.

	\begin{proof}[Proof of Theorem \ref{teo:gap}]
	      The proof follows by  Lemmas \ref{compacto} and \ref{convop} 
	      and using Proposition 7.4 and Remark 7.5 of \cite{Boffi}.
	\end{proof}

\section{Order of convergence: simple eigenfunctions}\label{sec.orden}

	In this section we assume that $\lambda^{(k)}$ is an eigenvalue of multiplicity one; in the next section we 
	extend the discussion in order to tackle the case of multiple eigenfunctions. 

	In order to simplify, we split the proof of Theorem \ref{teo:ordensimple} into three lemmas.
	We first estimate the rate of convergence of $u_h^{(k)}$ towards $u^{(k)}$ in the $L^2$-norm. 

	\begin{lemma}\label{lema:ordenl2}
		For any $\varepsilon>0$ there 
		is a positive constant $C$ independent of $h$ such that
		\begin{equation}
				\|u^{(k)} - u_h^{(k)} \|_{L^2(\Omega)} \leq C h^{\alpha + \nicefrac12 -\varepsilon},
					\label{eq:convL2}
		\end{equation} 
		where  $u_h^{(k)}$ is chosen in such a way that 
		$\left(\Pi_h u^{(k)}, u_h^{(k)}\right) \ge 0$ and $\alpha$ is given by Proposition \ref{estimacionlinf0}. 
	\end{lemma}
	\begin{proof}
		In first place, we define the $L^2$-projection of $\Pi_h u^{(k)}$ into the space spanned by $u_h^{(k)}$,
		$$
			\omega_h^{(k)} \coloneqq \left( \Pi_h u^{(k)}, u_h^{(k)} \right) u_h^{(k)} ,
		$$
		and the quantity
		$$
			\rho_h^{(k)} \coloneqq \max_{i \neq k} \frac{\lambda^{(k)}}{|\lambda^{(k)} - \lambda_h^{(i)}|}.
		$$

		Then
		\begin{equation}
			\|u^{(k)} - u_h^{(k)} \|_{L^2(\Omega)} \leq \|u^{(k)} - \Pi_h u^{(k)} \|_{L^2(\Omega)} + \| \Pi_h u^{(k)} - \omega_h^{(k)} 
			\|_{L^2(\Omega)} + \| \omega_h^{(k)} - u_h^{(k)} \|_{L^2(\Omega)} .
			\label{eq:triangular}
		\end{equation}
		We are going to estimate the terms in the right hand side separately. 
		
		Given $\varepsilon>0,$ it follows from our regularity estimate \eqref{eq:aproxL2} that
		there exists $C>0$ independent of $h$ such that 
		\begin{equation}\label{eq:des1}
				\| u^{(k)} - \Pi_h u^{(k)} \|_{L^2(\Omega)} \leq C h^{\alpha+\nicefrac12-\varepsilon} .
		\end{equation}
		Moreover, since 
		$$
			\left( \Pi_h u^{(k)}, u_h^{(i)} \right) = \frac{1}{\lambda_h^{(i)}} \left\langle \Pi_h u^{(k)}, u_h^{(i)} \right\rangle = 
			\frac{1}{\lambda_h^{(i)}} \left\langle u^{(k)}, u_h^{(i)} \right\rangle = \frac{\lambda^{(k)}}{\lambda_h^{(i)}} \left( u^{(k)}, 
			u_h^{(i)} \right), 
		$$
		we have
		\[
			\left| \left( \Pi_h u^{(k)}, u_h^{(i)} \right) \right| \leq \rho_h^{(k)} \left| \left( u^{(k)} - \Pi_h u^{(k)}, 
			u_h^{(i)} \right) \right|.
		\]

		So,  
		\begin{equation}\label{eq:des2}
			\begin{split}
				\| \Pi_h u^{(k)} - \omega_h^{(k)} \|_{L^2(\Omega)}^2 & = \sum_{i\neq k} \left( \Pi_h u^{(k)}, u_h^{(i)} \right)^2 \leq \left[\rho_h^{(k)}\right]^{2} \sum_{i\neq k} \left( u^{(k)} - \Pi_h u^{(k)}, u_h^{(i)} \right)^2  \\
				& \leq \left[\rho_h^{(k)}\right]^2 \| u^{(k)} - \Pi_h u^{(k)} \|^2_{L^2(\Omega)} \leq C h^{\alpha+\nicefrac12 -\varepsilon} .
			\end{split}
		\end{equation}


		Finally, let us show that 
		\begin{equation}
			\| \omega_h^{(k)} - u_h^{(k)} \|_{L^2(\Omega)} \leq \| \omega_h^{(k)} - u^{(k)} \|_{L^2(\Omega)} ,
			\label{eq:estimacion_proyeccion}
		\end{equation}
		so that
		$$ 
			\| \omega_h^{(k)} - u_h^{(k)} \|_{L^2(\Omega)} \leq \| u^{(k)} - \Pi_h u^{(k)} \|_{L^2(\Omega)} + \| \Pi_h u^{(k)} - \omega_h^{(k)} \|_{L^2(\Omega)} .
		$$
		Indeed, on one hand we have
		$$
			u_h^{(k)}-\omega_h^{(k)} = \left[1 - \left(\Pi_h u^{(k)}, u_h^{(k)}\right)  \right] u_h^{(k)}.
		$$

		On the other hand, due to the normalizations $\|u^{(k)}\|_{L^2(\Omega)}=\|u_h^{(k)}\|_{L^2(\Omega)}=1,$ we have 
		$$
			\left| 1 - \|\omega_h^{(k)} \|_{L^2(\Omega)} \right| \leq \|u^{(k)}- \omega_h^{(k)} \|_{L^2(\Omega)}
		$$
		and
		$$
			\|\omega_h^{(k)}\|_{L^2(\Omega)} = \left|\left(\Pi_h u^{(k)}, u_h^{(k)}\right)\right|.
		$$

		Therefore, if we choose the sign of $u_h^{(k)}$ in such a way that $\left(\Pi_h u^{(k)}, u_h^{(k)}\right) \ge 0$, 
		we deduce
		$$
			\|u_h^{(k)} - \omega_h^{(k)} \|_{L^2(\Omega)} = \left|1 - \left(\Pi_h u^{(k)}, u_h^{(k)}\right)  \right| = 
			\left| 1 -  \left|\left(\Pi_h u^{(k)}, u_h^{(k)}\right)\right| \, \right| \le 
			\|u^{(k)} - \omega_h^{(k)} \|_{L^2(\Omega)},
		$$
		as stated in \eqref{eq:estimacion_proyeccion}.


		Thus, estimate \eqref{eq:convL2} is obtained by combining \eqref{eq:triangular}, \eqref{eq:des1}, \eqref{eq:des2} and \eqref{eq:estimacion_proyeccion}.
	\end{proof}

	Our next step is to prove the rate of convergence for the eigenvalues.
	
	\begin{lemma}
		For any $\varepsilon>0$ there is
		a positive constant $C$ independent of $h$ such that
		\begin{equation}\label{ordenaut}
			0\leq \lambda_h^{(k)}-\lambda^{(k)}\leq C h^{1-\varepsilon} .
		\end{equation}
	\end{lemma}
	\begin{proof}
		By Remark \ref{re:compara} we only have to prove the second inequality. 
		Developing the inner product and taking into account that $u^{(k)}$ is an eigenfunction and $u^{(k)}_h$ is its 
		corresponding discrete eigenfunction, we obtain
		\begin{align*}
			\| u^{(k)}-u^{(k)}_h\|^2_{{\mathbb {V}}} & = \left\langle u^{(k)}-u^{(k)}_h, \, u^{(k)}-u^{(k)}_h \right\rangle\\
  			& = \lambda^{(k)} \left(u^{(k)}, u^{(k)}-u^{(k)}_h\right) - \lambda^{(k)} \left(u^{(k)}, u^{(k)}_h\right) + 
  			\lambda^{(k)}_h \left(u^{(k)}_h, u^{(k)}_h\right)  \\
			& = \lambda^{(k)} \left(u^{(k)}-u^{(k)}_h, u^{(k)}-u_h^{(k)}\right) + \left(\lambda^{(k)}_h - \lambda^{(k)}\right) 
			\left(u^{(k)}_h, u^{(k)}_h\right) .
		\end{align*}
		That is,
		\begin{equation}\label{igualdadauto}
			\left(\lambda^{(k)}_h-\lambda^{(k)}\right)\|u^{(k)}_h\|^2_{L^2(\Omega)}=\|u^{(k)}-u^{(k)}_h\|^2_{{\mathbb {V}}}-\lambda^{(k)} 
			\|u^{(k)}-u^{(k)}_h\|^2_{L^2(\Omega)}.
		\end{equation}
		On the other hand,
		\begin{align*}
 			\|u^{(k)}_h-\Pi_h u^{(k)}\|^2_{{\mathbb {V}}}&=
 			\left(\lambda^{(k)}_h u^{(k)}_h-\lambda^{(k)} u^{(k)},u^{(k)}_h-\Pi_h u^{(k)}\right)\\
			 &\leq  \|\lambda^{(k)}_h u^{(k)}_h-\lambda^{(k)} u^{(k)}\|_{L^2(\Omega)} \|u^{(k)}_h-\Pi_h u^{(k)} \|_{L^2(\Omega)}\\
 			&\leq  \left(\lambda^{(k)}_h \| u^{(k)}_h-u^{(k)}\|_{L^2(\Omega)}+|\lambda^{(k)}_h-\lambda^{(k)}| 
 			\|u^{(k)}\|_{L^2(\Omega)}\right) \|u^{(k)}_h-\Pi_h u^{(k)} \|_{L^2(\Omega)}.
 		\end{align*}
		Dividing by $ \|u^{(k)}_h-\Pi_h u^{(k)} \|_{L^2(\Omega)}$ and using the fact that 
		\[
			\|u^{(k)}_h-\Pi_h u^{(k)}\|_{{\mathbb {V}}}^2\geq \lambda^{(1)} \|u^{(k)}_h-\Pi_h u^{(k)} \|_{L^2(\Omega)}^2,
		\] 
		we get
		\begin{equation}\label{cota2auto}
			\|u^{(k)}_h-\Pi_h u^{(k)}\|_{{\mathbb {V}}}\leq 
			C \left( \lambda^{(k)}_h \| u^{(k)}_h-u^{(k)}\|_{L^2(\Omega)}+|\lambda^{(k)}_h-\lambda^{(k)}| \|u^{(k)}\|_{L^2(\Omega)} 
			\right).
		\end{equation}
 
		Therefore, by  \eqref{igualdadauto}, \eqref{cota2auto}, \eqref{estimacionlinf1} and taking into account that 
		$\|u_h^{(k)}\|_{L^2(\Omega)}= \|u^{(k)}\|_{L^2(\Omega)}=1,$ we arrive at
 		\begin{align*} 
 			\lambda^{(k)}_h-\lambda^{(k)} 
 			& \leq C\left(\| u^{(k)}-\Pi_h u^{(k)}\|^2_{{\mathbb {V}}}+\|u^{(k)}_h-\Pi_h u^{(k)}\|^2_{{\mathbb {V}}}\right)-\lambda^{(k)} 
 			\|u^{(k)}-u^{(k)}_h\|^2_{L^2(\Omega)}\\
 			& \leq C\left(h^{1-\varepsilon}+ \left[\lambda^{(k)}_h\right]^2 \| 
 			u^{(k)}_h-u^{(k)}\|^2_{L^2(\Omega)}+|\lambda^{(k)}_h-\lambda^{(k)}|^2 \right)- \lambda^{(k)} 
 			\|u^{(k)}-u^{(k)}_h\|^2_{L^2(\Omega)}\\
 			&= C\left(h^{1-\varepsilon}+|\lambda^{(k)}_h-\lambda^{(k)}|^2 
 			\right)+\left(C\left[\lambda^{(k)}_h\right]^2-\lambda^{(k)}\right)\|u^{(k)}_h-u^{(k)}\|^2_{L^2(\Omega)}.
 		\end{align*}
 
 		Finally, using that  $\lambda^{(k)}_h\to \lambda^{(k)}$ and that 
 		$\|u^{(k)}-u^{(k)}_h\|^2_{L^2(\Omega)}\leq C h^{2\alpha+1}$, we conclude
		$$
			\lambda^{(k)}_h-\lambda^{(k)}\leq C h^{1-\varepsilon}+O(h^{2\alpha+1})\leq C h^{1-\varepsilon}.
		$$ 
	\end{proof}

	In last place, we have the following estimate for the error in the energy norm. 

	\begin{lemma} \label{conv_V}
		 For any $\varepsilon>0$ there is a positive constant $C$ independent of $h$ 
		such that
		\[
			\|u^{(k)}-u^{(k)}_h\|_{{\mathbb {V}}}\leq C h^{\nicefrac12-\varepsilon},
		\]
		where  $u_h^{(k)}$ is chosen in such a way that $(\Pi_h u^{(k)}, u_h^{(k)}) \ge 0$. 
	\end{lemma}
	\begin{proof}
		By \eqref{igualdadauto} and   since $\|u^{(k)}_h\|_{L^2(\Omega)} = \| u^{(k)} \|_{L^2(\Omega)}=1$, we have that
		\[
			\|u^{(k)} - u^{(k)}_h\|_{{\mathbb {V}}} \leq C \left( \|u^{(k)}-u^{(k)}_h\|_{L^2(\Omega)} + \sqrt{\lambda^{(k)}_h - 
			\lambda^{(k)}} \right) \leq C h^{\nicefrac12 - \varepsilon}.
		\]
	\end{proof}

\section{Order of convergence: Multiple eigenfunctions}
\label{sec.orden_multiple}

	In this section we show the order of convergence in the case of multiple eigenspaces. For simplicity we prove the case of 
	double eigenvalues, but the argument we carry may be easily generalized to any multiplicity.
	We write $\lambda=\lambda^{(k)}=\lambda^{(k+1)}$ and take $u^{(k)}$ and $u^{(k+1)}$ to be two corresponding orthogonal eigenfunctions. 	
	Throughout this section we also employ the notation
	\[ 
		\alpha_h = \left( \Pi_h u^{(k)}, u_h^{(k)} \right), 
		\beta_h = \left( \Pi_h u^{(k)}, u_h^{(k+1)} \right),
		\delta_h = \left( \Pi_h u^{(k+1)}, u_h^{(k)} \right), 
		\gamma_h = \left( \Pi_h u^{(k+1)}, u_h^{(k+1)} \right)
	\] 
	\[
		w_h^{(k)}= \alpha_h u_h^{(k)} + \beta_h u_h^{(k+1)},\quad
		\text{ and } \quad
		w_h^{(k+1)}= \delta_h u_h^{(k)} + \gamma_h u_h^{(k+1)}.
	\]	
	
	As in the previous section, we divide the proof of Theorem \ref{teo:ordenmultiple} into three lemmas.	
	The proof of the first one is obtained as in \cite{Boffi}, following the lines  of Lemma \ref{lema:ordenl2}.
		

	\begin{lemma}\label{lema:ordenmultiple1}
		For any $\varepsilon>0,$ there is a positive constant $C$ independent of $h$ 
		such that
		\begin{equation}\label{ordenmultautovfun}
			\|u^{(j)}-w^{(j)}_h\|_{L^2(\Omega)}\leq C h^{\alpha+\nicefrac12-\varepsilon},
			 \quad j=k,k+1.
		\end{equation}
		where $\alpha$ is given by Proposition \ref{estimacionlinf0}.
	\end{lemma}
	
	Our next step is to prove the rate of convergence for the eigenvalues.

	\begin{lemma}\label{lem:ordenmultiple2}
		For any $\varepsilon>0,$ there is a positive constant $C$ independent of $h$ 
		such that
		\begin{equation}\label{ordenmultiaut}
				0\leq \lambda_h^{(j)}-\lambda\leq C h^{1-\varepsilon}, \quad j=k,k+1.
		\end{equation}
	\end{lemma}
	\begin{proof}
		
		By Remark \ref{re:compara}, in order to prove \eqref{ordenmultiaut}  we only have to show the second inequality. First observe 
		that an identity analogous to \eqref{igualdadauto} is feasible.
		More precisely, in this case we have
		\begin{equation}\label{igualdadauto2}
			\alpha_h^2 \left(\lambda_h^{(k)}-\lambda\right)+ 
			\beta_h^2 \left(\lambda_h^{(k+1)}-\lambda\right) =
			\|u^{(k)}-w^{(k)}_h\|^2_{{\mathbb {V}}}- 
			\lambda \|u^{(k)}-w^{(k)}_h\|^2_{L^2(\Omega)}.
		\end{equation}
		On the other hand,
		\begin{align*}
			 \|w^{(k)}_h-\Pi_h u^{(k)}\|^2_{{\mathbb {V}}}
		 		&=\left\langle u^{(k)}-w_h^{(k)}, \Pi_h u^{(k)}-w^{(k)}_h \right\rangle  \\
		 		&=\left(\lambda u^{(k)},\Pi_h u^{(k)}-w^{(k)}_h\right)-
		 		\left(\alpha_h \lambda_h^{(k)} u_h^{(k)}+\beta_h \lambda_h^{(k+1)} 
		 		u_h^{(k+1)}, \Pi_h u^{(k)}-w^{(k)}_h\right)\\
				& \leq \|\lambda u^{(k)}-\alpha_h \lambda_h^{(k)} u_h^{(k)}-\beta_h \lambda_h^{(k+1)} u_h^{(k+1)}\|_{L^2(\Omega)} 
				\|\Pi_h u^{(k)}-w^{(k)}_h\|_{L^2(\Omega)}. 
		\end{align*}
		Thus, by the Poincar\'e inequality 
		\begin{align*}
			\|w^{(k)}_h-\Pi_h u^{(k)}\|_{{\mathbb {V}}}&\leq
			C \left[ \|\lambda u^{(k)}-\alpha_h \lambda_h^{(k)} u_h^{(k)}-\beta_h 
				\lambda_h^{(k+1)} u_h^{(k+1)}\|_{L^2(\Omega)} \right] \\ 
			&\leq C \Big[ \lambda^{(k)}_h \| w^{(k)}_h-u^{(k)}\|_{L^2(\Omega)}+|\lambda^{(k)}_h-\lambda| 
				\|u^{(k)}\|_{L^2(\Omega)}\\
			&\hspace{1cm}+|\beta_h| \|u_h^{(k+1)}\|_{L^2(\Omega)}|\lambda^{(k)}_h-\lambda^{(k+1)}_h| \Big].
		\end{align*}
		Then, since $\|u^{(k)}\|^2_{L^2(\Omega)}= \|u_h^{(k+1)}\|^2_{L^2(\Omega)}=1,$ we obtain
		\begin{equation}\label{cota2auto2}
			\begin{aligned}
				\|w^{(k)}_h-\Pi_h u^{(k)}\|_{{\mathbb {V}}}&\leq  
				C \Big[ \lambda^{(k)}_h 
				\|w^{(k)}_h-u^{(k)}\|_{L^2(\Omega)}+
				|\lambda^{(k)}_h-\lambda|\\
				&\hspace{1.5cm}+|\beta_h|\left(|\lambda^{(k)}_h-\lambda|+|\lambda^{(k+1)}_h-\lambda|
				\right) \Big].
			\end{aligned}
		\end{equation}
		Therefore, by  \eqref{igualdadauto2}, \eqref{cota2auto2} and \eqref{estimacionlinf1} we arrive at
		 \begin{align*}
		 	\alpha_h^2 &\left(\lambda_h^{(k)}-\lambda\right)+ 
		 		\beta_h^2\left(\lambda_h^{(k+1)}-\lambda\right)\\ 
		 	& \leq C\left(
		 		\|u^{(k)}-\Pi_h u^{(k)}\|^2_{{\mathbb {V}}}+\|w^{(k)}_h-\Pi_h u^{(k)}\|^2_{{\mathbb {V}}}\right)-\lambda 
		 		\|u^{(k)}-w^{(k)}_h\|^2_{L^2(\Omega)}\\
		 	& \leq C\left(h^{1-\varepsilon}+ \left[\lambda^{(k)}_h\right]^2 \| 	
		 		w^{(k)}_h-u^{(k)}\|^2_{L^2(\Omega)}+|\lambda^{(k)}_h-\lambda|^2 
		 		+\beta_h^2\left(|\lambda^{(k)}_h-\lambda|+|\lambda^{(k+1)}_h-\lambda|\right)^2\right)\\
		 	&\hspace{1cm}- \lambda \|u^{(k)}-w^{(k)}_h\|^2_{L^2(\Omega)}. 
		 \end{align*}
		
		By Lemma \ref{lema:ordenmultiple1}, 
		$\|u^{(k)}-w^{(k)}_h\|^2_{L^2(\Omega)}\leq C h^{2\alpha+1-\varepsilon}\leq C h^{1-\varepsilon},$
		and thus
		\begin{equation}\label{multival} 
			\alpha_h^2 \left(\lambda_h^{(k)}-\lambda\right)+ \beta_h^2 \left(\lambda_h^{(k+1)}-\lambda\right)  
			\leq C \left(h^{1-\varepsilon}+  |\lambda^{(k)}_h-\lambda|^2+ |\lambda^{(k+1)}_h-\lambda|^2\right).
		\end{equation} 
		
		Taking into account that $\alpha_h^2+\beta_h^2=\|w^{(k)}_h\|^2_{L^2(\Omega)}\to \|u^{(k)}\|^2_{L^2(\Omega)}=1$, we split 
		the rest of the proof in three cases.
		
		\medskip
		
		{\noindent{\it Case 1}.} There is a positive constant $c$ such that $|\alpha_h|\geq c$ and $ |\beta_h|\geq c$.  
		 
		 Since  $\lambda^{(k)}_h\to \lambda$, $\lambda^{(k+1)}_h\to \lambda$ we conclude, by \eqref{multival}, that 
		 \eqref{ordenmultiaut} follows.
	
		\medskip
	
		{\noindent{\it Case 2}.} There exist subsequences satisfying $\alpha_{h_n}\to 0$ and $|\beta_{h_n}|\to 1$. 
			
			Observe that
			\begin{equation}\label{igualdadalphabeta}
				\begin{aligned}
					\left(\alpha_h^2+\beta_h^2\right)
							\left(\left(\lambda_h^{(k)}-\lambda\right)+ \left(\lambda_h^{(k+1)}-\lambda\right)\right)
							&=2\left(\alpha_h^2 \left(\lambda_h^{(k)}-\lambda\right)+ 
									\beta_h^2 \left(\lambda_h^{(k+1)}-\lambda\right) \right)\\
							&\hspace{1cm}+ \left(\alpha_h^2-\beta_h^2\right) \left(\lambda_h^{(k+1)}-\lambda_h^{(k)}\right).
				\end{aligned}
			\end{equation}
		Since $\alpha_{h_n}\to 0$ and $|\beta_{h_n}|\to 1$, for $n$ great enough $|\alpha_{h_n}|\leq |\beta_{h_n}|$, and  
		therefore by \eqref{igualdadalphabeta},  \eqref{multival} and the fact that $\alpha_h^2+\beta_h^2 \to 1$ we have that 
		\eqref{ordenmultiaut} is true for a subsequence. But since $\{\lambda^{(j)}_h \}$ is monotone in $h$ ($j = k, k+1$), the desired conclusion 
		follows.
		
		\medskip
			
		{\noindent{\it Case 3}.} There exist subsequences satisfying $\beta_{h_n}\to 0$ and $|\alpha_{h_n}|\to 1$.
		
		First, observe that
		$$
			u_h^{(k)}=w_h^{(k)}+(1-\alpha_h)u_h^{(k)}-\beta_h u_h^{(k+1)}.
		$$
		Then, since $\|u_h^{(k+1)}\|_{L^2(\Omega)}=\|u_h^{(k)}\|_{L^2(\Omega)}=1,$ $w_h^{(k)} \to u^{(k)},$ 
		$\alpha_{h_n}\to 1$ and $\beta_{h_n}\to 0,$ we get
		\begin{equation}\label{convato}
			u_{h_n}^{(k)}\to u^{(k)} \mbox{ in } L^2(\Omega).
		\end{equation}
		
		On the other hand, it follows form the definition of $\Pi_h$ that
		\[ 
			\alpha_h =\frac{\lambda}{\lambda_h^k} \left(u^{(k)},u_h^{(k)}\right),\ \
			\beta_h =\frac{\lambda}{\lambda_h^{k+1}} \left(u^{(k)},u_h^{(k+1)}\right),
		\]
		\[
			\delta_h =\frac{\lambda}{\lambda_h^{(k)}} \left(u^{(k+1)},u_h^{(k)}\right),\quad
			\gamma_h = \frac{\lambda}{\lambda_h^{(k+1)}} \left(u^{(k+1)},u_h^{(k+1)}\right).
		\] 
		
		By Lemma \ref{lema:ordenmultiple1}, 
		\[
			\|u^{(k+1)}-w^{(k+1)}_h\|_{L^2(\Omega)}\leq C h^{\alpha+\nicefrac12-\varepsilon}
		\]
		and then
		\begin{equation}\label{gamadelta} 
		 \delta_h^2+\gamma_h^2 \to 1.\end{equation} 
		
		Finally, by \eqref{convato},  using the orthogonality of the eigenfunctions and \eqref{gamadelta} we have 
		\[
			\delta_{h_n}=\frac{\lambda}{\lambda_{h_n}^{(k)}} \left(u^{(k+1)},u_{h_n}^{(k)}\right)\to 
			\left(u^{(k+1)},u^{(k)}\right)=0\quad\text{ and } \quad|\gamma_{h_n}|\to 1.
		\]
		Observe that changing $u^{(k)}$ by $u^{(k+1)}$, this is just {\it{Case 2}}. This concludes the proof.	
	\end{proof}

	Our last result follows from \eqref{igualdadauto2}, \eqref{ordenmultautovfun} and \eqref{ordenmultiaut}, in the same fashion as Lemma \ref{conv_V}.
	
	\begin{lemma}
		For any $\varepsilon>0,$ there is a positive constant $C$ independent of $h$ 
		such that
		\[
			\|u^{(j)}-w^{(j)}_h\|_{{\mathbb {V}}}\leq C h^{\nicefrac12-\varepsilon},
			\quad  j=k,k+1.			
		\]
	\end{lemma} 

\section{Numerical results}
\label{sec:numerico}
In order to illustrate the convergence estimates obtained in Sections \ref{sec.orden} and \ref{sec.orden_multiple}, we present the results of numerical tests for finite element discretizations of one and two-dimensional eigenvalue problems. Moreover, in the latter case, some examples in domains that do not satisfy the hypotheses of Proposition 
\ref{prop:regintro} are displayed. These examples provide numerical evidence that the assertion of this proposition still holds true under weaker assumptions about the domain.


Since in general no closed formula for the eigenvalues of the fractional Laplacian is available, we have estimated the order of convergence 
by means of a least-squares fitting of the model 
\[
\lambda^{(k)}_h = \lambda^{(k)} + Ch^\alpha .
\]
This allows to extrapolate approximations of the eigenvalues as well (in the tables we denote this extrapolated value of $\lambda^{(k)}$ as $\lambda_{ext}^{(k)}$).

Throughout this section, the results are compared with those available in the literature. In first place we consider one-dimensional problems, which have been widely studied both theoretically and from the numerical point of view. Next, we show some examples in two-dimensional domains: the unit ball, a square and an $L$-shaped domain.
As for the ball, the deep results of \cite{DydaKuznetsovKwasnicki} allow to obtain sharp estimates on the eigenvalues, and thus provide a point of comparison for the validity of the FE implementation. Regarding the square, some estimates for the eigenvalues are found in \cite{Kwasnicki}. The main interest of the $L$-shaped domain is that, although it does not satisfy the ``standard'' requirements to regularity of eigenfunctions to hold, the numerical order of convergence is the same as in problems posed on smooth, convex domains.


General estimates for eigenvalues, valid for a class of domains, are obtained by Chen and Song \cite{ChenSong}. In that paper, the authors state upper and lower bounds for eigenvalues of the fractional Laplacian on domains satisfying the exterior cone condition. Calling $\mu^{(k)}$ the $k$-th eigenvalue of the Laplacian with Dirichlet boundary conditione domain $\Omega$, they prove that there exists a constant $C=C(\Omega)$ such that
\begin{equation} \label{eq:chensong}
C \left(\mu^{(k)}\right)^s \le \lambda^{(k)} \le \left(\mu^{(k)}\right)^s .
\end{equation}
If $\Omega$ is a bounded convex domain, then $C$ can be taken as $\nicefrac12$. It is noteworthy that, due to the scaling property of the fractional Laplacian, eigenvalues for the dilations of a domain $\Omega$ are obtained by means of
 $\lambda^{(k)}(\gamma\Omega) = \gamma^{-2s} \lambda^{(k)}(\Omega).$ 


Since we are working with conforming methods, as well as providing approximations, the discrete eigenvalues yield upper bounds for the corresponding continuous eigenvalues. This is of special interest in those cases in which theoretical estimates are not sharp, or the non-symmetry of the domain precludes the possibility of developing arguments such as the ones in \cite{DydaKuznetsovKwasnicki}.



\subsection{One-dimensional intervals}
Eigenvalues for the fractional Laplacian in intervals have been studied by other authors previously. 
In \cite{ZoiaRossoKardar}, a discretized model of the fractional Laplacian is developed, and a numerical study of 
eigenfunctions and eigenvalues is implemented for different boundary conditions. 
In \cite{MR2679702}, the authors deal with one dimensional problems for $s=\nicefrac12$, and provide asymptotic expansion for eigenvalues.
Later, Kwa{\'s}nicki \cite{Kwasnicki} extended this work to the whole range $s\in (0,1)$. Namely, he showed the following identity 
for the $k$-th eigenvalue in the interval $(-1,1)$:
\begin{equation}
\lambda^{(k)} = \left( \frac{k\pi}{2} - \frac{(1-s)\pi}{4}\right)^{2s} + \frac{1-s}{\sqrt{s}} \, O\left(k^{-1}\right) .
\label{eq:kwasnicki}
\end{equation}
Moreover, in that work a method to obtain lower bounds  in arbitrary bounded domains is developed, and it is proved that, on one spacial dimension,
 eigenvalues are simple if $s\ge \nicefrac12$.
As eigenvalues are simple and we are working in one dimension, it is not difficult to numerically estimate the order of convergence of eigenfunctions 
in the $L^2$-norm. Indeed, normalizing the discrete eigenfunctions so that $\| u_h^{(k)} \|_{L^2(-1,1)} = 1$ and choosing their sign adequately (recall Section \ref{sec.orden}), 
these are then compared with a solution on a very fine grid.


On the other hand, in \cite{DuoZhang} it is performed a numerical study of the fractional Schr\"odinger equation in an infinite potential well 
in one spacial dimension. The authors find numerically the ground and first excited states and their corresponding eigenvalues for the stationary linear problem, 
which corresponds to the first two eigenpairs of our equation \eqref{eq:int.autovalores}. 


In Table \ref{tab:ordenes1d}, our results for the first $2$ eigenpairs are displayed. The extrapolated numerical values 
are compared with the estimates from \cite{DuoZhang, Kwasnicki}; the orders of convergence are in good 
agreement with those predicted correspondingly by \eqref{ordenaut} and \eqref{eq:convL2}. 

\begin{table}[htbp]
\label{tab:ordenes1d}	
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|} \cline{2-11}
\multicolumn{1}{c}{}		& \multicolumn{6}{|c|}{Numerical values} & \multicolumn{4}{c|}{Orders} \\ \hline
 $s$   & $\lambda_{ext}^{(1)}$ & $\lambda^{(1)}$ \cite{DuoZhang} & $\lambda^{(1)}$ \cite{Kwasnicki}  & $\lambda_{ext}^{(2)}$ & $\lambda^{(2)}$ \cite{DuoZhang} & $\lambda^{(2)}$ \cite{Kwasnicki}&  $\lambda^{(1)}$ & $\lambda^{(2)}$ & $u^{(1)}$ & $u^{(2)}$ \\ \hline
$0.05$&	$0.9726$								&	$0.9726$								&		$0.9809$								&				$1.0922$								&	$1.0922$								    &	$1.0913$										& $1.1082$										  &    	$1.1488$	 &     	$0.5507$	   &  $0.5676$  \\ 
$0.1$&			$0.9575$						& $0.9575$								&		$0.9712$								&			$1.1965$								&	$1.1966$									  &	$1.1948$																		  & $1.0706$   		 & $1.1017$     & $0.6117$  & $0.6250$   \\
$0.25$&	$0.9702$								& $0.9702$								&		$0.9908$								&				$1.6015$						&	$1.6016$									  &	$1.5977$										&		$1.0210$								  &  $1.0375$  		 &  $0.7616$   		   & $0.7823$   \\  
$0.5$&	$1.1577$								& $1.1578$									&	$1.1781$								&						$2.7548$				&	$2.7549$										&	$2.7488$																				  & $1.0005$   		 & $0.9793$		  	 & $0.9605$  & $0.9691$	  \\
$0.75$&	$1.5975$								& $1.5976$								&		$1.6114$								&								$5.0598$		&	$5.0600$									  &	$5.0545$										&			$0.9983$							  & $0.9990$   		 & $0.9983$    		 &  $0.9984$  \\ 
$0.9$&	$2.0487$								& ---					 						& $2.0555$									&		$7.5031$									&	---										&				$7.5003$							& $1.0035$   		 & $1.0214$    		 & $0.9989$  & $0.9989$   \\
$0.95$&	$2.2481$								& $2.2441$								&		$2.2477$								&		$8.5958$									&	$8.5959$									  &	$8.5942$										&		$1.0352$									  & $1.1418$   		 &  $0.9989$   		   & $0.9989$   \\ 
				\hline	 
\end{tabular}
\caption{First $2$ eigenpairs in the interval $(-1,1)$. On the left, the extrapolated numerical values are compared with the results from \cite{DuoZhang} and with approximation \eqref{eq:kwasnicki}, obtained in \cite{Kwasnicki}. On the right, orders of convergence for eigenvalues and eigenfunctions in the $L^2$-norm (obtained by a least-square fitting) are displayed.}
\end{table}


\subsection{Two-dimensional experiments} \label{sub:2d}
The theoretical order of convergence for eigenvalues is also attained in the following examples in ${\mathbb {R}}^2$. The implementation for these experiments is based on the code from \cite{AcostaBorthagaray}.


\subsubsection*{Unit ball} 
Let us consider the fractional eigenvalue problem on the two-dimensional unit ball. 
In \cite{DKKMeijer}, the weighted operator $u \mapsto (-\Delta)^s (\omega \, u)$  is studied, where $\omega (x) = (1-|x|^2)^s_+.$ In particular, explicit formulas for eigenvalues and eigenfunctions are established. 
Furthermore, in \cite{DydaKuznetsovKwasnicki} the same authors exploit these expressions to obtain two-sided bounds 
for the eigenvalues of the fractional Laplacian in the unit ball in any dimension.
This method provides sharp estimates; however, it depends on the decomposition of the fractional Laplacian as a weighted operator, and the weight $\omega$ is only explicitly known for the unit ball.

In Table \ref{tab:ordenes_bola}, our results for the first eigenvalue are compared with those of \cite{DydaKuznetsovKwasnicki} for different values of $s$. This comparison serves as a test for the validity of the code we are employing. As well as the extrapolated value of $\lambda^{(1)}$ and the numerical order of convergence, for every $s$ considered we exhibit an upper bound for the first eigenvalue. These outcomes are consistent with those from \cite{DydaKuznetsovKwasnicki} and the theoretical order of convergence given by \eqref{ordenaut}.


\begin{table}[htbp]	
\begin{tabular}{|c|c|c|c|c|} \hline  
		 $s$    &  ${\lambda^{(1)}}$& $\lambda_{ext}^{(1)}$ & $\lambda_h^{(1)}$ (UB)	& Order		\\ \hline 
		$0.005$ &  $1.00475$	 			& $1.00475$							& $1.00480$	  						& $0.9462$	\\
		$0.05$  &  $1.05095$ 	 			& $1.05094$							& $1.05145$								& $0.9455$	\\
		$0.25$  &  $1.34373$	 			& $1.34367$							& $1.34626$								& $0.9497$	\\
		$0.5$   &  $2.00612$	 			& $2.00607$							& $2.01060$								& $0.9686$	\\
		$0.75$  &  $3.27594$	 			& $3.27632$							& $3.28043$								& $1.0092$	\\
		\hline	 
		\end{tabular}
\caption{First eigenvalue in the unit ball in  ${\mathbb {R}}^2$. Estimate from \cite{DydaKuznetsovKwasnicki};  extrapolated value of $\lambda^{(1)}$; upper bound obtained by FEM with a meshsize $h\sim 0.02$; numerical order of convergence.} \label{tab:ordenes_bola}
\end{table}

\medskip

\subsubsection*{Square} Eigenvalue estimates for the case in which the domain $\Omega$ is a square in ${\mathbb {R}}^2$ were also addressed in \cite{Kwasnicki}. 
However, in order to obtain upper bounds, the method proposed in that work depends on 
having pointwise bounds of the Green function for the fractional Laplacian on $\Omega$. 
The estimates from \cite{ChenSong, Kwasnicki} are compared with our results in Table \ref{tab:ordenes_cuadrado}, where numerical orders of convergence are also displayed.

\begin{table}[htbp]	
\begin{tabular}{|c|c|c|c|c|c|} \hline
		 $s$   &  $\lambda^{(1)}$ (LB)& $\lambda^{(1)}$ (UB)  & $\lambda_h^{(1)}$ (UB)& $\lambda_{ext}^{(1)}$ & Order		\\ \hline 
		$0.05$ &  $1.0308^b$	      	& $1.0831^a$		   			& $1.0412$			   			& $1.0405$			 				& $0.9229$ 	\\
		$0.1$  & 	$1.0506^b$	      	& $1.1731^a$		 			  & $1.0895$			   			& $1.0882$			 				& $0.9230$	\\
		$0.25$ &	$1.1587^b$	      	& $1.4905^a$		  			& $1.2844$			   			& $1.2813$			 				& $0.9283$	\\
		$0.5$  & 	$1.3844^b$	      	& $2.2214^a$		   			& $1.8395$			   			& $1.8344$			 				& $0.9622$	\\
		$0.75$ &	$1.6555^a$	      	& $3.3109^a$		   			& $2.8921$			   			& $2.8872$			 				& $0.9940$	\\
		$0.9$  & 	$2.1034^a$	     		& $4.2067^a$		   			& $3.9492$			   			& $3.9467$			 				& $1.0654$	\\
		$0.95$ & 	$2.2781^a$	     		& $4.5562^a$		   			& $4.4083$			   			& $4.4062$			 				& $1.1496$	\\ \hline	 
		\multicolumn{6}{l}{\footnotesize \, $^a$See \cite{ChenSong}. $^b$See \cite{Kwasnicki}.} \\
		\end{tabular}
\caption{First eigenvalue in the square $[-1,1]^2$. Best lower (LB) and upper (UB) bounds known before; upper bound obtained by FEM with a meshsize $h\sim 0.04$; extrapolated value of $\lambda^{(1)}$; numerical order of convergence.}
\label{tab:ordenes_cuadrado}
\end{table}

\medskip


\subsubsection*{$L$-shaped domain}
To the authors knowledge, there is no \emph{efficient} method to estimate eigenvalues of the fractional Laplacian if the domain $\Omega$ lacks symmetry.
The bound \eqref{eq:chensong} remains valid as long as $\Omega$ satisfies the assumptions required, but the range that estimate provides is quite wide. 

The main advantage of employing the finite element method is that it is flexible enough to cope with a variety of domains.
Moreover, as we are working with conforming approximations, sharp upper bounds for the eigenvalues may be obtained by considering 
discrete solutions on refined meshes.

In Proposition \ref{prop:regintro}, which states that eigenfunctions belong to  $\widetilde{H}^{s+\nicefrac12-\varepsilon}(\Omega)$, it was assumed that the domain $\Omega$ was smooth. 
For the Laplacian, in order to prove regularity of solutions, it is customary to assume that $\Omega$ is either smooth or at least convex. In those cases, it is  well known that if $f \in H^{r}(\Omega)$ for some $r$, then the solutions of the Dirichlet problem with datum $f$ belong to $H^{r+2}(\Omega)$. However, if the domain has a re-entrant corner, solutions are less regular. This also applies to eigenfunctions: in the $L$-shaped domain $\Omega = [-1,1]^2\setminus [0,1]^2$, the first eigenvalue of the Laplacian is known not to belong to $H^{3/2}(\Omega)$.

Surprisingly, numerical evidence indicates that eigenvalues of the fractional Laplacian on this $L$-shaped domain converge with the same order as in the previous examples.
This motivates us to conjecture that eigenfunctions and solutions to the Dirichlet equation \eqref{eq:fuente} have the same Sobolev regularity than in smooth domains.

Our findings for the first eigenvalue are summarized in Table \ref{tab:ordenes_L}.


\begin{table}[htbp] 
\begin{tabular}{|c|c|c|c|} \hline
$s$  &  $\lambda_h^{(1)}$ (UB)& $\lambda_{ext}^{(1)}$ & Order $\lambda^{(1)}$ \\ \hline 
$0.1$ &  $1.1434$		 					& $1.1413$							& $0.9085$		\\
$0.2$ &  $1.3386$		 					& $1.3342$							& $0.9103$		\\
$0.3$ &  $1.6025$		 					& $1.5956$							& $0.9160$		\\
$0.4$ &  $1.9593$		 					& $1.9499$							& $0.9267$		\\
$0.5$ &  $2.4440$		 					& $2.4322$							& $0.9459$		\\
$0.6$ &  $3.1072$		 					& $3.0936$							& $0.9812$		\\
$0.7$ &  $4.0228$		 					& $4.0069$							& $0.9822$		\\
$0.8$ &  $5.2994$		 					& $5.2831$							& $1.0609$		\\
$0.9$ &  $7.0975$							& $7.0790$							& $1.1891$		\\ \hline	 
\end{tabular}
\caption{First eigenvalues in the $L-$shaped domain $[-1,1]^2\setminus [0,1]^2$. 
Upper bound obtained by FEM with a meshsize $h\sim 0.04$; extrapolated value of $\lambda^{(1)}$; numerical order of convergence.}
\label{tab:ordenes_L}
\end{table}



2016-03-02T02:13:01.000
1603
1603.00317
en
https://arxiv.org/abs/1603.00317
The purpose of this work is to study a finite element method for finding solutions to the eigenvalue problem for the fractional Laplacian. We prove that the discrete eigenvalue problem converges to the continuous one and we show the order of such convergence. Finally, we perform some numerical experiments and compare our results with previous work by other authors.
Numerical Analysis (math.NA)
Finite element approximation for the fractional eigenvalue problem
Qwen/Qwen-72B
1. YES
2. YES