In many lectures, it is stated that the maximal value $B_{\text{max}}$ of the magnetic field at the surface of a star can be found in Newton's gravitation theory by equating the gravitational potential energy with the magnetic field energy. For a sphere of mass $M$ and radius $R$, of uniform density and uniformly magnetized : \begin{equation}\tag{1} |\, U_{\text{grav}}| = \frac{3 G M^2}{5 R} = E_{\text{magn}} = \frac{\mu_0 \, \mu^2}{4 \pi R^3}, \end{equation} where $\mu$ is the sphere's dipolar magnetic moment. The right hand part is the total energy stored in the magnetic field of a dipole : \begin{align} E_{\text{magn}} = \int \frac{B^2}{2 \mu_0} \, d^3 x &= \int_0^R \frac{B_{\text{int}}^2}{2 \mu_0} \, 4 \pi r^2 \, dr + \int_R^{\infty} \frac{B_{\text{ext}}^2(r, \theta)}{2 \mu_0} \, r^2 \, dr \, \sin{\theta} \, d\theta \, d\varphi \\[12pt] &= \frac{\mu_0 \, \mu^2}{4 \pi R^3}. \tag{2} \end{align} Since the sphere is uniformly magnetized in its volume, the internal magnetic field is a constant : \begin{equation}\tag{3} B_{\text{int}} = \frac{2}{3} \, \mu_0 \, M = \frac{\mu_0 \, \mu}{2 \pi R^3} \quad \Rightarrow \quad \mu = \frac{2 \pi B_{\text{int}} \, R^3}{\mu_0}. \end{equation} Substituting this magnetic moment into equ. (1) gives the maximal field strength inside and at the surface of the sphere : \begin{equation}\tag{4} B_{\text{int max}} = \sqrt{\frac{3 \mu_0 \, G}{5 \pi}} \, \frac{M}{R^2}. \end{equation} So for a star of mass $M = 0.6 \, M_{\odot}$ and radius $R = 10^4 \, \mathrm{km}$ (a typical white dwarf), this give \begin{equation} B_{\text{int max}} \approx 5 \times 10^7 \, \mathrm{tesla} = 5 \times 10^{11} \, \mathrm{gauss}. \end{equation}
But how can we justify equation (1) ? Can it be made more rigorous ? Why should we have $E_{\text{magn}} + U_{\text{grav}} = 0$ for the maximal field strength ?
EDIT : In the case of a canonical neutron star of radius $R \approx 10 \, \mathrm{km}$ and mass $M \approx 1,44 \, M_{\odot}$, equation (4) gives \begin{equation}\tag{5} B_{\text{int max NS}} \approx 10^{14} \, \mathrm{tesla} = 10^{18} \, \mathrm{gauss}, \end{equation} which is pretty excessive (AFAIK). The strongest known magnetars have at most a field of about $10^{15} \, \mathrm{gauss}$. So is there a theoretical way in reducing the value (5) ?