  We present a novel numerical scheme to approximate the solution map $s\mapsto
u(s) := \mathcal{L}^{-s}f$ to partial differential equations involving
fractional elliptic operators. Reinterpreting $\mathcal{L}^{-s}$ as
interpolation operator allows us to derive an integral representation of $u(s)$
which includes solutions to parametrized reaction-diffusion problems. We
propose a reduced basis strategy on top of a finite element method to
approximate its integrand. Unlike prior works, we deduce the choice of
snapshots for the reduced basis procedure analytically. Avoiding further
discretization, the integral is interpreted in a spectral setting to evaluate
the surrogate directly. Its computation boils down to a matrix approximation
$L$ of the operator whose inverse is projected to a low-dimensional space,
where explicit diagonalization is feasible. The universal character of the
underlying $s$-independent reduced space allows the approximation of
$(u(s))_{s\in(0,1)}$ in its entirety. We prove exponential convergence rates
and confirm the analysis with a variety of numerical examples.
  Further improvements are proposed in the second part of this investigation to
avoid inversion of $L$. Instead, we directly project the matrix to the reduced
space, where its negative fractional power is evaluated. A numerical comparison
with the predecessor highlights its competitive performance.
