Analytical expression for step response of digital bessel filter
Yes, there is an analytical response, though it might not look like you would like to have it:

Your Bessel Filter can be written in Z-Domain as 

$$H(z)=K\frac{\prod_{i=0}^7(z-z_i)}{\prod_{i=0}^7(z-p_i)}$$

where $z_i$ are the zeros and $p_i$ are the poles of the z-Transform. Now, to get the step response, you "simply" need to calculate 

$$h[n]=\mathcal{Z}^{-1}\left(\frac{z}{z-1}H(z)\right).$$

i.e. you multiply the Z-Transform of the filter with the Z-Transform of a unit-step (which is $z/(z-1)$) and calculate the inverse Z-Transform. If you have access to Mathematica etc. it can do it for you. If not, you would need to perform partial fraction decomposition of the expression and then perform inverse Z-Transform of the several partial fractions.