  We study solution techniques for an evolution equation involving second order
derivative in time and the spectral fractional powers, of order $s \in (0,1)$,
of symmetric, coercive, linear, elliptic, second-order operators in bounded
domains $\Omega$. We realize fractional diffusion as the Dirichlet-to-Neumann
map for a nonuniformly elliptic problem posed on the semi-infinite cylinder
$\mathcal{C} = \Omega \times (0,\infty)$. We thus rewrite our evolution problem
as a quasi-stationary elliptic problem with a dynamic boundary condition and
derive space, time, and space-time regularity estimates for its solution. The
latter problem exhibits an exponential decay in the extended dimension and thus
suggests a truncation that is suitable for numerical approximation. We propose
and analyze two fully discrete schemes. The discretization in time is based on
finite difference discretization techniques: trapezoidal and leapfrog schemes.
The discretization in space relies on the tensorization of a first-degree FEM
in $\Omega$ with a suitable $hp$-FEM in the extended variable. For both schemes
we derive stability and error estimates.
