We analyze an expansion of the generalized block Krylov subspace framework of [Electron.\ Trans.\ Numer.\ Anal., 47 (2017), pp. 100-126]. This expansion allows the use of low-rank modifications of the matrix projected onto the block Krylov subspace and contains, as special cases, the block GMRES method and the new block Radau-Arnoldi method. Within this general setting, we present results that extend the interpolation property from the non-block case to a matrix polynomial interpolation property for the block case, and we relate the eigenvalues of the projected matrix to the latent roots of these matrix polynomials. Some convergence results for these modified block FOM methods for solving linear system are presented. We then show how {\em cospatial} residuals can be preserved in the case of families of shifted linear block systems. This result is used to derive computationally practical restarted algorithms for block Krylov approximations that compute the action of a matrix function on a set of several vectors simultaneously. We prove some convergence results and present numerical results showing that two modifications of FOM, the block harmonic and the block Radau-Arnoldi methods for matrix functions, can significantly improve the convergence behavior.