Magnetic energy in thermodynamics The potential energy per unit volume of a dipole in a magnetic field is given by:
$$ V = - \vec{M} \cdot \vec{H} $$
After taking the differential in $V$ we will have two terms. My question is why we only concern with one of the term ($-\vec{M} \cdot d\vec{H}$) when we represent the differential internal energy ($dU = T dS -\vec{M} \cdot d\vec{H}$) in thermodynamics, especially when we are dealing with magnetic cooling in paramagnetism?
 A: Actually, the infinitesimal magnetic work is $\delta w = HdB$ and when substituting $B=\mu_0 (H+M)$ it can be rewritten as  $\delta w = \mu_0 HdH + \mu_0 HdM$. The first term is a total differential, hence it is usually ignored when the energy change of the material is investigated. The 2nd term can be written as $HdM = d(HM)-MdH$ so the internal energy change $dU=TdS+\delta w$ can be written as $dU = TdS+\mu_0 HdH + \mu_0 HdM$ and then $$dU'=TdS -\mu_0 MdH $$ where I write $U' = U-\frac{1}{2}\mu_o HdH -\mu_0 HM$ 
Now with some more manipulation of the energy integral, see Stratton, one can show that in the term $-MdH$ the actual field $H$ can replaced with the external field of the biasing source (coil or magnet) $H_0$ that would be when the magnetic material is absent, that is $$dU^* = TdS - \mu_0 MdH_0$$ is a total differential and $U^*$ is a kind of magnetic enthalpy, details of the derivation using Stratton's formula are in Heine:"The Thermodynamics of bodies in static electromagnetic fields". 
A: In thermodynamics, when dealing with magnetic systems, $\vec{\bf M}$ is not the magnetic momentum of an isolated dipole but the average over the volume of the sample of the microscopic  magnetic moment. This has two main consequences. The first is that in general the magnetization density is not constant but it  a function of the point $\vec{\bf M}(\bf{r})$. The second is that it is not independent on the magnetic field. The two things together imply that, again in the general case, the magnetic energy of the system is not $-\vec{\bf M} \cdot \vec{\bf H}$. Such a statement it is not astonishing. Let's think for example to the gravitational potential energy. Only in the very special case of constant gravitational force, we can write the gravitational energy as $\vec{\bf F}\cdot \vec{\bf L}$, where $\vec{\bf L}$ is a vector representing the displacement from a reference position. What we do in general is to start from the formula for the infinitesimal work done by at each point by the force, integrating it over a path. The same procedure is required for the magnetic work to be used in the expression of the first principle of thermodynamics. Therefore, the starting point is
$$
dw_{magnetic} = \int d{\bf r} ~\vec{\bf H}({\bf r})\cdot \delta\vec{\bf B}({\bf r})
$$ 
where the space integration is intended over the whole volume where the fields $\vec{\bf H}$ and $\delta \vec{\bf B}$ are simultaneously different from zero. In general, such domain is not convenient because it is wider than the volume of the sample.
Notice that, taking into account the way such expression is derived,  $\delta \vec{\bf B}({\bf r})$ is the variation of the field $\vec{\bf B}({\bf r})$ and such a field has to be intended as the independent variable. The field $\vec{\bf H}({\bf r})$ instead, has to be considered as a function (well, in general a functional) of $\vec{\bf B}$.
It may sound weird, but it is not different from the the fact that in a fluid system from the expression of the work as $-pdV$ alone, one cannot conclude that $-pV$ is the contribution of the pressure to the internal energy. In that case again, $p$ has to be considered as a function of $V$ and the differential $-p(V)dV$ has to be integrated over $V$.
This clarification,  is basically answering a question close to the original: why cannot we write the internal enegy as  $\vec{\bf H} \cdot \vec{\bf B}$?
The case of the original question requires some additional work to rewrite the differential work 
$$
dw_{magnetic} = \int d{\bf r} ~\vec{\bf H}({\bf r})\cdot \delta\vec{\bf B}({\bf r})
$$ 
in an equivalent form 
$$
dw_{magnetic} = \int d{\bf r} ~\vec{\bf H_0}({\bf r})\cdot \delta\vec{\bf B_0}({\bf r}) + \int d{\bf r} ~\vec{\bf M}({\bf r})\cdot \delta\vec{\bf B_0}({\bf r})
$$
where the fields with subscript are the fields originating from the same external sources  but in the absence of the material sample. 
Notice that in general $\vec{\bf H}({\bf r}) \neq \vec{\bf H_0}({\bf r})$ and $\vec{\bf B}({\bf r}) \neq \vec{\bf B_0}({\bf r})$. However the last formula for $dw_{magnetic}$ is much more convenient for two reasons. The first is because the first integral has to be evaluated only once (it does not depend on the material system) and the second integral is over the volume of the sample. The second is that by using a special geometry, the external fields in absence of the sample can be taken as uniform. Therefore the integrals can be simplified a lot.
A final word about the role of independent variables. In the final formula, the magnetization density $\vec{\bf M}({\bf r})$ should be considered as a function of the magnetic field $\vec{\bf B_0}$ in the absence of the sample.
