I have a question on the terminology being used to compute results on the Ising Model. Every explanation I've seen defines the magnetisation of the system as some variation of:
$$ M = \frac{1}{N}\sum_i \langle\sigma_i\rangle = \frac{\langle \sum_i \sigma_i \rangle}{N} $$
Where $\langle . \rangle$ is the ensemble average. This formulation is then used to derive important results, such as the mean field approximations that derives a self-consistency:
$$ M = \tanh(\beta(h+qJM)) $$
Where $\beta$ is the inverse temperature, $h$ the external field, $q$ is the number of neighbours per site, $J$ the coupling constant. (Taken from here).
However this is leading to two questions for me which I haven't seen an obvious answer for:
How is $M$ ever non-zero for $h=0$? By the symmetry of the possible Hamiltonians, any positive value has an equal weight that is negative. The terms will always cancel, so surely the ensemble average must be zero.
Relatedly, in the ferromagnetic phase ($T<T_c$), how does one interpret the non-zero steady state solutions? When such results are discussed as observed magnetisations I don't quite follow as $M$ is defined as an ensemble average so it is never "observed" in that sense. Or are they referring to the a specific configuration quantity $M = \frac{\sum_i \sigma_i}{N}$ and the terminology is just a bit inconsistent?
I should note I'm not looking for a strict physical interpretation (it makes perfect sense, as the answers for this question lay out), my confusion is with the precise mathematics of the situation. I'm finding it uncomfortable to manipulate $M$ as the central quantity in such derivations when I don't have a clear idea of what it is consistently referring to. For the purposes of calculation, what is the appropriate definition of $M$?