  We design and analyze several Finite Element Methods (FEMs) applied to the
Caffarelli-Silvestre extension that localizes the fractional powers of
symmetric, coercive, linear elliptic operators in bounded domains with
Dirichlet boundary conditions. We consider open, bounded, polytopal but not
necessarily convex domains $\Omega \subset \mathbb{R}^d$ with $d=1,2$. For the
solution to the extension problem, we establish analytic regularity with
respect to the extended variable $y\in (0,\infty)$. We prove that the solution
belongs to countably normed, power-exponentially weighted Bochner spaces of
analytic functions with respect to $y$, taking values in corner-weighted
Kondat'ev type Sobolev spaces in $\Omega$. In $\Omega\subset \mathbb{R}^d$, we
discretize with continuous, piecewise linear, Lagrangian FEM ($P_1$-FEM) with
mesh refinement near corners, and prove that first order convergence rate is
attained for compatible data $f\in \mathbb{H}^{1-s}(\Omega)$.
  We also prove that tensorization of a $P_1$-FEM in $\Omega$ with a suitable
$hp$-FEM in the extended variable achieves log-linear complexity with respect
to $\mathcal{N}_\Omega$, the number of degrees of freedom in the domain
$\Omega$. In addition, we propose a novel, sparse tensor product FEM based on a
multilevel $P_1$-FEM in $\Omega$ and on a $P_1$-FEM on radical-geometric meshes
in the extended variable. We prove that this approach also achieves log-linear
complexity with respect to $\mathcal{N}_\Omega$. Finally, under the stronger
assumption that the data is analytic in $\overline{\Omega}$, and without
compatibility at $\partial \Omega$, we establish exponential rates of
convergence of $hp$-FEM for spectral, fractional diffusion operators. We also
report numerical experiments for model problems which confirm the theoretical
results. We indicate several extensions and generalizations of the proposed
methods.
