I have a question relating the internal energy U with the average of the Hamiltonian of a system, when external fields are involved. My example here is probably not physical and should only be seen as a Toy model, which shows my problems.
To my current understanding the first law of thermodynamics demands that the internal energy is a function of extensive quantities only, therefore not of an external field.
Let us consider electrons in a magnetic field which can interchange their spin. The Hamiltonian is given by
$$ H = \sum_p \left[ \epsilon_{p,\uparrow} c^\dagger_{p,\uparrow} c_{p,\uparrow} + \epsilon_{p,\downarrow} c^\dagger_{p,\downarrow} c_{p,\downarrow} + t( c^\dagger_{p,\uparrow} c_{p,\downarrow} +h.c.) \right] \qquad \text{with} \qquad \epsilon_{p,\sigma} = \frac{p^2}{2m} - \mu - \sigma B $$
So the chemical potential $\mu$ couples to the total particle number and is conserved. The magnetic field $B$ couples to the spin difference, which is not conserved.
The Hamiltonian can be diagonalised and in the new basis we find $$ H = \lambda_{p,+} n_{p,+} + \lambda_{p,-} n_{p,-} \qquad \lambda_{\sigma} = \frac{p^2}{2m} - \mu \pm \Delta \qquad \Delta=\sqrt{B^2+t^2} $$
From here on we can start calculating thermodynamics quantities in the new basis. The grand canonical potential is given by
$$ \Phi = - k_B T \sum_{\sigma} \int d\epsilon \rho(\epsilon) \ln\left(1+e^{-\beta(\epsilon -\mu + \sigma \Delta)} \right) $$
The magnetisation can now be calculated as $$ M = - \left( \frac{\partial \Phi}{\partial B}\right)_{\mu,T} \propto \frac{B}{\Delta} $$
Now my question is should we consider $\Phi$ as a thermodynamic function of $B$ or should we take $B$ just as a parameter?
If we take it as a thermodynamic function of $B$ and $U$ as a function of extensive quantities we can obtain $U$ with a Legendre transformation from $\Phi$. So we can calculate $U$ as
$$ U = \Phi + TS + \mu N + B M \propto \sum_{\sigma=\pm} \int \frac{ d\epsilon \rho(\epsilon)}{e^{\beta( \epsilon - \mu + \sigma \Delta)}+1} \left( \epsilon + \sigma \frac{t^2}{\Delta} \right) $$
This results is not the average of the overall Hamiltonian or the Hamiltonian of internal dynamics.
Further calculating the specific heat from the entropy suggests that $U$ should be the average of the given Hamiltonian.
So therefore a set of questions:
Should one really consider $U$ as a function of extensive quantities only or include external fields?
Should one take the fields only as parameters, which determine correlation function?
Is the internal energy generally at all related to the average of the Hamiltonian?
Is this problem particular, because the magnetization is not conserved?
EDIT:
The answer of tepsilon let me revisit the construction of the density matrix by the method of Lagrange multipliers.
Let the system without external coupling be described by the Hamiltonian $H$, in the Hamiltonian above $B=0$, $\mu=0$.
Let us try to impose that the operator $M$ has a certain average with $[M,H] \neq 0$.
We introduce a Lagrange function with multipliers to also enforce a certain average energy.
$$ L = Tr( \rho \log(\rho)) + \lambda (1-Tr(\rho)) + \alpha (M - Tr( \rho M) ) + \beta (E - Tr(\rho H )) $$
To get the density matrix we take the variation w.r.t. to the matrix entries of $\rho$ and use that we know it is diagonal in the eigenbasis of the Hamiltonian.
$$ \delta_{\rho_n} L = \delta \rho_n \left[ \ln(\rho_n) +1 -\lambda - \alpha M_{n,n} - \beta H_{n,n}\right] =0 $$
We used here that
$$ \delta_{\rho_k} Tr ( \rho M) = \delta_{\rho_k} \sum_n <n| \rho M |n> = \delta_{\rho_k} \sum_n \rho_n M_{n,n} = M_{k,k} \delta \rho_k $$
So the density matrix will be
$$ \rho_n \propto e^{\beta H_{n,n} + \alpha M_{n,n}} $$
Therefore the density matrix will not contain the operator $M$, but only the diagonal elements in the eigenbasis of the Hamiltonian.
$$ \rho \propto e^{\alpha \bar{M}} \quad \text{with} \quad \bar{M} = \sum_n <n| M |n> |n> <n| $$
The operators $M$ and $\bar{M}$ have the same average wrt. to this density matrix.
This coupling seems to be physically different from the direct coupling of the spin to the magnetic field. So i would suspect that generally the dynamics are different.