Starting with the case where $n = m + 1$, if I understand correctly you can take a point sampled from the uniform distribution on the surface of the unit sphere in $p$-space and obtain $d(x,S)^2$ by ignoring the contribution from the first $m$ coordinates. So using the representation of the uniform distribution on the sphere as a vector of $p$ standard normal random variables:
$$x = (X_1, X_2, \dots , X_p)/\sqrt{\sum\limits_{i=1}^n {X_i^2}}$$
we get
$$d(x,S)^2 = \sum\limits_{i=m+1}^n{X_i^2}/\sum\limits_{i=1}^n{X_i^2}$$
which is I think, after looking at Wikipedia for $\chi^2$ distributions, a Beta distribution with parameters determined by $m$ and $p$.
It looks like the distribution for $\sum\limits_i d(x_i,S)^2$ should then be a sum of iid Beta distributions.