Classical Heisenberg Model Using Mean Field Approximation Suppose we have the Classical Heisenberg Model of N spins $\vec S_i$ (unit vectors), external field $\vec H // \hat z $
$$ {\cal{H}} = -2J \sum_{<i,j>}  \vec S_i \cdot \vec S_j - gμ_B \sum_i \vec H \cdot \vec S_i $$
where $<i,j>$ indicates that we sum over nearest neighboors (let each site have $Z$ neighboors). In the molecular field approximation we replace the surrounding spins of each $S_i$ by the average $<\vec S>$ so the energy has the form:
$${\cal{H}}= -gμ_B \sum_i \vec H_{eff} \cdot \vec S_i = \sum_i  {\cal{H}}_i$$
where $\vec H_{eff}=\vec H  +\frac{2JZ}{gμ_B} <\vec S>$. Therefore every single spin sees an effective field and we can go ahead and write the partition function as $Z= Z_1 ^N$ with 
$$Z_1 = \int_{ |S|=1} e^{βgμ_B \vec H_{eff} \cdot \vec S } dS $$
If we calculate that we can go ahead and get the free energy, then the magnetization $\vec M = -\frac{ \partial F}{\partial {\vec H_{eff}}}$ and show that indeed it is different than $0$ below a certain temperature, like we do in the ising model.
However I'm having trouble working out the partition function. In (not very detailed) solutions I've seen that we simplify it so as to get $H_{eff}$ only in the $z$ direction, which makes things far simpler. I don't see the justification for this, so what do we do?
 A: I think that the justification for $\vec{H}_{eff}$ lying in the $z$ direction is that the system (i.e. the hamiltonian) is rotationally symmetric about the $z$ axis.
Looking at this in a bit more detail: the spin-spin part of the (full) hamiltonian is rotationally invariant (or, if you like, the spin-spin energy doesn't change if you globally rotate the spins by the same angle about a common, arbitrary, axis). This same symmetry applies to the term involving $\langle \vec{S}\rangle \cdot \vec{S}_i$ in the effective hamiltonian: if you rotate both $\langle \vec{S}\rangle$ and $\vec{S}_i$ in the same way, this term doesn't change.
So the symmetry of the  spin distribution function (which dictates $\langle \vec{S}\rangle$) must be  determined by that of the external field: cylindrically symmetric about $z$. Hence $\langle \vec{S}\rangle$ lies along $z$. Or, in other words, the $x$ and $y$ components of this vector must vanish, by symmetry.
As I read your question, you're not actually asking for the full solution, since you can proceed straightforwardly once this assumption is justified?
