Interpretation of Ising model simulations

If you set the inverse temperature to a value $beta>beta_{rm c}$, and let your Markov chain runs for an "infinitely long time" (or use a perfect sampling algorithm), then you obtain a configuration exactly distributed according to the corresponding finite-volume Gibbs measure (with periodic boundary condition, I assume form your pictures).

In that case, one can prove that, indeed, with large probability (going to $1$ when the system's size grow), the resulting configuration will have an average magnetization $frac1Nsum_{i} sigma_i$ very close to $+m^*(beta)$ or $-m^*(beta)$ (both alternatives occurring with probability $1/2$), where $m^*(beta)>0$ is the spontaneous magnetization.

One can also show that if you do the same with a positive field $H>0$, then the resulting configuration will have (with high probability) an average magnetization close to $m(beta,H)>m^*(beta)$. Moreover, $lim_{Hdownarrow 0} m(beta,H) = m^*(beta)$.

So the spontaneous magnetization corresponds both to the average magnetization in the $H=0$ system (up to an overall sign, which can be $+$ or $-$ with probability $1/2$), and to the limiting average magnetization as you let $Hdownarrow 0$.

All this was for the finite systems in your simulations. In infinite systems (the only ones undergoing a genuine phase transition), even after having specified the parameters (say, $beta$ and $H$ for the Ising model), there may still be infinitely many different Gibbs states that can describe it. This is due to the fact that several phases can occur for specific choices of $beta$ and $H$ (namely, when $beta>beta_{rm c}$ and $H=0$). Among all these Gibbs states, there are two that are physically most relevant; one will have an average magnetization equal to $+m^*(beta)$ with probability $1$, while the other one has an average magnetization equal to $-m^*(beta)$ with probability $1$.

The latter is completely analogous to what happens if you have a system of water molecules at 1 atmosphere and $100°$C. Then you don't know in which phase the system is from only this information: it may be liquid, it may be gaseous, or a mixture of both. You thus need different Gibbs states to describe these various possibilities.

The Hamiltonian has two global minima, corresponding to the spins all up or down, so you are equally likely to find the spins all up or all down below the Curie temperature, if you start with the checkerboard pattern and iterate the Metropolis algorithm.

It is not entirely clear from your question, but I gather this model does not include an external field in the Hamiltonian, in which case it is hardly surprising that you do see remnant magnetic field effects. Does that answer your question? If not you might want to clarify what you want to know.