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?