We study the use of the hybridizable discontinuous Galerkin (HDG) method for numerically solving fractional diffusion equations of order (Formula presented.) with (Formula presented.). For exact time-marching, we derive optimal algebraic error estimates assuming that the exact solution is sufficiently regular. Thus, if for each time (Formula presented.)∈[0,T] the approximations are taken to be piecewise polynomials of degree (Formula presented.) on the spatial domain (Formula presented.), the approximations to u in the (Formula presented.)-norm and to ∇u in the (Formula presented.)-norm are proven to converge with the rate (Formula presented.), where h is the maximum diameter of the elements of the mesh. Moreover, for k≥1 and quasi-uniform meshes, we obtain a superconvergence result which allows us to compute, in an elementwise manner, a new approximation for u converging with a rate of (Formula presented.).