First of all, equilibration can take a ``long'' time. Second, what the Mermin-Wagner theorem implies is that there's no spontaneous magnetization at any finite temperature, no matter how low, for the Heisenberg ferromagnet; it doesn't exclude the appearance of the Berezinskii-Kosterlitz-Thouless transition, for the O(2) Heisenberg ferromagnet, since, in both phases, the magnetization vanishes and the order parameter is non-local, not local. For the O(N) Heisenberg ferromagnet, for N>2, no such order parameter can be defined, that's why there's no Berezinkii-Kosterlitz-Thouless transition, when N>2. This doesn't imply that vortices can't appear, whose lifetime can be very ``long'' in Monte Carlo time, especially at ``low'' temperatures. Hence the usefulness of knowing the theory, before doing the simulations.

So the ground state configuration of the D=2, O(N) Heisenberg ferromagnet, when N>2, is characterized by the correlation functions of the spins as = 0, a=1,2,...,N, ~δabexp(-lx-y|m), with m finite for any finite temperature-there's a mass gap for any finite temperature. The other correlation functions satisfy, as for any field theory, the so-called Schwinger-Dyson identities.

In fact the same statements hold, also, for N=2, except, in that case, there's a phase transition at the Berezinskii-Kosterlitz-Thouless temperature, below which the vortices, which are non-local in terms of the spins, have power-law correlation functions, with exponent that is a function of the temperature.