Phase transition in a 2D classical Heisenberg model Monte Carlo simulation

Question

I have performed a Monte Carlo simulation of the classical Heisenberg model for a simple 2D square lattice. But I obtain some strange results. In fact, there is clearly a phase transition and a magnetized phase for low temperatures.

I attach some plots of the results.

I know that the symmetry can be broken only for $T=0$. I've checked multiple times the code I wrote and all seems ok. I have run some simulation of the 3D Heisenberg model and all went fine. I have also estimated some critical exponents for the 3D model with good accuracy so I don't think my code is bugged. What could be the reason of this behavior?

I have already read Phase transition in 2D Heisenberg model and I use the correct way to sample point on the sphere.

EDIT:

For each sample, I compute and save the average magnetization per spin $\vec m_i$ and the system energy $E_i$

In the analysis, I simply compute the mean of the magnitude $\langle | \vec m | \rangle$ for the first plot.

In the second one, the heat capacity per spin is estimated by

$$ c_v = \frac{1}{L^2} \frac{\langle E^2 \rangle - \langle E \rangle^2}{T^2} $$

Answer 1

This model (the classical Heisenberg model, with 3-component spins, on a 2D $L\times L$ square lattice) is one of the challenging ones. I believe that the current thinking is that there is not a true phase transition, but instead it exhibits pseudocritical behaviour: the correlation length becomes extremely large, but not infinite. The latest publication that I found is by Y Tomita, Phys Rev E, 90, 032109 (2014), and there a finite-size-scaling analysis was employed on systems up to $1024\times1024$. Also, the Monte Carlo simulations employed special techniques (Swendsen-Wang cluster updating). I'm afraid that I couldn't find an open-access version of the paper cited above.

---

[EDIT following OP comment]

In the paper I cited, the XY model is also studied, and compared with the Heisenberg model. The conclusion seems to be that they are different, but it is very difficult to show the difference unless considerable computational effort is spent, and a complicated finite size analysis is conducted. In both cases a peak in the heat capacity curve is seen, and a "critical temperature" can be identified. So, there will appear to be a phase transition. The next question is whether this is a true transition, with thermodynamic properties and the correlation length scaling in a systematic way (with critical exponents) as the system size becomes very large. For the XY model, this seems to be the case; for the Heisenberg model, it seems not.

On the face of it, while a BKT phase transition for the XY model can be explained in terms of topological defects in the configurations of spins confined to the plane, the Heisenberg model doesn't have these well-defined defects. So a BKT transition is not expected for the Heisenberg model. Nonetheless, the system has been known for years to exhibit behaviour looking like a phase transition. It has been an area of active research, trying to prove it one way or the other. This paper I found is just the latest (I think) of several papers on the topic.