В общем, сделал я вычисления для функции Ландау и среднего спина в случае четырёх спинов в кластере. Здесь можно посмотреть, какие состояния и с какими энергиями при этом рассматриваются. Здесь - программа и результат вычислений. Если хотите, можете сами проделать эти вычисления, текст программы для Математики внизу.

Что можно отметить. В случае аналитического решения Вакса критическая температура практически совпадает с параметром ε и практически не зависит от энергии заряженных конфигураций ω, которая, как пишет Вакс, раз 8-10 больше, чем ε. Когда я находил восприимчивость из аналитического решения Вакса у меня получалось, что экспериментальным даным лучьше соответсвует ω=3ε, но надо перепроверить. Здесь же получается, что температура перехода сильно зависит от обоих параметров и не совпадает с ε. Может, где-то в формулах ошибся. Опять же надо всё перепроверить.

Что касается, второй производной от функции Ландау по спину, которая является обратным временем релаксации. А не равна ли она просто обратной воприимчивости? Это видно, насколько я понимаю, если продифференцировать функцию Ландау сначала по спину, потом по внешнему полю. Поэтому эта величина может быть получена и из функции, которую пишет Вакс. Другое дело, что самосогласование Вакса несколько отличается от того, что используется при прямом счёте, если правильно понимаю.

В случае двух спинов, я вычислял обратное время релаксации и оно действительно почти совпадало (качественно, с точностью до множителя) с обратной восприимчивостью.


(* Средние значения четырёх спинов кластера, выраженные через \
вероятности 16 возможных состояний кластера p[i] *)
s1 = p[1] - p[2] - p[3] + p[4] + p[5] - p[6] - p[7] - p[8] + p[9] -
p[10] + p[11] + p[12] - p[13] + p[14] - p[15] + p[16];
s2 = p[1] - p[2] + p[3] - p[4] + p[5] - p[6] + p[7] - p[8] + p[9] +
p[10] - p[11] + p[12] - p[13] - p[14] + p[15] - p[16];
s3 = p[1] - p[2] + p[3] - p[4] - p[5] + p[6] + p[7] - p[8] - p[9] -
p[10] - p[11] + p[12] + p[13] + p[14] - p[15] + p[16];
s4 = p[1] - p[2] - p[3] + p[4] - p[5] + p[6] + p[7] + p[8] + p[9] -
p[10] - p[11] - p[12] - p[13] + p[14] + p[15] - p[16];

(* Энергии взаимодействия спинов v_ij, выраженные через энергии \
отдельных состояний кластера \[Epsilon] и \[Omega] (оба выражения \
поделил на два, чтобы не делать этого в определении самой функции) *)
v12 = v23 = v34 = v41 = (\[Omega]/2 - \[Epsilon]/4)/2;
v13 = v24 = (\[Epsilon]/2 - \[Omega]/2)/2;

(* Средняя энергия кластера, нуль энергии соответствует состояниям \
p[7]-p[14] *)
f1 = -\[Omega]*Sum[p[i], {i, 1, 2}] + (\[Epsilon] - \[Omega])*
Sum[p[i], {i, 3, 6}] + (3 \[Omega] - 2 \[Epsilon])*
Sum[p[i], {i, 15, 16}];
(* Энергия взаимодействия со средним полем и внешним полем e *)
f2 = -s1*(v12*x2 + v13*x3 + v41*x4 + e) -
s2*(v12*x1 + v23*x3 + v24*x4 + e) -
s3*(v13*x1 + v23*x2 + v34*x4 + e) -
s4*(v41*x1 + v24*x2 + v34*x3 + e);
(* Слагаемое -T*S, где S - энтропия кластера *)
f3 = t*Sum[p[i] Log[p[i]], {i, 1, 16}];
(* Лагрнажевы слагаемые, соответсвующие равенству единицы сыммарной \
вероятности и определению среднего значения спина по кластеру *)
f4 = a (Sum[p[i], {i, 1, 16}] - 1) +
b (p[1] -
p[2] + (1/2) (p[7] - p[8] + p[9] - p[10] - p[11] + p[12] -
p[13] + p[14]) - s);
(* Полная функция Ландау (свободной энергии) *)
f = f1 + f2 + f3 + f4;

(* Задание энергетических постоянных, соответствующих энергии \
различных состояний кластера *)
\[Epsilon]0 = 80;
k = 3;
правило = {\[Epsilon] -> \[Epsilon]0, \[Omega] -> k*\[Epsilon]0,
x1 -> s1, x2 -> s2, x3 -> s3, x4 -> s4};

(* Определение функции Ландау от среднего спина, температуры и \
внешнего поля в точке минимума по другим параметрам: вероятностям \
нахождения кластера в различных состояниях p[i] и множителям Лагранжа \
*)
p0 = 0.1;(* вспомогательный параметр для нахождении решения \
итерационным методом (начальное значение), может влиять на результат *)
h[s1_, t1_, e1_] :=
f /. правило /. s -> s1 /. t -> t1 /. e -> e1 /.
FindRoot[Join[
Table[D[f, p[i]] == 0, {i, 1, 16}], {D[f, a] == 0,
D[f, b] == 0}] /. правило /. s -> s1 /. t -> t1 /. e -> e1,
Join[Table[{p[i], p0}, {i, 1, 16}], {{a, p0}, {b, p0}}]];


e0 = 0.0;(* Будем интересоваться нулевым внешним полем *)
tmin = 50; tmax = 300;(* Рассматриваемый интервал температур *)

(* Фнкуция Ландау от среднего спина и температуры*)
time = AbsoluteTime[];
рисунок1 =
Plot3D[h[s, t, e0], {t, tmin, tmax}, {s, -1.1, 1.1},
AxesLabel -> {"температура", "спин", "ФЛ"}, Mesh -> 25,
ColorFunction -> Hue]
Print["Время вычислений = ", time = AbsoluteTime[] - time, " с"]

(* Определение первой производной от функции Ландау по спину, иначе \
ND[expr,x,Subscript[x, 0]] *)
eps = 0.000001;
h1[s1_, t1_, e1_] := (h[s1 + eps/2, t1, e1] - h[s1 - eps/2, t1, e1])/
eps;

(* Нахождение равновесных значение среднего спина при различных \
температурах - соответсвуют минимуму функции Ландау по спину или \
нулям производной функции Ландау по спину *)
time = AbsoluteTime[];
данные3 =
Table[sp /. FindRoot[h1[sp, t, e0] == 0, {sp, 1}], {t, tmin, tmax,
5}];

(* Отображение зависимости среднего спина от температуры *)
рисунок3 =
ListPlot[Transpose[{Table[t, {t, tmin, tmax, 5}], данные3}],
AxesOrigin -> {0, 0}, AxesLabel -> {"температура", "спин"}]
Print["Время вычислений = ", time = AbsoluteTime[] - time, " с"]