library(ggplot2)
library(dplyr)
library(tidyr)
library(gridExtra)
library(viridis)

theme_set(theme_minimal())



1 Abstract

On 6 October 2026 OpenAI released 722 manuscripts produced by an unreleased internal model (github.com/openai/math). Two of the result families bear directly on this repository. Family 215 claims a proof of the exact low-temperature mass gap of the two-dimensional O(4) lattice model, \(m_{\mathrm{lat}}(\beta)\sim 32\,e^{\pi/4-1/2}\sqrt{\beta}\,e^{-\pi\beta}\), the standard toy model for asymptotic freedom and dynamical mass generation in QCD. Family 271 claims a proof of Bloch’s \(T^{3/2}\) law with its exact coefficient for the 3D quantum Heisenberg ferromagnet, together with its first lattice correction. Neither result has a Lean formalization of these specific statements, and neither has been refereed.

This document benchmarks both claims numerically in R. The benchmarks test what can be tested at finite size and finite temperature. Neither benchmark proves the theorems.

For the O(4) model, three results are obtained. First, the claimed constant is shown to equal, to machine precision, the product of the Hasenfratz–Maggiore–Niedermayer Bethe-ansatz mass gap and the lattice \(\Lambda\)-parameter ratio, so the formula itself is the long-standing prediction, now claimed as a theorem. Second, a Wolff cluster Monte Carlo, validated against Onsager’s exact Ising gap, finds a local decay rate \(-d\ln m/d\beta\) that converges toward the asymptotic value \(\pi - 1/(2\beta)\). Third, the measured mass sits roughly 11–15 % above the asymptotic formula at correlation lengths up to about 19. An offset of this size is plausible at these correlation lengths, but a single \(1/\beta\) correction fits it poorly, and the data cannot independently confirm the limiting prefactor.

For Bloch’s law, the ideal-magnon density is evaluated exactly as a rapidly convergent Bessel-function lattice sum. Its low-temperature expansion reproduces the Bloch coefficient and the manuscript’s lattice correction \(3\zeta(5/2)(\beta S)^{-5/2}/(128\pi^{3/2})\) exactly. The residuals scale with the predicted powers over two decades in temperature. An exact diagonalization of the two-magnon sector shows that the zero-momentum two-magnon level is unshifted, as SU(2) symmetry requires. The lowest scattering level shifts only as \(L^{-5}\), the finite-size signature of the vanishing low-momentum magnon–magnon interaction behind Dyson’s \(T^4\) result. This illustrates why interactions should not disturb the \(T^{3/2}\) and \(T^{5/2}\) terms, but it does not test the interacting theorem itself.

Runtime. The Monte Carlo chunks take roughly 8 minutes on one core and are cached after the first knit. Set n_bins <- 20 in the parameter chunk for a quick draft run.



2 Scope and Epistemic Status

The flags used below are as follows.

  • [Established]: textbook or long-verified literature results.
  • [Claimed: unrefereed]: statements from the OpenAI manuscripts, as yet unrefereed and, for these statements, not Lean-formalized.
  • [Verified here]: checked numerically or algebraically in this document.
  • [Not tested here]: outside the reach of these numerics.

Two distinctions matter throughout. An asymptotic statement (\(\beta\to\infty\) or \(T\to 0\)) cannot be falsified by any finite-\(\beta\) or finite-\(T\) number; a numerical benchmark can only check consistency of rates, coefficients and the size of corrections. And for Bloch’s law, the free-magnon expansion is [Established] physics (Dyson, 1956); what the manuscripts claim is a rigorous proof that the interacting quantum model obeys it. The numerics below verify the reference function the theorem is stated against, not the theorem.

The conventions were checked against the manuscripts themselves. The O(4) Gibbs weight is \(\exp(\beta\sum_{\langle xy\rangle}\sigma_x\cdot\sigma_y)\) with unit spins on \(S^3\), and \(m_{\mathrm{lat}} = -\log\lVert T_\beta|_{\Omega^\perp}\rVert\) is the gap of the full transfer operator per lattice time step. The Heisenberg Hamiltonian is \(H=-\sum_{\langle ij\rangle}\mathbf S_i\cdot\mathbf S_j\) with every edge counted once, with one-magnon dispersion \(\varepsilon_S(\mathbf k)=2S\sum_{j=1}^3(1-\cos k_j)\).

# Monte Carlo statistics. n_bins = 100 is the production setting; 20 gives a quick draft.
n_bins         <- 100
sweeps_per_bin <- 10
n_therm        <- 100



3 Part I: The O(4) Mass Gap

3.1 The Claim

[Claimed: unrefereed] For the nearest-neighbour O(4) model on the square lattice,

\[ m_{\mathrm{lat}}(\beta)\;\sim\;32\,e^{\pi/4-1/2}\,\sqrt{\beta}\,e^{-\pi\beta}\qquad(\beta\to\infty). \]

A companion manuscript in the same family proves the weaker statement that \(m_{\mathrm{lat}}(\beta)\) is bounded above and below by constant multiples of \(\sqrt\beta\,e^{-\pi\beta}\). Another proves Polyakov’s conjecture: exponential decay at every temperature for all \(n\ge 3\). Only the Polyakov paper is Lean-formalized.

3.2 Where the Constant Comes From

