There are two problems at low temperatures. One is intrinsic to the Ising model, the other to the Metropolis algorithm. The Ising one is a bit more serious.
Problem with the Ising Model at low temperatures
The Ising model of ferromagnetism posits that the electrons in a model can be modeled as little magnetic dipoles, confined to only two quantized spin states, up and down, and that all the electrons in the metal have the same two possible orientations. The Metropolis algorithm counts up the quanta of energy to be gained if a particular dipole were flipped over by counting up the states of the dipole's neighbors and comparing it to the dipole itself. It compares the number it gets with the thermal energy from the Boltzmann distribution to determine whether or not the metal is cool enough that flipping to join the neighbors' alignments is energetically favorable.
Even if the dipoles have more possible orientations than this (ie if their spins don't really all align on the same plane), that's not the big problem at low temperatures. The real issue is that quantum mechanics is more subtle than this. In reality, it isn't the orientation of each individual which is quantized, but the sum of their magnetic moments which is quantized. At very low temperatures, this becomes an important failure in the Ising model, and the relevant states of a true ferromagnet are "magnons" in which all dipoles are nearly parallel and a single unit of opposite alignment is spread across many dipoles. So, we really shouldn't be comparing the dipole to its neighbors individually at very low temperatures.
Problem with the Metropolis algorithm at low temperatures
The problem here is also subtle. The Metropolis algorithm works by randomly picking a single dipole out of the whole magnet, checking it as described above, then checking another - it's just a Monte Carlo method. Now, theoretically, there should be a "critical temperature," below which it is energetically favorable for the dipoles to side with magnetism instead of temperature. But check out the following plots, which came from a simulation I did of a Metropolis-algorithm Ising magnet in C a while ago. In the units I used, the critical temperature happened to be $T_c=2.27$:
The lowest temperature systems converge much more quickly to a homogeneous magnetization, but even at $T=1.5$ you still see "meta-stable states." In these cases, large sections of the magnet are commonly aligned, but the entire magnet isn't. Change can only occur along the border region - and since you are using a random choice of dipoles, the chance that you randomly select one of those dipoles on the border region is pretty small, so you need a very long timeframe for the system to reach equilibrium using this Monte-Carlo scheme.
Since the efficiency is worse close to the critical point of temperature, you could say that the algorithm's efficiency is worse near the critical point of the acceptance ratio, since the acceptance is related to the temperature.
You can think of the whole thing as a "fluid of systems" and each one of them can be in any of the states $i$ available.
$\pi_{ij}$ tells you what is the speed at which a system in state $i$ will go to state $j$
$p_i$ tells you how likely it is for a system (in this fluid or ensemble of systems) to be in state $i$ and is proportional to the number of systems in this state. If your system is in contact with a thermostat $p_i \propto \exp(-\beta E_i)$
Now, the total flux that goes from $i$ to $j$ is the velocity to go from $i$ to $j$ for a single system times the number of systems in state $i$; hence $p_i \pi_{ij}$. Conversely, the total flux to go from $j$ to $i$ is $p_j \pi_{ji}$.
- The detailed balance criterion (which is a sufficient criterion for equilibrium) tells you that the flux from $i$ to $j$ has to be balanced by the flux from $j$ to $i$ so that state $i$ (stationary) population cannot suffer any loss/gain by moving to to state $j$ as it is exactly balanced by systems from $j$ moving to $i$ and vice versa.
Best Answer
What you are observing is the critical slowdown - getting to the equilibrium near to the critical temperature is very hard because of huge fluctuations in the system.
Therefore, as you are not averaging over enough runs, it is no wonder that you are getting very large deviations from the (expected) average values. This problem is less severe for more sophisticated algorithms, e.g. the Wolff algorithm, as they allow for large (for example, nonlocal) differences of states in adjacent time steps.
There are many references on critical slowdown online; Googling for "critical slowdown ising" returns videos of lectures, lecture notes, as well as research papers.