The numerical treatment of fractional differential operators is challenging due to their non-local nature. Additional numerical challenges arise in particular in the case of bounded domains from strong singularities of the solution at the boundary. In this talk we present recent results for high order finite element discretizations (hp-FEMs) of the spectral fractional Laplacian in bounded domains, in particular on polygonal domains. In this situation the solution has strong boundary singularities as well as corner singularities. We present mesh design principles that are based on geometric refinement towards the corners and anisotropic geometric refinement towards the boundary. We show that hp-FEM on such meshes can deliver exponential convergence. We discuss in more detail two high order discretization schemes. The first one is based on the Caffarelli-Silvestre extension, which realizes the non-local fractional Laplacian as a Dirichlet-to-Neumann map of a (degenerate) elliptic boundary value problem (BVP). This BVP is amenable to a discretization by high order finite element method (hp-FEM). The second discretization is based on the so-called Balakrishnan'' formula, an integral representation of the inverse of the spectral fractional Laplacian. The discretization of the integral leads to a collection of BVPs, which can be discretization by hp-FEM. For both discretization schemes, exponential convergence of hp-FEM is established. Extensions to time-dependent problems will also be given.