The formula is the combination of three ingredients.

  1. The two-loop renormalization group fixes the lattice \(\Lambda\)-parameter [Established]. With \(\beta = 1/g^2\), \(b_0=(N-2)/2\pi\) and \(b_1=(N-2)/(2\pi)^2\), this gives \(a\Lambda_L = \left(2\pi\beta/(N-2)\right)^{1/(N-2)}e^{-2\pi\beta/(N-2)}\). For \(N=4\) this is \(\sqrt{\pi\beta}\,e^{-\pi\beta}\).
  2. Lattice perturbation theory gives \(\Lambda_{\overline{\mathrm{MS}}}/\Lambda_L = 2^{5/2}e^{\pi/(2(N-2))}\) [Established].
  3. The Hasenfratz–Maggiore–Niedermayer Bethe-ansatz result gives \(m/\Lambda_{\overline{\mathrm{MS}}} = (8/e)^{1/(N-2)}/\Gamma\!\left(1+\tfrac{1}{N-2}\right)\) [Established, exact but non-rigorous].
N <- 4
m_over_LamMS    <- (8 / exp(1))^(1/(N-2)) / gamma(1 + 1/(N-2))
LamMS_over_LamL <- 2^(5/2) * exp(pi / (2*(N-2)))
# a * Lambda_L = sqrt(pi * beta) * exp(-pi * beta) for N = 4, so the sqrt(beta) e^{-pi beta}
# coefficient of the predicted mass is:
pref_hmn   <- m_over_LamMS * LamMS_over_LamL * sqrt(pi)
pref_claim <- 32 * exp(pi/4 - 1/2)

data.frame(quantity = c("m / Lambda_MSbar (Bethe ansatz)", "Lambda_MSbar / Lambda_L",
                        "HMN prefactor", "Claimed prefactor 32 e^(pi/4 - 1/2)", "Relative difference"),
           value = c(m_over_LamMS, LamMS_over_LamL, pref_hmn, pref_claim,
                     (pref_hmn - pref_claim) / pref_claim))
quantity value
m / Lambda_MSbar (Bethe ansatz) 1.935766
Lambda_MSbar / Lambda_L 12.407066
HMN prefactor 42.569331
Claimed prefactor 32 e^(pi/4 - 1/2) 42.569331
Relative difference 0.000000 
m_asym <- function(b) pref_claim * sqrt(b) * exp(-pi * b)

[Verified here] The claimed prefactor equals the HMN prediction to machine precision. The manuscript’s contribution is therefore not a new number. It claims a rigorous proof of a prediction that, until now, rested on the Bethe ansatz and perturbation theory. That is consistent with how family 215 is described.

3.3 Method: Wolff Clusters with Improved Estimators

The simulation uses the Wolff single-cluster algorithm, which suffers little critical slowing-down for O(N) models. Each update draws a random unit vector \(\mathbf r\), embeds an Ising model in the projections \(\mathbf r\cdot\sigma_x\), grows one cluster, and reflects it. Cluster growth is vectorized layer by layer, so that each frontier’s bonds are tested in one call. This keeps pure R fast enough for lattices up to \(176^2\).

Observables use the cluster (“improved”) estimators. Writing \(P_t=\sum_{x\in C,\,x_0=t}\mathbf r\cdot\sigma_x\) for the cluster’s slice sums,

