What is the correct expression for the magnetic energy density inside matter? I'll use a magnetized sphere as an example, of radius $R$, with a magnetization density $\vec{M}$.  The magnetic moment of the sphere is $\vec{\mu} = \vec{M} \, V$.  The magnetic field inside and outside is well known :
\begin{align}\tag{1}
\vec{B}_{\text{inside}} &= \frac{\mu_0 \, \vec{\mu}}{2 \pi R^3}, \\[12pt]
\vec{B}_{\text{outside}} &= \frac{\mu_0}{4 \pi} \Big( \frac{3 (\vec{\mu} \cdot \vec{r}) \vec{r}}{r^5} - \frac{\vec{\mu}}{r^3} \Big). \tag{2}
\end{align}
Then what is the correct way of calculating the total energy of that sphere ?
\begin{align}\tag{3}
U_{\text{magn}} &= \int \frac{B^2}{2 \mu_0} \, d^3 x \\[12pt]
&= \frac{1}{2 \mu_0} \int_{\mathcal{V}_{\text{inside}}} B_{\text{inside}}^2 \, d^3 x + \frac{1}{2 \mu_0} \int_{\mathcal{V}_{\text{outside}}} B_{\text{outside}}^2 \, d^3 x, \tag{4}
\end{align}
or this ?
\begin{align}\tag{5}
U_2 &= \frac{1}{2}\int \vec{B} \cdot \vec{H} \, d^3 x,
\end{align}
Or another expression involving the material properties of the sphere ?
What is the proper interpretation of this equation :
\begin{equation}\tag{6}
\int_{\text{all space}} \vec{B} \cdot \vec{H} \, d^3 x = 0,
\end{equation}
as shown in Jackson (second edition), page 207 (exercise 5.13).
 A: The energy of the magnetic field is the work required to establish a general steady-state distribution of currents and fields. This work is, in infinitesimal form, 
$$
\label{0}\tag{0}
\delta W =  \frac 1 2 \int (\delta \mathbf A \cdot \mathbf J) \ d^3 x
$$
where $\mathbf J$ is the current density.
If we are interested in work done on the free (macroscopic) currents, we have (a):
$$
\begin{align}
\delta W & =  \frac 1 2 \int (\delta \mathbf A \cdot \mathbf J_f) \ d^3 x\\
&=\int \delta \mathbf A \cdot (\nabla \times \mathbf H) \ d^3x\\
&=\int [\mathbf H \cdot (\nabla \times \delta \mathbf A) +\nabla \cdot(\mathbf H \times \delta \mathbf A)] \ d^3x
\end{align}
$$
Where $\mathbf A$ is the vector potential and $\mathbf H$ is the magnetic field (b). Assuming a localized field distribution, the second term in the integral vanishes, and using $\mathbf B = \nabla \times \mathbf A$ we get
$$
\delta W = \int \mathbf H \cdot \delta \mathbf B \ d^3x
$$
If we assume that the material is linear, i.e. that $\mathbf B = \mu \mathbf H$, we have
$$
\mathbf H \cdot \delta \mathbf B = \frac 1 2 \delta( \mathbf H \cdot \mathbf B)
$$
Therefore we finally get the following expression for the energy of the magnetic field in the presence of linear materials:
$$
U = \frac 1 2 \int \mathbf H \cdot \mathbf B \ d^3 x = \frac 1 {2 \mu} \int B^2 \ d^3 x
$$
where the magnetic energy density is written in the form
$$
\tag{1}\label{1}
u = \frac 1 2 \mathbf H \cdot \mathbf B = \frac {B^2} {2 \mu}  
$$
To derive \ref{1}, we made use of the macroscopic form of Maxwell's equations. In particular, we are assuming using for the 4th equation the form (neglecting the displacement current):
$$
\nabla \times \mathbf H =  \mathbf J_f
$$
where $\mathbf J_f$ is the free (macroscopic) current.
This means that \ref{1} represents the work done on free currents when establishing the steady-state current distribution of currents and fields.
It would also be possible to use the microscopic form of Maxwell's equations, in particular
$$
\nabla \times \mathbf B = \mu_0  \mathbf J
$$
where $\mathbf J$ is the total current, i.e. the sum of free and bound currents:
$$\mathbf J = \mathbf J_f + \mathbf J_b$$
In this case, there is no $\mathbf H$ vector and the magnetic energy density is (c)
$$
\tag{2}\label{2}
u' = \frac{B^2}{ 2 \mu_0}
$$
This represents the work done on every current when establishing the magnetic field, including the bound currents, i.e.
$$
\begin{split}
u & =  u_f\\
u'& =  u_f + u_b\\
\end{split}
$$
This means that the energy required to establish the bound currents can be calculated for a linear material as
$$
\tag{3}\label{3}
u_b = u'-u = \frac{B^2}2 \left( \frac 1 {\mu_0} - \frac 1 {\mu} \right) = \frac{B^2} 2 \left( \frac{\mu-\mu_0}{\mu \mu_0} \right)
$$
Since (for a linear material) we have
$$
\mathbf M = \left( \frac{\mu-\mu_0}{\mu \mu_0} \right) \mathbf B
$$
where $\mathbf M$ is the magnetization, we can write \ref{3} as
$$
\tag{4}\label{4}
u_b = \frac 1 2 \mathbf M \cdot \mathbf B
$$
And indeed this expression is valid for every material, not just for linear ones (d).

As for the expression
$$
\int \mathbf H \cdot \mathbf B \ d^3 x = 0
$$
on Jackson's book the full text says:

A magnetostatic field is due entrirely to a localized distribution of
  permanent magnetization. Show that $$ \int \mathbf H \cdot \mathbf B\ d^3 x = 0 $$ provided the integral is taken over all space.

So this is a statement that (in the presence of linear materials) the work done on free currents when establishing a magnetostatic field is 0. I suspect this to be more general, i.e. valid for any kind of relation between $\mathbf H$ and $\mathbf B$, but at the moment I cannot prove it.


(a) J.D. Jackson, Classical Electrodynamics (1962) 6.2
(b) I adopt the nomenclature in which $\mathbf H$ is the magnetic field and $\mathbf B$ the magnetic-flux density (or magnetic induction), which is the one used by Jackson. 
(c) D.J. Griffiths, Introduction to Electrodynamics, 3rd ed. (1999), 7.2.4 and 8.1.2. It is especially interesting to read the footnote at page 348.
(d) B.D. Popovic Evaluation of magnetic energy density in
magnetised matter,  PROC. IEE, Vol. 113, No. 7, July 1966.
