偉そうなことばっかり言って自分では何もしていない奴だと思われるのも嫌なので、僕もシミュレーションの結果を貼る。言語はC++, サンプル更新にはWolff algorithm を用いた。2次元正方格子でのIsing模型のBinder ratio のグラフ。横軸は逆温度

https://mathtod.online/media/fxqtfZepUWugM-AB6UI

各曲線の意味は、例えば4x4なら、
「4x4の格子でサンプル数100万、逆温度0.4から0.5の範囲をプロット」という意味。32x32はあまりにも時間がかかるのでサンプル数半分で妥協(数時間かかった)

グラフを見ると、確かに0.44付近で相転移していることがわかる

Binder ratio については次のトゥートで詳説

Binder ratio とは:
格子上の模型について相転移というとき、普通無限系での系の性質の著しい変化をさす。しかし当然ながら、コンピュータ上では有限系しか扱えない。そこで有限系の情報から無限系の転移点の情報を得るのに便利な値がBinder ratio である。Ising模型の場合は、磁化
\begin{align}
M=\frac1N\sum_{i=1}^N\sigma_i
\end{align}
に対して、Binder ratioを
\begin{align}
B=\frac{\langle M^4\rangle}{\langle M^2\rangle^2}
\end{align}
と定義する($\langle\cdot\rangle$はアンサンブル平均)。この値は、高温では磁化がガウス分布になるため3に近づき、低温では磁化が$\pm1$に分布するため1に近づく。理想的な無限系では臨界温度でステップ関数的に変化するが、有限系ではどうなるか?

実は、有限系であっても転移点でのBinder ratio は系の大きさに依存しないことが知られている。系を大きくするほど無限系でのグラフに近づくので、有限系の大きさを変えてBinder ratio をプロットすると、転移点で全てのグラフが交わることになる。

より詳しい説明は、
『別冊数理科学 演習形式で学ぶ相転移・臨界現象』(宮下精二・轟木義一)にある