\[ G(t) \;=\; N\,L\,\Big\langle \frac{1}{|C|}\sum_{t'}P_{t'}P_{t'+t}\Big\rangle_{\!W},\qquad \chi = N\Big\langle\frac{(\sum_{x\in C}\mathbf r\cdot\sigma_x)^2}{|C|}\Big\rangle_{\!W}, \]

where \(N = 4\) is the number of spin components and \(L\) the lattice side, and the subscript \(W\) denotes an average over Wolff clusters. \(F\) is defined like \(\chi\), at the smallest non-zero momentum. The mass gap is read from the effective mass \(m_{\mathrm{eff}}(t)=\operatorname{arccosh}\!\big[(G(t-1)+G(t+1))/2G(t)\big]\) of the zero-momentum time-slice correlator. Its plateau is averaged over \(t\in[\xi_2,\,2.5\,\xi_2]\). The second-moment mass \(m_2 = 1/\xi_2\), with \(\xi_2=\sqrt{\chi/F-1}\,/\,(2\sin(\pi/L))\), is reported as a cross-check.

Two physical facts make the plateau clean.

  • The spin field carries the vector representation, so the first state above the one-particle level in this channel is a three-particle state near \(3m\). Contamination therefore decays as \(e^{-2mt}\).
  • For the same reason \(\xi_{\mathrm{exp}}\) and \(\xi_2\) are known to agree closely in O(N≥3) models.

Errors are jackknife over bins of 10 sweep-equivalents. Lattices are chosen with \(L\approx 8/m_{\mathrm{asym}}\), which gives \(L/\xi\approx 9\).

make_nb <- function(L) {
  x <- rep(0:(L-1), times = L); y <- rep(0:(L-1), each = L)
  idx <- function(x, y) (x %% L) + L*(y %% L) + 1L
  cbind(idx(x+1, y), idx(x-1, y), idx(x, y+1), idx(x, y-1))
}

# One Wolff cluster update. Cluster growth is breadth-first and vectorized per frontier layer;
# each bond is tested once, from the cluster side, with p = 1 - exp(-2 beta (r.s_i)(r.s_j)).
wolff_step <- function(S, nb, beta) {
  N <- nrow(S); r <- rnorm(ncol(S)); r <- r / sqrt(sum(r^2))
  proj <- as.vector(S %*% r)
  inC <- logical(N); seed <- sample.int(N, 1L); inC[seed] <- TRUE; front <- seed
  while (length(front)) {
    i <- rep(front, 4L); j <- as.vector(nb[front, , drop = FALSE])
    k <- !inC[j]; i <- i[k]; j <- j[k]
    pp <- proj[i] * proj[j]
    act <- pp > 0
    act[act] <- runif(sum(act)) < -expm1(-2 * beta * pp[act])
    front <- unique(j[act]); inC[front] <- TRUE
  }
  id <- which(inC)
  S[id, ] <- S[id, ] - 2 * outer(proj[id], r)     # reflect the cluster through the plane normal to r
  list(S = S, id = id, proj = proj[id])
}

# O(ncomp) model on an L x L torus; ncomp = 1 is the Ising model (used as a control).
run_on <- function(L, beta, ncomp = 4, n_therm = 100, n_bins = 100, sweeps_per_bin = 10, seed = 1) {
  set.seed(seed)
  nb <- make_nb(L); N <- L*L
  xs <- (seq_len(N) - 1L) %% L + 1L; ys <- (seq_len(N) - 1L) %/% L + 1L
  S <- matrix(c(1, rep(0, ncomp - 1)), N, ncomp, byrow = TRUE)      # cold start
  flipped <- 0
  while (flipped < n_therm * N) { w <- wolff_step(S, nb, beta); S <- w$S; flipped <- flipped + length(w$id) }
  phase <- exp(2i * pi * (0:(L-1)) / L)
  Gb <- matrix(0, n_bins, L); chib <- Fb <- numeric(n_bins)
  for (b in seq_len(n_bins)) {
    flipped <- 0; nc <- 0; Gacc <- numeric(L); chiacc <- Facc <- 0
    while (flipped < sweeps_per_bin * N) {
      w <- wolff_step(S, nb, beta); S <- w$S; C <- length(w$id)
      flipped <- flipped + C; nc <- nc + 1; v <- w$proj
      Px <- Py <- numeric(L)
      ax <- rowsum(v, xs[w$id]); Px[as.integer(rownames(ax))] <- ax
      ay <- rowsum(v, ys[w$id]); Py[as.integer(rownames(ay))] <- ay
      cx <- Re(fft(Mod(fft(Px))^2, inverse = TRUE)) / L    # sum_t' P_t' P_t'+t, both axes
      cy <- Re(fft(Mod(fft(Py))^2, inverse = TRUE)) / L
      Gacc   <- Gacc + ncomp * L * (cx + cy) / (2 * C)
      chiacc <- chiacc + ncomp * sum(v)^2 / C
      Facc   <- Facc + ncomp * (Mod(sum(v * phase[xs[w$id]]))^2 + Mod(sum(v * phase[ys[w$id]]))^2) / (2 * C)
    }
    Gb[b, ] <- Gacc / nc; chib[b] <- chiacc / nc; Fb[b] <- Facc / nc
  }
  list(L = L, beta = beta, G = Gb, chi = chib, F = Fb)
}

analyse <- function(r) {
  L <- r$L; nb <- nrow(r$G)
  est <- function(G, chi, F) {
    list(xi2 = sqrt(chi / F - 1) / (2 * sin(pi / L)),
         meff = acosh((G[1:(L/2-1)] + G[3:(L/2+1)]) / (2 * G[2:(L/2)])))   # t = 1 .. L/2 - 1
  }
  full <- est(colMeans(r$G), mean(r$chi), mean(r$F))
  tmin <- max(1, ceiling(full$xi2)); tmax <- min(L/2 - 1, floor(2.5 * full$xi2))
  jk <- t(sapply(seq_len(nb), function(b) {
    e <- est(colMeans(r$G[-b, , drop = FALSE]), mean(r$chi[-b]), mean(r$F[-b]))
    c(1 / e$xi2, mean(e$meff[tmin:tmax]), e$meff)
  }))
  jkerr <- function(x) sqrt((nb - 1) / nb * sum((x - mean(x))^2))
  list(beta = r$beta, L = L, xi2 = full$xi2, m2 = 1 / full$xi2, m2_err = jkerr(jk[, 1]),
       m_exp = mean(full$meff[tmin:tmax]), m_exp_err = jkerr(jk[, 2]), tmin = tmin, tmax = tmax,
       meff = full$meff, meff_err = apply(jk[, -(1:2), drop = FALSE], 2, jkerr))
}

3.4 Control: The Ising Model Against Onsager

The same code with one spin component is the 2D Ising model, whose axial inverse correlation length above \(T_c\) is known exactly: \(m(\beta) = -\ln\tanh\beta - 2\beta\) [Established]. This tests the cluster update, the improved estimators and the plateau procedure against an exact answer.

ising <- lapply(c(0.38, 0.40), function(b)
  analyse(run_on(48, b, ncomp = 1, n_therm = n_therm, n_bins = 50,
                 sweeps_per_bin = sweeps_per_bin, seed = 7)))
ising_df <- do.call(rbind, lapply(ising, function(a) {
  exact <- -log(tanh(a$beta)) - 2 * a$beta
  data.frame(beta = a$beta, L = a$L, m_exp = a$m_exp, err = a$m_exp_err, m_exact = exact,
             pull_sigma = (a$m_exp - exact) / a$m_exp_err)
}))
ising_df
beta L m_exp err m_exact pull_sigma
0.38 48 0.2523791 0.0020572 0.2541586 -0.8650290
0.40 48 0.1666967 0.0018464 0.1677183 -0.5533124

[Verified here] Both measured gaps agree with Onsager’s exact values within about one standard error (see pull_sigma).

3.5 O(4) Simulations

betas <- c(1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3)
o4 <- lapply(betas, function(b) {
  L <- 8 * ceiling(8 / m_asym(b) / 8)          # L ~ 8 / m, a multiple of 8
  analyse(run_on(L, b, n_therm = n_therm, n_bins = n_bins,
                 sweeps_per_bin = sweeps_per_bin, seed = round(100 * b)))
})
# Finite-size control: the same beta on a lattice ~1.7x larger
fse <- analyse(run_on(96, 1.9, n_therm = n_therm, n_bins = n_bins,
                      sweeps_per_bin = sweeps_per_bin, seed = 1901))
res_df <- do.call(rbind, lapply(o4, function(a) data.frame(
  beta = a$beta, L = a$L, xi2 = a$xi2, L_over_xi = a$L / a$xi2,
  m_exp = a$m_exp, m_exp_err = a$m_exp_err, m2 = a$m2, m2_err = a$m2_err,
  m_asym = m_asym(a$beta))))
res_df <- res_df %>% mutate(R = m_exp / m_asym, R_err = m_exp_err / m_asym)
knitr::kable(res_df, digits = c(2, 0, 2, 1, 5, 5, 5, 5, 5, 4, 4),
             caption = "O(4) mass gap: Monte Carlo (exponential and second-moment) versus the asymptotic formula")
O(4) mass gap: Monte Carlo (exponential and second-moment) versus the asymptotic formula
beta L xi2 L_over_xi m_exp m_exp_err m2 m2_err m_asym R R_err
1.7 32 3.74 8.6 0.26817 0.00235 0.26744 0.00260 0.26600 1.0082 0.0088
1.8 48 4.77 10.1 0.21076 0.00176 0.20962 0.00184 0.19992 1.0542 0.0088
1.9 56 6.10 9.2 0.16409 0.00130 0.16404 0.00136 0.15002 1.0937 0.0087
2.0 72 7.97 9.0 0.12513 0.00104 0.12553 0.00106 0.11242 1.1130 0.0092
2.1 96 10.34 9.3 0.09639 0.00069 0.09673 0.00069 0.08414 1.1456 0.0082
2.2 128 14.15 9.0 0.07063 0.00050 0.07067 0.00053 0.06290 1.1229 0.0080
2.3 176 18.98 9.3 0.05302 0.00035 0.05269 0.00035 0.04698 1.1286 0.0074 
a19 <- o4[[which(betas == 1.9)]]
fse_pull <- (fse$m_exp - a19$m_exp) / sqrt(fse$m_exp_err^2 + a19$m_exp_err^2)
data.frame(L = c(a19$L, fse$L), m_exp = c(a19$m_exp, fse$m_exp), err = c(a19$m_exp_err, fse$m_exp_err),
           row.names = c("L ~ 9 xi", "L ~ 16 xi"))
L m_exp err
L ~ 9 xi 56 0.1640869 0.0013024
L ~ 16 xi 96 0.1659983 0.0006949

In the finite-size control at \(\beta=1.9\), the larger lattice differs by 1.3 combined standard errors, so \(L/\xi\approx 9\) is adequate at this precision. The exponential and second-moment masses agree within errors at every \(\beta\), as expected from the three-particle threshold.

a22 <- o4[[which(betas == 2.2)]]
meff_df <- data.frame(t = seq_along(a22$meff), m = a22$meff, err = a22$meff_err) %>% filter(t <= 50)
ggplot(meff_df, aes(t, m)) +
  annotate("rect", xmin = a22$tmin, xmax = a22$tmax, ymin = -Inf, ymax = Inf, alpha = 0.12, fill = "steelblue") +
  geom_hline(yintercept = a22$m_exp, colour = "steelblue", linewidth = 0.8) +
  geom_hline(yintercept = m_asym(2.2), colour = "firebrick", linetype = "dashed", linewidth = 0.8) +
  geom_errorbar(aes(ymin = m - err, ymax = m + err), width = 0.4, colour = "grey40") +
  geom_point(colour = "black", size = 1.6) +
  annotate("text", x = 44, y = m_asym(2.2) - 0.0015, label = "asymptotic formula", colour = "firebrick") +
  annotate("text", x = 25, y = a22$m_exp + 0.0045, label = "plateau fit (blue line)", colour = "steelblue") +
  coord_cartesian(ylim = c(0.055, 0.085)) +
  labs(title = expression("Effective mass at " * beta * " = 2.2 (L = 128)"),
       subtitle = "Shaded band: plateau window [xi_2, 2.5 xi_2]",
       x = "time separation t", y = expression(m[eff](t)))

3.6 Results

fit_df <- res_df %>% filter(beta >= 2.0)
w <- 1 / fit_df$R_err^2; u <- 1 / fit_df$beta
# Constrained form R = 1 + c / beta (the expected leading correction)
c_hat <- sum(w * (fit_df$R - 1) * u) / sum(w * u^2)
c_err <- 1 / sqrt(sum(w * u^2))
chi2_c <- sum(w * (fit_df$R - 1 - c_hat * u)^2)
# Unconstrained form R = a + c / beta (tests whether the intercept is 1)
X <- cbind(1, u); A <- solve(t(X) %*% (w * X)); ab <- A %*% t(X) %*% (w * fit_df$R)
a_hat <- ab[1]; a_err <- sqrt(A[1, 1])

# Local decay rate between neighbouring betas versus the asymptotic finite difference
slope_df <- data.frame(
  beta_mid = head(res_df$beta, -1) + diff(res_df$beta) / 2,
  rate = -diff(log(res_df$m_exp)) / diff(res_df$beta),
  rate_err = sqrt(head(res_df$m_exp_err / res_df$m_exp, -1)^2 + tail(res_df$m_exp_err / res_df$m_exp, -1)^2) / diff(res_df$beta),
  rate_asym = -diff(log(m_asym(res_df$beta))) / diff(res_df$beta))

# Widest weak-coupling span, beta = 2.0 to 2.3: less noisy than adjacent pairs
i0 <- which(res_df$beta == 2.0); i1 <- which(res_df$beta == 2.3); db <- res_df$beta[i1] - res_df$beta[i0]
wide_rate <- -(log(res_df$m_exp[i1]) - log(res_df$m_exp[i0])) / db
wide_err  <- sqrt((res_df$m_exp_err[i0] / res_df$m_exp[i0])^2 + (res_df$m_exp_err[i1] / res_df$m_exp[i1])^2) / db
wide_asym <- pi - log(res_df$beta[i1] / res_df$beta[i0]) / (2 * db)
p_c <- pchisq(chi2_c, df = nrow(fit_df) - 1, lower.tail = FALSE)
bb <- seq(1.65, 2.35, length.out = 200)
p1 <- ggplot() +
  geom_line(data = data.frame(beta = bb, m = m_asym(bb)), aes(beta, m),
            colour = "firebrick", linetype = "dashed", linewidth = 0.8) +
  geom_errorbar(data = res_df, aes(beta, ymin = m_exp - m_exp_err, ymax = m_exp + m_exp_err),
                width = 0.015, colour = "black") +
  geom_point(data = res_df, aes(beta, m_exp), colour = "black", size = 2) +
  scale_y_log10() +
  labs(title = "Mass gap: Monte Carlo (points) vs 32 exp(pi/4 - 1/2) sqrt(beta) exp(-pi beta) (dashed)",
       x = expression(beta), y = expression(m[lat]))

p2 <- ggplot(res_df, aes(beta, R)) +
  geom_hline(yintercept = 1, colour = "firebrick", linetype = "dashed") +
  geom_line(data = data.frame(beta = seq(2, 2.35, length.out = 50)) %>% mutate(R = 1 + c_hat / beta),
            aes(beta, R), colour = "steelblue", linewidth = 0.8) +
  geom_errorbar(aes(ymin = R - R_err, ymax = R + R_err), width = 0.015, colour = "black") +
  geom_point(colour = "black", size = 2) +
  labs(title = "Ratio to the asymptotic formula",
       subtitle = sprintf("Blue: fit 1 + c/beta, c = %.3f, chi2 = %.1f / %d dof", c_hat, chi2_c, nrow(fit_df) - 1),
       x = expression(beta), y = expression(m[MC] / m[asym]))

p3 <- ggplot(slope_df, aes(beta_mid, rate)) +
  geom_line(data = data.frame(beta_mid = bb, rate = pi - 1 / (2 * bb)), aes(beta_mid, rate),
            colour = "firebrick", linetype = "dashed", linewidth = 0.8) +
  geom_errorbar(aes(ymin = rate - rate_err, ymax = rate + rate_err), width = 0.015, colour = "black") +
  geom_point(colour = "black", size = 2) +
  annotate("segment", x = 2.0, xend = 2.3, y = wide_rate, yend = wide_rate, colour = "steelblue", linewidth = 0.8) +
  annotate("errorbar", x = 2.15, ymin = wide_rate - wide_err, ymax = wide_rate + wide_err,
           width = 0.02, colour = "steelblue", linewidth = 0.8) +
  annotate("point", x = 2.15, y = wide_rate, colour = "steelblue", shape = 18, size = 4) +
  labs(title = "Local decay rate -d ln m / d beta",
       subtitle = "Black: adjacent pairs. Blue: span 2.0-2.3. Dashed: pi - 1/(2 beta)",
       x = expression(beta), y = "rate")

grid.arrange(p1, arrangeGrob(p2, p3, ncol = 2), nrow = 2)

Five findings follow from the table and the figure.

  1. The exponential rate matches the theorem’s form. Over the widest weak-coupling span, \(\beta = 2.0\) to \(2.3\), the measured rate \(-\Delta\ln m/\Delta\beta\) is \(2.86 \pm 0.04\), against 2.91 for the asymptotic formula, a 1.6 % difference. Adjacent-pair rates are noisier, but they climb from about 2.4 at \(\beta\approx 1.75\) toward the asymptotic curve. [Verified here, within the accessible range]
  2. The prefactor is not reached at these correlation lengths. For \(\beta\ge 2.0\) the ratio \(m_{\mathrm{MC}}/m_{\mathrm{asym}}\) is 1.113–1.146.
  3. A single correction term does not describe the offset well. The constrained fit \(1+c/\beta\) gives \(c = 0.276 \pm 0.009\), but with \(\chi^2 =\) 11.8 for 3 degrees of freedom (\(p =\) 0.008). An unconstrained fit \(a + c/\beta\) gives \(a = 1.16 \pm 0.08\), which is too imprecise to confirm or refute the intercept \(a = 1\). [Not decisively tested]
  4. An offset of this size is plausible at these correlation lengths. Large-scale O(3) studies found deviations from asymptotic scaling of about 25 % at \(\xi\sim 10^2\), falling to about 4 % only at \(\xi\sim 10^5\), and O(4) is reported to approach asymptotic scaling sooner. A percent-level test of the prefactor would require correlation lengths of \(10^2\)–\(10^3\) and finite-size-scaling extrapolation, which is beyond a pure-R calculation.
  5. The strong-coupling points (\(\beta\le 1.9\)) do not follow the \(1/\beta\) form. Including them would bias any extrapolation, which is why the fits use \(\beta\ge 2.0\).



4 Part II: Bloch’s Law and Its Lattice Correction

4.1 The Claim

[Claimed: unrefereed] For the 3D nearest-neighbour quantum Heisenberg ferromagnet at every spin \(S\), the spontaneous magnetization \(m_S(\beta)\) satisfies

\[ S - m_S(\beta) \;=\; n^{\mathrm{sw}}_S(\beta) + o\!\left((\beta S)^{-5/2}\right), \]

where \(n^{\mathrm{sw}}_S\) is the ideal-magnon density built on the lattice dispersion \(\varepsilon_S(\mathbf k)=2S\sum_j(1-\cos k_j)\). This gives

\[ S - m_S(\beta) \;=\; \frac{\zeta(3/2)}{8\pi^{3/2}}(\beta S)^{-3/2} + \frac{3\,\zeta(5/2)}{128\,\pi^{3/2}}(\beta S)^{-5/2} + o\!\left((\beta S)^{-5/2}\right). \]

The manuscript’s own introduction notes that the \(T^{5/2}\) and \(T^{7/2}\) terms already arise from the lattice dispersion of independent magnons [Established]. The new content is the rigorous statement that the interacting quantum model inherits them.

4.2 The Ideal-Magnon Density as an Exact Lattice Sum

Expanding the Bose factor gives \(n_B(x)=\sum_{n\ge1}e^{-nx}\). Each Brillouin-zone direction then factorizes through \(\int_{-\pi}^{\pi}\frac{dk}{2\pi}e^{z\cos k}=I_0(z)\), which gives the exact representation

\[ n^{\mathrm{sw}}(t) \;=\; \sum_{n=1}^{\infty}\Big[e^{-x_n}I_0(x_n)\Big]^3,\qquad x_n = \frac{2n}{t},\quad t\equiv\frac{k_BT}{JS}. \]

The sum converges like \(n^{-3/2}\). The first 4000 terms are summed exactly and the tail is added from the large-\(x\) expansion of \(I_0\) via Euler–Maclaurin Hurwitz-zeta sums. R’s besselI underflows to zero for arguments above about \(10^5\), so arguments above 500 use the asymptotic series. Neglecting that underflow produces a spurious 2 % error at the lowest temperatures.

Inserting \(e^{-x}I_0(x)=(2\pi x)^{-1/2}\left[1+\tfrac{1}{8x}+\tfrac{9}{128x^2}+\dots\right]\) term by term gives Dyson’s free-magnon series [Established], with \(\tau = t/4\pi\):

\[ n^{\mathrm{sw}} = \zeta(\tfrac32)\,\tau^{3/2} + \tfrac{3\pi}{4}\,\zeta(\tfrac52)\,\tau^{5/2} + \tfrac{33\pi^2}{32}\,\zeta(\tfrac72)\,\tau^{7/2} + O(\tau^{9/2}). \]

# Hurwitz zeta sum_{n >= a} n^-s by Euler-Maclaurin (base R has no zeta function)
zeta_em <- function(s, a = 1, M = 60) {
  n <- a:(a + M - 1); N <- a + M
  sum(n^-s) + N^(1-s)/(s-1) + N^-s/2 + s*N^(-s-1)/12 - s*(s+1)*(s+2)*N^(-s-3)/720 +
    s*(s+1)*(s+2)*(s+3)*(s+4)*N^(-s-5)/30240
}
# exp(-x) I0(x). besselI() returns 0 for x >~ 1e5, so large arguments use the asymptotic series.
I0s <- function(x) {
  out <- numeric(length(x)); big <- x > 500
  out[!big] <- besselI(x[!big], 0, expon.scaled = TRUE)
  xb <- x[big]
  out[big] <- (1 + 1/(8*xb) + 9/(128*xb^2) + 75/(1024*xb^3) + 3675/(32768*xb^4)) / sqrt(2*pi*xb)
  out
}
# Exact ideal-magnon density per site, simple cubic, t = k_B T / (J S)
n_sw <- function(t, N0 = 4000) {
  head <- sum(I0s(2 * (1:N0) / t)^3)
  tail <- (4*pi/t)^-1.5 * (zeta_em(1.5, N0+1) + (3/8)*(t/2)*zeta_em(2.5, N0+1) +
                           (33/128)*(t/2)^2*zeta_em(3.5, N0+1))
  head + tail
}
z32 <- zeta_em(1.5); z52 <- zeta_em(2.5); z72 <- zeta_em(3.5)
dyson <- function(t, k) {
  tau <- t / (4*pi)
  sum(c(z32 * tau^1.5, 3*pi/4 * z52 * tau^2.5, 33*pi^2/32 * z72 * tau^3.5)[1:k])
}

data.frame(quantity = c("zeta(3/2)", "zeta(5/2)", "zeta(7/2)"),
           computed = c(z32, z52, z72),
           reference = c(2.612375348685488, 1.341487257250917, 1.126733867317057)) %>%
  mutate(abs_diff = abs(computed - reference))
quantity computed reference abs_diff
zeta(3/2) 2.612375 2.612375 0
zeta(5/2) 1.341487 1.341487 0
zeta(7/2) 1.126734 1.126734 0

4.3 The Coefficients Against the Manuscript

data.frame(
  term = c("(beta S)^-3/2 coefficient", "(beta S)^-5/2 coefficient"),
  manuscript = c(z32 / (8 * pi^1.5), 3 * z52 / (128 * pi^1.5)),
  from_dyson_series = c(z32 / (4*pi)^1.5, (3*pi/4) * z52 / (4*pi)^2.5)) %>%
  mutate(relative_difference = (manuscript - from_dyson_series) / from_dyson_series)
term manuscript from_dyson_series relative_difference
(beta S)^-3/2 coefficient 0.0586436 0.0586436 0
(beta S)^-5/2 coefficient 0.0056464 0.0056464 0

[Verified here] Under the manuscript’s normalization (each edge counted once, \(J=1\)), both claimed coefficients are exactly the first two terms of the ideal-magnon series.

4.4 Residual Scaling

If the coefficients are right, the relative residual after \(k\) terms should scale as \(\tau^{k}\). Concretely, the residual after Bloch alone should scale as \(\tau^1\), after the lattice correction as \(\tau^2\), and after the \(\tau^{7/2}\) term as \(\tau^3\).

ts <- 10^seq(-2, 0, by = 0.1)
bloch_df <- data.frame(t = ts) %>% rowwise() %>%
  mutate(exact = n_sw(t),
         `Bloch only` = abs(exact - dyson(t, 1)) / exact,
         `+ lattice correction` = abs(exact - dyson(t, 2)) / exact,
         `+ tau^(7/2) term` = abs(exact - dyson(t, 3)) / exact) %>% ungroup()

slopes <- sapply(c("Bloch only", "+ lattice correction", "+ tau^(7/2) term"), function(col) {
  d <- bloch_df %>% filter(t <= 0.1); coef(lm(log(d[[col]]) ~ log(d$t)))[2] })

bloch_long <- bloch_df %>% select(-exact) %>%
  pivot_longer(-t, names_to = "truncation", values_to = "rel_residual") %>%
  mutate(truncation = factor(truncation, levels = c("Bloch only", "+ lattice correction", "+ tau^(7/2) term")))
cols <- c("Bloch only" = "#440154", "+ lattice correction" = "#21908C", "+ tau^(7/2) term" = "#E3A600")
ggplot(bloch_long, aes(t, rel_residual)) +
  geom_line(aes(colour = truncation), linewidth = 0.9) +
  geom_point(aes(colour = truncation), size = 1.4) +
  scale_colour_manual(values = cols) +
  scale_x_log10() + scale_y_log10() +
  labs(title = "Relative residual of truncated Dyson series against the exact ideal-magnon density",
       subtitle = sprintf("Fitted log-log slopes for t <= 0.1: %.2f, %.2f, %.2f (expected 1, 2, 3)",
                          slopes[1], slopes[2], slopes[3]),
       x = expression(t == k[B] * T / (J * S)), y = "relative residual", colour = NULL)

[Verified here] The residuals fall with slopes 1, 2 and 3 over two decades of temperature. The \(\tau^{5/2}\) coefficient is therefore exactly the lattice correction, and nothing at order \(\tau^{5/2}\) is missing from the ideal-magnon reference function.

4.5 Why Interactions Should Not Enter at This Order: Two-Magnon Exact Diagonalization

The theorem’s non-trivial content is that magnon–magnon interactions do not alter the \(T^{3/2}\) or \(T^{5/2}\) terms. Dyson showed perturbatively that interactions first enter at \(T^4\), because the low-momentum scattering amplitude vanishes like \(\mathbf k_1\cdot\mathbf k_2\) [Established]. A finite-size signature of this can be computed exactly.

The construction is as follows. For \(S=\tfrac12\) on an \(L^3\) torus, the two-magnon sector at total momentum \(\mathbf K = 0\) reduces to a hopping problem in the relative coordinate \(\mathbf r\neq 0\). The diagonal term is \(6 - [\mathbf r \text{ nearest neighbour}]\) and the hopping amplitude is \(-1\) to each neighbour, with even parity under \(\mathbf r\to-\mathbf r\). The free reference energies are \(\varepsilon(\mathbf k)+\varepsilon(-\mathbf k)\).

For a generic short-range interaction with scattering length \(a\), the lowest \(\mathbf K=0\) level shifts by an amount \(\propto a/L^3\). For the Heisenberg ferromagnet, two predictions follow.

  • The zero-momentum level is exactly \(E=0\), because it is \((S^-_{\mathrm{tot}})^2|\Uparrow\rangle\).
  • The lowest scattering level, the pair \((\mathbf k,-\mathbf k)\) with \(|\mathbf k|=2\pi/L\), should shift as \(k^2/N\propto L^{-5}\).
two_magnon_K0 <- function(L) {
  N <- L^3
  co <- expand.grid(x = 0:(L-1), y = 0:(L-1), z = 0:(L-1))
  idx <- function(x, y, z) (x %% L) + L*((y %% L) + L*(z %% L)) + 1L
  r_sites <- 2:N; pos <- integer(N); pos[r_sites] <- seq_along(r_sites)
  H <- matrix(0, N-1, N-1)
  sh <- rbind(c(1,0,0), c(-1,0,0), c(0,1,0), c(0,-1,0), c(0,0,1), c(0,0,-1))
  for (a in seq_along(r_sites)) {
    s <- r_sites[a]; nn0 <- 0
    for (k in 1:6) {
      t <- idx(co$x[s] + sh[k,1], co$y[s] + sh[k,2], co$z[s] + sh[k,3])
      if (t == 1L) nn0 <- 1 else H[a, pos[t]] <- H[a, pos[t]] - 1   # hard core at r = 0
    }
    H[a, a] <- 6 - nn0                                             # (J/2)(12 - 2[r nn]), J = 1
  }
  minus <- sapply(r_sites, function(s) idx(-co$x[s], -co$y[s], -co$z[s]))
  P <- matrix(0, N-1, N-1); P[cbind(seq_along(r_sites), pos[minus])] <- 1
  e <- eigen(H, symmetric = TRUE)
  parity <- colSums(e$vectors * (P %*% e$vectors))
  sort(e$values[parity > 0.5])                                     # bosonic (even) states
}

tm <- do.call(rbind, lapply(c(6, 8, 10, 12), function(L) {
  ev <- two_magnon_K0(L); free <- 2 * (1 - cos(2*pi/L))
  lv <- ev[2:4]; d <- abs(outer(lv, lv, "-")); diag(d) <- Inf
  pair <- which(d == min(d), arr.ind = TRUE)[1, ]; single <- setdiff(1:3, pair)
  data.frame(L = L, E_zero_mode = sprintf("%.1e", ev[1]), E_free = free,
             E_doublet = mean(lv[pair]), E_singlet = lv[single],
             shift_doublet = mean(lv[pair]) - free, shift_singlet = lv[single] - free)
}))
knitr::kable(tm, digits = c(0, 0, 6, 6, 6, 7, 7), caption = "Two-magnon K = 0 spectrum, S = 1/2")
Two-magnon K = 0 spectrum, S = 1/2
L E_zero_mode E_free E_doublet E_singlet shift_doublet shift_singlet
6 7.1e-15 1.000000 0.993987 1.022536 -0.0060134 0.0225365
8 -5.6e-15 0.585786 0.584937 0.591816 -0.0008497 0.0060300
10 -4.3e-15 0.381966 0.381783 0.384066 -0.0001828 0.0021003
12 -3.0e-15 0.267949 0.267897 0.268823 -0.0000517 0.0008740 
loc_slope <- diff(log(abs(tm$shift_singlet))) / diff(log(tm$L))
tm_long <- tm %>% select(L, shift_doublet, shift_singlet) %>%
  pivot_longer(-L, names_to = "level", values_to = "shift") %>% mutate(abs_shift = abs(shift))
ref <- data.frame(L = c(6, 12)) %>% mutate(abs_shift = abs(tm$shift_singlet[tm$L == 12]) * (L/12)^-5)
ggplot(tm_long, aes(L, abs_shift)) +
  geom_line(data = ref, aes(L, abs_shift), colour = "grey50", linetype = "dashed") +
  geom_line(aes(colour = level), linewidth = 0.9) + geom_point(aes(colour = level), size = 2) +
  scale_colour_manual(values = c(shift_doublet = "#21908C", shift_singlet = "#440154"),
                      labels = c(shift_doublet = "doublet", shift_singlet = "singlet")) +
  scale_x_log10(breaks = c(6, 8, 10, 12)) + scale_y_log10() +
  labs(title = "Interaction shift of the lowest two-magnon scattering level",
       subtitle = sprintf("Dashed: L^-5 reference. Local slopes (singlet): %s",
                          paste(sprintf("%.2f", loc_slope), collapse = ", ")),
       x = "L", y = "|E - E_free|", colour = NULL)

[Verified here] The zero-momentum two-magnon level is zero to machine precision at every \(L\). The singlet’s interaction shift falls off with a local slope steepening toward \(-5\), and the doublet’s falls off faster still. A generic short-range interaction would instead shift the lowest level by an amount \(\propto L^{-3}\).

[Not tested here] This illustrates the mechanism that makes the claimed theorem plausible. It is not a test of the theorem, which concerns the interacting many-magnon thermodynamics in the infinite volume, with the volume limit taken before the field derivative and \(T\to 0\) last. A direct numerical test would need quantum Monte Carlo at large \(L\) and very low \(T\).



5 Status Summary

Claim Test in this document Outcome Flag
O(4) prefactor equals \(32e^{\pi/4-1/2}\) Exact identity with HMN + \(\Lambda\) ratio Equal to machine precision [Verified here] (input formula [Established])
O(4) mass \(\propto\sqrt\beta\,e^{-\pi\beta}\) Wolff MC, \(\beta = 1.7\)–\(2.3\), Ising-validated Rate over \(\beta = 2.0\)–\(2.3\) within 1.6 % of the asymptotic rate [Verified here] within range
O(4) limiting prefactor (the theorem) Ratio fits 11–15 % offset at \(\xi\le 19\); a single \(1/\beta\) term fits poorly [Not decisively tested]
Bloch coefficient \(\zeta(3/2)/(8\pi^{3/2})\) Exact lattice sum + residual scaling Exact [Verified here] (ideal magnons)
Lattice correction \(3\zeta(5/2)/(128\pi^{3/2})\) Same Exact [Verified here] (ideal magnons)
Interacting deficit = ideal + \(o(T^{5/2})\) Two-magnon ED Mechanism illustrated (\(E_0 = 0\), \(L^{-5}\) shifts) [Not tested here]



6 References

  • OpenAI, openai/math repository (6 October 2026): family 215, Canonical O(3) continuum limit and exact O(4) mass asymptotics; family 271, Bloch’s law, its lattice correction, and the spherical magnetization law. Unrefereed manuscripts.
  • P. Hasenfratz, M. Maggiore, F. Niedermayer, Phys. Lett. B 245, 522 (1990); P. Hasenfratz, F. Niedermayer, Phys. Lett. B 245, 529 (1990). Exact mass gap of the O(N) σ-model.
  • R. G. Edwards, S. J. Ferreira, J. Goodman, A. D. Sokal, Multi-grid Monte Carlo III: Two-dimensional O(4)-symmetric nonlinear σ-model, arXiv:hep-lat/9112002. Conventions, the \(\Lambda\) ratio, and MC asymptotic-scaling tests.
  • S. Caracciolo, R. G. Edwards, A. Pelissetto, A. D. Sokal, Asymptotic scaling in the two-dimensional O(3) σ-model at correlation length \(10^5\), arXiv:hep-lat/9411009.
  • U. Wolff, Phys. Rev. Lett. 62, 361 (1989). Single-cluster algorithm.
  • L. Onsager, Phys. Rev. 65, 117 (1944).
  • F. Bloch, Z. Phys. 61, 206 (1930).
  • F. J. Dyson, Phys. Rev. 102, 1217 and 1230 (1956). Spin-wave interactions; \(T^4\) onset of interaction corrections.