Запустив программу, с помощью команды 7 находим $$ L = 100;\qquad T = 10.0; \qquad h =0. $$
На каждом шаге алгоритма совершается $L^2$ операций. Поэтому время пропорционально $\tau \sim I L^2$. При $L = 100$ $I = 100$ операций выполняется за 2 минуты. Тогда в условиях задачи $$ \tau = 2~мин \times 5 \times 10^2 = 1000~мин \approx 17~часов. $$
При больших температурах магнитные моменты расположены хаотично. При низких температурах практически все моменты направлены в одну сторону. В точке фазового перехода соседние моменты начинают группироваться в "капли" сонаправленных магнитных моментов. Это происходит при температуре порядка $2.5$.
Запустив моделирование и посмотрев на график намагниченности от числа шагов, находим, что намагниченнось сходится к среднему значению за $I \sim 50 $ шагов.
Вся погрешность определения $m^2$ возникает из-за случайных отклонений. Это означает, что погрешность спадает с числом усредняемых слагаемых как $1/\sqrt{I}$. Тогда требуемое число шагов можно оценить как $I\sim 10^4$. Используя такое число шагов находим $m^2 = 0.205 \pm 0.010$. Запустив программу насколько раз, убеждаемся, что наша оценка погрешности действительно справедлива.
При $T \to 0$ практически все моменты сонаправлены и $\langle m^4 \rangle = \langle m^2 \rangle =1$. Тогда $U_0 = 2/3$. При больших температурах направления всех магнитных моментов практически независимы. Поэтому намагниченность представляет собой сумму большого числа независимых величин, а значит распределена по Гауссу. Тогда $U_\infty = 0$. Эти же результаты можно получить из численного моделирования при достаточно малых и достаточно больших температурах.
Снимем зависимость $m^2 (T)$, используя для усреднения $I = 10^4$ шагов алгоритма в интервале температур $T \in [1.5,\,4.5]$. Получим график
Построим график зависимости параметра Биндера от температуры для двух различных размеров решетки, например $L =5$ и $L = 10$. Вблизи точки фазового перехода эти зависимости можно аппроксимировать линейными. Находя точку пересечения соответствующих прямых, получаем $T_c \approx 2.21$. Для того, чтобы добиться хорошей точности, каждую точку требуется измерять, усредняя по $I \sim 10^4$ итерациям алгоритма. У этой задачи также есть точное аналитическое решение. Из этого решения следует $T_c = 2/ \ln(1+ \sqrt{2}) \approx 2.23$ .
Будем рассматривать температуры $T>3$. Для измерения магнитной восприимчивости будем включать небольшое магнитное поле $h \sim 0.1$. Величина поля должна быть не слишком большой, чтобы зависимость намагниченности линейно зависела от поля. С другой стороны, она должна быть не слишком маленькой, чтобы изменение намагниченности превышало погрешность, с которой она измеряется. Тогда магнитную восприимчивость можно определить как $\chi = \Delta m/h$. Построив график $1/\chi(T)$ получим линейную зависимость. Таким образом $$ \chi (T) = \frac{\alpha}{T-T_0}. $$ Константы $\alpha$ и $T_0$ зависят от выбранного размера решетки.
Возьмем $L = 100$, $T = T_c = 2.23$, $h=1$. Проделаем $I = 200$ итераций алгоритма, чтобы система пришла в тепловое равновесие. Снимем зависимость $C(l)$ при $l \in [1,\,30]$ и построим график $\ln C$ от $l$. Из графика находим, что зависимость приближенно линейна при $l \in [3, \,15]$, из наклона графика находим корреляционную длину $\chi \approx 8$. Отметим, что значение корреляционное длины сильно зависит от температуры.