We study solution techniques for an evolution equation involving second order derivative in time and the spectral fractional powers, of order s ∈ (0, 1), of symmetric, coercive, linear, elliptic, second-order operators in bounded domains Ω. We realize fractional diffusion as the Dirichlet-to-Neumann map for a nonuniformly elliptic problem posed on the semi-infinite cylinder C = Ω × (0, ∞). 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 Ω with a suitable hp-FEM in the extended variable. For both schemes we derive stability and error estimates.