I found this to be an interesting question, a solenoid when supplied with a current is a form of electromagnet, and your question led me to consider a rectangular electromagnet as a counter example.
I think, if my calculations are valid, that this result is just a coincidence.
At any rate, my logic:
Define the $x$ axis to be along the principle axis of the solenoid.
Cylinder
Using the Biot-Savart Law, we can calculate the magnetic field on axis for a current loop:
$$ d\textbf{B} = \frac{\mu_{0}I}{4\pi|\textbf{r-r'}|^{3}}d\textbf{l}\times(\textbf{r-r'})
$$
We can ignore the contributions in the $y$ and $z$ directions, as they will cancel out around the loop. Integrating around the loop:
Let $\textbf{r-r'} = \textbf{R}$
$$
B_{x} = \int_{0}^{2\pi}{ \frac{\mu_{0}I}{4\pi|\textbf{R}|^{2}}}asin\theta d\theta = \frac{\mu_{0}Iasin\theta}{2|\textbf{R}|^{2}}
$$
where $\theta$ is the angle between the $x$ axis and $\textbf{R}$.
Now we consider the cyinder, which we think of as an infinite number of current loops, each with current $Indx$, where n is the density of current turns, $\frac{N}{2a}$. Integrating this from one end of the finite cylinder to the other, changing the integral to be in terms of $\theta$ for convenience:
$$
B_{x} = \int_{-a}^{a}{\frac{\mu_{0}Inasin\theta}{2R^{2}}dx}
$$
$R = \frac{a}{sin\theta}$, $x = -Rcos\theta$ => $dx=\frac{R^{2}}{a}d\theta$
$$
B_{x} = \int_{0}^{\pi}{\frac{\mu_{0}Insin\theta}{2}d\theta} = \mu_{0}In = \frac{\mu_{0}IN}{2a}
$$
This corresponds to your first result, as the field is along the x axis.
Spherical
This is a similar calculation, but the radius of the current loop changes with x, while the vector $\textbf{R}$ has a constant magnitude $a$.
Conducting the same finite integral as before, but changing to suit the spherical system, $\rho$ is the radius of the current loop:
$$
B_{x} = \int_{-a}^{a}{\frac{\mu_{0}In\rho sin\theta}{2a^{2}}dx}
$$
$\rho = asin\theta$, $x = -acos\theta$, => $dx = asin\theta d\theta$
$$
B_{x} = \int_{0}^{\pi}{\frac{\mu_{0}Insin^{3}\theta}{2}d\theta} = \frac{2\mu_{0}n}{3} = \frac{\mu_{0}IN}{3a}
$$
Your second result.
Square
Consider a square current loop, of side length $2a$. By using the Biot-Savart law from above and some nifty cancellations, we can calculate the on axis field for one side of the loop. A similar procedure to the loop gives the x component of the field contributed by an infinitesimal length on the side as:
$$
dB_{x} = \frac{a\mu_{0}I}{4\pi R^{3}}dl
$$
If we consider the edge to be aligned with the $y$ direction, and expand $R$ before integrating across the edge:
$$
B_{xside} = \int_{-a}^{a}\frac{a\mu_{0}I}{4\pi(x^2+y^2+a^2)^{\frac{3}{2}}}dy = \frac{a^2\mu_{0}I}{2\pi(2a^{2}+x^{2})^{\frac{1}{2}}(a^2+x^2)}
$$
From the symmetry of the situation, it can be seen that each side of the square current loop contributes equally, so the total magnetic field in the $x$ direction is 4 times the prior value.
Again, we consider the infinite number of current loops:
$$
B_{x} = \int_{-a}^{a}\frac{2a^2\mu_{0}nI}{\pi(2a^{2}+x^{2})^{\frac{1}{2}}(a^2+x^2)}dx = \frac{2a^{2}\mu_{0}nI}{\pi}[\frac{tan^{-1}(\frac{x}{(2a^{2}+x^{2})^\frac{1}{2}})}{a^{2}}]_{-a}^{a} = \frac{2\mu_{0}In}{3} = \frac{\mu_{0}IN}{3}
$$
Now, comparing the volumes, the volume of this cylinder is $8a^3$, which is $\frac{4}{\pi}$ times the cylindrical volume, but the field is $\frac{2a}{3}$ times the cylindrical field.
So, the neat consistency between the fields and volumes for the spherical and cylindrical case doesn't hold for a rectangular electromagnet. Unfortunately, I currently can't give you any ideas as to how to approach this from a more theoretical and thought based approach, rather than a bunch of nasty integrals, whether the relation between the spherical and cylindrical systems is due to the similarities between the two coordinate systems, or is something more fundamental.
Any, hope this answer was at least interesting to you, as I'm not sure it's too much help, but I had fun.
Best Answer
The magnetization, $\mathbf{M}$, is related to the magnetic moment, $\mathbf{m}$, via $$ \mathbf{m}=\int\mathbf{M}dV $$ For a uniform, linear material with magnetic permeability $\mu$, we have $$ \mathbf{B}=\mu_0\left(\mathbf{H}+\mathbf{M}\right)=\mu_0\left(\frac{\mathbf{B}}{\mu}+\mathbf{M}\right) $$ or
$$ \mathbf{M}=\left(\frac{1}{\mu_0}-\frac{1}{\mu}\right)\mathbf{B}. $$ Inserting the magnetic field of your solenoid,
$$ \int\mathbf{M}dV=NI\left(1-\frac{\mu_0}{\mu}\right)\mathbf{S}, $$ where $\mathbf{S}$ is as you defined it. This is the contribution from just the core which must be added to the moment you calculated above which is due to the current in the wires.
Note that this neglects fringe fields, variation of $\mathbf{B}$ within the solenoid, nonlinear material behavior, and changes to the magnetization due to time dependence of the fields in which you are working. note that if $\mu_0/\mu\ll1$, much of these points won't matter much. Note also that if the material is subject to hysteresis (i.e. ferromagnetic materials), $\mathbf{B}\neq\mu\mathbf{H}$, or rather, $\mathbf{M}$ cannot in general be calculated a priori. You'd have to measure it, which I suggest you do anyway to account for what these calculations may miss. Put it in a Helmholtz coil and measure the torque in different orientations and hope it doesn't change much from here to space.